Source code for lsst.afw.geom.testUtils

#
# LSST Data Management System
# Copyright 2016 LSST Corporation.
#
# This product includes software developed by the
# LSST Project (http://www.lsst.org/).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the LSST License Statement and
# the GNU General Public License along with this program.  If not,
# see <http://www.lsstcorp.org/LegalNotices/>.
#

__all__ = ["BoxGrid", "makeSipIwcToPixel", "makeSipPixelToIwc"]

import itertools
import math
import os
import pickle

import astshim as ast
import numpy as np
from numpy.testing import assert_allclose, assert_array_equal
from astshim.test import makeForwardPolyMap, makeTwoWayPolyMap
from lsst.afw.geom.wcsUtils import getCdMatrixFromMetadata

import lsst.geom
import lsst.afw.geom as afwGeom
from lsst.pex.exceptions import InvalidParameterError
import lsst.utils
import lsst.utils.tests


[docs]class BoxGrid: """Divide a box into nx by ny sub-boxes that tile the region The sub-boxes will be of the same type as `box` and will exactly tile `box`; they will also all be the same size, to the extent possible (some variation is inevitable for integer boxes that cannot be evenly divided. Parameters ---------- box : `lsst.geom.Box2I` or `lsst.geom.Box2D` the box to subdivide; the boxes in the grid will be of the same type numColRow : pair of `int` number of columns and rows """ def __init__(self, box, numColRow): if len(numColRow) != 2: raise RuntimeError( "numColRow=%r; must be a sequence of two integers" % (numColRow,)) self._numColRow = tuple(int(val) for val in numColRow) if isinstance(box, lsst.geom.Box2I): stopDelta = 1 elif isinstance(box, lsst.geom.Box2D): stopDelta = 0 else: raise RuntimeError("Unknown class %s of box %s" % (type(box), box)) self.boxClass = type(box) self.stopDelta = stopDelta minPoint = box.getMin() self.pointClass = type(minPoint) dtype = np.array(minPoint).dtype self._divList = [np.linspace(start=box.getMin()[i], stop=box.getMax()[i] + self.stopDelta, num=self._numColRow[i] + 1, endpoint=True, dtype=dtype) for i in range(2)] @property def numColRow(self): return self._numColRow def __getitem__(self, indXY): """Return the box at the specified x,y index Parameters ---------- indXY : pair of `ints` the x,y index to return Returns ------- subBox : `lsst.geom.Box2I` or `lsst.geom.Box2D` """ beg = self.pointClass(*[self._divList[i][indXY[i]] for i in range(2)]) end = self.pointClass( *[self._divList[i][indXY[i] + 1] - self.stopDelta for i in range(2)]) return self.boxClass(beg, end) def __len__(self): return self.shape[0]*self.shape[1] def __iter__(self): """Return an iterator over all boxes, where column varies most quickly """ for row in range(self.numColRow[1]): for col in range(self.numColRow[0]): yield self[col, row]
class FrameSetInfo: """Information about a FrameSet Parameters ---------- frameSet : `ast.FrameSet` The FrameSet about which you want information Notes ----- **Fields** baseInd : `int` Index of base frame currInd : `int` Index of current frame isBaseSkyFrame : `bool` Is the base frame an `ast.SkyFrame`? isCurrSkyFrame : `bool` Is the current frame an `ast.SkyFrame`? """ def __init__(self, frameSet): self.baseInd = frameSet.base self.currInd = frameSet.current self.isBaseSkyFrame = frameSet.getFrame(self.baseInd).className == "SkyFrame" self.isCurrSkyFrame = frameSet.getFrame(self.currInd).className == "SkyFrame" def makeSipPolyMapCoeffs(metadata, name): """Return a list of ast.PolyMap coefficients for the specified SIP matrix The returned list of coefficients for an ast.PolyMap that computes the following function: f(dxy) = dxy + sipPolynomial(dxy)) where dxy = pixelPosition - pixelOrigin and sipPolynomial is a polynomial with terms `<name>n_m for x^n y^m` (e.g. A2_0 is the coefficient for x^2 y^0) Parameters ---------- metadata : lsst.daf.base.PropertySet FITS metadata describing a WCS with the specified SIP coefficients name : str The desired SIP terms: one of A, B, AP, BP Returns ------- list A list of coefficients for an ast.PolyMap that computes the specified SIP polynomial, including a term for out = in Note ---- This is an internal function for use by makeSipIwcToPixel and makeSipPixelToIwc """ outAxisDict = dict(A=1, B=2, AP=1, BP=2) outAxis = outAxisDict.get(name) if outAxis is None: raise RuntimeError("%s not a supported SIP name" % (name,)) width = metadata.getAsInt("%s_ORDER" % (name,)) + 1 found = False # start with a term for out = in coeffs = [] if outAxis == 1: coeffs.append([1.0, outAxis, 1, 0]) else: coeffs.append([1.0, outAxis, 0, 1]) # add SIP distortion terms for xPower in range(width): for yPower in range(width): coeffName = "%s_%s_%s" % (name, xPower, yPower) if not metadata.exists(coeffName): continue found = True coeff = metadata.getAsDouble(coeffName) coeffs.append([coeff, outAxis, xPower, yPower]) if not found: raise RuntimeError("No %s coefficients found" % (name,)) return coeffs
[docs]def makeSipIwcToPixel(metadata): """Make an IWC to pixel transform with SIP distortion from FITS-WCS metadata This function is primarily intended for unit tests. IWC is intermediate world coordinates, as described in the FITS papers. Parameters ---------- metadata : lsst.daf.base.PropertySet FITS metadata describing a WCS with inverse SIP coefficients Returns ------- lsst.afw.geom.TransformPoint2ToPoint2 Transform from IWC position to pixel position (zero-based) in the forward direction. The inverse direction is not defined. Notes ----- The inverse SIP terms APn_m, BPn_m are polynomial coefficients x^n y^m for computing transformed x, y respectively. If we call the resulting polynomial inverseSipPolynomial, the returned transformation is: pixelPosition = pixel origin + uv + inverseSipPolynomial(uv) where uv = inverseCdMatrix * iwcPosition """ crpix = (metadata.getScalar("CRPIX1") - 1, metadata.getScalar("CRPIX2") - 1) pixelRelativeToAbsoluteMap = ast.ShiftMap(crpix) cdMatrix = getCdMatrixFromMetadata(metadata) cdMatrixMap = ast.MatrixMap(cdMatrix.copy()) coeffList = makeSipPolyMapCoeffs(metadata, "AP") + makeSipPolyMapCoeffs(metadata, "BP") coeffArr = np.array(coeffList, dtype=float) sipPolyMap = ast.PolyMap(coeffArr, 2, "IterInverse=0") iwcToPixelMap = cdMatrixMap.getInverse().then(sipPolyMap).then(pixelRelativeToAbsoluteMap) return afwGeom.TransformPoint2ToPoint2(iwcToPixelMap)
[docs]def makeSipPixelToIwc(metadata): """Make a pixel to IWC transform with SIP distortion from FITS-WCS metadata This function is primarily intended for unit tests. IWC is intermediate world coordinates, as described in the FITS papers. Parameters ---------- metadata : lsst.daf.base.PropertySet FITS metadata describing a WCS with forward SIP coefficients Returns ------- lsst.afw.geom.TransformPoint2ToPoint2 Transform from pixel position (zero-based) to IWC position in the forward direction. The inverse direction is not defined. Notes ----- The forward SIP terms An_m, Bn_m are polynomial coefficients x^n y^m for computing transformed x, y respectively. If we call the resulting polynomial sipPolynomial, the returned transformation is: iwcPosition = cdMatrix * (dxy + sipPolynomial(dxy)) where dxy = pixelPosition - pixelOrigin """ crpix = (metadata.getScalar("CRPIX1") - 1, metadata.getScalar("CRPIX2") - 1) pixelAbsoluteToRelativeMap = ast.ShiftMap(crpix).getInverse() cdMatrix = getCdMatrixFromMetadata(metadata) cdMatrixMap = ast.MatrixMap(cdMatrix.copy()) coeffList = makeSipPolyMapCoeffs(metadata, "A") + makeSipPolyMapCoeffs(metadata, "B") coeffArr = np.array(coeffList, dtype=float) sipPolyMap = ast.PolyMap(coeffArr, 2, "IterInverse=0") pixelToIwcMap = pixelAbsoluteToRelativeMap.then(sipPolyMap).then(cdMatrixMap) return afwGeom.TransformPoint2ToPoint2(pixelToIwcMap)
class PermutedFrameSet: """A FrameSet with base or current frame possibly permuted, with associated information Only two-axis frames will be permuted. Parameters ---------- frameSet : `ast.FrameSet` The FrameSet you wish to permute. A deep copy is made. permuteBase : `bool` Permute the base frame's axes? permuteCurr : `bool` Permute the current frame's axes? Raises ------ RuntimeError If you try to permute a frame that does not have 2 axes Notes ----- **Fields** frameSet : `ast.FrameSet` The FrameSet that may be permuted. A local copy is made. isBaseSkyFrame : `bool` Is the base frame an `ast.SkyFrame`? isCurrSkyFrame : `bool` Is the current frame an `ast.SkyFrame`? isBasePermuted : `bool` Are the base frame axes permuted? isCurrPermuted : `bool` Are the current frame axes permuted? """ def __init__(self, frameSet, permuteBase, permuteCurr): self.frameSet = frameSet.copy() fsInfo = FrameSetInfo(self.frameSet) self.isBaseSkyFrame = fsInfo.isBaseSkyFrame self.isCurrSkyFrame = fsInfo.isCurrSkyFrame if permuteBase: baseNAxes = self.frameSet.getFrame(fsInfo.baseInd).nAxes if baseNAxes != 2: raise RuntimeError("Base frame has {} axes; 2 required to permute".format(baseNAxes)) self.frameSet.current = fsInfo.baseInd self.frameSet.permAxes([2, 1]) self.frameSet.current = fsInfo.currInd if permuteCurr: currNAxes = self.frameSet.getFrame(fsInfo.currInd).nAxes if currNAxes != 2: raise RuntimeError("Current frame has {} axes; 2 required to permute".format(currNAxes)) assert self.frameSet.getFrame(fsInfo.currInd).nAxes == 2 self.frameSet.permAxes([2, 1]) self.isBasePermuted = permuteBase self.isCurrPermuted = permuteCurr class TransformTestBaseClass(lsst.utils.tests.TestCase): """Base class for unit tests of Transform<X>To<Y> Subclasses must call `TransformTestBaseClass.setUp(self)` if they provide their own version. If a package other than afw uses this class then it must override the `getTestDir` method to avoid writing into afw's test directory. """ def getTestDir(self): """Return a directory where temporary test files can be written The default implementation returns the test directory of the `afw` package. If this class is used by a test in a package other than `afw` then the subclass must override this method. """ return os.path.join(lsst.utils.getPackageDir("afw"), "tests") def setUp(self): """Set up a test Subclasses should call this method if they override setUp. """ # tell unittest to use the msg argument of asserts as a supplement # to the error message, rather than as the whole error message self.longMessage = True # list of endpoint class name prefixes; the full name is prefix + "Endpoint" self.endpointPrefixes = ("Generic", "Point2", "SpherePoint") # GoodNAxes is dict of endpoint class name prefix: # tuple containing 0 or more valid numbers of axes self.goodNAxes = { "Generic": (1, 2, 3, 4), # all numbers of axes are valid for GenericEndpoint "Point2": (2,), "SpherePoint": (2,), } # BadAxes is dict of endpoint class name prefix: # tuple containing 0 or more invalid numbers of axes self.badNAxes = { "Generic": (), # all numbers of axes are valid for GenericEndpoint "Point2": (1, 3, 4), "SpherePoint": (1, 3, 4), } # Dict of frame index: identity name for frames created by makeFrameSet self.frameIdentDict = { 1: "baseFrame", 2: "frame2", 3: "frame3", 4: "currFrame", } @staticmethod def makeRawArrayData(nPoints, nAxes, delta=0.123): """Make an array of generic point data The data will be suitable for spherical points Parameters ---------- nPoints : `int` Number of points in the array nAxes : `int` Number of axes in the point Returns ------- np.array of floats with shape (nAxes, nPoints) The values are as follows; if nAxes != 2: The first point has values `[0, delta, 2*delta, ..., (nAxes-1)*delta]` The Nth point has those values + N if nAxes == 2 then the data is scaled so that the max value of axis 1 is a bit less than pi/2 """ delta = 0.123 # oneAxis = [0, 1, 2, ...nPoints-1] oneAxis = np.arange(nPoints, dtype=float) # [0, 1, 2...] # rawData = [oneAxis, oneAxis + delta, oneAxis + 2 delta, ...] rawData = np.array([j * delta + oneAxis for j in range(nAxes)], dtype=float) if nAxes == 2: # scale rawData so that max value of 2nd axis is a bit less than pi/2, # thus making the data safe for SpherePoint maxLatitude = np.max(rawData[1]) rawData *= math.pi * 0.4999 / maxLatitude return rawData @staticmethod def makeRawPointData(nAxes, delta=0.123): """Make one generic point Parameters ---------- nAxes : `int` Number of axes in the point delta : `float` Increment between axis values Returns ------- A list of `nAxes` floats with values `[0, delta, ..., (nAxes-1)*delta] """ return [i*delta for i in range(nAxes)] @staticmethod def makeEndpoint(name, nAxes=None): """Make an endpoint Parameters ---------- name : `str` Endpoint class name prefix; the full class name is name + "Endpoint" nAxes : `int` or `None`, optional number of axes; an int is required if `name` == "Generic"; otherwise ignored Returns ------- subclass of `lsst.afw.geom.BaseEndpoint` The constructed endpoint Raises ------ TypeError If `name` == "Generic" and `nAxes` is None or <= 0 """ EndpointClassName = name + "Endpoint" EndpointClass = getattr(afwGeom, EndpointClassName) if name == "Generic": if nAxes is None: raise TypeError("nAxes must be an integer for GenericEndpoint") return EndpointClass(nAxes) return EndpointClass() @classmethod def makeGoodFrame(cls, name, nAxes=None): """Return the appropriate frame for the given name and nAxes Parameters ---------- name : `str` Endpoint class name prefix; the full class name is name + "Endpoint" nAxes : `int` or `None`, optional number of axes; an int is required if `name` == "Generic"; otherwise ignored Returns ------- `ast.Frame` The constructed frame Raises ------ TypeError If `name` == "Generic" and `nAxes` is `None` or <= 0 """ return cls.makeEndpoint(name, nAxes).makeFrame() @staticmethod def makeBadFrames(name): """Return a list of 0 or more frames that are not a valid match for the named endpoint Parameters ---------- name : `str` Endpoint class name prefix; the full class name is name + "Endpoint" Returns ------- Collection of `ast.Frame` A collection of 0 or more frames """ return { "Generic": [], "Point2": [ ast.SkyFrame(), ast.Frame(1), ast.Frame(3), ], "SpherePoint": [ ast.Frame(1), ast.Frame(2), ast.Frame(3), ], }[name] def makeFrameSet(self, baseFrame, currFrame): """Make a FrameSet The FrameSet will contain 4 frames and three transforms connecting them. The idenity of each frame is provided by self.frameIdentDict Frame Index Mapping from this frame to the next `baseFrame` 1 `ast.UnitMap(nIn)` Frame(nIn) 2 `polyMap` Frame(nOut) 3 `ast.UnitMap(nOut)` `currFrame` 4 where: - `nIn` = `baseFrame.nAxes` - `nOut` = `currFrame.nAxes` - `polyMap` = `makeTwoWayPolyMap(nIn, nOut)` Returns ------ `ast.FrameSet` The FrameSet as described above Parameters ---------- baseFrame : `ast.Frame` base frame currFrame : `ast.Frame` current frame """ nIn = baseFrame.nAxes nOut = currFrame.nAxes polyMap = makeTwoWayPolyMap(nIn, nOut) # The only way to set the Ident of a frame in a FrameSet is to set it in advance, # and I don't want to modify the inputs, so replace the input frames with copies baseFrame = baseFrame.copy() baseFrame.ident = self.frameIdentDict[1] currFrame = currFrame.copy() currFrame.ident = self.frameIdentDict[4] frameSet = ast.FrameSet(baseFrame) frame2 = ast.Frame(nIn) frame2.ident = self.frameIdentDict[2] frameSet.addFrame(ast.FrameSet.CURRENT, ast.UnitMap(nIn), frame2) frame3 = ast.Frame(nOut) frame3.ident = self.frameIdentDict[3] frameSet.addFrame(ast.FrameSet.CURRENT, polyMap, frame3) frameSet.addFrame(ast.FrameSet.CURRENT, ast.UnitMap(nOut), currFrame) return frameSet @staticmethod def permuteFrameSetIter(frameSet): """Iterator over 0 or more frameSets with SkyFrames axes permuted Only base and current SkyFrames are permuted. If neither the base nor the current frame is a SkyFrame then no frames are returned. Returns ------- iterator over `PermutedFrameSet` """ fsInfo = FrameSetInfo(frameSet) if not (fsInfo.isBaseSkyFrame or fsInfo.isCurrSkyFrame): return permuteBaseList = [False, True] if fsInfo.isBaseSkyFrame else [False] permuteCurrList = [False, True] if fsInfo.isCurrSkyFrame else [False] for permuteBase in permuteBaseList: for permuteCurr in permuteCurrList: yield PermutedFrameSet(frameSet, permuteBase, permuteCurr) @staticmethod def makeJacobian(nIn, nOut, inPoint): """Make a Jacobian matrix for the equation described by `makeTwoWayPolyMap`. Parameters ---------- nIn, nOut : `int` the dimensions of the input and output data; see makeTwoWayPolyMap inPoint : `numpy.ndarray` an array of size `nIn` representing the point at which the Jacobian is measured Returns ------- J : `numpy.ndarray` an `nOut` x `nIn` array of first derivatives """ basePolyMapCoeff = 0.001 # see makeTwoWayPolyMap baseCoeff = 2.0 * basePolyMapCoeff coeffs = np.empty((nOut, nIn)) for iOut in range(nOut): coeffOffset = baseCoeff * iOut for iIn in range(nIn): coeffs[iOut, iIn] = baseCoeff * (iIn + 1) + coeffOffset coeffs[iOut, iIn] *= inPoint[iIn] assert coeffs.ndim == 2 # Avoid spurious errors when comparing to a simplified array assert coeffs.shape == (nOut, nIn) return coeffs def checkTransformation(self, transform, mapping, msg=""): """Check applyForward and applyInverse for a transform Parameters ---------- transform : `lsst.afw.geom.Transform` The transform to check mapping : `ast.Mapping` The mapping the transform should use. This mapping must contain valid forward or inverse transformations, but they need not match if both present. Hence the mappings returned by make*PolyMap are acceptable. msg : `str` Error message suffix describing test parameters """ fromEndpoint = transform.fromEndpoint toEndpoint = transform.toEndpoint mappingFromTransform = transform.getMapping() nIn = mapping.nIn nOut = mapping.nOut self.assertEqual(nIn, fromEndpoint.nAxes, msg=msg) self.assertEqual(nOut, toEndpoint.nAxes, msg=msg) # forward transformation of one point rawInPoint = self.makeRawPointData(nIn) inPoint = fromEndpoint.pointFromData(rawInPoint) # forward transformation of an array of points nPoints = 7 # arbitrary rawInArray = self.makeRawArrayData(nPoints, nIn) inArray = fromEndpoint.arrayFromData(rawInArray) if mapping.hasForward: self.assertTrue(transform.hasForward) outPoint = transform.applyForward(inPoint) rawOutPoint = toEndpoint.dataFromPoint(outPoint) assert_allclose(rawOutPoint, mapping.applyForward(rawInPoint), err_msg=msg) assert_allclose(rawOutPoint, mappingFromTransform.applyForward(rawInPoint), err_msg=msg) outArray = transform.applyForward(inArray) rawOutArray = toEndpoint.dataFromArray(outArray) self.assertFloatsAlmostEqual(rawOutArray, mapping.applyForward(rawInArray), msg=msg) self.assertFloatsAlmostEqual(rawOutArray, mappingFromTransform.applyForward(rawInArray), msg=msg) else: # Need outPoint, but don't need it to be consistent with inPoint rawOutPoint = self.makeRawPointData(nOut) outPoint = toEndpoint.pointFromData(rawOutPoint) rawOutArray = self.makeRawArrayData(nPoints, nOut) outArray = toEndpoint.arrayFromData(rawOutArray) self.assertFalse(transform.hasForward) if mapping.hasInverse: self.assertTrue(transform.hasInverse) # inverse transformation of one point; # remember that the inverse need not give the original values # (see the description of the `mapping` parameter) inversePoint = transform.applyInverse(outPoint) rawInversePoint = fromEndpoint.dataFromPoint(inversePoint) assert_allclose(rawInversePoint, mapping.applyInverse(rawOutPoint), err_msg=msg) assert_allclose(rawInversePoint, mappingFromTransform.applyInverse(rawOutPoint), err_msg=msg) # inverse transformation of an array of points; # remember that the inverse will not give the original values # (see the description of the `mapping` parameter) inverseArray = transform.applyInverse(outArray) rawInverseArray = fromEndpoint.dataFromArray(inverseArray) self.assertFloatsAlmostEqual(rawInverseArray, mapping.applyInverse(rawOutArray), msg=msg) self.assertFloatsAlmostEqual(rawInverseArray, mappingFromTransform.applyInverse(rawOutArray), msg=msg) else: self.assertFalse(transform.hasInverse) def checkInverseTransformation(self, forward, inverse, msg=""): """Check that two Transforms are each others' inverses. Parameters ---------- forward : `lsst.afw.geom.Transform` the reference Transform to test inverse : `lsst.afw.geom.Transform` the transform that should be the inverse of `forward` msg : `str` error message suffix describing test parameters """ fromEndpoint = forward.fromEndpoint toEndpoint = forward.toEndpoint forwardMapping = forward.getMapping() inverseMapping = inverse.getMapping() # properties self.assertEqual(forward.fromEndpoint, inverse.toEndpoint, msg=msg) self.assertEqual(forward.toEndpoint, inverse.fromEndpoint, msg=msg) self.assertEqual(forward.hasForward, inverse.hasInverse, msg=msg) self.assertEqual(forward.hasInverse, inverse.hasForward, msg=msg) # transformations of one point # we don't care about whether the transformation itself is correct # (see checkTransformation), so inPoint/outPoint need not be related rawInPoint = self.makeRawPointData(fromEndpoint.nAxes) inPoint = fromEndpoint.pointFromData(rawInPoint) rawOutPoint = self.makeRawPointData(toEndpoint.nAxes) outPoint = toEndpoint.pointFromData(rawOutPoint) # transformations of arrays of points nPoints = 7 # arbitrary rawInArray = self.makeRawArrayData(nPoints, fromEndpoint.nAxes) inArray = fromEndpoint.arrayFromData(rawInArray) rawOutArray = self.makeRawArrayData(nPoints, toEndpoint.nAxes) outArray = toEndpoint.arrayFromData(rawOutArray) if forward.hasForward: self.assertEqual(forward.applyForward(inPoint), inverse.applyInverse(inPoint), msg=msg) self.assertEqual(forwardMapping.applyForward(rawInPoint), inverseMapping.applyInverse(rawInPoint), msg=msg) # Assertions must work with both lists and numpy arrays assert_array_equal(forward.applyForward(inArray), inverse.applyInverse(inArray), err_msg=msg) assert_array_equal(forwardMapping.applyForward(rawInArray), inverseMapping.applyInverse(rawInArray), err_msg=msg) if forward.hasInverse: self.assertEqual(forward.applyInverse(outPoint), inverse.applyForward(outPoint), msg=msg) self.assertEqual(forwardMapping.applyInverse(rawOutPoint), inverseMapping.applyForward(rawOutPoint), msg=msg) assert_array_equal(forward.applyInverse(outArray), inverse.applyForward(outArray), err_msg=msg) assert_array_equal(forwardMapping.applyInverse(rawOutArray), inverseMapping.applyForward(rawOutArray), err_msg=msg) def checkTransformFromMapping(self, fromName, toName): """Check Transform_<fromName>_<toName> using the Mapping constructor Parameters ---------- fromName, toName : `str` Endpoint name prefix for "from" and "to" endpoints, respectively, e.g. "Point2" for `lsst.afw.geom.Point2Endpoint` fromAxes, toAxes : `int` number of axes in fromFrame and toFrame, respectively """ transformClassName = "Transform{}To{}".format(fromName, toName) TransformClass = getattr(afwGeom, transformClassName) baseMsg = "TransformClass={}".format(TransformClass.__name__) # check valid numbers of inputs and outputs for nIn, nOut in itertools.product(self.goodNAxes[fromName], self.goodNAxes[toName]): msg = "{}, nIn={}, nOut={}".format(baseMsg, nIn, nOut) polyMap = makeTwoWayPolyMap(nIn, nOut) transform = TransformClass(polyMap) # desired output from `str(transform)` desStr = "{}[{}->{}]".format(transformClassName, nIn, nOut) self.assertEqual("{}".format(transform), desStr) self.assertEqual(repr(transform), "lsst.afw.geom." + desStr) self.checkTransformation(transform, polyMap, msg=msg) # Forward transform but no inverse polyMap = makeForwardPolyMap(nIn, nOut) transform = TransformClass(polyMap) self.checkTransformation(transform, polyMap, msg=msg) # Inverse transform but no forward polyMap = makeForwardPolyMap(nOut, nIn).getInverse() transform = TransformClass(polyMap) self.checkTransformation(transform, polyMap, msg=msg) # check invalid # of output against valid # of inputs for nIn, badNOut in itertools.product(self.goodNAxes[fromName], self.badNAxes[toName]): badPolyMap = makeTwoWayPolyMap(nIn, badNOut) msg = "{}, nIn={}, badNOut={}".format(baseMsg, nIn, badNOut) with self.assertRaises(InvalidParameterError, msg=msg): TransformClass(badPolyMap) # check invalid # of inputs against valid and invalid # of outputs for badNIn, nOut in itertools.product(self.badNAxes[fromName], self.goodNAxes[toName] + self.badNAxes[toName]): badPolyMap = makeTwoWayPolyMap(badNIn, nOut) msg = "{}, badNIn={}, nOut={}".format(baseMsg, nIn, nOut) with self.assertRaises(InvalidParameterError, msg=msg): TransformClass(badPolyMap) def checkTransformFromFrameSet(self, fromName, toName): """Check Transform_<fromName>_<toName> using the FrameSet constructor Parameters ---------- fromName, toName : `str` Endpoint name prefix for "from" and "to" endpoints, respectively, e.g. "Point2" for `lsst.afw.geom.Point2Endpoint` """ transformClassName = "Transform{}To{}".format(fromName, toName) TransformClass = getattr(afwGeom, transformClassName) baseMsg = "TransformClass={}".format(TransformClass.__name__) for nIn, nOut in itertools.product(self.goodNAxes[fromName], self.goodNAxes[toName]): msg = "{}, nIn={}, nOut={}".format(baseMsg, nIn, nOut) baseFrame = self.makeGoodFrame(fromName, nIn) currFrame = self.makeGoodFrame(toName, nOut) frameSet = self.makeFrameSet(baseFrame, currFrame) self.assertEqual(frameSet.nFrame, 4) # construct 0 or more frame sets that are invalid for this transform class for badBaseFrame in self.makeBadFrames(fromName): badFrameSet = self.makeFrameSet(badBaseFrame, currFrame) with self.assertRaises(InvalidParameterError): TransformClass(badFrameSet) for badCurrFrame in self.makeBadFrames(toName): reallyBadFrameSet = self.makeFrameSet(badBaseFrame, badCurrFrame) with self.assertRaises(InvalidParameterError): TransformClass(reallyBadFrameSet) for badCurrFrame in self.makeBadFrames(toName): badFrameSet = self.makeFrameSet(baseFrame, badCurrFrame) with self.assertRaises(InvalidParameterError): TransformClass(badFrameSet) transform = TransformClass(frameSet) desStr = "{}[{}->{}]".format(transformClassName, nIn, nOut) self.assertEqual("{}".format(transform), desStr) self.assertEqual(repr(transform), "lsst.afw.geom." + desStr) self.checkPersistence(transform) mappingFromTransform = transform.getMapping() transformCopy = TransformClass(mappingFromTransform) self.assertEqual(type(transform), type(transformCopy)) self.assertEqual(transform.getMapping(), mappingFromTransform) polyMap = makeTwoWayPolyMap(nIn, nOut) self.checkTransformation(transform, mapping=polyMap, msg=msg) # If the base and/or current frame of frameSet is a SkyFrame, # try permuting that frame (in place, so the connected mappings are # correctly updated). The Transform constructor should undo the permutation, # (via SpherePointEndpoint.normalizeFrame) in its internal copy of frameSet, # forcing the axes of the SkyFrame into standard (longitude, latitude) order for permutedFS in self.permuteFrameSetIter(frameSet): if permutedFS.isBaseSkyFrame: baseFrame = permutedFS.frameSet.getFrame(ast.FrameSet.BASE) # desired base longitude axis desBaseLonAxis = 2 if permutedFS.isBasePermuted else 1 self.assertEqual(baseFrame.lonAxis, desBaseLonAxis) if permutedFS.isCurrSkyFrame: currFrame = permutedFS.frameSet.getFrame(ast.FrameSet.CURRENT) # desired current base longitude axis desCurrLonAxis = 2 if permutedFS.isCurrPermuted else 1 self.assertEqual(currFrame.lonAxis, desCurrLonAxis) permTransform = TransformClass(permutedFS.frameSet) self.checkTransformation(permTransform, mapping=polyMap, msg=msg) def checkGetInverse(self, fromName, toName): """Test Transform<fromName>To<toName>.getInverse Parameters ---------- fromName, toName : `str` Endpoint name prefix for "from" and "to" endpoints, respectively, e.g. "Point2" for `lsst.afw.geom.Point2Endpoint` """ transformClassName = "Transform{}To{}".format(fromName, toName) TransformClass = getattr(afwGeom, transformClassName) baseMsg = "TransformClass={}".format(TransformClass.__name__) for nIn, nOut in itertools.product(self.goodNAxes[fromName], self.goodNAxes[toName]): msg = "{}, nIn={}, nOut={}".format(baseMsg, nIn, nOut) self.checkInverseMapping( TransformClass, makeTwoWayPolyMap(nIn, nOut), "{}, Map={}".format(msg, "TwoWay")) self.checkInverseMapping( TransformClass, makeForwardPolyMap(nIn, nOut), "{}, Map={}".format(msg, "Forward")) self.checkInverseMapping( TransformClass, makeForwardPolyMap(nOut, nIn).getInverse(), "{}, Map={}".format(msg, "Inverse")) def checkInverseMapping(self, TransformClass, mapping, msg): """Test Transform<fromName>To<toName>.getInverse for a specific mapping. Parameters ---------- TransformClass : `type` The class of transform to test, such as TransformPoint2ToPoint2 mapping : `ast.Mapping` The mapping to use for the transform msg : `str` Error message suffix """ transform = TransformClass(mapping) inverse = transform.getInverse() inverseInverse = inverse.getInverse() self.checkInverseTransformation(transform, inverse, msg=msg) self.checkInverseTransformation(inverse, inverseInverse, msg=msg) self.checkTransformation(inverseInverse, mapping, msg=msg) def checkGetJacobian(self, fromName, toName): """Test Transform<fromName>To<toName>.getJacobian Parameters ---------- fromName, toName : `str` Endpoint name prefix for "from" and "to" endpoints, respectively, e.g. "Point2" for `lsst.afw.geom.Point2Endpoint` """ transformClassName = "Transform{}To{}".format(fromName, toName) TransformClass = getattr(afwGeom, transformClassName) baseMsg = "TransformClass={}".format(TransformClass.__name__) for nIn, nOut in itertools.product(self.goodNAxes[fromName], self.goodNAxes[toName]): msg = "{}, nIn={}, nOut={}".format(baseMsg, nIn, nOut) polyMap = makeForwardPolyMap(nIn, nOut) transform = TransformClass(polyMap) fromEndpoint = transform.fromEndpoint # Test multiple points to ensure correct functional form rawInPoint = self.makeRawPointData(nIn) inPoint = fromEndpoint.pointFromData(rawInPoint) jacobian = transform.getJacobian(inPoint) assert_allclose(jacobian, self.makeJacobian(nIn, nOut, rawInPoint), err_msg=msg) rawInPoint = self.makeRawPointData(nIn, 0.111) inPoint = fromEndpoint.pointFromData(rawInPoint) jacobian = transform.getJacobian(inPoint) assert_allclose(jacobian, self.makeJacobian(nIn, nOut, rawInPoint), err_msg=msg) def checkThen(self, fromName, midName, toName): """Test Transform<fromName>To<midName>.then(Transform<midName>To<toName>) Parameters ---------- fromName : `str` the prefix of the starting endpoint (e.g., "Point2" for a Point2Endpoint) for the final, concatenated Transform midName : `str` the prefix for the shared endpoint where two Transforms will be concatenated toName : `str` the prefix of the ending endpoint for the final, concatenated Transform """ TransformClass1 = getattr(afwGeom, "Transform{}To{}".format(fromName, midName)) TransformClass2 = getattr(afwGeom, "Transform{}To{}".format(midName, toName)) baseMsg = "{}.then({})".format(TransformClass1.__name__, TransformClass2.__name__) for nIn, nMid, nOut in itertools.product(self.goodNAxes[fromName], self.goodNAxes[midName], self.goodNAxes[toName]): msg = "{}, nIn={}, nMid={}, nOut={}".format( baseMsg, nIn, nMid, nOut) polyMap1 = makeTwoWayPolyMap(nIn, nMid) transform1 = TransformClass1(polyMap1) polyMap2 = makeTwoWayPolyMap(nMid, nOut) transform2 = TransformClass2(polyMap2) transform = transform1.then(transform2) fromEndpoint = transform1.fromEndpoint toEndpoint = transform2.toEndpoint inPoint = fromEndpoint.pointFromData(self.makeRawPointData(nIn)) outPointMerged = transform.applyForward(inPoint) outPointSeparate = transform2.applyForward( transform1.applyForward(inPoint)) assert_allclose(toEndpoint.dataFromPoint(outPointMerged), toEndpoint.dataFromPoint(outPointSeparate), err_msg=msg) outPoint = toEndpoint.pointFromData(self.makeRawPointData(nOut)) inPointMerged = transform.applyInverse(outPoint) inPointSeparate = transform1.applyInverse( transform2.applyInverse(outPoint)) assert_allclose( fromEndpoint.dataFromPoint(inPointMerged), fromEndpoint.dataFromPoint(inPointSeparate), err_msg=msg) # Mismatched number of axes should fail if midName == "Generic": nIn = self.goodNAxes[fromName][0] nOut = self.goodNAxes[toName][0] polyMap = makeTwoWayPolyMap(nIn, 3) transform1 = TransformClass1(polyMap) polyMap = makeTwoWayPolyMap(2, nOut) transform2 = TransformClass2(polyMap) with self.assertRaises(InvalidParameterError): transform = transform1.then(transform2) # Mismatched types of endpoints should fail if fromName != midName: # Use TransformClass1 for both args to keep test logic simple outName = midName joinNAxes = set(self.goodNAxes[fromName]).intersection( self.goodNAxes[outName]) for nIn, nMid, nOut in itertools.product(self.goodNAxes[fromName], joinNAxes, self.goodNAxes[outName]): polyMap = makeTwoWayPolyMap(nIn, nMid) transform1 = TransformClass1(polyMap) polyMap = makeTwoWayPolyMap(nMid, nOut) transform2 = TransformClass1(polyMap) with self.assertRaises(InvalidParameterError): transform = transform1.then(transform2) def assertTransformsEqual(self, transform1, transform2): """Assert that two transforms are equal""" self.assertEqual(type(transform1), type(transform2)) self.assertEqual(transform1.fromEndpoint, transform2.fromEndpoint) self.assertEqual(transform1.toEndpoint, transform2.toEndpoint) self.assertEqual(transform1.getMapping(), transform2.getMapping()) fromEndpoint = transform1.fromEndpoint toEndpoint = transform1.toEndpoint mapping = transform1.getMapping() nIn = mapping.nIn nOut = mapping.nOut if mapping.hasForward: nPoints = 7 # arbitrary rawInArray = self.makeRawArrayData(nPoints, nIn) inArray = fromEndpoint.arrayFromData(rawInArray) outArray = transform1.applyForward(inArray) outData = toEndpoint.dataFromArray(outArray) outArrayRoundTrip = transform2.applyForward(inArray) outDataRoundTrip = toEndpoint.dataFromArray(outArrayRoundTrip) assert_allclose(outData, outDataRoundTrip) if mapping.hasInverse: nPoints = 7 # arbitrary rawOutArray = self.makeRawArrayData(nPoints, nOut) outArray = toEndpoint.arrayFromData(rawOutArray) inArray = transform1.applyInverse(outArray) inData = fromEndpoint.dataFromArray(inArray) inArrayRoundTrip = transform2.applyInverse(outArray) inDataRoundTrip = fromEndpoint.dataFromArray(inArrayRoundTrip) assert_allclose(inData, inDataRoundTrip) def checkPersistence(self, transform): """Check persistence of a transform """ className = type(transform).__name__ # check writeString and readString transformStr = transform.writeString() serialVersion, serialClassName, serialRest = transformStr.split(" ", 2) self.assertEqual(int(serialVersion), 1) self.assertEqual(serialClassName, className) badStr1 = " ".join(["2", serialClassName, serialRest]) with self.assertRaises(lsst.pex.exceptions.InvalidParameterError): transform.readString(badStr1) badClassName = "x" + serialClassName badStr2 = " ".join(["1", badClassName, serialRest]) with self.assertRaises(lsst.pex.exceptions.InvalidParameterError): transform.readString(badStr2) transformFromStr1 = transform.readString(transformStr) self.assertTransformsEqual(transform, transformFromStr1) # check transformFromString transformFromStr2 = afwGeom.transformFromString(transformStr) self.assertTransformsEqual(transform, transformFromStr2) # Check pickling self.assertTransformsEqual(transform, pickle.loads(pickle.dumps(transform))) # Check afw::table::io persistence round-trip with lsst.utils.tests.getTempFilePath(".fits") as filename: transform.writeFits(filename) self.assertTransformsEqual(transform, type(transform).readFits(filename))