lsst::meas.astrom.sip::CreateWcsWithSip< MatchT > Class Template Reference

Measure the distortions in an image plane and express them a SIP polynomials. More...

#include <CreateWcsWithSip.h>

Public Types
typedef boost::shared_ptr < CreateWcsWithSip >	Ptr
typedef boost::shared_ptr < CreateWcsWithSip const >	ConstPtr

Public Member Functions
	CreateWcsWithSip (std::vector< MatchT > const &matches, afw::image::Wcs const &linearWcs, int const order, afw::geom::Box2I const &bbox=afw::geom::Box2I(), int const ngrid=0)
	Constructor. More...
boost::shared_ptr < afw::image::TanWcs >	getNewWcs ()
double	getScatterInPixels () const
afw::geom::Angle	getScatterOnSky () const
double	getLinearScatterInPixels () const
afw::geom::Angle	getLinearScatterOnSky () const
int	getOrder () const
	Return the number of terms in the SIP matrix. More...
int	getNPoints () const
	Return the number of points in the catalogue. More...
int	getNGrid () const
	Return the number of grid points (on each axis) used in inverse SIP transform. More...

Private Member Functions
void	_calculateForwardMatrices ()
void	_calculateReverseMatrices ()
afw::geom::Point2D	_getCrvalAsGeomPoint () const

Private Attributes
lsst::pex::logging::Log	_log
std::vector< MatchT > const	_matches
afw::geom::Box2I	_bbox
int	_ngrid
boost::shared_ptr < afw::image::Wcs const >	_linearWcs
int const	_sipOrder
int const	_reverseSipOrder
Eigen::MatrixXd	_sipA
Eigen::MatrixXd	_sipB
Eigen::MatrixXd	_sipAp
Eigen::MatrixXd	_sipBp
boost::shared_ptr < afw::image::TanWcs >	_newWcs

Detailed Description

template<class MatchT>
class lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >

Measure the distortions in an image plane and express them a SIP polynomials.

Given a list of matching sources between a catalogue and an image, and a linear Wcs that describes the mapping from pixel space in the image and ra/dec space in the catalogue, calculate discrepancies between the two and compute SIP distortion polynomials to describe the discrepancy

SIP polynomials are defined in Shupe at al. (2005) ASPC 347 491.

Note that the SIP standard insists (although it is only mentioned obliquly between Eqns 3 and 4) that the lowest three terms in the distortion polynomials be zero (A00, A10, A01, B00, etc.). To achieve this, we need to adjust the values of CD and CRPIX from the input wcs. This may not be the behaviour you expect.

A Wcs may be created in a variety of ways (e.g. lsst::meas::astrom::net::GlobalAstrometrySolution ), and the list of matched sources (matches) can be generated with the matchRaDec function.

#Example usage 
matches = matchRaDec(catSet, srcSet, 1.0*afwGeom.arcseconds, true)
wcs = getWcsFromSomewhere()

maxScatter=0.1
maxOrder= 10
sipObject = CreateWcsWithSip(matches, wcs, maxScatter, maxOrder)
wcs = sipObject.getNewWcs()

Note that the matches must be one-to-one; this is ensured by passing closest=true to matchRaDec.

Definition at line 87 of file CreateWcsWithSip.h.

Member Typedef Documentation

template<class MatchT>

typedef boost::shared_ptr<CreateWcsWithSip const> lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::ConstPtr

Definition at line 91 of file CreateWcsWithSip.h.

template<class MatchT>

typedef boost::shared_ptr<CreateWcsWithSip> lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::Ptr

Definition at line 90 of file CreateWcsWithSip.h.

Constructor & Destructor Documentation

template<class MatchT >

lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::CreateWcsWithSip	(	std::vector< MatchT > const &	matches,
		afw::image::Wcs const &	linearWcs,
		int const	order,
		afw::geom::Box2I const &	bbox = afw::geom::Box2I(),
		int const	ngrid = 0
	)

Constructor.

Construct a CreateWcsWithSip

Parameters

[in]	matches	list of matches
[in]	linearWcs	initial WCS, typically pure TAN but need not be
[in]	order	SIP order for fit WCS
[in]	bbox	bounding box over which to compute the reverse SIP transform. If empty then a bounding box is computed based on the matches, extended a bit to allow for the fact that the sources will not necessarily reach to each edge of the image. Specifially the box is grown by dimensions/sqrt(number of matches).
[in]	ngrid	number of points along x or y for the grid of points on which the reverse SIP transform is computed

Definition at line 107 of file CreateWcsWithSip.cc.

 :
    _log(pexLog::Log(pexLog::Log::getDefaultLog(), "meas.astrom.sip")),
    _matches(matches),
    _bbox(bbox),
    _ngrid(ngrid),
    _linearWcs(linearWcs.clone()),
    _sipOrder(order+1),
    _reverseSipOrder(order+2), //Higher order for reverse transform
    _sipA(Eigen::MatrixXd::Zero(_sipOrder, _sipOrder)),
    _sipB(Eigen::MatrixXd::Zero(_sipOrder, _sipOrder)),
    _sipAp(Eigen::MatrixXd::Zero(_reverseSipOrder, _reverseSipOrder)),
    _sipBp(Eigen::MatrixXd::Zero(_reverseSipOrder, _reverseSipOrder)),
    _newWcs()
{
    if  (order < 2) {
        throw LSST_EXCEPT(except::OutOfRangeError, "SIP must be at least 2nd order");        
    }
    if (_sipOrder > 9) {
        throw LSST_EXCEPT(except::OutOfRangeError,
                          str(boost::format("SIP forward order %d exceeds the convention limit of 9") %
                              _sipOrder));
    }
    if (_reverseSipOrder > 9) {
        throw LSST_EXCEPT(except::OutOfRangeError,
                          str(boost::format("SIP reverse order %d exceeds the convention limit of 9") %
                              _reverseSipOrder));
    }
    
    if (_matches.size() < std::size_t(_sipOrder)) {
        throw LSST_EXCEPT(except::LengthError, "Number of matches less than requested sip order");
    }

    if (_ngrid <= 0) {
        _ngrid = 5*_sipOrder;           // should be plenty
    }

    /*
     * We need a bounding box to define the region over which:
     *    The forward transformation should be valid
     *    We calculate the reverse transformartion
     * If no BBox is provided, guess one from the input points (extrapolated a bit to allow for fact
     * that a finite number of points won't reach to the edge of the image)
     */
    if (_bbox.isEmpty() && !_matches.empty() > 0) {
        for (
            typename std::vector<MatchT>::const_iterator ptr = _matches.begin();
            ptr != _matches.end();
            ++ptr
        ) {
            afwTable::SourceRecord const & src = *ptr->second;
            _bbox.include(afwGeom::PointI(src.getX(), src.getY()));
        }
        float const borderFrac = 1/::sqrt(_matches.size()); // fractional border to add to exact BBox
        afwGeom::Extent2I border(borderFrac*_bbox.getWidth(), borderFrac*_bbox.getHeight());

        _bbox.grow(border);
    }
    // Calculate the forward part of the SIP distortion
    _calculateForwardMatrices();

    //Build a new wcs incorporating the forward sip matrix, it's all we know so far
    afwGeom::Point2D crval = _getCrvalAsGeomPoint();
    afwGeom::Point2D crpix = _linearWcs->getPixelOrigin();
    Eigen::MatrixXd CD = _linearWcs->getCDMatrix();
    
    _newWcs = PTR(afwImg::TanWcs)(new afwImg::TanWcs(crval, crpix, CD, _sipA, _sipB, _sipAp, _sipBp));
    // Use _newWcs to calculate the forward transformation on a grid, and derive the back transformation
    _calculateReverseMatrices();

    //Build a new wcs incorporating both of the sip matrices
    _newWcs = PTR(afwImg::TanWcs)(new afwImg::TanWcs(crval, crpix, CD, _sipA, _sipB, _sipAp, _sipBp));

}

Member Function Documentation

template<class MatchT >

void lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_calculateForwardMatrices ( )

private

Definition at line 189 of file CreateWcsWithSip.cc.

