Add Transform-based overload of warpImage

Add an overload of warpImage that uses a
Transform<Point2Endpoint, Point2Endpoint>
to describe the transform.

examples/timeWarpExposureUsingTransform.py

@@ -0,0 +1,153 @@
+#!/usr/bin/env python
+
+#
+# LSST Data Management System
+# Copyright 2017 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/>.
+#
+
+from __future__ import absolute_import, division, print_function
+from builtins import range
+import sys
+import os
+import time
+
+import lsst.utils
+import lsst.afw.geom as afwGeom
+import lsst.afw.image as afwImage
+import lsst.afw.math as afwMath
+
+MaxIter = 20
+MaxTime = 1.0  # seconds
+WarpSubregion = True  # set False to warp more pixels
+SaveImages = False
+
+afwdataDir = lsst.utils.getPackageDir("afwdata")
+
+InputExposurePath = os.path.join(
+    afwdataDir, "ImSim/calexp/v85408556-fr/R23/S11.fits")
+
+
+def timeWarp(destMaskedImage, srcMaskedImage, destToSrc, warpingControl):
+    """Time warpImage
+
+    Parameters
+    ----------
+    destMaskedImage : `lsst.afw.image.MaskedImage`
+        Destination (output) masked image
+    srcMaskedImage : `lsst.afw.image.MaskedImage`
+        Source (input) masked image
+    destToSrc : `lsst.afw.geom.TransformPoint2ToPoint2`
+        Transform from source pixels to destination pixels
+    warpingControl : `lsst.afw.geom.WarpingControl`
+        Warning control parameters
+
+    Returns
+    -------
+    `tuple(float, int, int)`
+        - Elapsed time in seconds
+        - Number of iterations
+        - Number of good pixels
+    """
+    startTime = time.time()
+    for nIter in range(1, MaxIter + 1):
+        goodPix = afwMath.warpImage(destMaskedImage, srcMaskedImage, destToSrc, warpingControl)
+        endTime = time.time()
+        if endTime - startTime > MaxTime:
+            break
+
+    return (endTime - startTime, nIter, goodPix)
+
+
+def run():
+    if len(sys.argv) < 2:
+        srcExposure = afwImage.ExposureF(InputExposurePath)
+        if WarpSubregion:
+            bbox = afwGeom.Box2I(afwGeom.Point2I(
+                0, 0), afwGeom.Extent2I(2000, 2000))
+            srcExposure = afwImage.ExposureF(
+                srcExposure, bbox, afwImage.LOCAL, False)
+    else:
+        srcExposure = afwImage.ExposureF(sys.argv[1])
+    srcMaskedImage = srcExposure.getMaskedImage()
+    srcMetadata = afwImage.DecoratedImageF(InputExposurePath).getMetadata()
+    srcWcs = afwGeom.SkyWcs(srcMetadata)
+    srcDim = srcExposure.getDimensions()
+    srcCtrPos = afwGeom.Box2D(srcMaskedImage.getBBox()).getCenter()
+    srcCtrSky = srcWcs.tranForward(srcCtrPos)
+    srcScale = srcWcs.getPixelScale()
+
+    # make the destination exposure small enough that even after rotation and offset
+    # (by reasonable amounts) there are no edge pixels
+    destDim = afwGeom.Extent2I(*[int(sd * 0.5) for sd in srcDim])
+    destMaskedImage = afwImage.MaskedImageF(destDim)
+    destCrPix = afwGeom.Box2D(destMaskedImage.getBBox()).getCenter()
+
+    maskKernelName = ""
+    cacheSize = 0
+
+    print("Warping", InputExposurePath)
+    print("Source (sub)image size:", srcDim)
+    print("Destination image size:", destDim)
+    print()
+
+    print("test# interp  scaleFac    destSkyOff    rotAng   kernel   goodPix time/iter")
+    print('       (pix)              (bear°, len")    (deg)                      (sec)')
+    testNum = 1
+    for interpLength in (0, 1, 5, 10):
+        for scaleFac in (1.2,):
+            destScale = srcScale / scaleFac
+            for offsetOrientDegLenArcsec in ((0.0, 0.0),):  # ((0.0, 0.0), (-35.0, 10.5)):
+                # offset (bearing, length) from sky at center of source to sky at center of dest
+                offset = (offsetOrientDegLenArcsec[0] * afwGeom.degrees,
+                          offsetOrientDegLenArcsec[1] * afwGeom.arcseconds)
+                destCtrSky = srcCtrSky.offset(*offset)
+                for rotAngDeg, kernelName in (
+                    (0.0, "bilinear"),
+                    (0.0, "lanczos2"),
+                    (0.0, "lanczos3"),
+                    (45.0, "lanczos3"),
+                ):
+                    warpingControl = afwMath.WarpingControl(
+                        kernelName,
+                        maskKernelName,
+                        cacheSize,
+                        interpLength,
+                    )
+                    destWcs = afwGeom.SkyWcs(
+                        crpix = destCrPix,
+                        crval = destCtrSky,
+                        cdMatrix = afwGeom.makeCdMatrix(scale = destScale,
+                                                        orientation = rotAngDeg * afwGeom.degrees,
+                                                        flipX = False))
+                    destToSrc = srcWcs.getInverse().of(destWcs)
+                    dTime, nIter, goodPix = timeWarp(
+                        destMaskedImage, srcMaskedImage, destToSrc, warpingControl)
+                    print("%4d  %5d  %8.1f  %6.1f, %6.1f  %7.1f %10s %8d %6.2f" % (
+                        testNum, interpLength, scaleFac, offsetOrientDegLenArcsec[0],
+                        offsetOrientDegLenArcsec[1], rotAngDeg, kernelName, goodPix, dTime/float(nIter)))
+
+                    if SaveImages:
+                        destMaskedImage.writeFits(
+                            "warpedMaskedImage%03d.fits" % (testNum,))
+                    testNum += 1
+
+
+if __name__ == "__main__":
+    run()

