AstrometryTransform.h

Go to the documentation of this file.

// -*- LSST-C++ -*-
/*
 * This file is part of jointcal.
 *
 * 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_JOINTCAL_ASTROMETRY_TRANSFORM_H
#define LSST_JOINTCAL_ASTROMETRY_TRANSFORM_H
 
#include <iostream>
#include <memory>
#include <string>
#include <sstream>
#include <vector>
 
#include "Eigen/Core"
 
#include "lsst/pex/exceptions.h"
#include "lsst/afw/geom/SkyWcs.h"
#include "lsst/jointcal/FatPoint.h"
#include "lsst/jointcal/Frame.h"
 
namespace pexExcept = lsst::pex::exceptions;
 
namespace lsst {
namespace jointcal {
 
class StarMatchList;
class Frame;
class AstrometryTransformLinear;
 
 
class AstrometryTransform {
public:
    virtual void apply(double xIn, double yIn, double &xOut, double &yOut) const = 0;
 
    void apply(Point const &in, Point &out) const { apply(in.x, in.y, out.x, out.y); }
 
    Point apply(Point const &in) const {
        double xout, yout;
        apply(in.x, in.y, xout, yout);
        return Point(xout, yout);
    }
 
    Frame apply(Frame const &inputframe, bool inscribed) const;
 
    virtual void print(std::ostream &out) const = 0;
 
    std::string __str__() const {
        std::stringstream s;
        print(s);
        return s.str();
    }
 
 
    virtual double fit(StarMatchList const &starMatchList) = 0;
 
    void transformStar(FatPoint &in) const { transformPosAndErrors(in, in); }
 
    virtual double getJacobian(Point const &point) const { return getJacobian(point.x, point.y); }
 
    virtual std::unique_ptr<AstrometryTransform> clone() const = 0;
 
    virtual std::unique_ptr<AstrometryTransform> composeAndReduce(AstrometryTransform const &right) const;
 
    virtual double getJacobian(double x, double y) const;
 
    virtual void computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
                                   double step = 0.01) const;
 
    virtual AstrometryTransformLinear linearApproximation(Point const &where, double step = 0.01) const;
 
    virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const;
 
    virtual void transformErrors(Point const &where, const double *vIn, double *vOut) const;
 
 
    virtual std::unique_ptr<AstrometryTransform> inverseTransform(double precision,
                                                                  const Frame &region) const;
 
    void getParams(double *params) const;
 
    void offsetParams(Eigen::VectorXd const &delta);
 
    virtual double paramRef(Eigen::Index i) const;
 
    virtual double &paramRef(Eigen::Index i);
 
    virtual void paramDerivatives(Point const &where, double *dx, double *dy) const;
 
 
    virtual std::unique_ptr<AstrometryTransform> roughInverse(const Frame &region) const;
 
    virtual std::size_t getNpar() const { return 0; }
 
    virtual std::shared_ptr<ast::Mapping> toAstMap(jointcal::Frame const &domain) const {
        throw std::logic_error("toAstMap is not implemented for this class.");
    }
 
    void write(const std::string &fileName) const;
 
    virtual void write(std::ostream &stream) const;
 
    virtual ~AstrometryTransform() = default;;
};
 
std::ostream &operator<<(std::ostream &stream, AstrometryTransform const &transform);
 
std::unique_ptr<AstrometryTransform> compose(AstrometryTransform const &left,
                                             AstrometryTransform const &right);
 
/*=============================================================*/
 
class AstrometryTransformIdentity : public AstrometryTransform {
public:
    AstrometryTransformIdentity() = default;
 
    void apply(const double xIn, const double yIn, double &xOut, double &yOut) const override {
        xOut = xIn;
        yOut = yIn;
    }  // to speed up
 
    double fit(StarMatchList const &starMatchList) override {
        throw pexExcept::TypeError(
                "AstrometryTransformIdentity is the identity transformation: it cannot be fit to anything.");
    }
 
