lsst::jointcal::GtransfoPoly Class Reference

Polynomial transformation class. More...

#include <Gtransfo.h>

Inheritance diagram for lsst::jointcal::GtransfoPoly:

Public Member Functions
	GtransfoPoly (const unsigned order=1)
	Default transfo : identity for all orders (>=1 ). More...
	GtransfoPoly (const Gtransfo *gtransfo, const Frame &frame, unsigned order, unsigned nPoint=1000)
	Constructs a "polynomial image" from an existing transfo, over a specified domain. More...
	GtransfoPoly (std::shared_ptr< afw::geom::TransformPoint2ToPoint2 > transform, jointcal::Frame const &domain, unsigned const order, unsigned const nSteps=50)
	Constructs a polynomial approximation to an afw::geom::TransformPoint2ToPoint2. More...
void	setOrder (const unsigned order)
	Sets the polynomial order (the highest sum of exponents of the largest monomial). More...
unsigned	getOrder () const
	Returns the polynomial order. More...
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
	specialised analytic routine More...
virtual void	transformPosAndErrors (const FatPoint &in, FatPoint &out) const override
	a mix of apply and Derivative More...
int	getNpar () const override
	total number of parameters More...
void	dump (std::ostream &stream=std::cout) const override
	print out of coefficients in a readable form. More...
double	fit (StarMatchList const &starMatchList) override
	guess what More...
GtransfoPoly	operator* (GtransfoPoly const &right) const
	Composition (internal stuff in quadruple precision) More...
GtransfoPoly	operator+ (GtransfoPoly const &right) const
	Addition. More...
GtransfoPoly	operator- (GtransfoPoly const &right) const
	Subtraction. More...
std::unique_ptr< Gtransfo >	composeAndReduce (GtransfoPoly const &right) const
	Return a reduced composition of newTransfo = this(right()), or nullptr if it cannot be reduced. More...
std::unique_ptr< Gtransfo >	clone () const override
	returns a copy (allocated by new) of the transformation. More...
double	coeff (const unsigned powX, const unsigned powY, const unsigned whichCoord) const
	access to coefficients (read only) More...
double &	coeff (const unsigned powX, const unsigned powY, const unsigned whichCoord)
	write access More...
double	coeffOrZero (const unsigned powX, const unsigned powY, const unsigned whichCoord) const
	read access, zero if beyond order More...
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
	Derivative w.r.t parameters. More...
std::shared_ptr< ast::Mapping >	toAstMap (jointcal::Frame const &domain) const override
	Create an equivalent AST mapping for this transformation, including an analytic inverse if possible. More...
void	write (std::ostream &s) const override
void	read (std::istream &s)
void	apply (Point const &in, Point &out) const
	applies the tranfo to in and writes into out. Is indeed virtual. More...
Point	apply (Point const &in) const
	All these apply(..) shadow the virtual one in derived classes, unless one writes "using Gtransfo::apply". More...
Frame	apply (Frame const &inputframe, bool inscribed) const
	Transform a bounding box, taking either the inscribed or circumscribed box. More...
std::string	__str__ ()
void	transformStar (FatPoint &in) const
	allows to write MyTransfo(MyStar) More...
virtual double	getJacobian (Point const &point) const
	returns the local jacobian. More...
virtual double	getJacobian (const double x, const double y) const
	returns the local jacobian. More...
virtual std::unique_ptr< Gtransfo >	composeAndReduce (Gtransfo const &right) const
	Return a reduced composition of newTransfo = this(right()), or nullptr if it cannot be reduced. More...
virtual GtransfoLin	linearApproximation (Point const &where, const double step=0.01) const
	linear (local) approximation. More...
virtual void	transformErrors (Point const &where, const double *vIn, double *vOut) const
	transform errors (represented as double[3] in order V(xx),V(yy),Cov(xy)) More...
virtual std::unique_ptr< Gtransfo >	inverseTransfo (const double precision, const Frame &region) const
	returns an inverse transfo. Numerical if not overloaded. More...
void	getParams (double *params) const
	params should be at least Npar() long More...
void	offsetParams (Eigen::VectorXd const &delta)
virtual std::unique_ptr< Gtransfo >	roughInverse (const Frame &region) const
	Rough inverse. More...
void	write (const std::string &fileName) const

Detailed Description

Polynomial transformation class.

Definition at line 276 of file Gtransfo.h.

Constructor & Destructor Documentation

GtransfoPoly() [1/3]

lsst::jointcal::GtransfoPoly::GtransfoPoly