include/lsst/afw/math/warpExposure.h

@@ -460,6 +460,30 @@ int warpImage(DestImageT &destImage,                            ///< remapped %i
               ///< use this value for undefined (edge) pixels
               );
 
+/**
+ * @brief A variant of warpImage that uses a Transform<Point2Endpoint, Point2Endpoint>
+ * instead of a pair of WCS to describe the transformation.
+ *
+ * @param[in,out] destImage  Destination image; all pixels are set
+ * @param[in] srcImage  Source image
+ * @param[in] destToSrc  Transformation from destination to source pixels in parent coordinates.
+ *    This can be computed as srcWcs.getInverse().of(destWcs)
+ *    because WCS computes pixels to sky in the forward direction
+ * @param[in] control  Warning control parameters
+ * @param[in] padValue  Value used for pixels in the destination image that are outside
+ *   the region of pixels that can be computed from the source image
+ * @return the number of good pixels
+ */
+template <typename DestImageT, typename SrcImageT>
+int warpImage(DestImageT &destImage,
+              SrcImageT const &srcImage,
+              geom::Transform<geom::Point2Endpoint, geom::Point2Endpoint> const & destToSrc,
+              WarpingControl const &control,
+              typename DestImageT::SinglePixel padValue = lsst::afw::math::edgePixel<DestImageT>(
+                      typename lsst::afw::image::detail::image_traits<DestImageT>::image_category())
+              );
+
+
 /**
  * Warp an image with a LinearTranform about a specified point.
  *

python/lsst/afw/math/warpExposure.cc

@@ -97,6 +97,13 @@ void declareImageWarpingFunctions(py::module &mod) {
                                  warpImage<DestImageT, SrcImageT>,
             "destImage"_a, "srcImage"_a, "xyTransform"_a, "control"_a, "padValue"_a = EdgePixel);
 
+    mod.def("warpImage",
+            (int (*)(DestImageT &, SrcImageT const &,
+                     geom::Transform<geom::Point2Endpoint, geom::Point2Endpoint> const &,
+                     WarpingControl const &, typename DestImageT::SinglePixel)) &
+                    warpImage<DestImageT, SrcImageT>,
+            "destImage"_a, "srcImage"_a, "destToSrc"_a, "control"_a, "padValue"_a = EdgePixel);
+
     mod.def("warpCenteredImage", &warpCenteredImage<DestImageT, SrcImageT>, "destImage"_a, "srcImage"_a,
             "linearTransform"_a, "centerPoint"_a, "control"_a, "padValue"_a = EdgePixel);
 }