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

#include <CreateWcsWithSip.h>

Public Types
typedef std::shared_ptr< CreateWcsWithSip >	Ptr

typedef std::shared_ptr< CreateWcsWithSip const >	ConstPtr

Public Member Functions
	CreateWcsWithSip (std::vector< MatchT > const &matches, afw::geom::SkyWcs const &linearWcs, int const order, geom::Box2I const &bbox=geom::Box2I(), int const ngrid=0)
	Construct a CreateWcsWithSip. More...

std::shared_ptr< afw::geom::SkyWcs >	getNewWcs ()

double	getScatterInPixels () const
	Compute the median separation, in pixels, between items in this object's match list. More...

geom::Angle	getScatterOnSky () const
	Compute the median on-sky separation between items in this object's match list. More...

double	getLinearScatterInPixels () const
	Compute the median radial separation between items in this object's match list. More...

geom::Angle	getLinearScatterOnSky () const
	Compute the median on-sky separation between items in this object's match list,. More...

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...

Eigen::MatrixXd const	getSipA ()

Eigen::MatrixXd const	getSipB ()

Eigen::MatrixXd const	getSipAp ()

Eigen::MatrixXd const	getSipBp ()

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 76 of file CreateWcsWithSip.h.

Member Typedef Documentation

◆ ConstPtr

template<class MatchT >

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

Definition at line 79 of file CreateWcsWithSip.h.

◆ Ptr

template<class MatchT >

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

Definition at line 78 of file CreateWcsWithSip.h.

Constructor & Destructor Documentation

◆ CreateWcsWithSip()

template<class MatchT >

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

Construct a CreateWcsWithSip.

Constructor.

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 102 of file CreateWcsWithSip.cc.

         : _matches(matches),
           _bbox(bbox),
           _ngrid(ngrid),
           _linearWcs(std::make_shared<afw::geom::SkyWcs>(linearWcs)),
           _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(pex::exceptions::OutOfRangeError, "SIP must be at least 2nd order");
     }
     if (_sipOrder > 9) {
         throw LSST_EXCEPT(
                 pex::exceptions::OutOfRangeError,
                 str(boost::format("SIP forward order %d exceeds the convention limit of 9") % _sipOrder));
     }
     if (_reverseSipOrder > 9) {
         throw LSST_EXCEPT(pex::exceptions::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(pex::exceptions::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()) {
         for (typename std::vector<MatchT>::const_iterator ptr = _matches.begin(); ptr != _matches.end();
              ++ptr) {
             afw::table::SourceRecord const& src = *ptr->second;
             _bbox.include(geom::PointI(src.getX(), src.getY()));
         }
         float const borderFrac = 1 / ::sqrt(_matches.size());  // fractional border to add to exact BBox
         geom::Extent2I border(borderFrac * _bbox.getWidth(), borderFrac * _bbox.getHeight());
  
         _bbox.grow(border);
     }
  
     // If crpix is too far from the center of the fit bbox, move it to the center to improve the fit
     auto const initialCrpix = _linearWcs->getPixelOrigin();
     auto const bboxCenter = geom::Box2D(_bbox).getCenter();
     if (std::hypot(initialCrpix[0] - bboxCenter[0], initialCrpix[1] - bboxCenter[1]) >
         MAX_DISTANCE_CRPIX_TO_BBOXCTR) {
         LOGL_DEBUG(_log,
                    "_linearWcs crpix = %d, %d farther than %f from bbox center; shifting to center %d, %d",
                    initialCrpix[0], initialCrpix[1], MAX_DISTANCE_CRPIX_TO_BBOXCTR, bboxCenter[0],
                    bboxCenter[1]);
         _linearWcs = _linearWcs->getTanWcs(bboxCenter);
     }
  
     // 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
     auto const crval = _linearWcs->getSkyOrigin();
     auto const crpix = _linearWcs->getPixelOrigin();
     Eigen::MatrixXd cdMatrix = _linearWcs->getCdMatrix();
     _newWcs = afw::geom::makeTanSipWcs(crpix, crval, cdMatrix, _sipA, _sipB);
  
     // Use _newWcs to calculate the forward transformation on a grid, and derive the back transformation
     _calculateReverseMatrices();
  
     // Build a new wcs incorporating all SIP matrices
     _newWcs = afw::geom::makeTanSipWcs(crpix, crval, cdMatrix, _sipA, _sipB, _sipAp, _sipBp);
 }

Member Function Documentation

◆ getLinearScatterInPixels()

template<class MatchT >

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

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

◆ getLinearScatterOnSky()

template<class MatchT >

geom::Angle lsst::meas::astrom::sip::CreateWcsWithSip< MatchT >::getLinearScatterOnSky

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 370 of file CreateWcsWithSip.cc.

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

◆ getNewWcs()

template<class MatchT >

std::shared_ptr<afw::geom::SkyWcs> lsst::meas::astrom::sip::CreateWcsWithSip< MatchT >::getNewWcs ( )

inline

Definition at line 98 of file CreateWcsWithSip.h.

98 { return _newWcs; }

◆ getNGrid()

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 137 of file CreateWcsWithSip.h.

137 { return _ngrid; }

◆ getNPoints()

template<class MatchT >

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

inline

Return the number of points in the catalogue.

Definition at line 135 of file CreateWcsWithSip.h.

135 { return _matches.size(); }

◆ getOrder()

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 133 of file CreateWcsWithSip.h.

133 { return _sipA.rows(); }

◆ getScatterInPixels()

template<class MatchT >

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

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

◆ getScatterOnSky()

template<class MatchT >

geom::Angle lsst::meas::astrom::sip::CreateWcsWithSip< MatchT >::getScatterOnSky

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() * geom::radians;
 }

◆ getSipA()

template<class MatchT >

Eigen::MatrixXd const lsst::meas::astrom::sip::CreateWcsWithSip< MatchT >::getSipA ( )

inline

Definition at line 140 of file CreateWcsWithSip.h.

140 { return _sipA; }

◆ getSipAp()

template<class MatchT >

Eigen::MatrixXd const lsst::meas::astrom::sip::CreateWcsWithSip< MatchT >::getSipAp ( )

inline

Definition at line 144 of file CreateWcsWithSip.h.

144 { return _sipAp; }

◆ getSipB()

template<class MatchT >

Eigen::MatrixXd const lsst::meas::astrom::sip::CreateWcsWithSip< MatchT >::getSipB ( )

inline

Definition at line 142 of file CreateWcsWithSip.h.

142 { return _sipB; }

◆ getSipBp()

template<class MatchT >

Eigen::MatrixXd const lsst::meas::astrom::sip::CreateWcsWithSip< MatchT >::getSipBp ( )

inline

Definition at line 146 of file CreateWcsWithSip.h.

146 { return _sipBp; }

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

/j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/meas_astrom/22.0.1-12-g44dc1dc+2a20fdde0d/include/lsst/meas/astrom/sip/CreateWcsWithSip.h
/j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/meas_astrom/22.0.1-12-g44dc1dc+2a20fdde0d/src/sip/CreateWcsWithSip.cc

Public Types

Public Member Functions

Detailed Description

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

Member Typedef Documentation

◆ ConstPtr

◆ Ptr

Constructor & Destructor Documentation

◆ CreateWcsWithSip()

Member Function Documentation

◆ getLinearScatterInPixels()

◆ getLinearScatterOnSky()

◆ getNewWcs()

◆ getNGrid()

◆ getNPoints()

◆ getOrder()

◆ getScatterInPixels()

◆ getScatterOnSky()

◆ getSipA()

◆ getSipAp()

◆ getSipB()

◆ getSipBp()

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