    std::unique_ptr<AstrometryTransform> composeAndReduce(AstrometryTransform const &right) const override {
        return right.clone();
    }
 
    void print(std::ostream &out) const override { out << "x' = x\ny' = y" << std::endl; }
 
    std::size_t getNpar() const override { return 0; }
 
    std::unique_ptr<AstrometryTransform> clone() const override {
        return std::unique_ptr<AstrometryTransform>(new AstrometryTransformIdentity);
    }
 
    void computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
                           double step = 0.01) const override;
 
    virtual AstrometryTransformLinear linearApproximation(Point const &where,
                                                          double step = 0.01) const override;
 
    std::shared_ptr<ast::Mapping> toAstMap(jointcal::Frame const &domain) const override;
 
    void write(std::ostream &s) const override;
 
    void read(std::istream &s);
 
    //    ClassDef(AstrometryTransformIdentity,1)
};
 
std::unique_ptr<AstrometryTransform> compose(AstrometryTransform const &left,
                                             AstrometryTransformIdentity const &right);
 
bool isIntegerShift(const AstrometryTransform *transform);
 
/*====================   AstrometryTransformPolynomial  =======================*/
 
class AstrometryTransformPolynomial : public AstrometryTransform {
public:
    AstrometryTransformPolynomial(std::size_t order = 1);
 
    AstrometryTransformPolynomial(const AstrometryTransform *transform, const Frame &frame, std::size_t order,
                                  std::size_t nPoint = 1000);
 
    AstrometryTransformPolynomial(const std::shared_ptr<afw::geom::TransformPoint2ToPoint2>& transform,
                                  jointcal::Frame const &domain, std::size_t order, std::size_t nSteps = 50);
 
    void setOrder(std::size_t order);
    std::size_t getOrder() const { return _order; }
 
    using AstrometryTransform::apply;  // to unhide AstrometryTransform::apply(Point const &)
 
    void apply(double xIn, double yIn, double &xOut, double &yOut) const override;
 
    void computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
                           double step = 0.01) const override;
 
    virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const override;
 
    std::size_t getNpar() const override { return 2 * _nterms; }
 
    void print(std::ostream &out) const override;
 
    double fit(StarMatchList const &starMatchList) override;
 
    AstrometryTransformPolynomial operator*(AstrometryTransformPolynomial const &right) const;
 
    AstrometryTransformPolynomial operator+(AstrometryTransformPolynomial const &right) const;
 
    AstrometryTransformPolynomial operator-(AstrometryTransformPolynomial const &right) const;
 
    using AstrometryTransform::composeAndReduce;  // to unhide
                                                  // AstrometryTransform::composeAndReduce(AstrometryTransform
                                                  // const &)
    std::unique_ptr<AstrometryTransform> composeAndReduce(AstrometryTransformPolynomial const &right) const;
 
    std::unique_ptr<AstrometryTransform> clone() const override {
        return std::unique_ptr<AstrometryTransform>(new AstrometryTransformPolynomial(*this));
    }
 
    double getCoefficient(std::size_t powX, std::size_t powY, std::size_t whichCoord) const;
    double &getCoefficient(std::size_t powX, std::size_t powY, std::size_t whichCoord);
    double coeffOrZero(std::size_t powX, std::size_t powY, std::size_t whichCoord) const;
 
    double determinant() const;
 
    double paramRef(Eigen::Index i) const override;
 
    double &paramRef(Eigen::Index i) override;
 
    void paramDerivatives(Point const &where, double *dx, double *dy) const override;
 
    std::shared_ptr<ast::Mapping> toAstMap(jointcal::Frame const &domain) const override;
 
    void write(std::ostream &s) const override;
    void read(std::istream &s);
 
