linearfit

A module that provides algorithms for performing linear fit between sets of 2D points.

Authors:

Mihai Cara, Warren Hack

License:

LICENSE

tweakwcs.linearfit.build_fit_matrix(rot, scale=1)

Create an affine transformation matrix (2x2) from the provided rotation angle(s) and scale(s):

\[\begin{split}M = \begin{bmatrix} s_x \cos(\theta_x) & s_y \sin(\theta_y) \\ -s_x \sin(\theta_x) & s_y \cos(\theta_y) \end{bmatrix}\end{split}\]

Parameters:

rot: tuple, float, optional

Rotation angle in degrees. Two values (one for each axis) can be provided as a tuple.

scale: tuple, float, optional

Scale of the liniar transformation. Two values (one for each axis) can be provided as a tuple.

Returns:

matrix: numpy.ndarray

A 2x2 numpy.ndarray containing coefficients of a liniear transformation.

tweakwcs.linearfit.iter_linear_fit(xy, uv, wxy=None, wuv=None, fitgeom='general', center=None, nclip=3, sigma=(3.0, 'rmse'), clip_accum=False)

Compute linear transformation parameters that “best” (in the sense of minimizing residuals) transform uv source position to xy sources iteratively using sigma-clipping.

More precisely, this functions attempts to find a 2x2 matrix F and a shift vector s that minimize the residuals between the transformed reference source coordinates uv

(1)\[\mathbf{xy}'_k = \mathbf{F}\cdot(\mathbf{uv}_k-\mathbf{c})+\ \mathbf{s} + \mathbf{c}\]

and the “observed” source positions xy:

(2)\[\epsilon^2 = \Sigma_k w_k \|\mathbf{xy}_k-\mathbf{xy}'_k\|^2.\]

In the above equations, \(\mathbf{F}\) is a 2x2 matrix while \(\mathbf{xy}_k\) and \(\mathbf{uv}_k\) are the position coordinates of the k-th source (row in input xy and uv arrays).

One of the two catalogs (xy or uv) contains what we refer to as “image” source positions and the other one as “reference” source positions. The meaning assigned to xy and uv parameters are up to the caller of this function.

Parameters:

xy: numpy.ndarray

A (N, 2)-shaped array of source positions (one 2-coordinate position per line).

uv: numpy.ndarray

A (N, 2)-shaped array of source positions (one 2-coordinate position per line). This array must have the same length (shape) as the xy array.

wxy: numpy.ndarray, None, optional

A 1-dimensional array of weights of the same length (N) as xy array indicating how much a given coordinate should be weighted in the fit. If not provided or set to None, all positions will be contribute equally to the fit if wuv is also set to None. See Notes section for more details.

wuv: numpy.ndarray, None, optional

A 1-dimensional array of weights of the same length (N) as xy array indicating how much a given coordinate should be weighted in the fit. If not provided or set to None, all positions will be contribute equally to the fit if wxy is also set to None. See Notes section for more details.

fitgeom: {‘shift’, ‘rshift’, ‘rscale’, ‘general’}, optional

The fitting geometry to be used in fitting the matched object lists. This parameter is used in fitting the shifts (offsets), rotations and/or scale changes from the matched object lists. The ‘general’ fit geometry allows for independent scale and rotation for each axis.

center: tuple, list, numpy.ndarray, None, optional

A list-like container with two X- and Y-positions of the center (origin) of rotations in the uv and xy coordinate frames. If not provided, center is estimated as a (weighted) mean position in the uv frame.

nclip: int, None, optional

Number (a non-negative integer) of clipping iterations in fit. Clipping will be turned off if nclip is either None or 0.

sigma: float, tuple of the form (float, str), optional

When a tuple is provided, first value (a positive number) indicates the number of “fit error estimates” to use for clipping. The second value (a string) indicates the statistic to be used for “fit error estimate”. Currently the following values are supported: 'rmse', 'mae', and 'std' - see Notes section for more details.

When sigma is a single number, it must be a positive number and the default error estimate 'rmse' is assumed.

This parameter is ignored when nclip is either None or 0.

clip_accum: bool, optional

Indicates whether or not to reset the list of “bad” (clipped out) sources after each clipping iteration. When set to True the list only grows with each iteration as “bad” positions never re-enter the pool of available position for the fit. By default the list of “bad” source positions is purged at each iteration. This parameter is ignored when nclip is either None or 0.

Returns:

fit: dict

• 'shift': A numpy.ndarray with two components of the computed shift.

• 'shift_ld': A numpy.ndarray with two components of the computed shift of type numpy.longdouble.

• 'matrix': A 2x2 numpy.ndarray with the computed generalized rotation matrix.

• 'matrix_ld': A 2x2 numpy.ndarray with the computed generalized rotation matrix of type numpy.longdouble.

• 'proper_rot': Rotation angle (degree) as if the rotation is proper.

• 'rot': A tuple of (rotx, roty) - the rotation angles with regard to the X and Y axes.

• '<rot>': Arithmetic mean of the angles of rotation around X and Y axes.

• 'scale': A tuple of (sx, sy) - scale change in the direction of the X and Y axes.

• '<scale>': Geometric mean of scales sx and sy.

• 'skew': Computed skew.

• 'proper': a boolean indicating whether the rotation is proper.

• 'fitgeom': Fit geometry (allowed transformations) used for fitting data (to minimize residuals). This is copy of the input argument fitgeom.

• 'center': Center of rotation

• 'center_ld': Center of rotation as a numpy.longdouble.

• 'fitmask': A boolean array indicating which source positions where used for fitting (True) and which were clipped out (False). NOTE For weighted fits, positions with zero weights are automatically excluded from the fits.

• 'eff_nclip': Effective number of clipping iterations

• 'rmse': Root-Mean-Square Error

• 'mae': Mean Absolute Error

• 'std': Standard Deviation of the residuals

• 'resids': An array of residuals of the fit. NOTE: Only the residuals for the “valid” points are reported here. Therefore the length of this array may be smaller than the length of input arrays of positions.

Notes

Weights

Weights can be provided for both “image” source positions and “reference” source positions. When no weights are given, all positions are weighted equally. When only one set of positions have weights (i.e., either wxy or wuv is not None) then weights in (2) are set to be equal to the provided set of weights. When weights for both “image” source positions and “reference” source positions are provided, then the combined weight that is used in (2) is computed as:

\[1/w = 1/w_{xy} + 1/w_{uv}.\]

Statistics for clipping

Several statistics are available for clipping iterations and all of them are reported in the returned fit dictionary regardless of the setting in sigma:

\[\mathrm{RMSE} = \sqrt{\Sigma_k w_k \|\mathbf{r}_k\|^2}\]

\[\mathrm{MAE} = \sqrt{\Sigma_k w_k \|\mathbf{r}_k\|}\]

\[\mathrm{STD} = \sqrt{\Sigma_k w_k \|\mathbf{r}_k - \ \mathbf{\overline{r}}\|^2}/(1-V_2)\]

where \(\mathbf{r}_k=\mathbf{xy}_k-\mathbf{xy}'_k\), \(\Sigma_k w_k = 1\), and \(V_2=\Sigma_k w_k^2\).