(

const unsigned

order = 1

)

Default transfo : identity for all orders (>=1 ).

Default transfo : identity for all orders (>=1 )

Parameters

order

The highest total power (x+y) of monomials of this polynomial.

Definition at line 447 of file Gtransfo.cc.

                                               : _order(order) {
    _nterms = (order + 1) * (order + 2) / 2;

    // allocate and fill coefficients
    _coeffs.resize(2 * _nterms, 0.);
    // the default is supposed to be the identity, (for order>=1).
    if (_order >= 1) {
        coeff(1, 0, 0) = 1;
        coeff(0, 1, 1) = 1;
    }
}

GtransfoPoly() [2/3]

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

Constructs a "polynomial image" from an existing transfo, over a specified domain.

Definition at line 460 of file Gtransfo.cc.

                                                                                                        {
    StarMatchList sm;

    double step = std::sqrt(fabs(frame.getArea()) / double(nPoint));
    for (double x = frame.xMin + step / 2; x <= frame.xMax; x += step)
        for (double y = frame.yMin + step / 2; y <= frame.yMax; y += step) {
            auto pix = std::make_shared<BaseStar>(x, y, 0, 0);
            double xtr, ytr;
            gtransfo->apply(x, y, xtr, ytr);
            auto tp = std::make_shared<BaseStar>(xtr, ytr, 0, 0);
            /* These are fake stars so no need to transform fake errors.
               all errors (and weights) will be equal : */
            sm.push_back(StarMatch(*pix, *tp, pix, tp));
        }
    GtransfoPoly ret(order);
    ret.fit(sm);
    *this = ret;
}

GtransfoPoly() [3/3]

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

Constructs a polynomial approximation to an afw::geom::TransformPoint2ToPoint2.

Parameters

[in]	transform	The transform to be approximated.
[in]	domain	The valid domain of the transform.
[in]	order	The polynomial order to use when approximating.
[in]	nSteps	The number of sample points per axis (nSteps^2 total points).

Definition at line 480 of file Gtransfo.cc.

                                                                                                     {
    jointcal::StarMatchList starMatchList;
    double xStart = domain.xMin;
    double yStart = domain.yMin;
    double xStep = domain.getWidth() / (nSteps + 1);
    double yStep = domain.getHeight() / (nSteps + 1);
    for (unsigned i = 0; i < nSteps; ++i) {
        for (unsigned j = 0; j < nSteps; ++j) {
            // TODO: once DM-4044 is done, we can remove the redundancy in `Point`/`Point2D` here
            jointcal::Point in(xStart + i * xStep, yStart + j * yStep);
            afw::geom::Point2D inAfw(in.x, in.y);
            afw::geom::Point2D outAfw = transform->applyForward(inAfw);
            jointcal::Point out(outAfw.getX(), outAfw.getY());
            starMatchList.emplace_back(in, out, nullptr, nullptr);
        }
    }
    GtransfoPoly poly(order);
    poly.fit(starMatchList);
    *this = poly;
}

Member Function Documentation

__str__()

std::string lsst::jointcal::Gtransfo::__str__ ( )

inlineinherited

Definition at line 94 of file Gtransfo.h.

                        {
        std::stringstream s;
        dump(s);
        return s.str();
    }

apply() [1/4]

void lsst::jointcal::Gtransfo::apply ( Point const &  in, Point &  out  ) const

inlineinherited

applies the tranfo to in and writes into out. Is indeed virtual.

Definition at line 71 of file Gtransfo.h.

{ apply(in.x, in.y, out.x, out.y); }

apply() [2/4]

Point lsst::jointcal::Gtransfo::apply ( Point const &  in ) const

inlineinherited

All these apply(..) shadow the virtual one in derived classes, unless one writes "using Gtransfo::apply".

Definition at line 75 of file Gtransfo.h.

                                       {
        double xout, yout;
        apply(in.x, in.y, xout, yout);
        return Point(xout, yout);
    }

apply() [3/4]

Frame lsst::jointcal::Gtransfo::apply ( Frame const &  inputframe, bool  inscribed  ) const

inherited

Transform a bounding box, taking either the inscribed or circumscribed box.

Parameters

[in]	inputframe	The frame to be transformed.
[in]	inscribed	Return the inscribed (true) or circumscribed (false) box.

Returns

The transformed frame.