private:
    double computeFit(StarMatchList const &starMatchList, AstrometryTransform const &shiftToCenter,
                      bool useErrors);
 
    std::size_t _order{};           // The highest sum of exponents of the largest monomial.
    std::size_t _nterms{};          // number of parameters per coordinate
    std::vector<double> _coeffs;  // the actual coefficients
                                  // both polynomials in a single vector to speed up allocation and copies
 
    /* use std::vector rather than double * to avoid
       writing copy constructor and "operator =".
       Vect would work as well but introduces a dependence
       that can be avoided */
 
    /* This routine take a double * for the vector because the array can
       then be allocated on the execution stack, which speeds thing
       up. However this uses Variable Length Array (VLA) which is not
       part of C++, but gcc implements it. */
    void computeMonomials(double xIn, double yIn, double *monomial) const;
 
    ndarray::Array<double, 2, 2> toAstPolyMapCoefficients() const;
};
 
std::shared_ptr<AstrometryTransformPolynomial> inversePolyTransform(AstrometryTransform const &forward,
                                                                    Frame const &domain,
                                                                    double precision,
                                                                    std::size_t maxOrder = 9,
                                                                    std::size_t nSteps = 50);
 
AstrometryTransformLinear normalizeCoordinatesTransform(const Frame &frame);
 
/*=============================================================*/
class AstrometryTransformLinear : public AstrometryTransformPolynomial {
public:
    using AstrometryTransformPolynomial::apply;  // to unhide AstrometryTransform::apply(Point const &)
 
    AstrometryTransformLinear() : AstrometryTransformPolynomial(1){};
 
    explicit AstrometryTransformLinear(AstrometryTransformPolynomial const &transform);
 
    AstrometryTransformLinear operator*(AstrometryTransformLinear const &right) const;
 
    AstrometryTransformLinear inverted() const;
 
    void print(std::ostream &out) const override;
 
    // useful?    double jacobian(const double x, const double y) const { return determinant();}
 
    void computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
                           double step = 0.01) const override;
    AstrometryTransformLinear linearApproximation(Point const &where,
                                                  double step = 0.01) const override;
 
    //  void print(std::ostream &out) const;
 
    // double fit(StarMatchList const &starMatchList);
 
    AstrometryTransformLinear(double ox, double oy, double aa11, double aa12,
                              double aa21, double aa22);
 
    AstrometryTransformLinear(AstrometryTransformIdentity const &) : AstrometryTransformPolynomial(1){};
 
    std::unique_ptr<AstrometryTransform> clone() const override {
        return std::unique_ptr<AstrometryTransform>(new AstrometryTransformLinear(*this));
    }
 
    std::unique_ptr<AstrometryTransform> inverseTransform(double precision,
                                                          const Frame &region) const override;
 
    double A11() const { return getCoefficient(1, 0, 0); }
    double A12() const { return getCoefficient(0, 1, 0); }
    double A21() const { return getCoefficient(1, 0, 1); }
    double A22() const { return getCoefficient(0, 1, 1); }
    double Dx() const { return getCoefficient(0, 0, 0); }
    double Dy() const { return getCoefficient(0, 0, 1); }
 
protected:
    double &a11() { return getCoefficient(1, 0, 0); }
    double &a12() { return getCoefficient(0, 1, 0); }
    double &a21() { return getCoefficient(1, 0, 1); }
    double &a22() { return getCoefficient(0, 1, 1); }
    double &dx() { return getCoefficient(0, 0, 0); }
    double &dy() { return getCoefficient(0, 0, 1); }
 
    friend class AstrometryTransform;
    friend class AstrometryTransformIdentity;    // for AstrometryTransform::Derivative
    friend class AstrometryTransformPolynomial;  // // for AstrometryTransform::Derivative
 
private:
    void setOrder(std::size_t order);  // to hide AstrometryTransformPolynomial::setOrder
};
 
/*=============================================================*/
 
