Add function for linearizing a Transform.

The function has been placed in a new transformFactory header that
will soon have other Transform-creating functions.

include/lsst/afw/geom.h

@@ -46,5 +46,6 @@
 #include "lsst/afw/geom/TransformMap.h"
 #include "lsst/afw/geom/Endpoint.h"
 #include "lsst/afw/geom/Transform.h"
+#include "lsst/afw/geom/transformFactory.h"
 
 #endif  // LSST_GEOM_H

include/lsst/afw/geom/transformFactory.h

@@ -0,0 +1,77 @@
+// -*- lsst-c++ -*-
+
+#ifndef LSST_AFW_GEOM_TRANSFORMFACTORY_H
+#define LSST_AFW_GEOM_TRANSFORMFACTORY_H
+
+/*
+ * LSST Data Management System
+ * Copyright 2008-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/>.
+ */
+
+/*
+ * Functions for producing Transforms with commonly desired properties.
+ */
+
+#include "lsst/afw/geom/Transform.h"
+
+namespace lsst {
+namespace afw {
+namespace geom {
+
+/**
+ * Approximate a Transform by its local linearization.
+ *
+ * @tparam FromEndpoint, ToEndpoint The endpoints of the transform.
+ *
+ * @param original the Transform to linearize
+ * @param point the point at which a linear approximation is desired
+ * @returns a linear Transform whose value and Jacobian at `point` match those
+ *         of `original`.
+ *
+ * @throws pex::exceptions::InvalidParameterError Thrown if `original` does not
+ *             have a well-defined value and Jacobian at `point`
+ * @exceptsafe Provides basic exception safety.
+ */
+template <class FromEndpoint, class ToEndpoint>
+Transform<FromEndpoint, ToEndpoint> linearizeTransform(
+        Transform<FromEndpoint, ToEndpoint> const &original,
+        typename Transform<FromEndpoint, ToEndpoint>::FromPoint const &point);
+
+/*
+ * The correct behavior for linearization is unclear where SpherePoints are involved (see discussion on
+ * DM-10542). Forbid usage until somebody needs it. Note to maintainers: the template specializations MUST
+ * be deleted in the header for compilers to complain correctly.
+ */
+#define DISABLE(From, To)                                                         \
+    template <>                                                                   \
+    Transform<From, To> linearizeTransform<From, To>(Transform<From, To> const &, \
+                                                     Transform<From, To>::FromPoint const &) = delete;
+DISABLE(GenericEndpoint, SpherePointEndpoint);
+DISABLE(Point2Endpoint, SpherePointEndpoint);
+DISABLE(SpherePointEndpoint, GenericEndpoint);
+DISABLE(SpherePointEndpoint, Point2Endpoint);
+DISABLE(SpherePointEndpoint, SpherePointEndpoint);
+#undef DISABLE
+
+}  // geom
+}  // afw
+}  // lsst
+
+#endif  // LSST_AFW_GEOM_TRANSFORMFACTORY_H

python/lsst/afw/geom/SConscript

@@ -16,6 +16,7 @@ scripts.BasicSConscript.pybind11(
         'xyTransform',
         'endpoint',
         'transform/transform',
+        'transformFactory',
         'skyWcs/skyWcs',
     ],
     addUnderscore=False

python/lsst/afw/geom/__init__.py

@@ -42,5 +42,6 @@
 from .utils import *
 from .endpoint import *
 from .transform import *
+from .transformFactory import *
 from .skyWcs import *
 from .readTransform import *

python/lsst/afw/geom/transformFactory.cc

@@ -0,0 +1,70 @@
+/*
+* LSST Data Management System
+* See COPYRIGHT file at the top of the source tree.
+*
+* 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/>.
+*/
+#include "pybind11/pybind11.h"
+#include "pybind11/stl.h"  // Needed because some Endpoints are vectors
+
+#include "numpy/arrayobject.h"
+#include "ndarray/pybind11.h"
+
+#include "lsst/afw/geom/Endpoint.h"
+#include "lsst/afw/geom/transformFactory.h"
+
+namespace py = pybind11;
+using namespace py::literals;
+
+namespace lsst {
+namespace afw {
+namespace geom {
+namespace {
+
+// Declare methods on Transform<FromEndpoint, ToEndpoint>
+template <class FromEndpoint, class ToEndpoint>
+void declareFunctionTemplates(py::module &mod) {
+    using Class = Transform<FromEndpoint, ToEndpoint>;
+
+    mod.def("linearizeTransform",
+            (Class(*)(Class const &, typename Class::FromPoint const &)) & linearizeTransform);
+}
+
+PYBIND11_PLUGIN(transformFactory) {
+    py::module mod("transformFactory");
+
+    py::module::import("lsst.afw.geom.transform");
+
+    // Need to import numpy for ndarray and eigen conversions
+    if (_import_array() < 0) {
+        PyErr_SetString(PyExc_ImportError, "numpy.core.multiarray failed to import");
+        return nullptr;
+    }
+
+    declareFunctionTemplates<GenericEndpoint, GenericEndpoint>(mod);
+    declareFunctionTemplates<GenericEndpoint, Point2Endpoint>(mod);
+    declareFunctionTemplates<Point2Endpoint, GenericEndpoint>(mod);
+    declareFunctionTemplates<Point2Endpoint, Point2Endpoint>(mod);
+
+    return mod.ptr();
+}
+
+}  // <anonymous>
+}  // geom
+}  // afw
+}  // lsst

src/geom/transformFactory.cc

@@ -0,0 +1,125 @@
+// -*- LSST-C++ -*-
+/*
+ * LSST Data Management System
+ * See COPYRIGHT file at the top of the source tree.
+ *
+ * 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/>.
+ */
+
+#include <sstream>
+#include <cmath>
+
+#include "astshim.h"
+
+#include "Eigen/Core"
+
+#include "lsst/afw/geom/Endpoint.h"
+#include "lsst/afw/geom/transformFactory.h"
+#include "lsst/pex/exceptions.h"
+
+#include "ndarray.h"
+#include "ndarray/eigen.h"
+
+namespace lsst {
+namespace afw {
+namespace geom {
+
+namespace {
+/**
+ * Print a vector to a stream.
+ *
+ * The exact details of the representation are unspecified and subject to
+ * change, but the following may be regarded as typical:
+ *
+ *     [1.0, -3.560, 42.0]
+ *
+ * @tparam T the element type. Must support stream output.
+ */
+template <typename T>
+std::ostream &operator<<(std::ostream &os, std::vector<T> const &v) {
+    os << '[';
+    bool first = true;
+    for (T element : v) {
+        if (first) {
+            first = false;
+        } else {
+            os << ", ";
+        }
+        os << element;
+    }
+    os << ']';
+    return os;
+}
+
+/**
+ * Convert a Matrix to the equivalent ndarray.
+ *
+ * @param matrix The matrix to convert.
+ * @returns an ndarray containing a copy of the data in `matrix`
+ */
+template <class T, int Rows, int Cols>
+ndarray::Array<T, 2, 2> toNdArray(Eigen::Matrix<T, Rows, Cols> const &matrix) {
+    ndarray::Array<T, 2, 2> array = ndarray::allocate(ndarray::makeVector(matrix.rows(), matrix.cols()));
+    array.asEigen() = matrix;
+    return array;
+}
+}  // namespace {}
+
+template <class From, class To>
+Transform<From, To> linearizeTransform(Transform<From, To> const &original,
+                                       typename Transform<From, To>::FromPoint const &inPoint) {
+    using Transform = Transform<From, To>;
+    auto fromEndpoint = original.getFromEndpoint();
+    auto toEndpoint = original.getToEndpoint();
+
+    auto outPoint = original.applyForward(inPoint);
+    auto jacobian = original.getJacobian(inPoint);
+    for (int i = 0; i < toEndpoint.getNAxes(); ++i) {
+        if (!std::isfinite(outPoint[i])) {
+            std::ostringstream buffer;
+            buffer << "Transform ill-defined: " << inPoint << " -> " << outPoint;
+            throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, buffer.str());
+        }
+    }
+    if (!jacobian.allFinite()) {
+        std::ostringstream buffer;
+        buffer << "Transform not continuous at " << inPoint << ": J = " << jacobian;
+        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, buffer.str());
+    }
+
+    // y(x) = J (x - x0) + y0
+    auto map = ast::ShiftMap(fromEndpoint.dataFromPoint(inPoint))
+                       .getInverse()
+                       ->then(ast::MatrixMap(toNdArray(jacobian)))
+                       .then(ast::ShiftMap(toEndpoint.dataFromPoint(outPoint)));
+    // TODO: remove false flag as part of DM-10947
+    return Transform(map, false);
+}
+
+#define INSTANTIATE_FACTORIES(From, To)                                                             \
+    template Transform<From, To> linearizeTransform<From, To>(Transform<From, To> const &transform, \
+                                                              Transform<From, To>::FromPoint const &point);
+
+// explicit instantiations
+INSTANTIATE_FACTORIES(GenericEndpoint, GenericEndpoint);
+INSTANTIATE_FACTORIES(GenericEndpoint, Point2Endpoint);
+INSTANTIATE_FACTORIES(Point2Endpoint, GenericEndpoint);
+INSTANTIATE_FACTORIES(Point2Endpoint, Point2Endpoint);
+}
+}
+} /* namespace lsst::afw::geom */