DM-9599: Support concatenation of Transforms

include/lsst/afw/geom/Transform.h

@@ -60,7 +60,7 @@ since the base and current frames in the FrameSet can be checked against by the
 because data must be copied when converting from LSST data types to the type used by astshim,
 so it didn't seem worth the bother.
 */
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 class Transform {
 public:
     using FromArray = typename FromEndpoint::Array;
@@ -189,6 +189,29 @@ class Transform {
      */
     Eigen::MatrixXd getJacobian(FromPoint const &x) const;
 
+    /**
+     * Concatenate two Transforms.
+     *
+     * @tparam FirstFromEndpoint the starting Endpoint of `first`
+     * @param first the Transform to apply before this one
+     * @returns a Transform that first applies `first` to its input, and then
+     *          this Transform to the result. Its inverse shall first apply the
+     *          inverse of this Transform, then the inverse of `first`.
+     *
+     * @throws InvalidParameterErrror Thrown if `first.getToEndpoint()` and
+     *                                `this->getFromEndpoint()` do not have
+     *                                the same number of axes.
+     * @exceptsafe Provides basic exception safety.
+     *
+     * More than two Transforms can be combined in series. For example:
+     *
+     *     auto skyFromPixels = skyFromPupil.of(pupilFromFp)
+     *                                      .of(fpFromPixels);
+     */
+    template <class FirstFromEndpoint>
+    Transform<FirstFromEndpoint, ToEndpoint> of(
+            Transform<FirstFromEndpoint, FromEndpoint> const &first) const;
+
 private:
     FromEndpoint const _fromEndpoint;
     std::shared_ptr<const ast::FrameSet> _frameSet;
@@ -202,7 +225,7 @@ The format is "Transform<_fromEndpoint_, _toEndpoint_>"
 where _fromEndpoint_ and _toEndpoint_ are the appropriate endpoint printed to the ostream;
 for example "Transform<GenericEndpoint(4), Point3Endpoint()>"
 */
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 std::ostream &operator<<(std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform);
 
 }  // geom

python/lsst/afw/geom/transform/transform.cc

@@ -42,7 +42,7 @@ namespace {
 
 // Return a string consisting of "_pythonClassName_[_fromNAxes_->_toNAxes_]",
 // for example "TransformGenericToPoint3[4->3]"
-template <typename Class>
+template <class Class>
 std::string formatStr(Class const &self, std::string const &pyClassName) {
     std::ostringstream os;
     os << pyClassName;
@@ -51,10 +51,21 @@ std::string formatStr(Class const &self, std::string const &pyClassName) {
     return os.str();
 }
 
+template <class ExtraEndpoint, class FromEndpoint, class ToEndpoint, class PyClass>
+void declareMethodTemplates(PyClass &cls) {
+    using FirstTransform = Transform<ExtraEndpoint, FromEndpoint>;
+    using SecondTransform = Transform<FromEndpoint, ToEndpoint>;
+    using FinalTransform = Transform<ExtraEndpoint, ToEndpoint>;
+    // Need Python-specific logic to give sensible errors for mismatched Transform types
+    cls.def("_of", (FinalTransform (SecondTransform::*)(FirstTransform const &) const) &
+                           SecondTransform::template of<ExtraEndpoint>,
+            "first"_a);
+}
+
 // Declare Transform<FromEndpoint, ToEndpoint> using python class name TransformFrom<X>To<Y>
 // where <X> and <Y> are the name of the from endpoint and to endpoint class, respectively,
 // for example TransformFromGenericToPoint3
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 void declareTransform(py::module &mod, std::string const &fromName, std::string const &toName) {
     using Class = Transform<FromEndpoint, ToEndpoint>;
     using ToPoint = typename ToEndpoint::Point;
@@ -81,7 +92,15 @@ void declareTransform(py::module &mod, std::string const &fromName, std::string
     cls.def("tranInverse", (FromArray (Class::*)(ToArray const &) const) & Class::tranInverse, "array"_a);
     cls.def("tranInverse", (FromPoint (Class::*)(ToPoint const &) const) & Class::tranInverse, "point"_a);
     cls.def("getInverse", &Class::getInverse);
+    /* Need some extra handling of ndarray return type in Python to prevent dimensions
+     * of length 1 from being deleted */
     cls.def("_getJacobian", &Class::getJacobian);
+
+    declareMethodTemplates<GenericEndpoint, FromEndpoint, ToEndpoint>(cls);
+    declareMethodTemplates<Point2Endpoint, FromEndpoint, ToEndpoint>(cls);
+    declareMethodTemplates<Point3Endpoint, FromEndpoint, ToEndpoint>(cls);
+    declareMethodTemplates<SpherePointEndpoint, FromEndpoint, ToEndpoint>(cls);
+
     // str(self) = "<Python class name>[<nIn>-><nOut>]"
     cls.def("__str__", [pyClassName](Class const &self) { return formatStr(self, pyClassName); });
     // repr(self) = "lsst.afw.geom.<Python class name>[<nIn>-><nOut>]"

python/lsst/afw/geom/transform/transformContinued.py

@@ -21,6 +21,8 @@
 #
 from __future__ import absolute_import, division, print_function
 
+from lsst.pex.exceptions import InvalidParameterError
+
 from . import transform
 
 __all__ = []
@@ -33,10 +35,20 @@ def getJacobian(self, x):
                     self.getFromEndpoint().getNAxes())
     return matrix
 
+
+def of(self, first):
+    if first.getToEndpoint() == self.getFromEndpoint():
+        return self._of(first)
+    else:
+        raise InvalidParameterError(
+            "Cannot concatenate %r and %r: endpoints do not match."
+            % (first, self))
+
 endpoints = ("Generic", "Point2", "Point3", "SpherePoint")
 
 for fromPoint in endpoints:
     for toPoint in endpoints:
         name = "Transform" + fromPoint + "To" + toPoint
         cls = getattr(transform, name)
         cls.getJacobian = getJacobian
+        cls.of = of

src/geom/Transform.cc

@@ -37,7 +37,7 @@ namespace lsst {
 namespace afw {
 namespace geom {
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 Transform<FromEndpoint, ToEndpoint>::Transform(ast::Mapping const &mapping, bool simplify)
         : _fromEndpoint(mapping.getNin()), _frameSet(), _toEndpoint(mapping.getNout()) {
     auto fromFrame = _fromEndpoint.makeFrame();
@@ -49,7 +49,7 @@ Transform<FromEndpoint, ToEndpoint>::Transform(ast::Mapping const &mapping, bool
     }
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 Transform<FromEndpoint, ToEndpoint>::Transform(ast::FrameSet const &frameSet, bool simplify)
         : _fromEndpoint(frameSet.getNin()), _frameSet(), _toEndpoint(frameSet.getNout()) {
     // Normalize the base and current frame in a way that affects its behavior as a mapping.
@@ -75,39 +75,39 @@ Transform<FromEndpoint, ToEndpoint>::Transform(ast::FrameSet const &frameSet, bo
     _frameSet = frameSetCopy;
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 typename ToEndpoint::Point Transform<FromEndpoint, ToEndpoint>::tranForward(
         typename FromEndpoint::Point const &point) const {
     auto const rawFromData = _fromEndpoint.dataFromPoint(point);
     auto rawToData = _frameSet->tranForward(rawFromData);
     return _toEndpoint.pointFromData(rawToData);
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 typename ToEndpoint::Array Transform<FromEndpoint, ToEndpoint>::tranForward(
         typename FromEndpoint::Array const &array) const {
     auto const rawFromData = _fromEndpoint.dataFromArray(array);
     auto rawToData = _frameSet->tranForward(rawFromData);
     return _toEndpoint.arrayFromData(rawToData);
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 typename FromEndpoint::Point Transform<FromEndpoint, ToEndpoint>::tranInverse(
         typename ToEndpoint::Point const &point) const {
     auto const rawFromData = _toEndpoint.dataFromPoint(point);
     auto rawToData = _frameSet->tranInverse(rawFromData);
     return _fromEndpoint.pointFromData(rawToData);
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 typename FromEndpoint::Array Transform<FromEndpoint, ToEndpoint>::tranInverse(
         typename ToEndpoint::Array const &array) const {
     auto const rawFromData = _toEndpoint.dataFromArray(array);
     auto rawToData = _frameSet->tranInverse(rawFromData);
     return _fromEndpoint.arrayFromData(rawToData);
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 Transform<ToEndpoint, FromEndpoint> Transform<FromEndpoint, ToEndpoint>::getInverse() const {
     auto inverse = std::dynamic_pointer_cast<ast::FrameSet>(_frameSet->getInverse());
     if (!inverse) {
@@ -119,7 +119,7 @@ Transform<ToEndpoint, FromEndpoint> Transform<FromEndpoint, ToEndpoint>::getInve
     return Transform<ToEndpoint, FromEndpoint>(*inverse);
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
 Eigen::MatrixXd Transform<FromEndpoint, ToEndpoint>::getJacobian(FromPoint const &x) const {
     try {
         int const nIn = _fromEndpoint.getNAxes();
@@ -138,16 +138,36 @@ Eigen::MatrixXd Transform<FromEndpoint, ToEndpoint>::getJacobian(FromPoint const
     }
 }
 
-template <typename FromEndpoint, typename ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
+template <class FirstFromEndpoint>
+Transform<FirstFromEndpoint, ToEndpoint> Transform<FromEndpoint, ToEndpoint>::of(
+        Transform<FirstFromEndpoint, FromEndpoint> const &first) const {
+    if (_fromEndpoint.getNAxes() == first.getToEndpoint().getNAxes()) {
+        return Transform<FirstFromEndpoint, ToEndpoint>(*ast::prepend(*_frameSet, *(first.getFrameSet())));
+    } else {
+        auto message = "Cannot match " + std::to_string(first.getToEndpoint().getNAxes()) +
+                       "-D to-endpoint to " + std::to_string(_fromEndpoint.getNAxes()) + "-D from-endpoint.";
+        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, message);
+    }
+}
+
+template <class FromEndpoint, class ToEndpoint>
 std::ostream &operator<<(std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform) {
     auto const frameSet = transform.getFrameSet();
     os << "Transform<" << transform.getFromEndpoint() << ", " << transform.getToEndpoint() << ">";
     return os;
 };
 
-#define INSTANTIATE_TRANSFORM(FromEndpoint, ToEndpoint)          \
-    template class Transform<FromEndpoint, ToEndpoint>;          \
-    template std::ostream &operator<<<FromEndpoint, ToEndpoint>( \
+#define INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, ExtraEndpoint)                                    \
+    template Transform<ExtraEndpoint, ToEndpoint> Transform<FromEndpoint, ToEndpoint>::of<ExtraEndpoint>( \
+            Transform<ExtraEndpoint, FromEndpoint> const &) const;
+#define INSTANTIATE_TRANSFORM(FromEndpoint, ToEndpoint)                  \
+    template class Transform<FromEndpoint, ToEndpoint>;                  \
+    INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, GenericEndpoint)     \
+    INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, Point2Endpoint)      \
+    INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, Point3Endpoint)      \
+    INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, SpherePointEndpoint) \
+    template std::ostream &operator<<<FromEndpoint, ToEndpoint>(         \
             std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform);
 
 // explicit instantiations

tests/test_transform.cc

@@ -76,7 +76,7 @@ ast::PolyMap makeForwardPolyMap(size_t nIn, size_t nOut) {
 }
 
 /**
- * Tests whether the result of SpherePoint::getJacobian(FromPoint const&)
+ * Tests whether the result of Transform::getJacobian(FromPoint const&)
  * has the specified dimensions.
  *
  * The Python version of this method follows a slightly different spec to