diff --git a/include/astshim.h b/include/astshim.h index a228d5f..a7d95d4 100644 --- a/include/astshim.h +++ b/include/astshim.h @@ -41,6 +41,7 @@ #include "astshim/ShiftMap.h" #include "astshim/SlaMap.h" #include "astshim/SphMap.h" +#include "astshim/SplineMap.h" #include "astshim/TimeMap.h" #include "astshim/TranMap.h" #include "astshim/UnitMap.h" diff --git a/include/astshim/SplineMap.h b/include/astshim/SplineMap.h new file mode 100644 index 0000000..588961d --- /dev/null +++ b/include/astshim/SplineMap.h @@ -0,0 +1,118 @@ +/* + * This file is part of astshim. + * + * Developed for the LSST Data Management System. + * This product includes software developed by the LSST Project + * (https://www.lsst.org). + * See the COPYRIGHT file at the top-level directory of this distribution + * for details of code ownership. + * + * 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 GNU General Public License + * along with this program. If not, see . + */ +#ifndef ASTSHIM_SPLINEMAP_H +#define ASTSHIM_SPLINEMAP_H + +#include +#include + +#include "astshim/base.h" +#include "astshim/Mapping.h" + +namespace ast { + +/** +SplineMap is a @ref Mapping which performs a 2-in, 2-out Mapping defined by a pair of two- +dimensional B-spline surfaces. The inverse transformation is fit iteratively, and requires that +the output axes are a small (i.e. much less than the range of the input values) perturbation of the +input axes. +### Attributes + +All those of @ref Mapping plus: + +- "InvNiter": Maximum number of iterations for iterative inverse. +- "InvTol": Target relative error for iterative inverse. +- "OutUnit": Determines how out-of-bounds input positions are handled. If +true, input points outside the bounding box are copied unchanged to output. If false, the output +for these points will be nan. +*/ +class SplineMap : public Mapping { + friend class Object; + +public: + /** + Construct a @ref SplineMap with the forward transform specified by the supplied b-splines and + the inverse transform fit iteratively. + + @param[in] kx Polynomial order in X input direction. + @param[in] ky Polynomial order in Y input direction. + @param[in] nx Number of coefficients in X input direction. + @param[in] ny Number of coefficients in Y input direction + @param[in] tx Array of knot positions for X input. + @param[in] ty Array of knot positions for Y input. + @param[in] cu Array of coefficients defining U output. + @param[in] cv Array of coefficients defining V output. + @param[in] options Comma-separated list of attribute assignments. + + @throws std::invalid_argument if the polynomial order of number of coefficients in either + direction are negative. + @throws std::invalid_argument if the polynomial order plus the number of coefficients does + not equal the number of knot positions for a given direction. + @throws std::invalid_argument if the product of the number of coefficients in each direction + does not equal the number of coefficients provided. + + The coefficients cu and cv are 1-dimensional arrays. + */ + explicit SplineMap(int kx, int ky, int nx, int ny, std::vector const &tx, std::vector const &ty, + std::vector const &cu, std::vector const &cv, std::string const &options = "") + : Mapping(reinterpret_cast( + _makeRawSplineMap(kx, ky, nx, ny, tx, ty, cu, cv, options))) {} + + ~SplineMap() override = default; + + /// Copy constructor: make a deep copy + SplineMap(SplineMap const &) = default; + SplineMap(SplineMap &&) = default; + SplineMap &operator=(SplineMap const &) = delete; + SplineMap &operator=(SplineMap &&) = default; + + /// Return a deep copy of this object. + std::shared_ptr copy() const { return std::static_pointer_cast(copyPolymorphic()); } + + /// Get "InvNiter": Get maximum number of iterations for iterative inverse. + int getInvNiter() const { return getI("InvNiter"); } + + /// Get "OutUnit": Is the out-of-bounds position transform an identity map? + int getOutUnit() const { return getI("OutUnit"); } + + /// Get "InvTol": Get target relative error for iterative inverse. + double getInvTol() const { return getD("InvTol"); } + + +protected: + std::shared_ptr copyPolymorphic() const override { + return copyImpl(); + } + + /// Construct a SplineMap from an raw AST pointer + explicit SplineMap(AstSplineMap *map); + +private: + /// Make a raw AstSplineMap with a specified forward transform and an iterative inverse. + AstSplineMap *_makeRawSplineMap(int kx, int ky, int nx, int ny, std::vector const &tx, std::vector const &ty, + std::vector const &cu, std::vector const &cv, std::string const &options = "") const; +}; + +} // namespace ast + +#endif diff --git a/python/astshim/SConscript b/python/astshim/SConscript index 31d135b..2e656de 100755 --- a/python/astshim/SConscript +++ b/python/astshim/SConscript @@ -33,6 +33,7 @@ scripts.BasicSConscript.python(['_astshimLib'], [ "slaMap.cc", "specFrame.cc", "sphMap.cc", + "splineMap.cc", "stream.cc", "table.cc", "timeFrame.cc", diff --git a/python/astshim/_astshimLib.cc b/python/astshim/_astshimLib.cc index f9cd2d6..7a99b62 100644 --- a/python/astshim/_astshimLib.cc +++ b/python/astshim/_astshimLib.cc @@ -61,6 +61,7 @@ void wrapSkyFrame(WrapperCollection&); void wrapSlaMap(WrapperCollection&); void wrapSpecFrame(WrapperCollection&); void wrapSphMap(WrapperCollection&); +void wrapSplineMap(WrapperCollection&); void wrapStream(WrapperCollection&); void wrapTable(WrapperCollection&); void wrapTimeFrame(WrapperCollection&); @@ -110,6 +111,7 @@ PYBIND11_MODULE(_astshimLib, mod) { wrapSlaMap(wrappers); wrapSpecFrame(wrappers); wrapSphMap(wrappers); + wrapSplineMap(wrappers); wrapStream(wrappers); wrapTimeFrame(wrappers); wrapTimeMap(wrappers); diff --git a/python/astshim/splineMap.cc b/python/astshim/splineMap.cc new file mode 100644 index 0000000..0d84a27 --- /dev/null +++ b/python/astshim/splineMap.cc @@ -0,0 +1,52 @@ +/* + * This file is part of astshim. + * + * Developed for the LSST Data Management System. + * This product includes software developed by the LSST Project + * (https://www.lsst.org). + * See the COPYRIGHT file at the top-level directory of this distribution + * for details of code ownership. + * + * 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 GNU General Public License + * along with this program. If not, see . + */ +#include +#include + +#include +#include +#include "lsst/cpputils/python.h" + +#include "astshim/Mapping.h" +#include "astshim/SplineMap.h" + +namespace py = pybind11; +using namespace pybind11::literals; + +namespace ast { +void wrapSplineMap(lsst::cpputils::python::WrapperCollection &wrappers) { + using PySplineMap=py::classh; + wrappers.wrapType(PySplineMap(wrappers.module, "SplineMap"), [](auto &mod, auto &cls) { + + cls.def(py::init const &, std::vector const &, + std::vector const &, std::vector const &, std::string const &>(), "kx"_a, + "ky"_a, "nx"_a, "ny"_a, "tx"_a, "ty"_a, "cu"_a, "cv"_a, "options"_a = ""); + cls.def(py::init()); + cls.def("copy", &SplineMap::copy); + cls.def_property_readonly("invNIter", &SplineMap::getInvNiter); + cls.def_property_readonly("outUnit", &SplineMap::getOutUnit); + cls.def_property_readonly("invTol", &SplineMap::getInvTol); + }); +} + +} // namespace ast diff --git a/src/Object.cc b/src/Object.cc index f7da7b3..58b5642 100644 --- a/src/Object.cc +++ b/src/Object.cc @@ -54,6 +54,7 @@ #include "astshim/SlaMap.h" #include "astshim/SpecFrame.h" #include "astshim/SphMap.h" +#include "astshim/SplineMap.h" #include "astshim/Table.h" #include "astshim/TimeFrame.h" #include "astshim/TimeMap.h" @@ -114,6 +115,7 @@ std::shared_ptr Object::_basicFromAstObject(AstObject *rawObj) { {"SlaMap", makeShim}, {"SpecFrame", makeShim}, {"SphMap", makeShim}, + {"SplineMap", makeShim}, {"TimeFrame", makeShim}, {"Table", makeShim}, {"TimeMap", makeShim}, diff --git a/src/SplineMap.cc b/src/SplineMap.cc new file mode 100644 index 0000000..5984228 --- /dev/null +++ b/src/SplineMap.cc @@ -0,0 +1,64 @@ +/* + * This file is part of astshim. + * + * Developed for the LSST Data Management System. + * This product includes software developed by the LSST Project + * (https://www.lsst.org). + * See the COPYRIGHT file at the top-level directory of this distribution + * for details of code ownership. + * + * 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 GNU General Public License + * along with this program. If not, see . + */ +#include +#include + +#include "astshim/SplineMap.h" + +namespace ast { + +SplineMap::SplineMap(AstSplineMap *map) : Mapping(reinterpret_cast(map)) { + if (!astIsASplineMap(getRawPtr())) { + std::ostringstream os; + os << "this is a " << getClassName() << ", which is not a SplineMap"; + throw std::invalid_argument(os.str()); + } +} + +/// Make a raw AstSplineMap with a specified forward transform. +AstSplineMap *SplineMap::_makeRawSplineMap(int kx, int ky, int nx, int ny, std::vector const &tx, + std::vector const &ty, std::vector const &cu, std::vector const &cv, + std::string const &options) const { + + const size_t nTx = tx.size(); + const size_t nTy = ty.size(); + const size_t nCoeffsU = cu.size(); + const size_t nCoeffsV = cv.size(); + + if ((kx < 0) || (ky < 0) || (nx < 0) || (ny < 0)) { + throw std::invalid_argument("The polynomial order and number of coefficients must not be negative."); + } + if (((size_t)(kx + nx) != nTx) || ((size_t)(ky + ny) != nTy)) { + throw std::invalid_argument("The length of the knot positions must equal the polynomial order plus " + "the number of coefficients."); + } + if (((size_t)(nx * ny) != nCoeffsU) || ((size_t)(nx * ny) != nCoeffsV)) { + throw std::invalid_argument("The length of the coefficients must equal product of the arguments nx " + "and ny."); + } + + return reinterpret_cast(astSplineMap(kx, ky, nx, ny, tx.data(), ty.data(), cu.data(), + cv.data(), options.c_str())); +} + +} // namespace ast diff --git a/tests/test_splineMap.py b/tests/test_splineMap.py new file mode 100644 index 0000000..e1865a1 --- /dev/null +++ b/tests/test_splineMap.py @@ -0,0 +1,99 @@ +import unittest + +import numpy as np +import numpy.testing as npt +from scipy.interpolate import RectBivariateSpline + +import astshim as ast +from astshim.test import MappingTestCase + + +class TestSplineMap(MappingTestCase): + + def setUp(self): + self.nx = 150 + self.ny = 100 + self.k = 4 + a = 4.0 + b = 0.2 + c = 3.0 + + self.xs = np.linspace(1, self.nx, self.nx) + self.ys = np.linspace(1, self.ny, self.ny) + + fcx = np.zeros((self.nx, self.ny)) + fcy = np.zeros((self.nx, self.ny)) + for j in range(self.ny): + y = self.ys[j] + for i in range(self.nx): + x = self.xs[i] + fcx[i][j] = x + a * np.sin(b * x) * np.cos(b * y) + fcy[i][j] = y + c * np.cos(b * x) * np.sin(b * y) + + self.splinex = RectBivariateSpline(self.xs, self.ys, fcx, s=0) + (self.tx, self.ty) = self.splinex.get_knots() + self.arx = self.splinex.get_coeffs() + + self.spliney = RectBivariateSpline(self.xs, self.ys, fcy, s=0) + self.ary = self.spliney.get_coeffs() + + self.splineMap = ast.SplineMap( + self.k, self.k, self.nx, self.ny, self.tx, self.ty, self.arx, self.ary + ) + + def test_SplineMap(self): + """Test that the forward and inverse transforms match + scipy.interpolate.RectBivariateSpline. + """ + xval = np.array([0.0, 12.5, 12.0, 1.0, 15.0, 1.0, 95.0, 149.5, 151.2, 77.77]) + yval = np.array([-1.0, 8.8, 8.0, 1.0, 15.0, 76.0, 100.0, 99.8, 82.3, 54.3]) + + u = self.splinex.ev(xval, yval) + v = self.spliney.ev(xval, yval) + + outData = self.splineMap.applyForward(np.array([xval, yval])) + + outOfBounds = ( + (xval > self.xs.max()) + | (xval < self.xs.min()) + | (yval < self.ys.min()) + | (yval > self.ys.max()) + ) + + npt.assert_equal(outData[:, outOfBounds], np.nan) + npt.assert_almost_equal(outData[0, ~outOfBounds], u[~outOfBounds]) + npt.assert_almost_equal(outData[1, ~outOfBounds], v[~outOfBounds]) + + reverseTrip = self.splineMap.applyInverse(outData[:, ~outOfBounds].copy()) + + npt.assert_almost_equal(reverseTrip[0], xval[~outOfBounds]) + npt.assert_almost_equal(reverseTrip[1], yval[~outOfBounds]) + + def test_splineMapAttributes(self): + """Check that SplineMap attributes can be accessed and return expected + values.""" + self.assertEqual(self.splineMap.invTol, 1e-6) + self.assertEqual(self.splineMap.invNIter, 6) + self.assertEqual(self.splineMap.outUnit, False) + + # Initialize SplineMap with other attributes and check they are set + # correctly. + options = "InvNIter=5,InvTol=1e-7,OutUnit=1" + splineMap = ast.SplineMap( + self.k, + self.k, + self.nx, + self.ny, + self.tx, + self.ty, + self.arx, + self.ary, + options=options, + ) + self.assertEqual(splineMap.invTol, 1e-7) + self.assertEqual(splineMap.invNIter, 5) + self.assertEqual(splineMap.outUnit, 1) + + +if __name__ == "__main__": + unittest.main()