{
    // Assumes FITS (1-indexed) coordinates.
    afwGeom::Point2D crpix = _linearWcs->getPixelOrigin();

    // Calculate u, v and intermediate world coordinates
    int const nPoints = _matches.size();
    Eigen::VectorXd u(nPoints), v(nPoints), iwc1(nPoints), iwc2(nPoints);

    int i = 0;
    for (
        typename std::vector<MatchT>::const_iterator ptr = _matches.begin();
        ptr != _matches.end();
        ++ptr, ++i
    ) {
        afwTable::ReferenceMatch const & match = *ptr;

        // iwc: intermediate world coordinate positions of catalogue objects
        afwCoord::IcrsCoord c = match.first->getCoord();
        afwGeom::Point2D p = _linearWcs->skyToIntermediateWorldCoord(c);
        iwc1[i] = p[0];
        iwc2[i] = p[1];
        // u and v are intermediate pixel coordinates of observed (distorted) positions
        u[i] = match.second->getX() - crpix[0];
        v[i] = match.second->getY() - crpix[1];
    }
    // Scale u and v down to [-1,,+1] in order to avoid too large numbers in the polynomials
    double uMax = u.cwiseAbs().maxCoeff();
    double vMax = v.cwiseAbs().maxCoeff();
    double norm = std::max(uMax, vMax);
    u = u/norm;
    v = v/norm;
   
    // Forward transform
    int ord = _sipOrder;
    Eigen::MatrixXd forwardC = calculateCMatrix(u, v, ord);
    Eigen::VectorXd mu = leastSquaresSolve(iwc1, forwardC);
    Eigen::VectorXd nu = leastSquaresSolve(iwc2, forwardC);

    // Use mu and nu to refine CD

    // Given the implementation of indexToPQ(), the refined values
    // of the elements of the CD matrices are in elements 1 and "_sipOrder" of mu and nu
    // If the implementation of indexToPQ() changes, these assertions
    // will catch that change.
    assert ((indexToPQ(0,   ord) == std::pair<int, int>(0, 0)));
    assert ((indexToPQ(1,   ord) == std::pair<int, int>(0, 1)));
    assert ((indexToPQ(ord, ord) == std::pair<int, int>(1, 0)));

    // Scale back CD matrix
    Eigen::Matrix2d CD;
    CD(1,0) = nu[ord]/norm;
    CD(1,1) = nu[1]/norm;
    CD(0,0) = mu[ord]/norm;
    CD(0,1) = mu[1]/norm;

    Eigen::Matrix2d CDinv = CD.inverse();   //Direct inverse OK for 2x2 matrix in Eigen

    // The zeroth elements correspond to a shift in crpix
    crpix[0] -= mu[0]*CDinv(0,0) + nu[0]*CDinv(0,1); 
    crpix[1] -= mu[0]*CDinv(1,0) + nu[0]*CDinv(1,1);

    afwGeom::Point2D crval = _getCrvalAsGeomPoint();

    _linearWcs = afwImg::Wcs::Ptr( new afwImg::Wcs(crval, crpix, CD));

    //Get Sip terms
    
    //The rest of the elements correspond to
    //mu[i] == CD11*Apq + CD12*Bpq and
    //nu[i] == CD21*Apq + CD22*Bpq and
    //
    //We solve for Apq and Bpq with the equation
    // (Apq)  = (CD11 CD12)-1  * (mu[i])  
    // (Bpq)    (CD21 CD22)      (nu[i])
    
    for(int i=1; i< mu.rows(); ++i) {
        std::pair<int, int> pq = indexToPQ(i, ord);
        int p = pq.first, q = pq.second;

        if(p + q > 1 && p + q < ord) {    
            Eigen::Vector2d munu(2,1);
            munu(0) = mu(i);
            munu(1) = nu(i);
            Eigen::Vector2d AB = CDinv*munu;
            // Scale back sip coefficients
            _sipA(p,q) = AB[0]/::pow(norm,p+q);
            _sipB(p,q) = AB[1]/::pow(norm,p+q);
        }
    }
}

template<class MatchT >

void lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_calculateReverseMatrices ( )

private