Definition at line 74 of file Gtransfo.cc.

                                                                   {
    // 2 opposite corners
    double xtmin1, xtmax1, ytmin1, ytmax1;
    apply(inputframe.xMin, inputframe.yMin, xtmin1, ytmin1);
    apply(inputframe.xMax, inputframe.yMax, xtmax1, ytmax1);
    Frame fr1(std::min(xtmin1, xtmax1), std::min(ytmin1, ytmax1), std::max(xtmin1, xtmax1),
              std::max(ytmin1, ytmax1));
    // 2 other corners
    double xtmin2, xtmax2, ytmin2, ytmax2;
    apply(inputframe.xMin, inputframe.yMax, xtmin2, ytmax2);
    apply(inputframe.xMax, inputframe.yMin, xtmax2, ytmin2);
    Frame fr2(std::min(xtmin2, xtmax2), std::min(ytmin2, ytmax2), std::max(xtmin2, xtmax2),
              std::max(ytmin2, ytmax2));

    if (inscribed) return fr1 * fr2;
    return fr1 + fr2;
}

apply() [4/4]

void lsst::jointcal::GtransfoPoly::apply ( const double  xIn, const double  yIn, double &  xOut, double &  yOut  ) const

overridevirtual

Implements lsst::jointcal::Gtransfo.

Definition at line 546 of file Gtransfo.cc.

                                                                                             {
    /*
      This routine computes the monomials only once for both
      polynomials.  This is why GtransfoPoly does not use an auxilary
      class (such as PolyXY) to handle each polynomial.

      The code works even if &xIn == &xOut (or &yIn == &yOut)
      It uses Variable Length Allocation (VLA) rather than a vector<double>
      because allocating the later costs about 50 ns. All VLA uses are tagged.
    */
    double monomials[_nterms];  // this is VLA, which is (perhaps) not casher C++
    computeMonomials(xIn, yIn, monomials);

    xOut = 0;
    yOut = 0;
    const double *c = &_coeffs[0];
    const double *pm = &monomials[0];
    // the ordering of the coefficients and the monomials are identical.
    for (int k = _nterms; k--;) xOut += (*(pm++)) * (*(c++));
    pm = &monomials[0];
    for (int k = _nterms; k--;) yOut += (*(pm++)) * (*(c++));
}

clone()

std::unique_ptr<Gtransfo> lsst::jointcal::GtransfoPoly::clone ( ) const

inlineoverridevirtual

returns a copy (allocated by new) of the transformation.

Implements lsst::jointcal::Gtransfo.

Reimplemented in lsst::jointcal::GtransfoLin.

Definition at line 337 of file Gtransfo.h.

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

coeff() [1/2]

double lsst::jointcal::GtransfoPoly::coeff	(	const unsigned	powX,
		const unsigned	powY,
		const unsigned	whichCoord
	)		const

access to coefficients (read only)

Definition at line 711 of file Gtransfo.cc.

                                                                                                    {
    assert((degX + degY <= _order) && whichCoord < 2);
    /* this assertion above is enough to ensure that the index used just
       below is within bounds since the reserved length is
       2*_nterms=(order+1)*(order+2) */
    return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
}

coeff() [2/2]

double & lsst::jointcal::GtransfoPoly::coeff	(	const unsigned	powX,
		const unsigned	powY,
		const unsigned	whichCoord
	)

write access

Definition at line 719 of file Gtransfo.cc.

                                                                                               {
    assert((degX + degY <= _order) && whichCoord < 2);
    return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
}

coeffOrZero()

double lsst::jointcal::GtransfoPoly::coeffOrZero	(	const unsigned	powX,
		const unsigned	powY,
		const unsigned	whichCoord
	)		const

read access, zero if beyond order

Definition at line 724 of file Gtransfo.cc.

                                                                                                          {
    //  assert((degX+degY<=order) && whichCoord<2);
    assert(whichCoord < 2);
    if (degX + degY <= _order)
        return _coeffs[(degX + degY) * (degX + degY + 1) / 2 + degY + whichCoord * _nterms];
    return 0;
}

composeAndReduce() [1/2]

std::unique_ptr< Gtransfo > lsst::jointcal::Gtransfo::composeAndReduce ( Gtransfo const &  right ) const

virtualinherited

Return a reduced composition of newTransfo = this(right()), or nullptr if it cannot be reduced.

"Reduced" in this context means that they are capable of being merged into a single transform, for example, for two polynomials:

\[ f(x) = 1 + x^2, g(x) = -1 + 3x \]

we would have h = f.composeAndReduce(g) == 2 - 6x + 9x^2.

To be overloaded by derived classes if they can properly reduce the composition.

Parameters

right

