Gtransfo.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_GTRANSFO_H
#define LSST_JOINTCAL_GTRANSFO_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 GtransfoLin;


class Gtransfo {
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<Gtransfo> clone() const = 0;

    virtual std::unique_ptr<Gtransfo> composeAndReduce(Gtransfo const &right) const;

    virtual double getJacobian(const double x, const double y) const;

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

    virtual GtransfoLin 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<Gtransfo> inverseTransfo(const double precision, const Frame &region) const;

    void getParams(double *params) const;

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

    virtual double paramRef(const int i) const;

    virtual double &paramRef(const int i);

    virtual void paramDerivatives(Point const &where, double *dx, double *dy) const;


    virtual std::unique_ptr<Gtransfo> roughInverse(const Frame &region) const;

    virtual int 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 ~Gtransfo(){};
};

std::ostream &operator<<(std::ostream &stream, Gtransfo const &gtransfo);

std::unique_ptr<Gtransfo> gtransfoCompose(Gtransfo const &left, Gtransfo const &right);

/*=============================================================*/

class GtransfoIdentity : public Gtransfo {
public:
    GtransfoIdentity() {}

    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(
                "GtransfoIdentity is the identity transformation: it cannot be fit to anything.");
    }

    std::unique_ptr<Gtransfo> composeAndReduce(Gtransfo const &right) const override { return right.clone(); }

    void dump(std::ostream &stream = std::cout) const override { stream << "x' = x\ny' = y" << std::endl; }

    int getNpar() const override { return 0; }

    std::unique_ptr<Gtransfo> clone() const override {
        return std::unique_ptr<Gtransfo>(new GtransfoIdentity);
    }

    void computeDerivative(Point const &where, GtransfoLin &derivative,
                           const double step = 0.01) const override;

    virtual GtransfoLin 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(GtransfoIdentity,1)
};

std::unique_ptr<Gtransfo> gtransfoCompose(Gtransfo const &left, GtransfoIdentity const &right);

bool isIntegerShift(const Gtransfo *gtransfo);

/*====================   GtransfoPoly  =======================*/

class GtransfoPoly : public Gtransfo {
public:
    GtransfoPoly(const unsigned order = 1);

    GtransfoPoly(const Gtransfo *gtransfo, const Frame &frame, unsigned order, unsigned nPoint = 1000);

    GtransfoPoly(std::shared_ptr<afw::geom::TransformPoint2ToPoint2> transform, jointcal::Frame const &domain,
                 unsigned const order, unsigned const nSteps = 50);

    void setOrder(const unsigned order);
    unsigned getOrder() const { return _order; }

    using Gtransfo::apply;  // to unhide Gtransfo::apply(Point const &)

    void apply(const double xIn, const double yIn, double &xOut, double &yOut) const override;

    void computeDerivative(Point const &where, GtransfoLin &derivative,
                           const double step = 0.01) const override;

    virtual void transformPosAndErrors(const FatPoint &in, FatPoint &out) const override;

    int getNpar() const override { return 2 * _nterms; }

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

    double fit(StarMatchList const &starMatchList) override;

    GtransfoPoly operator*(GtransfoPoly const &right) const;

    GtransfoPoly operator+(GtransfoPoly const &right) const;

    GtransfoPoly operator-(GtransfoPoly const &right) const;

    using Gtransfo::composeAndReduce;  // to unhide Gtransfo::composeAndReduce(Gtransfo const &)
    std::unique_ptr<Gtransfo> composeAndReduce(GtransfoPoly const &right) const;

    std::unique_ptr<Gtransfo> clone() const override {
        return std::unique_ptr<Gtransfo>(new GtransfoPoly(*this));
    }

    double coeff(const unsigned powX, const unsigned powY, const unsigned whichCoord) const;

    double &coeff(const unsigned powX, const unsigned powY, const unsigned whichCoord);

    double coeffOrZero(const unsigned powX, const unsigned powY, const unsigned whichCoord) const;

    double determinant() const;

    double paramRef(const int i) const override;

    double &paramRef(const int 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, Gtransfo const &InTransfo, const bool UseErrors);

    unsigned _order;              // The highest sum of exponents of the largest monomial.
    unsigned _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<GtransfoPoly> inversePolyTransfo(Gtransfo const &forward, Frame const &domain,
                                                 double const precision, int const maxOrder = 9,
                                                 unsigned const nSteps = 50);

GtransfoLin normalizeCoordinatesTransfo(const Frame &frame);

/*=============================================================*/
class GtransfoLin : public GtransfoPoly {
public:
    using GtransfoPoly::apply;  // to unhide Gtransfo::apply(Point const &)

    GtransfoLin() : GtransfoPoly(1){};

    explicit GtransfoLin(GtransfoPoly const &gtransfoPoly);

    GtransfoLin operator*(GtransfoLin const &right) const;

    GtransfoLin inverted() const;

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

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

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

    // double fit(StarMatchList const &starMatchList);

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

    GtransfoLin(GtransfoIdentity const &) : GtransfoPoly(1){};

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