Definition at line 282 of file CreateWcsWithSip.cc.

                                                         {
    int const ngrid2 = _ngrid*_ngrid;

    Eigen::VectorXd U(ngrid2), V(ngrid2);
    Eigen::VectorXd delta1(ngrid2), delta2(ngrid2);
    
    int const x0 = _bbox.getMinX();
    double const dx = _bbox.getWidth()/(double)(_ngrid - 1);
    int const y0 = _bbox.getMinY();
    double const dy = _bbox.getHeight()/(double)(_ngrid - 1);

    // wcs->getPixelOrigin() returns LSST-style (0-indexed) pixel coords.
    afwGeom::Point2D crpix = _newWcs->getPixelOrigin();

    _log.debugf("_calcReverseMatrices: x0,y0 %i,%i, W,H %i,%i, ngrid %i, dx,dy %g,%g, CRPIX %g,%g",
                x0, y0, _bbox.getWidth(), _bbox.getHeight(), _ngrid, dx, dy, crpix[0], crpix[1]);

    int k = 0;
    for (int i = 0; i < _ngrid; ++i) {
        double const y = y0 + i*dy;
        for (int j = 0; j < _ngrid; ++j, ++k) {
            double const x = x0 + j*dx;
            double u,v;
            // u and v are intermediate pixel coordinates on a grid of positions
            u = x - crpix[0];
            v = y - crpix[1];

            // U and V are the result of applying the "forward" (A,B) SIP coefficients
            // NOTE that the "undistortPixel()" function accepts 1-indexed (FITS-style)
            // coordinates, and here we are treating "x" and "y" as LSST-style.
            afwGeom::Point2D xy = _newWcs->undistortPixel(afwGeom::Point2D(x + 1, y + 1));
            // "crpix", on the other hand, is LSST-style 0-indexed, so we have to remove
            // the FITS-style 1-index from "xy"
            U[k] = xy[0] - 1 - crpix[0];
            V[k] = xy[1] - 1 - crpix[1];

            if ((i == 0 || i == (_ngrid-1) || i == (_ngrid/2)) &&
                (j == 0 || j == (_ngrid-1) || j == (_ngrid/2))) {
                _log.debugf("  x,y (%.1f, %.1f), u,v (%.1f, %.1f), U,V (%.1f, %.1f)", x, y, u, v, U[k], V[k]);
            }

            delta1[k] = u - U[k];
            delta2[k] = v - V[k];
        }
    }
    
    // Scale down U and V in order to avoid too large numbers in the polynomials
    double UMax = U.cwiseAbs().maxCoeff();
    double VMax = V.cwiseAbs().maxCoeff();
    double norm = (UMax > VMax) ? UMax : VMax;
    U = U/norm;
    V = V/norm;

    // Reverse transform
    int const ord = _reverseSipOrder;
    Eigen::MatrixXd reverseC = calculateCMatrix(U, V, ord); 
    Eigen::VectorXd tmpA = leastSquaresSolve(delta1, reverseC);
    Eigen::VectorXd tmpB = leastSquaresSolve(delta2, reverseC);
    
    assert(tmpA.rows() == tmpB.rows());
    for(int j=0; j< tmpA.rows(); ++j) {
        std::pair<int, int> pq = indexToPQ(j, ord);
        int p = pq.first, q = pq.second;
        // Scale back sip coefficients
        _sipAp(p, q) = tmpA[j]/::pow(norm,p+q);
        _sipBp(p, q) = tmpB[j]/::pow(norm,p+q);   
    } 
}

template<class MatchT >

afwGeom::Point2D lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_getCrvalAsGeomPoint ( ) const

private

Definition at line 378 of file CreateWcsWithSip.cc.

                                                                    {
    afwCoord::Fk5Coord coo = _linearWcs->getSkyOrigin()->toFk5();
    return coo.getPosition(afwGeom::degrees);
}

template<class MatchT >

double lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getLinearScatterInPixels

(

)

const

Compute the median radial separation between items in this object's match list

For each match, project the reference object coord to pixels using the initial "linearWcs" WCS, and measure the radial separation to the source centroid.

Definition at line 358 of file CreateWcsWithSip.cc.

                                                                {
    assert(_linearWcs.get());
    return makeMatchStatisticsInPixels(*_linearWcs, _matches, afw::math::MEDIAN).getValue();
}

template<class MatchT >

afwGeom::Angle lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getLinearScatterOnSky

(

)

const

Compute the median on-sky separation between items in this object's match list,