The transform to apply first.

Returns

The new reduced and composed gtransfo, or nullptr if no such reduction is possible.

Reimplemented in lsst::jointcal::GtransfoIdentity.

Definition at line 92 of file Gtransfo.cc.

                                {  // by default no way to compose
    return std::unique_ptr<Gtransfo>(nullptr);
}

composeAndReduce() [2/2]

std::unique_ptr< Gtransfo > lsst::jointcal::GtransfoPoly::composeAndReduce

(

GtransfoPoly const &

right

)

const

Return a reduced composition of newTransfo = this(right()), or nullptr if it cannot be reduced.

"Reduced" in this context means that they are capable of being merged into a single transform, for example, for two polynomials:

\[ f(x) = 1 + x^2, g(x) = -1 + 3x \]

we would have h = f.composeAndReduce(g) == 2 - 6x + 9x^2.

To be overloaded by derived classes if they can properly reduce the composition.

Parameters

right

The transform to apply first.

Returns

The new reduced and composed gtransfo, or nullptr if no such reduction is possible.

Definition at line 898 of file Gtransfo.cc.

                                                                                      {
    if (getOrder() == 1 && right.getOrder() == 1)
        return std::make_unique<GtransfoLin>((*this) * (right));  // does the composition
    else
        return std::make_unique<GtransfoPoly>((*this) * (right));  // does the composition
}

computeDerivative()

void lsst::jointcal::GtransfoPoly::computeDerivative ( Point const &  where, GtransfoLin &  derivative, const double  step = 0.01  ) const

overridevirtual

specialised analytic routine

Reimplemented from lsst::jointcal::Gtransfo.

Reimplemented in lsst::jointcal::GtransfoLin.

Definition at line 569 of file Gtransfo.cc.

              { /* routine checked against numerical derivatives from Gtransfo::Derivative */
    if (_order == 1) {
        derivative = GtransfoLin(*this);
        derivative.dx() = derivative.dy() = 0;
        return;
    }

    double dermx[2 * _nterms];  // VLA
    double *dermy = dermx + _nterms;
    double xin = where.x;
    double yin = where.y;

    double xx = 1;
    double xxm1 = 1;  // xx^(ix-1)
    for (unsigned ix = 0; ix <= _order; ++ix) {
        unsigned k = (ix) * (ix + 1) / 2;
        // iy = 0
        dermx[k] = ix * xxm1;
        dermy[k] = 0;
        k += ix + 2;
        double yym1 = 1;  // yy^(iy-1)
        for (unsigned iy = 1; iy <= _order - ix; ++iy) {
            dermx[k] = ix * xxm1 * yym1 * yin;
            dermy[k] = iy * xx * yym1;
            yym1 *= yin;
            k += ix + iy + 2;
        }
        xx *= xin;
        if (ix >= 1) xxm1 *= xin;
    }

    derivative.dx() = 0;
    derivative.dy() = 0;

    const double *mx = &dermx[0];
    const double *my = &dermy[0];
    const double *c = &_coeffs[0];
    // dx'
    double a11 = 0, a12 = 0;
    for (int k = _nterms; k--;) {
        a11 += (*(mx++)) * (*c);
        a12 += (*(my++)) * (*(c++));
    }
    derivative.a11() = a11;
    derivative.a12() = a12;
    // dy'
    double a21 = 0, a22 = 0;
    mx = &dermx[0];
    my = &dermy[0];
    for (int k = _nterms; k--;) {
        a21 += (*(mx++)) * (*c);
        a22 += (*(my++)) * (*(c++));
    }
    derivative.a21() = a21;
    derivative.a22() = a22;
}

determinant()

double lsst::jointcal::GtransfoPoly::determinant

(

)

const

Definition at line 782 of file Gtransfo.cc.

                                       {
    if (_order >= 1) return coeff(1, 0, 0) * coeff(0, 1, 1) - coeff(0, 1, 0) * coeff(1, 0, 1);
    return 0;
}

dump()

void lsst::jointcal::GtransfoPoly::dump ( std::ostream &  stream = std::cout ) const

overridevirtual

print out of coefficients in a readable form.

Implements lsst::jointcal::Gtransfo.

