LSST Applications g0265f82a02+093ff98f47,g02d81e74bb+10352d6f39,g1f3e9fa97e+40b0fc677d,g2079a07aa2+14824f138e,g2bbee38e9b+093ff98f47,g337abbeb29+093ff98f47,g3ddfee87b4+064c02c7ee,g487adcacf7+7e842ebf4b,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g568d43a26c+1d7ac31cb0,g591dd9f2cf+fb1f69e2ea,g858d7b2824+10352d6f39,g8a8a8dda67+a6fc98d2e7,g8cdfe0ae6a+66d966b544,g99cad8db69+7ce8a7c20a,g9ddcbc5298+d4bad12328,ga1e77700b3+246acaaf9c,ga2e4dd1c03+064c02c7ee,ga8c6da7877+04f6ba86dc,gae46bcf261+093ff98f47,gb0e22166c9+3863383f4c,gba4ed39666+9664299f35,gbb8dafda3b+db40f59a7d,gbeb006f7da+e6a448e96d,gbf5cecdb8a+10352d6f39,gc0f3af6251+10a3fd39cd,gc120e1dc64+5839e50a77,gc28159a63d+093ff98f47,gcf0d15dbbd+064c02c7ee,gd2a12a3803+0c2c227a2b,gdaeeff99f8+a38ce5ea23,ge79ae78c31+093ff98f47,gee10cc3b42+a6fc98d2e7,gf1cff7945b+10352d6f39,w.2024.15
LSST Data Management Base Package
Loading...
Searching...
No Matches
Public Types | Public Member Functions | List of all members
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 std::shared_ptr< CreateWcsWithSipPtr
 
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.
 
std::shared_ptr< afw::geom::SkyWcsgetNewWcs ()
 
double getScatterInPixels () const
 Compute the median separation, in pixels, between items in this object's match list.
 
geom::Angle getScatterOnSky () const
 Compute the median on-sky separation between items in this object's match list.
 
double getLinearScatterInPixels () const
 Compute the median radial separation between items in this object's match list.
 
geom::Angle getLinearScatterOnSky () const
 Compute the median on-sky separation between items in this object's match list,.
 
int getOrder () const
 Return the number of terms in the SIP matrix.
 
int getNPoints () const
 Return the number of points in the catalogue.
 
int getNGrid () const
 Return the number of grid points (on each axis) used in inverse SIP transform.
 
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()
Measure the distortions in an image plane and express them a SIP polynomials.
std::shared_ptr< afw::geom::SkyWcs > 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]matcheslist of matches
[in]linearWcsinitial WCS, typically pure TAN but need not be
[in]orderSIP order for fit WCS
[in]bboxbounding 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]ngridnumber 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.

