LSQ Image Intensity Matching¶
A module that provides main API for optimal (LSQ) “matching” of weighted N-dimensional image intensity data using (multivariate) polynomials.
- Author
Mihai Cara (contact: help@stsci.edu)
- jwst.wiimatch.match.match_lsq(images, masks=None, sigmas=None, degree=0, center=None, image2world=None, center_cs='image', ext_return=False, solver='RLU')[source]¶
Compute coefficients of (multivariate) polynomials that once subtracted from input images would provide image intensity matching in the least squares sense.
- Parameters
images (list of numpy.ndarray) – A list of 1D, 2D, etc.
numpy.ndarray
data array whose “intensities” must be “matched”. All arrays must have identical shapes.masks (list of numpy.ndarray, None) – A list of
numpy.ndarray
arrays of same length asimages
. Non-zero mask elements indicate valid data in the correspondingimages
array. Mask arrays must have identical shape to that of the arrays in inputimages
. Default value ofNone
indicates that all pixels in input images are valid.sigmas (list of numpy.ndarray, None) – A list of
numpy.ndarray
data array of same length asimages
representing the uncertainties of the data in the corresponding array inimages
. Uncertainty arrays must have identical shape to that of the arrays in inputimages
. The default value ofNone
indicates that all pixels will be assigned equal weights.degree (iterable, int) – A list of polynomial degrees for each dimension of data arrays in
images
. The length of the input list must match the dimensionality of the input images. When a single integer number is provided, it is assumed that the polynomial degree in each dimension is equal to that integer.center (iterable, None, optional) – An iterable of length equal to the number of dimensions in
image_shape
that indicates the center of the coordinate system in image coordinates whencenter_cs
is'image'
otherwise center is assumed to be in world coordinates (whencenter_cs
is'world'
). Whencenter
isNone
thencenter
is set to the middle of the “image” ascenter[i]=image_shape[i]//2
. Ifimage2world
is notNone
andcenter_cs
is'image'
, then supplied center will be converted to world coordinates.image2world (function, None, optional) – Image-to-world coordinates transformation function. This function must be of the form
f(x,y,z,...)
and accept a number of argumentsnumpy.ndarray
arguments equal to the dimensionality of images.center_cs ({'image', 'world'}, optional) – Indicates whether
center
is in image coordinates or in world coordinates. This parameter is ignored whencenter
is set toNone
: it is assumed to beFalse
.center_cs
cannot be'world'
whenimage2world
isNone
unlesscenter
isNone
.ext_return (bool, optional) – Indicates whether this function should return additional values besides optimal polynomial coefficients (see
bkg_poly_coeff
return value below) that match image intensities in the LSQ sense. See Returns section for more details.solver ({'RLU', 'PINV'}, optional) – Specifies method for solving the system of equations.
- Returns
bkg_poly_coeff (numpy.ndarray) – When
nimages
isNone
, this function returns a 1Dnumpy.ndarray
that holds the solution (polynomial coefficients) to the system.When
nimages
is notNone
, this function returns a 2Dnumpy.ndarray
that holds the solution (polynomial coefficients) to the system. The solution is grouped by image.a (numpy.ndarray) – A 2D
numpy.ndarray
that holds the coefficients of the linear system of equations. This value is returned only whenext_return
isTrue
.b (numpy.ndarray) – A 1D
numpy.ndarray
that holds the free terms of the linear system of equations. This value is returned only whenext_return
isTrue
.coord_arrays (list) – A list of
numpy.ndarray
coordinate arrays each ofimage_shape
shape. This value is returned only whenext_return
isTrue
.eff_center (tuple) – A tuple of coordinates of the effective center as used in generating coordinate arrays. This value is returned only when
ext_return
isTrue
.coord_system ({‘image’, ‘world’}) – Coordinate system of the coordinate arrays and returned
center
value. This value is returned only whenext_return
isTrue
.
Notes
match_lsq()
builds a system of linear equations\[a \cdot c = b\]whose solution \(c\) is a set of coefficients of (multivariate) polynomials that represent the “background” in each input image (these are polynomials that are “corrections” to intensities of input images) such that the following sum is minimized:
\[L = \sum^N_{n,m=1,n \neq m} \sum_k \frac{\left[I_n(k) - I_m(k) - P_n(k) + P_m(k)\right]^2} {\sigma^2_n(k) + \sigma^2_m(k)}.\]In the above equation, index \(k=(k_1,k_2,...)\) labels a position in input image’s pixel grid [NOTE: all input images share a common pixel grid].
“Background” polynomials \(P_n(k)\) are defined through the corresponding coefficients as:
\[P_n(k_1,k_2,...) = \sum_{d_1=0,d_2=0,...}^{D_1,D_2,...} c_{d_1,d_2,...}^n \cdot k_1^{d_1} \cdot k_2^{d_2} \cdot \ldots .\]Coefficients \(c_{d_1,d_2,...}^n\) are arranged in the vector \(c\) in the following order:
\[(c_{0,0,\ldots}^1,c_{1,0,\ldots}^1,\ldots,c_{0,0,\ldots}^2,c_{1,0,\ldots}^2,\ldots).\]match_lsq()
returns coefficients of the polynomials that minimize L.Examples
>>> import numpy as np >>> im1 = np.zeros((5, 5, 4), dtype=float) >>> cbg = 1.32 * np.ones_like(im1) >>> ind = np.indices(im1.shape, dtype=float) >>> im3 = cbg + 0.15 * ind[0] + 0.62 * ind[1] + 0.74 * ind[2] >>> mask = np.ones_like(im1, dtype=np.int8) >>> sigma = np.ones_like(im1, dtype=float) >>> match_lsq([im1, im3], [mask, mask], [sigma, sigma], ... degree=(1, 1, 1), center=(0, 0, 0)) array([[-6.60000000e-01, -7.50000000e-02, -3.10000000e-01, -1.19371180e-15, -3.70000000e-01, -1.62003744e-15, -1.10844667e-15, 5.11590770e-16], [ 6.60000000e-01, 7.50000000e-02, 3.10000000e-01, 1.19371180e-15, 3.70000000e-01, 1.62003744e-15, 1.10844667e-15, -5.11590770e-16]])