Definition at line 763 of file Gtransfo.cc.

                                             {
    auto oldPrecision = stream.precision();
    stream.precision(12);
    for (unsigned ic = 0; ic < 2; ++ic) {
        if (ic == 0)
            stream << "newx = ";
        else
            stream << "newy = ";
        for (unsigned p = 0; p <= _order; ++p)
            for (unsigned py = 0; py <= p; ++py) {
                if (p + py != 0) stream << " + ";
                stream << coeff(p - py, py, ic) << monomialString(p - py, py);
            }
        stream << endl;
    }
    if (_order > 0) stream << " Linear determinant = " << determinant() << endl;
    stream.precision(oldPrecision);
}

fit()

double lsst::jointcal::GtransfoPoly::fit ( StarMatchList const &  starMatchList )

overridevirtual

guess what

Implements lsst::jointcal::Gtransfo.

Reimplemented in lsst::jointcal::GtransfoLinRot, and lsst::jointcal::GtransfoLinShift.

Definition at line 880 of file Gtransfo.cc.

                                                           {
    if (starMatchList.size() < _nterms) {
        LOGLS_FATAL(_log, "GtransfoPoly::fit trying to fit a polynomial transfo of order "
                                  << _order << " with only " << starMatchList.size() << " matches.");
        return -1;
    }

    GtransfoPoly conditionner = shiftAndNormalize(starMatchList);

    computeFit(starMatchList, conditionner, false);               // get a rough solution
    computeFit(starMatchList, conditionner, true);                // weight with it
    double chi2 = computeFit(starMatchList, conditionner, true);  // once more

    (*this) = (*this) * conditionner;
    if (starMatchList.size() == _nterms) return 0;
    return chi2;
}

getJacobian() [1/2]

virtual double lsst::jointcal::Gtransfo::getJacobian ( Point const &  point ) const

inlinevirtualinherited

returns the local jacobian.

Definition at line 111 of file Gtransfo.h.

{ return getJacobian(point.x, point.y); }

getJacobian() [2/2]

double lsst::jointcal::Gtransfo::getJacobian ( const double  x, const double  y  ) const

virtualinherited

returns the local jacobian.

Definition at line 97 of file Gtransfo.cc.

                                                                 {
    double x2, y2;
    double eps = x * 0.01;
    if (eps == 0) eps = 0.01;
    apply(x, y, x2, y2);
    double dxdx, dydx;
    apply(x + eps, y, dxdx, dydx);
    dxdx -= x2;
    dydx -= y2;
    double dxdy, dydy;
    apply(x, y + eps, dxdy, dydy);
    dxdy -= x2;
    dydy -= y2;
    return ((dxdx * dydy - dxdy * dydx) / (eps * eps));
}

getNpar()

int lsst::jointcal::GtransfoPoly::getNpar ( ) const

inlineoverridevirtual

total number of parameters

Reimplemented from lsst::jointcal::Gtransfo.

Reimplemented in lsst::jointcal::GtransfoLinScale, lsst::jointcal::GtransfoLinRot, and lsst::jointcal::GtransfoLinShift.

Definition at line 316 of file Gtransfo.h.

{ return 2 * _nterms; }

getOrder()

unsigned lsst::jointcal::GtransfoPoly::getOrder ( ) const

inline

Returns the polynomial order.

Definition at line 302 of file Gtransfo.h.

{ return _order; }

getParams()

void lsst::jointcal::Gtransfo::getParams ( double *  params ) const

inherited

params should be at least Npar() long

Definition at line 210 of file Gtransfo.cc.

                                             {
    int npar = getNpar();
    for (int i = 0; i < npar; ++i) params[i] = paramRef(i);
}

inverseTransfo()

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

virtualinherited

returns an inverse transfo. Numerical if not overloaded.

precision and region refer to the "input" side of this, and hence to the output side of the returned Gtransfo.

Reimplemented in lsst::jointcal::TanRaDecToPixel, lsst::jointcal::TanSipPixelToRaDec, lsst::jointcal::TanPixelToRaDec, lsst::jointcal::GtransfoLin, and lsst::jointcal::GtransfoInverse.

Definition at line 292 of file Gtransfo.cc.

                                                                                                  {
    return std::unique_ptr<Gtransfo>(new GtransfoInverse(this, precision, region));
}

linearApproximation()

GtransfoLin lsst::jointcal::Gtransfo::linearApproximation ( Point const &  where, const double  step = 0.01  ) const

virtualinherited

linear (local) approximation.

Reimplemented in lsst::jointcal::GtransfoLin, and lsst::jointcal::GtransfoIdentity.

Definition at line 133 of file Gtransfo.cc.

                                                                                     {
    Point outwhere = apply(where);
    GtransfoLin der;
    computeDerivative(where, der, step);
    return GtransfoLinShift(outwhere.x, outwhere.y) * der * GtransfoLinShift(-where.x, -where.y);
}