105 : _matches(matches),
106 _bbox(bbox),
107 _ngrid(ngrid),
108 _linearWcs(std::make_shared<afw::geom::SkyWcs>(linearWcs)),
109 _sipOrder(order + 1),
110 _reverseSipOrder(order + 2), // Higher order for reverse transform
111 _sipA(Eigen::MatrixXd::Zero(_sipOrder, _sipOrder)),
112 _sipB(Eigen::MatrixXd::Zero(_sipOrder, _sipOrder)),
113 _sipAp(Eigen::MatrixXd::Zero(_reverseSipOrder, _reverseSipOrder)),
114 _sipBp(Eigen::MatrixXd::Zero(_reverseSipOrder, _reverseSipOrder)),
115 _newWcs() {
116 if (order < 2) {
117 throw LSST_EXCEPT(pex::exceptions::OutOfRangeError, "SIP must be at least 2nd order");
118 }
119 if (_sipOrder > 9) {
120 throw LSST_EXCEPT(
121 pex::exceptions::OutOfRangeError,
122 str(boost::format("SIP forward order %d exceeds the convention limit of 9") % _sipOrder));
123 }
124 if (_reverseSipOrder > 9) {
125 throw LSST_EXCEPT(pex::exceptions::OutOfRangeError,
126 str(boost::format("SIP reverse order %d exceeds the convention limit of 9") %
127 _reverseSipOrder));
128 }
129
130 if (_matches.size() < std::size_t(_sipOrder)) {
131 throw LSST_EXCEPT(pex::exceptions::LengthError, "Number of matches less than requested sip order");
132 }
133
134 if (_ngrid <= 0) {
135 _ngrid = 5 * _sipOrder; // should be plenty
136 }
137
138 /*
139 * We need a bounding box to define the region over which:
140 * The forward transformation should be valid
141 * We calculate the reverse transformartion
142 * If no BBox is provided, guess one from the input points (extrapolated a bit to allow for fact
143 * that a finite number of points won't reach to the edge of the image)
144 */
145 if (_bbox.isEmpty() && !_matches.empty()) {
146 for (typename std::vector<MatchT>::const_iterator ptr = _matches.begin(); ptr != _matches.end();
147 ++ptr) {
148 afw::table::SourceRecord const& src = *ptr->second;
149 _bbox.include(geom::PointI(src.getX(), src.getY()));
150 }
151 float const borderFrac = 1 / ::sqrt(_matches.size()); // fractional border to add to exact BBox
152 geom::Extent2I border(borderFrac * _bbox.getWidth(), borderFrac * _bbox.getHeight());
153
154 _bbox.grow(border);
155 }
156
157 // If crpix is too far from the center of the fit bbox, move it to the center to improve the fit
158 auto const initialCrpix = _linearWcs->getPixelOrigin();
159 auto const bboxCenter = geom::Box2D(_bbox).getCenter();
160 if (std::hypot(initialCrpix[0] - bboxCenter[0], initialCrpix[1] - bboxCenter[1]) >
161 MAX_DISTANCE_CRPIX_TO_BBOXCTR) {
162 LOGL_DEBUG(_log,
163 "_linearWcs crpix = %d, %d farther than %f from bbox center; shifting to center %d, %d",
164 initialCrpix[0], initialCrpix[1], MAX_DISTANCE_CRPIX_TO_BBOXCTR, bboxCenter[0],
165 bboxCenter[1]);
166 _linearWcs = _linearWcs->getTanWcs(bboxCenter);
167 }
168
169 // Calculate the forward part of the SIP distortion
170 _calculateForwardMatrices();
171
172 // Build a new wcs incorporating the forward SIP matrix, it's all we know so far
173 auto const crval = _linearWcs->getSkyOrigin();
174 auto const crpix = _linearWcs->getPixelOrigin();
175 Eigen::MatrixXd cdMatrix = _linearWcs->getCdMatrix();
176 _newWcs = afw::geom::makeTanSipWcs(crpix, crval, cdMatrix, _sipA, _sipB);
177
178 // Use _newWcs to calculate the forward transformation on a grid, and derive the back transformation
179 _calculateReverseMatrices();
180
181 // Build a new wcs incorporating all SIP matrices
182 _newWcs = afw::geom::makeTanSipWcs(crpix, crval, cdMatrix, _sipA, _sipB, _sipAp, _sipBp);
183}
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition Log.h:515
std::shared_ptr< RecordT > src
Definition Match.cc:48
table::PointKey< double > crpix
Definition OldWcs.cc:129
table::PointKey< double > crval
Definition OldWcs.cc:128
uint64_t * ptr
Definition RangeSet.cc:95
T begin(T... args)
A floating-point coordinate rectangle geometry.
Definition Box.h:413
Point2D const getCenter() const noexcept
Return true if the box contains no points.
Definition Box.h:549
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition Box.cc:152
int getHeight() const noexcept
Definition Box.h:188
bool isEmpty() const noexcept
Return true if the box contains no points.
Definition Box.h:213
void grow(int buffer)
Increase the size of the box by the given buffer amount in all directions.
Definition Box.h:249
int getWidth() const noexcept
Definition Box.h:187
T empty(T... args)
T end(T... args)
T hypot(T... args)
std::shared_ptr< SkyWcs > makeTanSipWcs(lsst::geom::Point2D const &crpix, lsst::geom::SpherePoint const &crval, Eigen::Matrix2d const &cdMatrix, Eigen::MatrixXd const &sipA, Eigen::MatrixXd const &sipB)
Construct a TAN-SIP SkyWcs with forward SIP distortion terms and an iterative inverse.
Definition SkyWcs.cc:538
T size(T... args)
table::Key< int > order

Member Function Documentation

◆ getLinearScatterInPixels()

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.

358 {
359 assert(_linearWcs.get());
360 return makeMatchStatisticsInPixels(*_linearWcs, _matches, afw::math::MEDIAN).getValue();
361}
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
@ MEDIAN
estimate sample median
Definition Statistics.h:60
afw::math::Statistics makeMatchStatisticsInPixels(afw::geom::SkyWcs const &wcs, std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
Compute statistics of on-detector radial separation for a match list, in pixels.

◆ getLinearScatterOnSky()

template<class MatchT >
geom::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 370 of file CreateWcsWithSip.cc.

370 {
371 assert(_linearWcs.get());
372 return makeMatchStatisticsInRadians(*_linearWcs, _matches, afw::math::MEDIAN).getValue() * geom::radians;
373}
AngleUnit constexpr radians
constant with units of radians
Definition Angle.h:109
afw::math::Statistics makeMatchStatisticsInRadians(afw::geom::SkyWcs const &wcs, std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
Compute statistics of on-sky radial separation for a match list, in 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 ( ) 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.

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

◆ getScatterOnSky()

template<class MatchT >
geom::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.

364 {
365 assert(_newWcs.get());
367}

◆ 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: