Add TransformBoundedField

include/lsst/afw/math/TransformBoundedField.h

@@ -0,0 +1,94 @@
+// -*- LSST-C++ -*-
+/*
+ * 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/>.
+ */
+
+#ifndef LSST_AFW_MATH_TransformBoundedField_h_INCLUDED
+#define LSST_AFW_MATH_TransformBoundedField_h_INCLUDED
+
+#include "ndarray.h"
+
+#include "lsst/pex/config.h"
+#include "lsst/afw/geom/Endpoint.h"
+#include "lsst/afw/geom/Transform.h"
+#include "lsst/afw/math/BoundedField.h"
+
+namespace lsst {
+namespace afw {
+namespace math {
+
+/**
+ * A BoundedField based on geom::Transform<Poin2Endpoint, GenericEndpoint<1>>.
+ *
+ * TransformBoundedField supports arbitrary transforms.
+ */
+class TransformBoundedField : public table::io::PersistableFacade<TransformBoundedField>,
+                              public BoundedField {
+public:
+    using Transform = geom::Transform<geom::Point2Endpoint, geom::GenericEndpoint>;
+
+    /**
+     *  Create a TransformBoundedField from a bounding box and transform.
+     */
+    TransformBoundedField(geom::Box2I const& bbox, Transform const& transform);
+
+    /// Get the contained Transform
+    Transform getTransform() const { return _transform; }    
+
+    /// @copydoc BoundedField::evaluate
+    virtual double evaluate(geom::Point2D const & position) const override;
+
+    /// @copydoc BoundedField::evaluate
+    virtual ndarray::Array<double, 1, 1> evaluate(ndarray::Array<double const, 1> const& x,
+                                                  ndarray::Array<double const, 1> const& y) const override;
+
+    using BoundedField::evaluate;
+
+    /// TransformBoundedField is always persistable.
+    virtual bool isPersistable() const override { return true; }
+
+    /// @copydoc BoundedField::operator*
+    virtual std::shared_ptr<BoundedField> operator*(double const scale) const;
+
+    /// @copydoc BoundedField::operator==
+    virtual bool operator==(BoundedField const& rhs) const;
+
+protected:
+    virtual std::string getPersistenceName() const;
+
+    virtual std::string getPythonModule() const;
+
+    virtual void write(OutputArchiveHandle & handle) const;
+
+private:
+    // Internal constructor for fit() routines: just initializes the transform,
+    // leaves coefficients empty.
+    explicit TransformBoundedField(afw::geom::Box2I const & bbox);
+
+    virtual std::string toString() const;
+
+    Transform _transform;
+};
+}
+}
+}  // namespace lsst::afw::math
+
+#endif // !LSST_AFW_MATH_TransformBoundedField_h_INCLUDED

python/lsst/afw/math/SConscript

@@ -7,6 +7,7 @@ scripts.BasicSConscript.pybind11(['function',
                                   'spatialCell',
                                   'boundedField',
                                   'chebyshevBoundedField',
+                                  'transformBoundedField',
                                   'leastSquares',
                                   'random',
                                   'convolveImage',

python/lsst/afw/math/mathLib.py

@@ -13,6 +13,7 @@
 from ._boundedField import *
 from ._chebyshevBoundedField import *
 from .chebyshevBoundedFieldConfig import ChebyshevBoundedFieldConfig
+from ._transformBoundedField import *
 from ._leastSquares import *
 from ._random import *
 from ._convolveImage import *

python/lsst/afw/math/transformBoundedField.cc

@@ -0,0 +1,85 @@
+/*
+ * LSST Data Management System
+ * Copyright 2017  AURA/LSST.
+ *
+ * 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 <https://www.lsstcorp.org/LegalNotices/>.
+ */
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "numpy/arrayobject.h"
+#include "ndarray/pybind11.h"
+#include "ndarray/converter.h"
+
+#include "lsst/pex/config/python.h"  // defines LSST_DECLARE_CONTROL_FIELD
+#include "lsst/afw/table/io/python.h"
+
+#include "lsst/afw/table/io/Persistable.h"
+#include "lsst/afw/math/BoundedField.h"
+#include "lsst/afw/math/TransformBoundedField.h"
+
+namespace py = pybind11;
+using namespace py::literals;
+
+namespace lsst {
+namespace afw {
+namespace math {
+namespace {
+
+using ClsField = py::class_<TransformBoundedField, std::shared_ptr<TransformBoundedField>, BoundedField,
+                            lsst::afw::table::io::PersistableFacade<TransformBoundedField>>;
+
+PYBIND11_PLUGIN(_transformBoundedField) {
+    py::module mod("_transformBoundedField", "Python wrapper for afw _transformBoundedField library");
+
+    if (_import_array() < 0) {
+        PyErr_SetString(PyExc_ImportError, "numpy.core.multiarray failed to import");
+        return nullptr;
+    };
+
+    /* Module level */
+    lsst::afw::table::io::python::declarePersistableFacade<TransformBoundedField>(mod,
+                                                                                  "TransformBoundedField");
+    ClsField cls(mod, "TransformBoundedField");
+
+    cls.def(py::init<lsst::afw::geom::Box2I const &, TransformBoundedField::Transform const &>(), "bbox"_a,
+            "transform"_a);
+
+    cls.def("__rmul__", [](TransformBoundedField &bf, double const scale) { return bf * scale; },
+            py::is_operator());
+    cls.def("__mul__", [](TransformBoundedField &bf, double const scale) { return bf * scale; },
+            py::is_operator());
+    cls.def("__eq__", &TransformBoundedField::operator==, py::is_operator());
+
+    cls.def("getTransform", &TransformBoundedField::getTransform);
+    cls.def("evaluate", (double (BoundedField::*)(double, double) const) & BoundedField::evaluate);
+    cls.def("evaluate",
+            (ndarray::Array<double, 1, 1>(TransformBoundedField::*)(
+                    ndarray::Array<double const, 1> const &, ndarray::Array<double const, 1> const &) const) &
+                    TransformBoundedField::evaluate);
+    cls.def("evaluate", (double (TransformBoundedField::*)(lsst::afw::geom::Point2D const &) const) &
+                                TransformBoundedField::evaluate);
+
+    return mod.ptr();
+}
+
+}  // namespace
+}  // namespace math
+}  // namespace afw
+}  // namespace lsst

src/math/TransformBoundedField.cc

@@ -0,0 +1,172 @@
+// -*- LSST-C++ -*-
+/*
+ * LSST Data Management System
+ * Copyright 2008-2014 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/>.
+ */
+
+#include <memory>
+#include <string>
+
+#include "ndarray/eigen.h"
+#include "astshim.h"
+#include "lsst/afw/math/TransformBoundedField.h"
+#include "lsst/afw/table/io/InputArchive.h"
+#include "lsst/afw/table/io/OutputArchive.h"
+#include "lsst/afw/table/io/CatalogVector.h"
+#include "lsst/afw/table/aggregates.h"
+
+namespace lsst {
+namespace afw {
+namespace math {
+
+TransformBoundedField::TransformBoundedField(
+        geom::Box2I const& bbox,
+        geom::Transform<geom::Point2Endpoint, geom::GenericEndpoint> const& transform)
+        : BoundedField(bbox), _transform(transform) {
+    if (transform.getToEndpoint().getNAxes() != 1) {
+        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
+                          "To EndPoint of transform must have 1 axis but has " +
+                                  std::to_string(transform.getToEndpoint().getNAxes()));
+    }
+}
+
+double TransformBoundedField::evaluate(geom::Point2D const& position) const {
+    return _transform.applyForward(position)[0];
+}
+
+ndarray::Array<double, 1, 1> TransformBoundedField::evaluate(ndarray::Array<double const, 1> const& x,
+                                                             ndarray::Array<double const, 1> const& y) const {
+    if (x.getSize<0>() != y.getSize<0>()) {
+        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
+                          "x length " + std::to_string(x.getSize<0>()) + "!= y length " +
+                                  std::to_string(x.getSize<0>()));
+    }
+
+    // TODO if Mapping.applyForward gains support for x, y then use that instead of copying data
+    int const nPoints = x.getSize<0>();
+    ndarray::Array<double, 2, 2> xy = ndarray::allocate(ndarray::makeVector(2, nPoints));
+    for (int col = 0; col < nPoints; ++col) {
+        xy[0][col] = x[col];
+        xy[1][col] = y[col];
+    }
+
+    auto res2D = _transform.getFrameSet()->applyForward(xy);
+
+    // res2D has shape 1 x N; return a 1-D view with the extra dimension stripped
+    auto resShape = ndarray::makeVector(nPoints);
+    auto resStrides = ndarray::makeVector(1);
+    return ndarray::external(res2D.getData(), resShape, resStrides, res2D);
+}
+
+// ------------------ persistence ---------------------------------------------------------------------------
+
+namespace {
+
+struct PersistenceHelper {
+    table::Schema schema;
+    table::Key<int> orderX;
+    table::PointKey<int> bboxMin;
+    table::PointKey<int> bboxMax;
+    table::Key<std::string> frameSet;
+
+    PersistenceHelper()
+            : schema(),
+              bboxMin(table::PointKey<int>::addFields(schema, "bbox_min", "lower-left corner of bounding box",
+                                                      "pixel")),
+              bboxMax(table::PointKey<int>::addFields(schema, "bbox_max",
+                                                      "upper-right corner of bounding box", "pixel")),
+              frameSet(schema.addField<std::string>("frameSet", "FrameSet contained in the Transform", 0,
+                                                    false)) {}
+
+    PersistenceHelper(table::Schema const& s)
+            : schema(s), bboxMin(s["bbox_min"]), bboxMax(s["bbox_max"]), frameSet(s["frameSet"]) {}
+};
+
+class TransformBoundedFieldFactory : public table::io::PersistableFactory {
+public:
+    virtual std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
+                                                         CatalogVector const& catalogs) const {
+        LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
+        LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
+        table::BaseRecord const& record = catalogs.front().front();
+        PersistenceHelper const keys(record.getSchema());
+        // NOTE: needed invert=false in case min=-1, max=0 (empty bbox). See RFC-324 and DM-10200
+        geom::Box2I bbox(record.get(keys.bboxMin), record.get(keys.bboxMax), false);
+        auto frameSetStr = record.get(keys.frameSet);
+        auto stream = ast::StringStream(frameSetStr);
+        auto astObjectPtr = ast::Channel(stream).read();
+        auto frameSet = std::dynamic_pointer_cast<ast::FrameSet>(astObjectPtr);
+        if (!frameSet) {
+            throw LSST_EXCEPT(pex::exceptions::RuntimeError,
+                              "The contents of the frameSet field could not be parsed as an ast::FrameSet");
+        }
+        auto transform = geom::Transform<geom::Point2Endpoint, geom::GenericEndpoint>(*frameSet);
+        return std::make_shared<TransformBoundedField>(bbox, transform);
+    }
+
+    TransformBoundedFieldFactory(std::string const& name) : afw::table::io::PersistableFactory(name) {}
+};
+
+std::string getTransformBoundedFieldPersistenceName() { return "TransformBoundedField"; }
+
+TransformBoundedFieldFactory registration(getTransformBoundedFieldPersistenceName());
+
+}  // namespace
+
+std::string TransformBoundedField::getPersistenceName() const {
+    return getTransformBoundedFieldPersistenceName();
+}
+
+std::string TransformBoundedField::getPythonModule() const { return "lsst.afw.math"; }
+
+void TransformBoundedField::write(OutputArchiveHandle& handle) const {
+    PersistenceHelper const keys;
+    table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
+    std::shared_ptr<table::BaseRecord> record = catalog.addNew();
+    record->set(keys.bboxMin, getBBox().getMin());
+    record->set(keys.bboxMax, getBBox().getMax());
+    record->set(keys.frameSet, getTransform().getFrameSet()->show());
+    handle.saveCatalog(catalog);
+}
+
+PTR(BoundedField) TransformBoundedField::operator*(double const scale) const {
+    auto zoomMap = ast::ZoomMap(1, scale);
+    auto newMapping = getTransform().getFrameSet()->then(zoomMap);
+    auto newTransform = Transform(newMapping);
+    return std::make_shared<TransformBoundedField>(getBBox(), newTransform);
+}
+
+bool TransformBoundedField::operator==(BoundedField const& rhs) const {
+    auto rhsCasted = dynamic_cast<TransformBoundedField const*>(&rhs);
+    if (!rhsCasted) return false;
+
+    return getBBox() == rhsCasted->getBBox() &&
+           *(getTransform().getFrameSet()) == *(rhsCasted->getTransform().getFrameSet());
+}
+
+std::string TransformBoundedField::toString() const {
+    std::ostringstream os;
+    os << "TransformBoundedField (containing " << _transform.getFrameSet()->getNFrame() << " frames)";
+    return os.str();
+}
+
+}  // namespace math
+}  // namespace afw
+}  // namespace lsst