For each match, project the source centroid to RA,Dec using the initial "linearWcs" WCS, and measure the on-sky angular separation to the reference source coord.

Definition at line 371 of file CreateWcsWithSip.cc.

                                                                   {
    assert(_linearWcs.get());
    return makeMatchStatisticsInRadians(
        *_linearWcs, _matches, afw::math::MEDIAN).getValue()*afw::geom::radians;
}

template<class MatchT>

boost::shared_ptr< afw::image::TanWcs > lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getNewWcs ( )

inline

Definition at line 115 of file CreateWcsWithSip.h.

{ return _newWcs; }

template<class MatchT>

int lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getNGrid ( ) const

inline

Return the number of grid points (on each axis) used in inverse SIP transform.

Definition at line 154 of file CreateWcsWithSip.h.

{ return _ngrid; }

template<class MatchT>

int lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getNPoints ( ) const

inline

Return the number of points in the catalogue.

Definition at line 152 of file CreateWcsWithSip.h.

{ return _matches.size(); }

template<class MatchT>

int lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getOrder ( ) const

inline

Return the number of terms in the SIP matrix.

Definition at line 150 of file CreateWcsWithSip.h.

{ return  _sipA.rows(); }

template<class MatchT >

double lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getScatterInPixels

(

)

const

Compute the median separation, in pixels, between items in this object's match list

For each match, project the reference object coord to pixels using the fit TAN-SIP WCS, and measure the radial separation to the source centroid

Definition at line 352 of file CreateWcsWithSip.cc.

                                                          {
    assert(_newWcs.get());
    return makeMatchStatisticsInPixels(*_newWcs, _matches, afw::math::MEDIAN).getValue();
}

template<class MatchT >

afwGeom::Angle lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::getScatterOnSky

(

)

const

Compute the median on-sky separation between items in this object's match list

For each match, project the source centroid to RA,Dec using the fit TAN-SIP WCS, and measure the on-sky angular separation to the reference source coord.

Definition at line 364 of file CreateWcsWithSip.cc.

                                                             {
    assert(_newWcs.get());
    return makeMatchStatisticsInRadians(
        *_newWcs, _matches, afw::math::MEDIAN).getValue()*afw::geom::radians;
}

Member Data Documentation

template<class MatchT>

afw::geom::Box2I lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_bbox

mutableprivate

Definition at line 161 of file CreateWcsWithSip.h.

template<class MatchT>

boost::shared_ptr< afw::image::Wcs const> lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_linearWcs

private

Definition at line 163 of file CreateWcsWithSip.h.

template<class MatchT>

lsst::pex::logging::Log lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_log

private

Definition at line 158 of file CreateWcsWithSip.h.

template<class MatchT>

std::vector<MatchT> const lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_matches

private

Definition at line 160 of file CreateWcsWithSip.h.

template<class MatchT>

boost::shared_ptr< afw::image::TanWcs > lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_newWcs

private

Definition at line 171 of file CreateWcsWithSip.h.

template<class MatchT>

int lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_ngrid

private

Definition at line 162 of file CreateWcsWithSip.h.

template<class MatchT>

int const lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_reverseSipOrder

private

Definition at line 166 of file CreateWcsWithSip.h.

template<class MatchT>

Eigen::MatrixXd lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_sipA

private

Definition at line 168 of file CreateWcsWithSip.h.

template<class MatchT>

Eigen::MatrixXd lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_sipAp

private

Definition at line 169 of file CreateWcsWithSip.h.

template<class MatchT>

Eigen::MatrixXd lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_sipB

private

Definition at line 168 of file CreateWcsWithSip.h.

template<class MatchT>

Eigen::MatrixXd lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_sipBp

private

Definition at line 169 of file CreateWcsWithSip.h.

template<class MatchT>

int const lsst::meas.astrom.sip::CreateWcsWithSip< MatchT >::_sipOrder

private

Definition at line 166 of file CreateWcsWithSip.h.

---

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

• /home/lsstsw/stack/Linux64/meas_astrom/11.0-1-g60db491/include/lsst/meas/astrom/sip/CreateWcsWithSip.h
• /home/lsstsw/stack/Linux64/meas_astrom/11.0-1-g60db491/src/sip/CreateWcsWithSip.cc