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(const double xIn, const 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 dump(std::ostream &stream = std::cout) const = 0;

    std::string __str__() {
        std::stringstream s;
        dump(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(const double x, const double y) const;

    virtual void computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
                                   const double step = 0.01) const;

    virtual AstrometryTransformLinear linearApproximation(Point const &where, const 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(const double precision,
                                                                  const Frame &region) const;

    void getParams(double *params) const;

    void offsetParams(Eigen::VectorXd const &delta);

    virtual double paramRef(Eigen::Index const i) const;

    virtual double &paramRef(Eigen::Index const 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(){};
};

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() {}

    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 dump(std::ostream &stream = std::cout) const override { stream << "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,
                           const double step = 0.01) const override;

    virtual AstrometryTransformLinear linearApproximation(Point const &where,
                                                          const 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(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(const double xIn, const double yIn, double &xOut, double &yOut) const override;

    void computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
                           const 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 dump(std::ostream &stream = std::cout) 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 coeff(std::size_t powX, std::size_t powY, std::size_t whichCoord) const;

    double &coeff(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 const i) const override;

    double &paramRef(Eigen::Index const 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,
                      const 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 const 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;

    // useful?    double jacobian(const double x, const double y) const { return determinant();}

    void computeDerivative(Point const &where, AstrometryTransformLinear &derivative,
                           const double step = 0.01) const;
    AstrometryTransformLinear linearApproximation(Point const &where, const double step = 0.01) const;

    //  void dump(std::ostream &stream = std::cout) const;

    // double fit(StarMatchList const &starMatchList);

    AstrometryTransformLinear(const double ox, const double oy, const double aa11, const double aa12,
                              const double aa21, const double aa22);

    AstrometryTransformLinear(AstrometryTransformIdentity const &) : AstrometryTransformPolynomial(1){};

    std::unique_ptr<AstrometryTransform> clone() const {
        return std::unique_ptr<AstrometryTransform>(new AstrometryTransformLinear(*this));
    }

    std::unique_ptr<AstrometryTransform> inverseTransform(const double precision, const Frame &region) const;

    double A11() const { return coeff(1, 0, 0); }
    double A12() const { return coeff(0, 1, 0); }
    double A21() const { return coeff(1, 0, 1); }
    double A22() const { return coeff(0, 1, 1); }
    double Dx() const { return coeff(0, 0, 0); }
    double Dy() const { return coeff(0, 0, 1); }

protected:
    double &a11() { return coeff(1, 0, 0); }
    double &a12() { return coeff(0, 1, 0); }
    double &a21() { return coeff(1, 0, 1); }
    double &a22() { return coeff(0, 1, 1); }
    double &dx() { return coeff(0, 0, 0); }
    double &dy() { return coeff(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(const double angleRad, const Point *center = nullptr,
                                 const 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(const double xIn, const double yIn, double &xOut, double &yOut) const override;

    void dump(std::ostream &stream = std::cout) 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);

    void operator=(const BaseTanWcs &original);

    void apply(const double xIn, const 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(const double precision, const Frame &region) const;

    std::unique_ptr<AstrometryTransform> clone() const;

    void dump(std::ostream &stream) 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(const double precision, const Frame &region) const;

    std::unique_ptr<AstrometryTransform> clone() const;

    void dump(std::ostream &stream) const;

    double fit(StarMatchList const &starMatchList);
};


class TanRaDecToPixel : public AstrometryTransform {
public:
    using AstrometryTransform::apply;  // to unhide apply(const Point&)

    TanRaDecToPixel(AstrometryTransformLinear const &tan2Pix, Point const &tangentPoint);

    TanRaDecToPixel();

    AstrometryTransformLinear getLinPart() const;

    void setTangentPoint(Point const &tangentPoint);

    Point getTangentPoint() const;

    void apply(const double xIn, const 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(const double precision, const Frame &region) const;

    void dump(std::ostream &stream) 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
};

typedef void(AstrometryTransformFun)(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(const double xIn, const double yIn, double &xOut, double &yOut) const;

    void dump(std::ostream &stream = std::cout) 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