DM-13065: SipApproximation

include/lsst/afw/geom/SipApproximation.h

@@ -0,0 +1,370 @@
+// -*- lsst-c++ -*-
+/*
+ * 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 <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef LSST_AFW_GEOM_SipApproximation_h_INCLUDED
+#define LSST_AFW_GEOM_SipApproximation_h_INCLUDED
+
+#include <memory>
+#include <vector>
+
+#include "Eigen/Core"
+
+#include "lsst/afw/geom/Transform.h"
+#include "lsst/afw/geom/Box.h"
+#include "lsst/afw/geom/LinearTransform.h"
+
+namespace lsst { namespace afw { namespace geom {
+
+/**
+ *  A fitter and results class for approximating a general Transform in a form compatible with
+ *  FITS WCS persistence.
+ *
+ *  The Simple Imaging Polynomial (SIP) convention (Shupe et al 2005) adds
+ *  forward and reverse polynomial mappings to a standard projection FITS WCS
+ *  projection (e.g. "TAN" for gnomonic) that relate Intermediate World
+ *  Coordinates (see Calabretta & Greisen 2002) to image pixel coordinates.
+ *  The SIP "forward" transform is defined by polynomial coeffients @f$A@f$
+ *  and @f$B@f$ that map pixel coordinates @f$(u, v)@f$ to Intermediate World
+ *  Coordinates @f$(x, y)@f$ via
+ *  @f[
+ *     \boldsymbol{S}\left[\begin{array}{c}
+ *        x \\
+ *        y
+ *     \end{array}\right]
+ *     \equiv
+ *     \left[\begin{array}{c}
+ *        x_s \\
+ *        y_s
+ *     \end{array}\right]
+ *     =
+ *     \left[\begin{array}{c}
+ *        (u - u_0) + \displaystyle\sum_{p,q}^{0 \le p + q \le N} \mathrm{A}_{p,q} (u - u_0)^p (v - v_0)^q \\
+ *        (v - v_0) + \displaystyle\sum_{p,q}^{0 \le p + q \le N} \mathrm{B}_{p,q} (u - u_0)^p (v - v_0)^q
+ *     \end{array}\right]
+ *  @f]
+ *  The reverse transform has essentially the same form:
+ *  @f[
+ *     \left[\begin{array}{c}
+ *        u - u_0 \\
+ *        v - v_0
+ *     \end{array}\right]
+ *     =
+ *     \left[\begin{array}{c}
+ *        x_s + \displaystyle\sum_{p,q}^{0 \le p + q \le N} \mathrm{AP}_{p,q} x_s^p y_s^q \\
+ *        y_s + \displaystyle\sum_{p,q}^{0 \le p + q \le N} \mathrm{BP}_{p,q} x_s^p y_s^q
+ *     \end{array}\right]
+ *  @f]
+ *  In both cases, @f$(u_0, v_0)@f$ is the pixel origin (CRPIX in FITS WCS) and
+ *  @f$\boldsymbol{S}@f$ is the *inverse* of the Jacobian "CD" matrix.  Both
+ *  CRPIX and CD are considered fixed inputs, and we do not attempt to null
+ *  the zeroth- and first-order terms of @f$A@f$ and @f$B@f$ (as some SIP
+ *  fitters do); together, these conventions make solving for the coefficients
+ *  a much simpler linear problem.
+ *
+ *  @note In the implementation, we typically refer to @f$(u-u_0, v-v_0)@f$ as
+ *  @c dpix (for "pixel delta"), and @f$(x_s, y_s)@f$ as @c siwc (for "scaled
+ *  intermediate world coordinates").
+ *
+ *  While LSST WCSs are in general too complex to be described exactly in FITS
+ *  WCS, they can generally be closely approximated by standard FITS WCS
+ *  projection with additional SIP distortions.  This class fits such an
+ *  approximation, given a TransformPoint2ToPoint2 object that represents the
+ *  exact mapping from pixels to Intermediate World Coordinates with a SIP
+ *  distortion.
+ */
+class SipApproximation {
+public:
+
+    /**
+     *  Construct a new approximation by fitting on a grid of points.
+     *
+     *  @param[in]  pixelToIwc      The true Transform to approximate.  Should go
+     *                              from pixels to Intermediate World Coordinates
+     *                              when applyForward is called.
+     *  @param[in]  crpix           Pixel origin, using the LSST 0-indexed
+     *                              convention rather than the FITS 1-indexed
+     *                              convention; equal to (CRPIX1 - 1, CRPIX2 - 1).
+     *  @param[in]  cd              Nominal Jacobian ("CD" in FITS WCS).
+     *  @param[in]  bbox            Pixel-coordinate bounding box over which the
+     *                              approximation should be valid.  Used to construct
+     *                              the grid of points to fit.
+     *  @param[in]  gridShape       Number of points in x and y for the grid of points.
+     *  @param[in]  order           Order of the polynomial (same for forward and
+     *                              reverse transforms).
+     *  @param[in]  useInverse      If true, the inverse SIP transform will be fit
+     *                              and compared to data points generated by calls to
+     *                              pixelToIwc.applyInverse instead of
+     *                              pixelToIwc.applyForward.
+     *  @param[in]  svdThreshold    Fraction of the largest singular value at which to
+     *                              declare smaller singular values zero in the least
+     *                              squares solution.  Negative values use Eigen's
+     *                              internal default.
+     *
+     *  @throws lsst::pex::exceptions::InvalidParameterError Thrown if order is negative
+     *      or gridShape is non-positive.
+     *
+     *  @exceptsafe strong
+     */
+    SipApproximation(
+        std::shared_ptr<TransformPoint2ToPoint2> pixelToIwc,
+        Point2D const & crpix,
+        Eigen::Matrix2d const & cd,
+        Box2D const & bbox,
+        Extent2I const & gridShape,
+        int order,
+        bool useInverse=true,
+        double svdThreshold=-1
+    );
+
+    /**
+     *  Construct from existing SIP coefficients.
+     *
+     *  This constructor is primarily intended for testing purposes.
+     *
+     *  @param[in]  pixelToIwc      The true Transform to approximate.  Should go
+     *                              from pixels to Intermediate World Coordinates
+     *                              when applyForward is called.
+     *  @param[in]  crpix           Pixel origin, using the LSST 0-indexed
+     *                              convention rather than the FITS 1-indexed
+     *                              convention; equal to (CRPIX1 - 1, CRPIX - 1).
+     *  @param[in]  cd              Nominal Jacobian ("CD" in FITS WCS).
+     *  @param[in]  bbox            Pixel-coordinate bounding box over which the
+     *                              approximation should be valid.  Used to construct
+     *                              the grid of points to fit.
+     *  @param[in]  gridShape       Number of points in x and y for the grid of points.
+     *  @param[in]  a               Matrix of A coefficients, with the first dimension
+     *                              corresponding to powers of @f$(u - u_0)@f$ and the
+     *                              second corresponding to powers of @f$(v - v_0)@f$.
+     *  @param[in]  b               Matrix of B coefficients, with the first dimension
+     *                              corresponding to powers of @f$(u - u_0)@f$ and the
+     *                              second corresponding to powers of @f$(v - v_0)@f$.
+     *  @param[in]  ap              Matrix of AP coefficients, with the first dimension
+     *                              corresponding to powers of @f$x_s@f$ and the
+     *                              second corresponding to powers of @f$y_s@f$.
+     *  @param[in]  bp              Matrix of BP coefficients, with the first dimension
+     *                              corresponding to powers of @f$x_s@f$ and the
+     *                              second corresponding to powers of @f$y_s@f$.
+     *  @param[in]  useInverse      If true, the inverse SIP transform will be compared
+     *                              to data points generated by calls to
+     *                              pixelToIwc.applyInverse instead of
+     *                              pixelToIwc.applyForward.
+     *
+     *  @throws lsst::pex::exceptions::InvalidParameterError Thrown if gridShape
+     *       is non-positive, or any matrix argument is non-square.
+     *
+     *  @exceptsafe strong
+     */
+    SipApproximation(
+        std::shared_ptr<TransformPoint2ToPoint2> pixelToIwc,
+        Point2D const & crpix,
+        Eigen::Matrix2d const & cd,
+        Box2D const & bbox,
+        Extent2I const & gridShape,
+        ndarray::Array<double const, 2> const & a,
+        ndarray::Array<double const, 2> const & b,
+        ndarray::Array<double const, 2> const & ap,
+        ndarray::Array<double const, 2> const & bp,
+        bool useInverse=true
+    );
+
+    // No copies just because they'd be a pain to implement and probably aren't necessary
+    SipApproximation(SipApproximation const &) = delete;
+    SipApproximation & operator=(SipApproximation const &) = delete;
+
+    // Moves are fine.
+    SipApproximation(SipApproximation &&) noexcept = default;
+    SipApproximation & operator=(SipApproximation &&) noexcept = default;
+
+    // Needs to be defined in .cc, where compiler can see the definitions of forward-declared
+    // private implementation classes.
+    ~SipApproximation() noexcept;
+
+    /// Return the polynomial order of the current solution (same for forward and reverse).
+    int getOrder() const noexcept;
+
+    /**
+     * Return a coefficient of the forward transform polynomial.
+     *
+     * Out-of-bounds arguments yields undefined behavior.
+     *
+     * @exceptsafe strong
+     */
+    double getA(int p, int q) const;
+
+    /**
+     * Return a coefficient of the forward transform polynomial.
+     *
+     * Out-of-bounds arguments yields undefined behavior.
+     *
+     * @exceptsafe strong
+     */
+    double getB(int p, int q) const;
+
+    /**
+     * Return a coefficient of the reverse transform polynomial.
+     *
+     * Out-of-bounds arguments yields undefined behavior.
+     *
+     * @exceptsafe strong
+     */
+    double getAP(int p, int q) const;
+
+    /**
+     * Return a coefficient of the reverse transform polynomial.
+      *
+     * Out-of-bounds arguments yields undefined behavior.
+     *
+     * @exceptsafe strong
+    */
+    double getBP(int p, int q) const;
+
+    /// Return the coefficients of the forward transform polynomial.
+    Eigen::MatrixXd getA() const noexcept;
+
+    /// Return the coefficients of the forward transform polynomial.
+    Eigen::MatrixXd getB() const noexcept;
+
+    /// Return the coefficients of the reverse transform polynomial.
+    Eigen::MatrixXd getAP() const noexcept;
+
+    /// Return the coefficients of the reverse transform polynomial.
+    Eigen::MatrixXd getBP() const noexcept;
+
+    /**
+     *  Convert a point from pixels to intermediate world coordinates.
+     *
+     *  This method is inefficient and should only be used for diagnostic purposes.
+     *
+     *  @exceptsafe strong
+     */
+    Point2D applyForward(Point2D const & pix) const;
+
+    /**
+     *  Convert an array of points from pixels to intermediate world coordinates.
+     *
+     *  @exceptsafe strong
+     */
+    std::vector<Point2D> applyForward(std::vector<Point2D> const & pix) const;
+
+    /**
+     *  Convert a point from intermediate world coordinates to pixels.
+     *
+     *  This method is inefficient and should only be used for diagnostic purposes.
+     *
+     *  @exceptsafe strong
+     */
+    Point2D applyInverse(Point2D const & iwcs) const;
+
+    /**
+     *  Convert an array of points from intermediate world coordinates to pixels.
+     *
+     *  @exceptsafe strong
+     */
+    std::vector<Point2D> applyInverse(std::vector<Point2D> const & iwcs) const;
+
+    /// Return the distance between grid points in pixels.
+    Extent2D getGridStep() const noexcept;
+
+    /// Return the number of grid points in x and y.
+    Extent2I getGridShape() const noexcept;
+
+    /// Return the pixel-coordinate bounding box over which the approximation should be valid.
+    Box2D getBBox() const noexcept { return _bbox; }
+
+    /// Return the pixel origin of the WCS being approximated.
+    Point2D getPixelOrigin() const noexcept { return Point2D(_crpix); }
+
+    /// Return the CD matrix of the WCS being approximated.
+    Eigen::Matrix2d getCdMatrix() const noexcept { return _cdInv.inverted().getMatrix(); }
+
+    /**
+     *  Update the grid to the given number of points in x and y.
+     *
+     *  This does not invalidate or modify the current solution; this allows
+     *  the user to fit with a coarse grid and then check whether the solution
+     *  still works well on a finer grid.
+     *
+     *  @throws lsst::pex::exceptions::InvalidParameterError Thrown if shape is
+     *     non-positive.
+     *
+     *  @exceptsafe strong
+     */
+    void updateGrid(Extent2I const & shape);
+
+    /**
+     *  Update the grid by making it finer by a given integer factor.
+     *
+     *  @throws lsst::pex::exceptions::InvalidParameterError Thrown if factor is
+     *     non-positive.
+     *
+     *  @exceptsafe strong
+     */
+    void refineGrid(int factor=2);
+
+    /**
+     *  Obtain a new solution at the given order with the current grid.
+     *
+     *  @param[in]  order           Polynomial order to fit.
+     *  @param[in]  svdThreshold    Fraction of the largest singular value at which to
+     *                              declare smaller singular values zero in the least
+     *                              squares solution.  Negative values use Eigen's
+     *                              internal default.
+     *
+     *  @throws pex::exceptions::LogicError Thrown if the number of free
+     *          parameters implied by order is larger than the number of
+     *          data points defined by the grid.
+     *
+     *  @exceptsafe strong
+     */
+    void fit(int order, double svdThreshold=-1);
+
+    /**
+     *  Return the maximum deviation of the solution from the exact transform
+     *  on the current grid.
+     *
+     *  The deviations are in scaled intermediate world coordinates
+     *  @f$\sqrt{\delta x_s^2 \delta y_s^2}@f$ for the forward transform and in
+     *  pixels @f$(\delta u^2, \delta v^2)@f$ for the reverse transform
+     *  (respectively).  Note that in the common case where the CD matrix
+     *  includes the scaling from angle units to pixel units, the *scaled*
+     *  intermediate world coordinate values are also in (nominal) pixel
+     *  units.
+     */
+    std::pair<double, double> computeMaxDeviation() const noexcept;
+
+private:
+
+    struct Grid;
+    struct Solution;
+
+    bool _useInverse;
+    std::shared_ptr<TransformPoint2ToPoint2> _pixelToIwc;
+    Box2D _bbox;
+    Extent2D _crpix;
+    LinearTransform _cdInv;
+    std::unique_ptr<Grid const> _grid;
+    std::unique_ptr<Solution const> _solution;
+};
+
+}}}  // namespace lsst::afw::geom
+
+#endif // !LSST_AFW_GEOM_SipApproximation_h_INCLUDED

python/lsst/afw/geom/SConscript

@@ -4,6 +4,7 @@ scripts.BasicSConscript.pybind11(
     [
         'detail/frameSetUtils',
         'wcsUtils/wcsUtils',
+        'sipApproximation',
         'span',
         'spanSet',
         'endpoint',

python/lsst/afw/geom/__init__.py

@@ -41,3 +41,4 @@
 from .skyWcs import *
 from .transformFromString import *
 from . import wcsUtils
+from .sipApproximation import *