LSQ Equation Construction and Solving

A module that provides core algorithm for optimal matching of backgrounds of N-dimensional images using (multi-variate) polynomials.

Author:Mihai Cara (contact: help@stsci.edu)
License:LICENSE
jwst.wiimatch.lsq_optimizer.build_lsq_eqs(images, masks, sigmas, degree, center=None, image2world=None, center_cs=u'image')[source]

Build system of linear equations whose solution would provide image intensity matching in the least squares sense.

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
A list of numpy.ndarray arrays of same length as images. Non-zero mask elements indicate valid data in the corresponding images array. Mask arrays must have identical shape to that of the arrays in input images.
sigmas : list of numpy.ndarray
A list of numpy.ndarray data array of same length as images representing the uncertainties of the data in the corresponding array in images. Uncertainty arrays must have identical shape to that of the arrays in input images.
degree : iterable
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.
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 when center_cs is 'image' otherwise center is assumed to be in world coordinates (when center_cs is 'world'). When center is None then center is set to the middle of the “image” as center[i]=image_shape[i]//2. If image2world is not None and center_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 arguments numpy.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 when center is set to None: it is assumed to be False. center_cs cannot be 'world' when image2world is None unless center is None.
a : numpy.ndarray
A 2D numpy.ndarray that holds the coefficients of the linear system of equations.
b : numpy.ndarray
A 1D numpy.ndarray that holds the free terms of the linear system of equations.
coord_arrays : list
A list of numpy.ndarray coordinate arrays each of image_shape shape.
eff_center : tuple
A tuple of coordinates of the effective center as used in generating coordinate arrays.
coord_system : {‘image’, ‘world’}
Coordinate system of the coordinate arrays and returned center value.

build_lsq_eqs() 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).

>>> import wiimatch
>>> import numpy as np
>>> im1 = np.zeros((5, 5, 4), dtype=np.float)
>>> cbg = 1.32 * np.ones_like(im1)
>>> ind = np.indices(im1.shape, dtype=np.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=np.float)
>>> a, b, ca, ef, cs = wiimatch.lsq_optimizer.build_lsq_eqs([im1, im3],
... [mask, mask], [sigma, sigma], degree=(1,1,1), center=(0,0,0))
>>> print(a)
[[   50.   100.   100.   200.    75.   150.   150.   300.   -50.  -100.
    -100.  -200.   -75.  -150.  -150.  -300.]
    [  100.   300.   200.   600.   150.   450.   300.   900.  -100.  -300.
    -200.  -600.  -150.  -450.  -300.  -900.]
    [  100.   200.   300.   600.   150.   300.   450.   900.  -100.  -200.
    -300.  -600.  -150.  -300.  -450.  -900.]
    [  200.   600.   600.  1800.   300.   900.   900.  2700.  -200.  -600.
    -600. -1800.  -300.  -900.  -900. -2700.]
    [   75.   150.   150.   300.   175.   350.   350.   700.   -75.  -150.
    -150.  -300.  -175.  -350.  -350.  -700.]
    [  150.   450.   300.   900.   350.  1050.   700.  2100.  -150.  -450.
    -300.  -900.  -350. -1050.  -700. -2100.]
    [  150.   300.   450.   900.   350.   700.  1050.  2100.  -150.  -300.
    -450.  -900.  -350.  -700. -1050. -2100.]
    [  300.   900.   900.  2700.   700.  2100.  2100.  6300.  -300.  -900.
    -900. -2700.  -700. -2100. -2100. -6300.]
    [  -50.  -100.  -100.  -200.   -75.  -150.  -150.  -300.    50.   100.
    100.   200.    75.   150.   150.   300.]
    [ -100.  -300.  -200.  -600.  -150.  -450.  -300.  -900.   100.   300.
    200.   600.   150.   450.   300.   900.]
    [ -100.  -200.  -300.  -600.  -150.  -300.  -450.  -900.   100.   200.
    300.   600.   150.   300.   450.   900.]
    [ -200.  -600.  -600. -1800.  -300.  -900.  -900. -2700.   200.   600.
    600.  1800.   300.   900.   900.  2700.]
    [  -75.  -150.  -150.  -300.  -175.  -350.  -350.  -700.    75.   150.
    150.   300.   175.   350.   350.   700.]
    [ -150.  -450.  -300.  -900.  -350. -1050.  -700. -2100.   150.   450.
    300.   900.   350.  1050.   700.  2100.]
    [ -150.  -300.  -450.  -900.  -350.  -700. -1050. -2100.   150.   300.
    450.   900.   350.   700.  1050.  2100.]
    [ -300.  -900.  -900. -2700.  -700. -2100. -2100. -6300.   300.   900.
    900.  2700.   700.  2100.  2100.  6300.]]
>>> print(b)
[ -198.5  -412.   -459.   -948.   -344.   -710.5  -781.  -1607.    198.5
    412.    459.    948.    344.    710.5   781.   1607. ]
jwst.wiimatch.lsq_optimizer.lsq_solve(matrix, free_term, nimages=None, tol=None)[source]

Computes least-square solution of a system of linear equations

a \cdot c = b.

This function uses Moore-Penrose pseudoinverse to solve the above system of equations.

matrix : numpy.ndarray
A 2D array containing coefficients of the system.
free_term : numpy.ndarray
A 1D array containing free terms of the system of the equations.
nimages : int, None, optional
Number of images for which the system is being solved.
tol : float, None, optional
Cutoff for small singular values for Moore-Penrose pseudoinverse. When provided, singular values smaller (in modulus) than tol * |largest_singular_value| are set to zero. When tol is None (default), cutoff value is determined based on the type of the input matrix argument.
bkg_poly_coeff : numpy.ndarray

When nimages is None, this function returns a 1D numpy.ndarray that holds the solution (polynomial coefficients) to the system.

When nimages is not None, this function returns a 2D numpy.ndarray that holds the solution (polynomial coefficients) to the system. The solution is grouped by image.

>>> import wiimatch
>>> import numpy as np
>>> im1 = np.zeros((5, 5, 4), dtype=np.float)
>>> cbg = 1.32 * np.ones_like(im1)
>>> ind = np.indices(im1.shape, dtype=np.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=np.float)
>>> a, b = wiimatch.lsq_optimizer.build_lsq_eqs([im1, im3], [mask, mask],
... [sigma, sigma], degree=(1,1,1), center=(0,0,0))
>>> wiimatch.lsq_optimizer.lsq_solve(a, b, 2)
array([[ -6.60000000e-01,  -7.50000000e-02,  -3.10000000e-01,
          3.33066907e-15,  -3.70000000e-01,   5.44009282e-15,
          7.88258347e-15,  -2.33146835e-15],
       [  6.60000000e-01,   7.50000000e-02,   3.10000000e-01,
         -4.44089210e-15,   3.70000000e-01,  -4.21884749e-15,
         -7.43849426e-15,   1.77635684e-15]])