    std::unique_ptr<Gtransfo> inverseTransfo(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 Gtransfo;
    friend class GtransfoIdentity;  // for Gtransfo::Derivative
    friend class GtransfoPoly;      // // for Gtransfo::Derivative

private:
    void setOrder(const unsigned order);  // to hide GtransfoPoly::setOrder
};

/*=============================================================*/

class GtransfoLinShift : public GtransfoLin {
public:
    using Gtransfo::apply;  // to unhide Gtransfo::apply(Point const &)
    GtransfoLinShift(double ox = 0., double oy = 0.) : GtransfoLin(ox, oy, 1., 0., 0., 1.) {}
    GtransfoLinShift(Point const &point) : GtransfoLin(point.x, point.y, 1., 0., 0., 1.){};
    double fit(StarMatchList const &starMatchList);

    int getNpar() const { return 2; }
};

/*=============================================================*/
class GtransfoLinRot : public GtransfoLin {
public:
    using Gtransfo::apply;  // to unhide apply(const Point&)

    GtransfoLinRot() : GtransfoLin(){};
    GtransfoLinRot(const double angleRad, const Point *center = nullptr, const double scaleFactor = 1.0);
    double fit(StarMatchList const &starMatchList);

    int getNpar() const { return 4; }
};

/*=============================================================*/

class GtransfoLinScale : public GtransfoLin {
public:
    using Gtransfo::apply;  // to unhide apply(const Point&)
    GtransfoLinScale(const double scale = 1) : GtransfoLin(0.0, 0.0, scale, 0., 0., scale){};
    GtransfoLinScale(const double scaleX, const double scaleY)
            : GtransfoLin(0.0, 0.0, scaleX, 0., 0., scaleY){};

    int getNpar() const { return 2; }
};

class GtransfoSkyWcs : public Gtransfo {
public:
    GtransfoSkyWcs(std::shared_ptr<afw::geom::SkyWcs> skyWcs);

    using Gtransfo::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<Gtransfo> clone() const override;

    std::shared_ptr<afw::geom::SkyWcs> getSkyWcs() const { return _skyWcs; }

private:
    std::shared_ptr<afw::geom::SkyWcs> _skyWcs;
};

/*==================WCS's transfo's =====================================*/

class BaseTanWcs : public Gtransfo {
public:
    using Gtransfo::apply;  // to unhide apply(const Point&)

    BaseTanWcs(GtransfoLin const &pixToTan, Point const &tangentPoint,
               const GtransfoPoly *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;

    GtransfoLin getLinPart() const;

    const GtransfoPoly *getCorr() const { return corr.get(); }

    void setCorrections(std::unique_ptr<GtransfoPoly> corrections);

    Point getCrPix() const;

    virtual GtransfoPoly getPixelToTangentPlane() const = 0;

    virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
                                   double &yTangentPlane) const = 0;

    ~BaseTanWcs();

protected:
    GtransfoLin linPixelToTan;  // transform from pixels to tangent plane (degrees)
                                // a linear approximation centered at the pixel and sky origins
    std::unique_ptr<GtransfoPoly> 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 Gtransfo::apply;  // to unhide apply(const Point&)
    TanPixelToRaDec(GtransfoLin const &pixToTan, Point const &tangentPoint,
                    const GtransfoPoly *corrections = nullptr);

    GtransfoPoly getPixelToTangentPlane() const;

    virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
                                   double &yTangentPlane) const;

    TanPixelToRaDec();

    TanPixelToRaDec operator*(GtransfoLin const &right) const;

    using Gtransfo::composeAndReduce;  // to unhide Gtransfo::composeAndReduce(Gtransfo const &)
    std::unique_ptr<Gtransfo> composeAndReduce(GtransfoLin const &right) const;

    TanRaDecToPixel inverted() const;

    std::unique_ptr<Gtransfo> roughInverse(const Frame &region) const;

    std::unique_ptr<Gtransfo> inverseTransfo(const double precision, const Frame &region) const;

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

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

    double fit(StarMatchList const &starMatchList);
};

class TanSipPixelToRaDec : public BaseTanWcs {
public:
    TanSipPixelToRaDec(GtransfoLin const &pixToTan, Point const &tangentPoint,
                       const GtransfoPoly *corrections = nullptr);

    GtransfoPoly getPixelToTangentPlane() const;

    virtual void pixToTangentPlane(double xPixel, double yPixel, double &xTangentPlane,
                                   double &yTangentPlane) const;

    TanSipPixelToRaDec();

    std::unique_ptr<Gtransfo> inverseTransfo(const double precision, const Frame &region) const;

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

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

    double fit(StarMatchList const &starMatchList);
};


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

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

    TanRaDecToPixel();

    GtransfoLin 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<Gtransfo> roughInverse(const Frame &region) const;

    std::unique_ptr<Gtransfo> inverseTransfo(const double precision, const Frame &region) const;

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

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

    double fit(StarMatchList const &starMatchList);

private:
    double ra0, dec0;  // tangent point (radians)
    double cos0, sin0;
    GtransfoLin linTan2Pix;  // tangent plane (probably degrees) to pixels
};

typedef void(GtransfoFun)(const double, const double, double &, double &, const void *);

class UserTransfo : public Gtransfo {
public:
    using Gtransfo::apply;  // to unhide apply(const Point&)

    UserTransfo(GtransfoFun &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<Gtransfo> clone() const;

private:
    GtransfoFun *_userFun;
    const void *_userData;
};

std::unique_ptr<Gtransfo> gtransfoRead(const std::string &fileName);
std::unique_ptr<Gtransfo> gtransfoRead(std::istream &s);
}  // namespace jointcal
}  // namespace lsst

#endif  // LSST_JOINTCAL_GTRANSFO_H