class AstrometryTransformLinearShift : public AstrometryTransformLinear {
public:
    using AstrometryTransform::apply;  // to unhide AstrometryTransform::apply(Point const &)
    AstrometryTransformLinearShift(double ox = 0., double oy = 0.)
            : AstrometryTransformLinear(ox, oy, 1., 0., 0., 1.) {}
    AstrometryTransformLinearShift(Point const &point)
            : AstrometryTransformLinear(point.x, point.y, 1., 0., 0., 1.){};
    double fit(StarMatchList const &starMatchList);
 
    std::size_t getNpar() const { return 2; }
};
 
/*=============================================================*/
class AstrometryTransformLinearRot : public AstrometryTransformLinear {
public:
    using AstrometryTransform::apply;  // to unhide apply(const Point&)
 
    AstrometryTransformLinearRot() : AstrometryTransformLinear(){};
    AstrometryTransformLinearRot(double angleRad, const Point *center = nullptr,
                                 double scaleFactor = 1.0);
    double fit(StarMatchList const &starMatchList);
 
    std::size_t getNpar() const { return 4; }
};
 
/*=============================================================*/
 
class AstrometryTransformLinearScale : public AstrometryTransformLinear {
public:
    using AstrometryTransform::apply;  // to unhide apply(const Point&)
    AstrometryTransformLinearScale(const double scale = 1)
            : AstrometryTransformLinear(0.0, 0.0, scale, 0., 0., scale){};
    AstrometryTransformLinearScale(const double scaleX, const double scaleY)
            : AstrometryTransformLinear(0.0, 0.0, scaleX, 0., 0., scaleY){};
 
    std::size_t getNpar() const { return 2; }
};
 
class AstrometryTransformSkyWcs : public AstrometryTransform {
public:
    AstrometryTransformSkyWcs(std::shared_ptr<afw::geom::SkyWcs> skyWcs);
 
    using AstrometryTransform::apply;
 
    // Input is x, y pixels; output is ICRS RA, Dec in degrees
    void apply(double xIn, double yIn, double &xOut, double &yOut) const override;
 
    void print(std::ostream &out) const override;
 
    double fit(const StarMatchList &starMatchList) override;
 
    std::unique_ptr<AstrometryTransform> clone() const override;
 
    std::shared_ptr<afw::geom::SkyWcs> getSkyWcs() const { return _skyWcs; }
 
private:
    std::shared_ptr<afw::geom::SkyWcs> _skyWcs;
};
 
/*==================WCS's transform's =====================================*/
 
class BaseTanWcs : public AstrometryTransform {
public:
    using AstrometryTransform::apply;  // to unhide apply(const Point&)
 
    BaseTanWcs(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint,
               const AstrometryTransformPolynomial *corrections = nullptr);
 
    BaseTanWcs(const BaseTanWcs &original);
 
    BaseTanWcs &operator=(const BaseTanWcs &original);
 
    void apply(double xIn, double yIn, double &xOut, double &yOut) const;
 
    Point getTangentPoint() const;
 
    AstrometryTransformLinear getLinPart() const;
 
    const AstrometryTransformPolynomial *getCorr() const { return corr.get(); }
 
    void setCorrections(std::unique_ptr<AstrometryTransformPolynomial> corrections);
 
    Point getCrPix() const;
 
    virtual AstrometryTransformPolynomial getPixelToTangentPlane() const = 0;
 
    virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
                                   double &yTangentPlane) const = 0;
 
    ~BaseTanWcs();
 
protected:
    AstrometryTransformLinear linPixelToTan;  // transform from pixels to tangent plane (degrees)
                                              // a linear approximation centered at the pixel and sky origins
    std::unique_ptr<AstrometryTransformPolynomial> corr;
    double ra0{}, dec0{};   // sky origin (radians)
    double cos0{}, sin0{};  // cos(dec0), sin(dec0)
};
 
class TanRaDecToPixel;  // the inverse of TanPixelToRaDec.
 
class TanPixelToRaDec : public BaseTanWcs {
public:
    using AstrometryTransform::apply;  // to unhide apply(const Point&)
    TanPixelToRaDec(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint,
                    const AstrometryTransformPolynomial *corrections = nullptr);
 
