From 6278c240dec7bd58c4454989ed3fb152e6736870 Mon Sep 17 00:00:00 2001 From: Nate Lust Date: Thu, 13 Jul 2023 12:58:12 -0400 Subject: [PATCH] Change the algorithm FitAffineWcsTask uses This commit changes FitAffineWcsTask to directly solve for affine and shift parameters by solving a matrix equation instead of creating wcs objects while iterating in a solver. This also absorbs linear shifts inside the linear transforms instead of solving for a sphere point offset in the initial wcs. --- python/lsst/meas/astrom/fitAffineWcs.py | 156 +++++++++--------------- 1 file changed, 60 insertions(+), 96 deletions(-) diff --git a/python/lsst/meas/astrom/fitAffineWcs.py b/python/lsst/meas/astrom/fitAffineWcs.py index 7d5031ad..3764c505 100644 --- a/python/lsst/meas/astrom/fitAffineWcs.py +++ b/python/lsst/meas/astrom/fitAffineWcs.py @@ -24,11 +24,11 @@ import astshim import numpy as np -from scipy.optimize import least_squares +from scipy.linalg import lstsq from lsst.afw.geom import makeSkyWcs, SkyWcs import lsst.afw.math -from lsst.geom import Point2D, degrees, arcseconds, radians +from lsst.geom import Point2D, radians import lsst.pex.config as pexConfig import lsst.pipe.base as pipeBase from lsst.utils.timer import timeMethod @@ -37,53 +37,6 @@ from .setMatchDistance import setMatchDistance -def _chiFunc(x, refPoints, srcPixels, wcsMaker): - """Function to minimize to fit the shift and rotation in the WCS. - - Parameters - ---------- - x : `numpy.ndarray` - Current fit values to test. Float values in array are: - - - ``bearingTo``: Direction to move the wcs coord in. - - ``separation``: Distance along sphere to move wcs coord in. - - ``affine0,0``: [0, 0] value of the 2x2 affine transform matrix. - - ``affine0,1``: [0, 1] value of the 2x2 affine transform matrix. - - ``affine1,0``: [1, 0] value of the 2x2 affine transform matrix. - - ``affine1,1``: [1, 1] value of the 2x2 affine transform matrix. - refPoints : `list` of `lsst.afw.geom.SpherePoint` - Reference object on Sky locations. - srcPixels : `list` of `lsst.geom.Point2D` - Source object positions on the pixels. - wcsMaker : `TransformedSkyWcsMaker` - Container class for producing the updated Wcs. - - Returns - ------- - outputSeparations : `list` of `float` - Separation between predicted source location and reference location in - radians. - """ - wcs = wcsMaker.makeWcs(x[:2], x[2:].reshape((2, 2))) - - outputSeparations = [] - # Fit both sky to pixel and pixel to sky to avoid any non-invertible - # affine matrices. - for ref, src in zip(refPoints, srcPixels): - skySep = ref.getTangentPlaneOffset(wcs.pixelToSky(src)) - outputSeparations.append(skySep[0].asArcseconds()) - outputSeparations.append(skySep[1].asArcseconds()) - xySep = src - wcs.skyToPixel(ref) - # Convert the pixel separations to units, arcseconds to match units - # of sky separation. - outputSeparations.append( - xySep[0] * wcs.getPixelScale(src).asArcseconds()) - outputSeparations.append( - xySep[1] * wcs.getPixelScale(src).asArcseconds()) - - return outputSeparations - - # Keeping this around for now in case any of the fit parameters need to be # configurable. Likely the maximum allowed shift magnitude (parameter 2 in the # fit.) @@ -163,52 +116,62 @@ def fitWcs(self, # appends the new transform. wcsMaker = TransformedSkyWcsMaker(initWcs) - refPoints = [] - srcPixels = [] - offsetLong = 0 - offsetLat = 0 - # Grab reference coordinates and source centroids. Compute the average - # direction and separation between the reference and the sources. - for match in matches: + # Grab the initial transformations going back from sky coordinates + # and forward from pixel coordinates. + back = wcsMaker.frameDict.getMapping(wcsMaker.frameMax, wcsMaker.frameMax-1) + forward = wcsMaker.lastMapBeforeSky + + # Create containers for the data going into the Affine fit. This will + # be done by approximating the solution to Ax=b where x will be the + # affine parameters and a linear shift. The approximate solution is + # calculated using a least squares minimization of the Ax=b equation. + # + # This is looking to find the affine transform of the following form: + # [x', y'] = [[a, b], [c, d]] [x, y] + [s, t] + # + # where a,b,c,d are the parameters of the affine transform, and s,t + # are linear shift parameters. + # + # To solve for these unknown parameters the unknown matrix x in the + # equation Ax=b will be of the form: + # [a, b, c, d, s, t]. + # + # This implies that each constraining point will correspond to two rows + # in the A matrix in the following form: + # [x_i, y_i, 0, 0, 1, 0] + # [0, 0, x_i, y_i, 0, 1]. + # + # The corresponding output points in the b vector will have the form: + # [x'_i, y'_i, x'_(i+i), y'_(i+1)....] + A = np.zeros((len(matches)*2, 6), dtype=float) + b = np.empty(len(matches)*2, dtype=float) + + # Constant terms related to the shift in x and and y parameters. + A[::2, 4] = 1 + A[1::2, 5] = 1 + + # loop over each of the matches and populate the matrices. + for i, match in enumerate(matches): refCoord = match.first.getCoord() - refPoints.append(refCoord) + b[i*2:i*2+2] = back.applyForward(refCoord) + srcCentroid = match.second.getCentroid() - srcPixels.append(srcCentroid) - srcCoord = initWcs.pixelToSky(srcCentroid) - deltaLong, deltaLat = srcCoord.getTangentPlaneOffset(refCoord) - offsetLong += deltaLong.asArcseconds() - offsetLat += deltaLat.asArcseconds() - offsetLong /= len(srcPixels) - offsetLat /= len(srcPixels) - offsetDist = np.sqrt(offsetLong ** 2 + offsetLat ** 2) - if offsetDist > 0.: - offsetDir = np.degrees(np.arccos(offsetLong / offsetDist)) - else: - offsetDir = 0. - offsetDir *= np.sign(offsetLat) - self.log.debug("Initial shift guess: Direction: %.3f, Dist %.3f...", - offsetDir, offsetDist) - - # Best performing fitter in scipy tried so far (vs. default settings in - # minimize). Exits early because of the xTol value which cannot be - # disabled in scipy1.2.1. Matrix starting values are non-zero as this - # results in better fit off-diagonal terms. - fit = least_squares( - _chiFunc, - x0=[offsetDir, offsetDist, 1., 1e-8, 1e-8, 1.], - args=(refPoints, srcPixels, wcsMaker), - method='dogbox', - bounds=[[-360, -np.inf, -np.inf, -np.inf, -np.inf, -np.inf], - [360, np.inf, np.inf, np.inf, np.inf, np.inf]], - ftol=2.3e-16, - gtol=2.31e-16, - xtol=2.3e-16) - self.log.debug("Best fit: Direction: %.3f, Dist: %.3f, " + val = forward.applyForward(srcCentroid) + A[i*2, :2] = val + A[i*2+1, 2:4] = val + + # solve for the affine and shift parameters + # The lapack_driver parameter is set to the quickest routine tested for + # this application at the time of writing. + fit = lstsq(A, b, lapack_driver='gelsy')[0] + + self.log.debug("Linear shift in x: %.3f, y: %.3f, " "Affine matrix: [[%.6f, %.6f], [%.6f, %.6f]]...", - fit.x[0], fit.x[1], - fit.x[2], fit.x[3], fit.x[4], fit.x[5]) + fit[4], fit[5], + fit[0], fit[1], fit[2], fit[3]) - wcs = wcsMaker.makeWcs(fit.x[:2], fit.x[2:].reshape((2, 2))) + # create the final wcs + wcs = wcsMaker.makeWcs(fit[4:], fit[:4].reshape((2, 2))) # Copied from other fit*WcsTasks. if refCat is not None: @@ -286,7 +249,7 @@ def __init__(self, inputSkyWcs): self.origin = inputSkyWcs.getSkyOrigin() - def makeWcs(self, crvalOffset, affMatrix): + def makeWcs(self, linearShift, affMatrix): """Apply a shift and affine transform to the WCS internal to this class. @@ -294,8 +257,9 @@ def makeWcs(self, crvalOffset, affMatrix): Parameters ---------- - crval_shift : `numpy.ndarray`, (2,) - Shift in radians to apply to the Wcs origin/crvals. + linearShift : `numpy.ndarray`, (2,) + A linear shift to apply at the same time as applying the affine + matrix transform. aff_matrix : 'numpy.ndarray', (3, 3) Affine matrix to apply to the mapping/transform to add to the WCS. @@ -311,14 +275,14 @@ def makeWcs(self, crvalOffset, affMatrix): # map and gives us our final shift position. iwcsToSkyWcs = makeSkyWcs( Point2D(0., 0.), - self.origin.offset(crvalOffset[0] * degrees, - crvalOffset[1] * arcseconds), + self.origin, np.array([[1., 0.], [0., 1.]])) iwcToSkyMap = iwcsToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY") # Append a simple affine Matrix transform to the current to the # second to last frame mapping. e.g. the one just before IWC to SKY. newMapping = self.lastMapBeforeSky.then(astshim.MatrixMap(affMatrix)) + newMapping = newMapping.then(astshim.ShiftMap(linearShift)) # Create a new frame dict starting from the input_sky_wcs's first # frame. Append the correct mapping created above and our new on