Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DM-22114: Remove deprecated makeDistortedTanWcs #541

Merged
merged 2 commits into from
Oct 29, 2020
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
98 changes: 1 addition & 97 deletions python/lsst/afw/geom/wcsUtils/wcsUtilsContinued.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,13 +20,10 @@
# see <http://www.lsstcorp.org/LegalNotices/>.
#

__all__ = ["getSipMatrixFromMetadata", "makeDistortedTanWcs", "computePixelToDistortedPixel"]

from deprecated.sphinx import deprecated
__all__ = ["getSipMatrixFromMetadata", "computePixelToDistortedPixel"]

import lsst.geom
from ..transformFactory import linearizeTransform, makeTransform
from ..skyWcs import makeModifiedWcs
from .wcsUtils import _getSipMatrixFromMetadata


Expand Down Expand Up @@ -61,99 +58,6 @@ def getSipMatrixFromMetadata(metadata, name):
return arr


@deprecated(reason="Camera geometry-based SkyWcs are now set when reading raws. To be removed after v20.",
category=FutureWarning)
def makeDistortedTanWcs(tanWcs, pixelToFocalPlane, focalPlaneToFieldAngle):
"""Compute a WCS that includes a model of optical distortion.

This is useful in the common case that the initial WCS entirely ignores
the effect of optical distortion.

Parameters
----------
tanWcs : `lsst.afw.geom.SkyWcs`
A pure TAN WCS, such as is usually provided in raw data.
This should have no existing compensation for optical distortion
(though it may include an ``ACTUAL_PIXELS`` frame to model pixel-level
distortions).
pixelToFocalPlane : `lsst.afw.geom.TransformPoint2ToPoint2`
Transform parent pixel coordinates to focal plane coordinates.
This models the location of the CCD on the focal plane
and is almost always an affine transformation.
This can be obtained from the detector of an exposure.
focalPlaneToFieldAngle : `lsst.afw.geom.TransformPoint2ToPoint2`
Transform focal plane coordinates to field angle coordinates.
This is a model for optical distortion, and is often a radial
polynomial. This can be obtained from the camera geometry.


Returns
-------
lsst.afw.geom.SkyWcs
A copy of `tanWcs` that includes the effect of optical distortion.

Raises
------
RuntimeError
If the current frame of `wcs` is not a SkyFrame;
LookupError
If 2-dimensional Frames with Domain "PIXELS" and "IWC"
are not all found.
"""
# The math is as follows:
#
# Our input TAN WCS is:
# tanWcs = PIXELS frame -> pixelToIwc -> IWC frame -> iwcToSky -> SkyFrame
# See lsst.afw.geom.SkyWcs for a description of these frames.
# tanWcs may also contain an ACTUAL_PIXELS frame before the PIXELS frame;
# if so it will be preserved, but it is irrelevant to the computation
# and so not discussed further.
#
# Our desired WCS must still contain the PIXELS and IWC frames.
# The distortion will be inserted just after the PIXELS frame,
# So the new WCS will be as follows:
# wcs = PIXELS frame -> pixelToDistortedPixel -> pixelToIwc -> IWC frame -> iwcToSky -> SkyFrame
#
# We compute pixelToDistortedPixel as follows...
#
# We will omit the frames from now on, for simplicity. Thus:
# tanWcs = pixelToIwc -> iwcToSksy
# and:
# wcs = pixelToDistortedPixel -> pixelToIwc -> iwcToSky
#
# We also know pixelToFocalPlane and focalPlaneToFieldAngle,
# and can use them as follows:
#
# The tan WCS can be expressed as:
# tanWcs = pixelToFocalPlane -> focalPlaneToTanFieldAngle -> fieldAngleToIwc -> iwcToSky
# where:
# - focalPlaneToTanFieldAngle is the linear approximation to
# focalPlaneToFieldAngle at the center of the focal plane
#
# The desired WCS can be expressed as:
# wcs = pixelToFocalPlane -> focalPlaneToFieldAngle -> fieldAngleToIwc -> iwcToSky
#
# By equating the two expressions for tanWcs, we get:
# pixelToIwc = pixelToFocalPlane -> focalPlaneToTanFieldAngle -> fieldAngleToIwc
# fieldAngleToIwc = tanFieldAngleToFocalPlane -> focalPlaneToPixel -> pixelToIwc
#
# By equating the two expressions for desired wcs we get:
# pixelToDistortedPixel -> pixelToIwc = pixelToFocalPlane -> focalPlaneToFieldAngle -> fieldAngleToIwc
#
# Substitute our expression for fieldAngleToIwc from tanWcs into the
# previous equation, we get:
# pixelToDistortedPixel -> pixelToIwc
# = pixelToFocalPlane -> focalPlaneToFieldAngle -> tanFieldAngleToFocalPlane -> focalPlaneToPixel
# -> pixelToIwc
#
# Thus:
# pixelToDistortedPixel
# = pixelToFocalPlane -> focalPlaneToFieldAngle -> tanFieldAngleToFocalPlane -> focalPlaneToPixel

