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:
- wiimatch.lsq_optimizer.build_lsq_eqs(images, masks, sigmas, degree, center=None, image2world=None, center_cs='image', container_cls=<class 'wiimatch.containers.WMInMemoryData'>)[source]¶
build_lsq_eqs(images, masks, sigmas, degree, center=None, image2world=None, center_cs=’image’, container_cls=WMInMemoryData): Build system of linear equations whose solution would provide image intensity matching in the least squares sense.
- Parameters:
- imageslist of WMData
A list of
WMData
to 1D, 2D, etc.numpy.ndarray
data arrays whose “intensities” must be “matched”. All arrays must have identical shapes. Whenimages
is a list ofnumpy.ndarray
, the container class specified by thedefault_container
will be used to convertnumpy.ndarray
toWMData
objects. Input list may mixWMData
,numpy.ndarray
, andNone
objects.- maskslist of WMData and/or None
A list of
WMData
of the same length asimages
. Non-zero mask elements in data arrays 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 (the corresponding) input images are valid. Whenmasks
is a list ofnumpy.ndarray
, the container class specified by thedefault_container
will be used to convertnumpy.ndarray
toWMData
objects. Input list may mixWMData
,numpy.ndarray
, andNone
objects.- sigmaslist of WMData, list of None
A list of
WMData
of the 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 in all images will be assigned equal weights of 1. Whensigmas
is a list ofnumpy.ndarray
, the container class specified by thedefault_container
will be used to convertnumpy.ndarray
toWMData
objects.- degreeiterable
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.- centeriterable, None, optional
An iterable of length equal to the number of dimensions of images in
images
parameter 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.- image2worldfunction, 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
.
- Returns:
- anumpy.ndarray
A 2D
numpy.ndarray
that holds the coefficients of the linear system of equations.- bnumpy.ndarray
A 1D
numpy.ndarray
that holds the free terms of the linear system of equations.- coord_arrayslist
A list of
numpy.ndarray
coordinate arrays each ofimages[0].shape
shape.- eff_centertuple
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.
Notes
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).\]Examples
>>> from wiimatch.lsq_optimizer import build_lsq_eqs >>> from wiimatch.containers import WMInMemoryData >>> 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) >>> a, b, ca, ef, cs = build_lsq_eqs( ... [WMInMemoryData(im1), WMInMemoryData(im3)], ... [WMInMemoryData(mask), WMInMemoryData(mask)], ... [WMInMemoryData(sigma), WMInMemoryData(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. ]
- wiimatch.lsq_optimizer.pinv_solve(matrix, free_term, nimages, tol=None)[source]¶
Solves a system of linear equations
\[a \cdot c = b.\]using Moore-Penrose pseudoinverse.
- Parameters:
- matrixnumpy.ndarray
A 2D array containing coefficients of the system.
- free_termnumpy.ndarray
A 1D array containing free terms of the system of the equations.
- nimagesint
Number of images for which the system is being solved.
- tolfloat, 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. Whentol
isNone
(default), cutoff value is determined based on the type of the inputmatrix
argument.
- Returns:
- bkg_poly_coeffnumpy.ndarray
A 2D
numpy.ndarray
that holds the solution (polynomial coefficients) to the system. The solution is grouped by image.
Examples
>>> from wiimatch.lsq_optimizer import build_lsq_eqs, pinv_solve >>> from wiimatch.containers import WMInMemoryData >>> 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) >>> a, b, _, _, _ = build_lsq_eqs( ... [WMInMemoryData(im1), WMInMemoryData(im3)], ... [WMInMemoryData(mask), WMInMemoryData(mask)], ... [WMInMemoryData(sigma), WMInMemoryData(sigma)], ... degree=(1, 1, 1), center=(0, 0, 0) ... ) >>> pinv_solve(a, b, 2) array([[-6.60000000e-01, -7.50000000e-02, -3.10000000e-01, -4.44089210e-15, -3.70000000e-01, -7.66053887e-15, 3.69704267e-14, 8.37108161e-14], [ 6.60000000e-01, 7.50000000e-02, 3.10000000e-01, 3.55271368e-15, 3.70000000e-01, 4.32986980e-15, 4.88498131e-14, 7.87148124e-14]])
- wiimatch.lsq_optimizer.rlu_solve(matrix, free_term, nimages)[source]¶
Computes solution of a “reduced” system of linear equations
\[a' \cdot c' = b'.\]using LU-decomposition. If the original system contained a set of linearly-dependent equations, then the “reduced” system is formed by dropping equations and unknowns related to the first image. The unknowns corresponding to the first image initially are assumed to be 0. Upon solving the reduced system, these unknowns are recomputed so that mean correction coefficients for all images are 0. This function uses
lu_solve
andlu_factor
functions.- Parameters:
- matrixnumpy.ndarray
A 2D array containing coefficients of the system.
- free_termnumpy.ndarray
A 1D array containing free terms of the system of the equations.
- nimagesint
Number of images for which the system is being solved.
- Returns:
- bkg_poly_coeffnumpy.ndarray
A 2D
numpy.ndarray
that holds the solution (polynomial coefficients) to the system. The solution is grouped by image.
Examples
>>> from wiimatch.lsq_optimizer import build_lsq_eqs, pinv_solve >>> from wiimatch.containers import WMInMemoryData >>> 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) >>> a, b, _, _, _ = build_lsq_eqs( ... [WMInMemoryData(im1), WMInMemoryData(im3)], ... [WMInMemoryData(mask), WMInMemoryData(mask)], ... [WMInMemoryData(sigma), WMInMemoryData(sigma)], ... degree=(1, 1, 1), center=(0, 0, 0) ... ) >>> rlu_solve(a, b, 2) array([[-6.60000000e-01, -7.50000000e-02, -3.10000000e-01, -6.96331881e-16, -3.70000000e-01, -1.02318154e-15, -5.96855898e-16, 2.98427949e-16], [ 6.60000000e-01, 7.50000000e-02, 3.10000000e-01, 6.96331881e-16, 3.70000000e-01, 1.02318154e-15, 5.96855898e-16, -2.98427949e-16]])