offsetParams()

void lsst::jointcal::Gtransfo::offsetParams ( Eigen::VectorXd const &  delta )

inherited

Definition at line 215 of file Gtransfo.cc.

                                                      {
    int npar = getNpar();
    for (int i = 0; i < npar; ++i) paramRef(i) += delta[i];
}

operator*()

GtransfoPoly lsst::jointcal::GtransfoPoly::operator*

(

GtransfoPoly const &

right

)

const

Composition (internal stuff in quadruple precision)

Definition at line 1003 of file Gtransfo.cc.

                                                                    {
    // split each transfo into 2d polynomials
    PolyXY plx(*this, 0);
    PolyXY ply(*this, 1);
    PolyXY prx(right, 0);
    PolyXY pry(right, 1);

    // compute the compositions
    PolyXY rx(composition(plx, prx, pry));
    PolyXY ry(composition(ply, prx, pry));

    // copy the results the hard way.
    GtransfoPoly result(_order * right._order);
    for (unsigned px = 0; px <= result._order; ++px)
        for (unsigned py = 0; py <= result._order - px; ++py) {
            result.coeff(px, py, 0) = rx.coeff(px, py);
            result.coeff(px, py, 1) = ry.coeff(px, py);
        }
    return result;
}

operator+()

GtransfoPoly lsst::jointcal::GtransfoPoly::operator+

(

GtransfoPoly const &

right

)

const

Addition.

Definition at line 1024 of file Gtransfo.cc.

                                                                    {
    if (_order >= right._order) {
        GtransfoPoly res(*this);
        for (unsigned i = 0; i <= right._order; ++i)
            for (unsigned j = 0; j <= right._order - i; ++j) {
                res.coeff(i, j, 0) += right.coeff(i, j, 0);
                res.coeff(i, j, 1) += right.coeff(i, j, 1);
            }
        return res;
    } else
        return (right + (*this));
}

operator-()

GtransfoPoly lsst::jointcal::GtransfoPoly::operator-

(

GtransfoPoly const &

right

)

const

Subtraction.

Definition at line 1037 of file Gtransfo.cc.

                                                                    {
    GtransfoPoly res(std::max(_order, right._order));
    for (unsigned i = 0; i <= res._order; ++i)
        for (unsigned j = 0; j <= res._order - i; ++j) {
            res.coeff(i, j, 0) = coeffOrZero(i, j, 0) - right.coeffOrZero(i, j, 0);
            res.coeff(i, j, 1) = coeffOrZero(i, j, 1) - right.coeffOrZero(i, j, 1);
        }
    return res;
}

paramDerivatives()

void lsst::jointcal::GtransfoPoly::paramDerivatives ( Point const &  where, double *  dx, double *  dy  ) const

overridevirtual

Derivative w.r.t parameters.

Derivatives should be al least 2*NPar long. first Npar, for x, last Npar for y.

Reimplemented from lsst::jointcal::Gtransfo.

Definition at line 743 of file Gtransfo.cc.

              { /* first half : dxout/dpar, second half : dyout/dpar */
    computeMonomials(where.x, where.y, dx);
    for (unsigned k = 0; k < _nterms; ++k) {
        dy[_nterms + k] = dx[k];
        dx[_nterms + k] = dy[k] = 0;
    }
}

paramRef() [1/2]

double lsst::jointcal::GtransfoPoly::paramRef ( const int  i ) const

overridevirtual

Reimplemented from lsst::jointcal::Gtransfo.

Definition at line 733 of file Gtransfo.cc.

                                               {
    assert(unsigned(i) < 2 * _nterms);
    return _coeffs[i];
}

paramRef() [2/2]

double & lsst::jointcal::GtransfoPoly::paramRef ( const int  i )

overridevirtual

Reimplemented from lsst::jointcal::Gtransfo.

Definition at line 738 of file Gtransfo.cc.

                                          {
    assert(unsigned(i) < 2 * _nterms);
    return _coeffs[i];
}

read()

void lsst::jointcal::GtransfoPoly::read

(

std::istream &

s

)

Definition at line 1062 of file Gtransfo.cc.

                                  {
    int format;
    s >> format;
    if (format != 1)
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, " GtransfoPoly::read : format is not 1 ");

    string order;
    s >> order >> _order;
    if (order != "order")
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          " GtransfoPoly::read : expecting \"order\" and found " + order);
    setOrder(_order);
    for (unsigned k = 0; k < 2 * _nterms; ++k) s >> _coeffs[k];
}