pixelToDistortedPixel = computePixelToDistortedPixel(pixelToFocalPlane, focalPlaneToFieldAngle)
return makeModifiedWcs(pixelTransform=pixelToDistortedPixel, wcs=tanWcs, modifyActualPixels=False)


def computePixelToDistortedPixel(pixelToFocalPlane, focalPlaneToFieldAngle):
"""Compute the transform ``pixelToDistortedPixel``, which applies optical
distortion specified by ``focalPlaneToFieldAngle``.
Expand Down
108 changes: 29 additions & 79 deletions tests/test_wcsUtils.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
from lsst.afw.geom.wcsUtils import createTrivialWcsMetadata, deleteBasicWcsMetadata, \
getCdMatrixFromMetadata, getSipMatrixFromMetadata, getImageXY0FromMetadata, \
hasSipMatrix, makeSipMatrixMetadata, makeTanSipMetadata, \
computePixelToDistortedPixel, makeDistortedTanWcs
computePixelToDistortedPixel


def makeRotationMatrix(angle, scale):
Expand All @@ -47,8 +47,7 @@ def makeRotationMatrix(angle, scale):


class BaseTestCase(lsst.utils.tests.TestCase):
"""Base class for testing makeDistortedTanWcs and
computePixelToDistortedPixel
"""Base class for testing computePixelToDistortedPixel
"""
def setUp(self):
# define the position and size of one CCD in the focal plane
Expand Down Expand Up @@ -85,133 +84,84 @@ def makeAffineTransform(self, offset=(0, 0), rotation=0*lsst.geom.degrees, scale
return afwGeom.makeTransform(rotScale*offset)


class MakeDistortedTanWcsTestCase(BaseTestCase):
"""Test lsst.afw.geom.makeDistortedTanWcs
class ComputePixelToDistortedPixelTestCase(BaseTestCase):
"""Test lsst.afw.geom.computePixelToDistortedPixel
"""

def testNoDistortion(self):
"""Test makeDistortedTanWcs using an affine transform for pixelToFocalPlane
"""Test computePixelToDistortedPixel without distortion

Construct pixelToFocalPlane to match the plate scale used to
generate self.tanWcs, the input to makeDistortedTanWcs. Thus the WCS
returned by makeDistortedTanWcs should match self.tanWcs.
Use an affine transform for pixelToFocalPlane; the transform returned
arunkannawadi marked this conversation as resolved.
Show resolved Hide resolved
by computePixelToDistortedPixel should be the identity transform
"""
focalPlaneToFieldAngle = self.makeAffineTransform(scale=self.radPerMm)
wcs = makeDistortedTanWcs(
tanWcs=self.tanWcs,
pixelToDistortedPixel = computePixelToDistortedPixel(
pixelToFocalPlane=self.pixelToFocalPlane,
focalPlaneToFieldAngle=focalPlaneToFieldAngle,
)
self.assertWcsAlmostEqualOverBBox(wcs, self.tanWcs, bbox=self.bbox)
bboxD = lsst.geom.Box2D(self.bbox)
pixelPoints = bboxD.getCorners()
pixelPoints.append(bboxD.getCenter())

assert_allclose(pixelToDistortedPixel.applyForward(pixelPoints), pixelPoints)
assert_allclose(pixelToDistortedPixel.applyInverse(pixelPoints), pixelPoints)

def testDistortion(self):
"""Test makeDistortedTanWcs using a non-affine transform for pixelToFocalPlane
"""Test computePixelToDistortedPixel with distortion

Compare the output of computePixelToDistortedPixel to a non-affine
transform for pixelToFocalPlane
"""
# Compute a distorted wcs that matches self.tanWcs at the center of the field;
# 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)
wcs = makeDistortedTanWcs(
tanWcs=self.tanWcs,
pixelToDistortedPixel = computePixelToDistortedPixel(
pixelToFocalPlane=self.pixelToFocalPlane,
focalPlaneToFieldAngle=focalPlaneToFieldAngle,
)

tanWcsTransform = afwGeom.TransformPoint2ToSpherePoint(self.tanWcs.getFrameDict())
pixelToDistortedSky = pixelToDistortedPixel.then(tanWcsTransform)

# 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 = wcs.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 -> wcs.skyToPixel -> pixelToFocalPlane -> focalPlaneToFieldAngle
# sky -> pixelToDistortedPixel.applyInverse -> pixelToFocalPlane -> focalPlaneToFieldAngle
# where focalPlaneToTanFieldAngle is the linear approximation to
# focalPlaneToFieldAngle at the center of the field (where tanWcs and wcs match),
# 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(wcs.skyToPixel(skyPoints)))
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 -> wcs.pixelToSky
# fieldAngle -> fieldAngleToFocalPlane -> focalPlaneToPixel -> pixelToDistortedPixel.applyForward
focalPlaneToPixel = self.pixelToFocalPlane.inverted()
fieldAngleToTanFocalPlane = focalPlaneToTanFieldAngle.inverted()
tanSkyPoints2 = self.tanWcs.pixelToSky(
focalPlaneToPixel.applyForward(
fieldAngleToTanFocalPlane.applyForward(fieldAnglePoints)))

skyPoints2 = wcs.pixelToSky(
skyPoints2 = pixelToDistortedSky.applyForward(
focalPlaneToPixel.applyForward(
fieldAngleToFocalPlane.applyForward(fieldAnglePoints)))

self.assertSpherePointListsAlmostEqual(tanSkyPoints2, skyPoints2)


class ComputePixelToDistortedPixelTestCase(BaseTestCase):
"""Test lsst.afw.geom.computePixelToDistortedPixel
"""

def testNoDistortion(self):
"""Test computePixelToDistortedPixel without distortion

Use an affine transform for pixelToFocalPlane; the transform
returned by computePixelToDistortedPixel should be the identity transform
"""
focalPlaneToFieldAngle = self.makeAffineTransform(scale=self.radPerMm)
pixelToDistortedPixel = computePixelToDistortedPixel(
pixelToFocalPlane=self.pixelToFocalPlane,
focalPlaneToFieldAngle=focalPlaneToFieldAngle,
)
bboxD = lsst.geom.Box2D(self.bbox)
pixelPoints = bboxD.getCorners()
pixelPoints.append(bboxD.getCenter())

assert_allclose(pixelToDistortedPixel.applyForward(pixelPoints), pixelPoints)
assert_allclose(pixelToDistortedPixel.applyInverse(pixelPoints), pixelPoints)

def testDistortion(self):
"""Test computePixelToDistortedPixel with distortion

pixelToDistortedPixel -> self.tanWcs should match a WCS
created with makeDistortedTanWcs
"""
focalPlaneToFieldAngle = afwGeom.makeRadialTransform([0.0, self.radPerMm, 0.0, 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)


class DetailTestCase(lsst.utils.tests.TestCase):
"""Test functions in the detail sub-namespace
"""
Expand Down