Expand the unit test for ComputePixelToDistortedPixel

lsst · Oct 26, 2020 · bf3c8df · bf3c8df
1 parent 89af182
commit bf3c8df
Showing 1 changed file with 45 additions and 24 deletions.
diff --git a/tests/test_wcsUtils.py b/tests/test_wcsUtils.py
@@ -91,8 +91,8 @@ class ComputePixelToDistortedPixelTestCase(BaseTestCase):
     def testNoDistortion(self):
         """Test computePixelToDistortedPixel without distortion
 
-        Use an affine transform for pixelToFocalPlane; the transform
-        returned by computePixelToDistortedPixel should be the identity transform
+        Use an affine transform for pixelToFocalPlane; the transform returned
+        by computePixelToDistortedPixel should be the identity transform
         """
         focalPlaneToFieldAngle = self.makeAffineTransform(scale=self.radPerMm)
         pixelToDistortedPixel = computePixelToDistortedPixel(
@@ -109,36 +109,57 @@ def testNoDistortion(self):
     def testDistortion(self):
         """Test computePixelToDistortedPixel with distortion
 
-        pixelToDistortedPixel -> self.tanWcs should match a WCS
-        created with makeDistortedTanWcs
+        Compare the output of computePixelToDistortedPixel to a non-affine
+        transform for pixelToFocalPlane
         """
-        focalPlaneToFieldAngle = afwGeom.makeRadialTransform([0.0, self.radPerMm, 0.0, self.radPerMm])
+        # Compute a pixelToDistortedSky that matches self.tanWCS at the center of the field;
+        # the amount of distortion is 10s of pixels over the detector
+        fieldAngleToFocalPlane = afwGeom.makeRadialTransform([0.0, 1/self.radPerMm, 0.0, 1000/self.radPerMm])
+        focalPlaneToFieldAngle = fieldAngleToFocalPlane.inverted()
+        focalPlaneToTanFieldAngle = self.makeAffineTransform(scale=self.radPerMm)
         pixelToDistortedPixel = computePixelToDistortedPixel(
             pixelToFocalPlane=self.pixelToFocalPlane,
             focalPlaneToFieldAngle=focalPlaneToFieldAngle,
         )
-        # Do not try to make pixelToDistortedPixel -> self.tanWcs into a WCS
-        # because the frame names will be wrong; use a TransformPoint2Tolsst.geom.SpherePoint instead
+
         tanWcsTransform = afwGeom.TransformPoint2ToSpherePoint(self.tanWcs.getFrameDict())
         pixelToDistortedSky = pixelToDistortedPixel.then(tanWcsTransform)
 
-        wcs = makeDistortedTanWcs(
-            tanWcs=self.tanWcs,
-            pixelToFocalPlane=self.pixelToFocalPlane,
-            focalPlaneToFieldAngle=focalPlaneToFieldAngle,
-        )
-
-        bboxD = lsst.geom.Box2D(self.bbox)
-        pixelPoints = bboxD.getCorners()
-        pixelPoints.append(bboxD.getCenter())
-
-        skyPoints1 = pixelToDistortedSky.applyForward(pixelPoints)
-        skyPoints2 = wcs.pixelToSky(pixelPoints)
-        self.assertSpherePointListsAlmostEqual(skyPoints1, skyPoints2)
-
-        pixelPoints1 = pixelToDistortedSky.applyInverse(skyPoints1)
-        pixelPoints2 = wcs.skyToPixel(skyPoints1)
-        assert_allclose(pixelPoints1, pixelPoints2)
+        # At the center of the focal plane both WCS should give the same sky position
+        pixelAtCtr = self.pixelToFocalPlane.applyInverse(lsst.geom.Point2D(0, 0))
+        tanSkyAtCtr = self.tanWcs.pixelToSky(pixelAtCtr)
+        skyAtCtr = pixelToDistortedSky.applyForward(pixelAtCtr)
+        self.assertPairsAlmostEqual(tanSkyAtCtr, skyAtCtr)
+
+        # At all reasonable sky points the following field angles should be almost equal:
+        #   sky -> tanWcs.skyToPixel -> pixelToFocalPlane -> focalPlaneToTanFieldAngle
+        #   sky -> pixelToDistortedPixel.applyInverse -> pixelToFocalPlane -> focalPlaneToFieldAngle
+        # where focalPlaneToTanFieldAngle is the linear approximation to
+        # focalPlaneToFieldAngle at the center of the field (where tanWcs and pixelToDistortedSky match),
+        # since for a given pointing, field angle gives position on the sky
+        skyPoints = self.tanWcs.pixelToSky(self.pixelPoints)
+
+        tanFieldAnglePoints = focalPlaneToTanFieldAngle.applyForward(
+            self.pixelToFocalPlane.applyForward(self.tanWcs.skyToPixel(skyPoints)))
+        fieldAnglePoints = focalPlaneToFieldAngle.applyForward(
+            self.pixelToFocalPlane.applyForward(pixelToDistortedSky.applyInverse(skyPoints)))
+        assert_allclose(tanFieldAnglePoints, fieldAnglePoints)
+
+        # The inverse should also be true: for a set of field angle points
+        # the following sky positions should be almost equal:
+        # fieldAngle -> fieldAngleToTanFocalPlane -> focalPlaneToPixel -> tanWcs.pixelToSky
+        # fieldAngle -> fieldAngleToFocalPlane -> focalPlaneToPixel -> pixelToDistortedPixel.applyForward
+        focalPlaneToPixel = self.pixelToFocalPlane.inverted()
+        fieldAngleToTanFocalPlane = focalPlaneToTanFieldAngle.inverted()
+        tanSkyPoints2 = self.tanWcs.pixelToSky(
+            focalPlaneToPixel.applyForward(
+                fieldAngleToTanFocalPlane.applyForward(fieldAnglePoints)))
+
+        skyPoints2 = pixelToDistortedSky.applyForward(
+            focalPlaneToPixel.applyForward(
+                fieldAngleToFocalPlane.applyForward(fieldAnglePoints)))
+
+        self.assertSpherePointListsAlmostEqual(tanSkyPoints2, skyPoints2)
 
 
 class DetailTestCase(lsst.utils.tests.TestCase):