    AstrometryTransformPolynomial getPixelToTangentPlane() const;
 
    virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
                                   double &yTangentPlane) const;
 
    TanPixelToRaDec();
 
    TanPixelToRaDec operator*(AstrometryTransformLinear const &right) const;
 
    using AstrometryTransform::composeAndReduce;  // to unhide
                                                  // AstrometryTransform::composeAndReduce(AstrometryTransform
                                                  // const &)
    std::unique_ptr<AstrometryTransform> composeAndReduce(AstrometryTransformLinear const &right) const;
 
    TanRaDecToPixel inverted() const;
 
    std::unique_ptr<AstrometryTransform> roughInverse(const Frame &region) const;
 
    std::unique_ptr<AstrometryTransform> inverseTransform(double precision, const Frame &region) const;
 
    std::unique_ptr<AstrometryTransform> clone() const;
 
    void print(std::ostream &out) const;
 
    double fit(StarMatchList const &starMatchList);
};
 
class TanSipPixelToRaDec : public BaseTanWcs {
public:
    TanSipPixelToRaDec(AstrometryTransformLinear const &pixToTan, Point const &tangentPoint,
                       const AstrometryTransformPolynomial *corrections = nullptr);
 
    AstrometryTransformPolynomial getPixelToTangentPlane() const;
 
    virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
                                   double &yTangentPlane) const;
 
    TanSipPixelToRaDec();
 
    std::unique_ptr<AstrometryTransform> inverseTransform(double precision, const Frame &region) const;
 
    std::unique_ptr<AstrometryTransform> clone() const;
 
    void print(std::ostream &out) const;
 
    double fit(StarMatchList const &starMatchList);
};
 
 
class TanRaDecToPixel : public AstrometryTransform {
public:
    using AstrometryTransform::apply;  // to unhide apply(const Point&)
 
    TanRaDecToPixel(AstrometryTransformLinear tan2Pix, Point const &tangentPoint);
 
    TanRaDecToPixel();
 
    AstrometryTransformLinear getLinPart() const;
 
    void setTangentPoint(Point const &tangentPoint);
 
    Point getTangentPoint() const;
 
    void apply(double xIn, double yIn, double &xOut, double &yOut) const;
 
    void transformPosAndErrors(const FatPoint &in, FatPoint &out) const;
 
    TanPixelToRaDec inverted() const;
 
    std::unique_ptr<AstrometryTransform> roughInverse(const Frame &region) const;
 
    std::unique_ptr<AstrometryTransform> inverseTransform(double precision, const Frame &region) const;
 
    void print(std::ostream &out) const;
 
    std::unique_ptr<AstrometryTransform> clone() const;
 
    double fit(StarMatchList const &starMatchList);
 
private:
    double ra0{}, dec0{};  // tangent point (radians)
    double cos0{}, sin0{};
    AstrometryTransformLinear linTan2Pix;  // tangent plane (probably degrees) to pixels
};
 
using AstrometryTransformFun = void (const double, const double, double &, double &, const void *);
 
class UserTransform : public AstrometryTransform {
public:
    using AstrometryTransform::apply;  // to unhide apply(const Point&)
 
    UserTransform(AstrometryTransformFun &fun, const void *userData);
 
    void apply(double xIn, double yIn, double &xOut, double &yOut) const;
 
    void print(std::ostream &out) const;
 
    double fit(StarMatchList const &starMatchList);
 
    std::unique_ptr<AstrometryTransform> clone() const;
 
private:
    AstrometryTransformFun *_userFun;
    const void *_userData;
};
 
std::unique_ptr<AstrometryTransform> astrometryTransformRead(const std::string &fileName);
std::unique_ptr<AstrometryTransform> astrometryTransformRead(std::istream &s);
}  // namespace jointcal
}  // namespace lsst
 
#endif  // LSST_JOINTCAL_ASTROMETRY_TRANSFORM_H