roughInverse()

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

virtualinherited

Rough inverse.

Stored by the numerical inverter to guess starting point for the trials. Just here to enable overloading.

Reimplemented in lsst::jointcal::TanRaDecToPixel, lsst::jointcal::TanPixelToRaDec, and lsst::jointcal::GtransfoInverse.

Definition at line 190 of file Gtransfo.cc.

                                                                        {
    // "in" and "out" refer to the inverse direction.
    Point centerOut = region.getCenter();
    Point centerIn = apply(centerOut);
    GtransfoLin der;
    computeDerivative(centerOut, der, std::sqrt(region.getArea()) / 5.);
    der = der.inverted();
    der = GtransfoLinShift(centerOut.x, centerOut.y) * der * GtransfoLinShift(-centerIn.x, -centerIn.y);
    return std::unique_ptr<Gtransfo>(new GtransfoLin(der));
}

setOrder()

void lsst::jointcal::GtransfoPoly::setOrder

(

const unsigned

order

)

Sets the polynomial order (the highest sum of exponents of the largest monomial).

Definition at line 525 of file Gtransfo.cc.

                                                {
    _order = order;
    unsigned old_nterms = _nterms;
    _nterms = (_order + 1) * (_order + 2) / 2;

    // temporarily save coefficients
    vector<double> old_coeffs = _coeffs;
    // reallocate enough size
    _coeffs.resize(2 * _nterms);
    // reassign to zero (this is necessary because ycoeffs
    // are after xcoeffs and so their meaning changes
    for (unsigned k = 0; k < _nterms; ++k) _coeffs[k] = 0;
    // put back what we had before
    unsigned kmax = min(old_nterms, _nterms);
    for (unsigned k = 0; k < kmax; ++k) {
        _coeffs[k] = old_coeffs[k];                         // x terms
        _coeffs[k + _nterms] = old_coeffs[k + old_nterms];  // y terms
    }
}

toAstMap()

std::shared_ptr< ast::Mapping > lsst::jointcal::GtransfoPoly::toAstMap ( jointcal::Frame const &  domain ) const

overridevirtual

Create an equivalent AST mapping for this transformation, including an analytic inverse if possible.

Parameters

domain

The domain of the transfo, to help find an inverse.

Returns

An AST Mapping that represents this transformation.

Reimplemented from lsst::jointcal::Gtransfo.

Definition at line 1047 of file Gtransfo.cc.

                                                                                  {
    auto inverse = inversePolyTransfo(*this, domain, 1e-7, _order + 2, 100);
    return std::make_shared<ast::PolyMap>(toAstPolyMapCoefficients(), inverse->toAstPolyMapCoefficients());
}

transformErrors()

void lsst::jointcal::Gtransfo::transformErrors ( Point const &  where, const double *  vIn, double *  vOut  ) const

virtualinherited

transform errors (represented as double[3] in order V(xx),V(yy),Cov(xy))

Definition at line 157 of file Gtransfo.cc.

                                                                                        {
    GtransfoLin der;
    computeDerivative(where, der, 0.01);
    double a11 = der.A11();
    double a22 = der.A22();
    double a21 = der.A21();
    double a12 = der.A12();

    /*      (a11 a12)          (vxx  vxy)
       M =  (       )  and V = (        )
            (a21 a22)          (xvy  vyy)

       Vxx = Vin[0], vyy = Vin[1], Vxy = Vin[2];
       we want to compute M*V*tp(M)
       A lin alg light package would be perfect...
    */
    int xx = 0;
    int yy = 1;
    int xy = 2;
    // M*V :

    double b11 = a11 * vIn[xx] + a12 * vIn[xy];
    double b22 = a21 * vIn[xy] + a22 * vIn[yy];
    double b12 = a11 * vIn[xy] + a12 * vIn[yy];
    double b21 = a21 * vIn[xx] + a22 * vIn[xy];

    // (M*V) * tp(M)

    vOut[xx] = b11 * a11 + b12 * a12;
    vOut[xy] = b11 * a21 + b12 * a22;
    vOut[yy] = b21 * a21 + b22 * a22;
}

transformPosAndErrors()

void lsst::jointcal::GtransfoPoly::transformPosAndErrors ( const FatPoint &  in, FatPoint &  out  ) const

overridevirtual

a mix of apply and Derivative

Reimplemented from lsst::jointcal::Gtransfo.

Definition at line 627 of file Gtransfo.cc.

                                                                                {
    /*
       The results from this routine were compared to what comes out
       from apply and transformErrors. The Derivative routine was
       checked against numerical derivatives from
       Gtransfo::Derivative. (P.A dec 2009).

       This routine could be made much simpler by calling apply and
       Derivative (i.e. you just suppress it, and the fallback is the
       generic version in Gtransfo).  BTW, I checked that both routines
       provide the same result. This version is however faster
       (monomials get recycled).
    */
    double monomials[_nterms];  // VLA

    FatPoint res;  // to store the result, because nothing forbids &in == &out.

    double dermx[2 * _nterms];        // monomials for derivative w.r.t. x (VLA)
    double *dermy = dermx + _nterms;  // same for y
    double xin = in.x;
    double yin = in.y;

    double xx = 1;
    double xxm1 = 1;  // xx^(ix-1)
    for (unsigned ix = 0; ix <= _order; ++ix) {
        unsigned k = (ix) * (ix + 1) / 2;
        // iy = 0
        dermx[k] = ix * xxm1;
        dermy[k] = 0;
        monomials[k] = xx;
        k += ix + 2;
        double yy = yin;
        double yym1 = 1;  // yy^(iy-1)
        for (unsigned iy = 1; iy <= _order - ix; ++iy) {
            monomials[k] = xx * yy;
            dermx[k] = ix * xxm1 * yy;
            dermy[k] = iy * xx * yym1;
            yym1 *= yin;
            yy *= yin;
            k += ix + iy + 2;
        }
        xx *= xin;
        if (ix >= 1) xxm1 *= xin;
    }

    // output position
    double xout = 0, yout = 0;
    const double *c = &_coeffs[0];
    const double *pm = &monomials[0];
    for (int k = _nterms; k--;) xout += (*(pm++)) * (*(c++));
    pm = &monomials[0];
    for (int k = _nterms; k--;) yout += (*(pm++)) * (*(c++));
    res.x = xout;
    res.y = yout;

    // derivatives
    c = &_coeffs[0];
    const double *mx = &dermx[0];
    const double *my = &dermy[0];
    double a11 = 0, a12 = 0;
    for (int k = _nterms; k--;) {
        a11 += (*(mx++)) * (*c);
        a12 += (*(my++)) * (*(c++));
    }

    double a21 = 0, a22 = 0;
    mx = &dermx[0];
    my = &dermy[0];
    for (int k = _nterms; k--;) {
        a21 += (*(mx++)) * (*c);
        a22 += (*(my++)) * (*(c++));
    }

    // output co-variance
    res.vx = a11 * (a11 * in.vx + 2 * a12 * in.vxy) + a12 * a12 * in.vy;
    res.vy = a21 * a21 * in.vx + a22 * a22 * in.vy + 2. * a21 * a22 * in.vxy;
    res.vxy = a21 * a11 * in.vx + a22 * a12 * in.vy + (a21 * a12 + a11 * a22) * in.vxy;
    out = res;
}

transformStar()

void lsst::jointcal::Gtransfo::transformStar ( FatPoint &  in ) const

inlineinherited

allows to write MyTransfo(MyStar)

Definition at line 108 of file Gtransfo.h.

{ transformPosAndErrors(in, in); }

write() [1/2]

void lsst::jointcal::Gtransfo::write ( const std::string &  fileName ) const

inherited

Definition at line 239 of file Gtransfo.cc.

                                                    {
    ofstream s(fileName.c_str());
    write(s);
    bool ok = !s.fail();
    s.close();
    if (!ok)
        throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                          "Gtransfo::write, something went wrong for file " + fileName);
}

write() [2/2]

void lsst::jointcal::GtransfoPoly::write ( std::ostream &  s ) const

overridevirtual

Reimplemented from lsst::jointcal::Gtransfo.

Definition at line 1052 of file Gtransfo.cc.

                                         {
    s << " GtransfoPoly 1" << endl;
    s << "order " << _order << endl;
    int oldprec = s.precision();
    s << setprecision(12);
    for (unsigned k = 0; k < 2 * _nterms; ++k) s << _coeffs[k] << ' ';
    s << endl;
    s << setprecision(oldprec);
}

---

The documentation for this class was generated from the following files:

• /j/snowflake/release/lsstsw/stack/Linux64/jointcal/16.0-25-g07af9f2+1/include/lsst/jointcal/Gtransfo.h
• /j/snowflake/release/lsstsw/stack/Linux64/jointcal/16.0-25-g07af9f2+1/src/Gtransfo.cc