LSSTApplications  19.0.0-11-g2ce9f25+4,20.0.0+1,20.0.0+10,20.0.0+11,20.0.0+13,20.0.0+2,20.0.0+3,20.0.0+4,20.0.0+6,20.0.0-1-g10df615+10,20.0.0-1-g253301a+5,20.0.0-1-g596936a+11,20.0.0-1-g8a53f90+1,20.0.0-1-gc96f8cb+12,20.0.0-1-gd1c87d7+1,20.0.0-17-g41c5faf,20.0.0-2-g04cfba9+4,20.0.0-2-gd11eeda,20.0.0-2-gec03fae+3,20.0.0-3-g082faa5+1,20.0.0-3-gbdbfa727+3,20.0.0-3-gc53c7b6,20.0.0-4-gde602ef96+4,20.0.0-4-ge48a6ca+6,20.0.0-4-ge987224+1,20.0.0-8-g7eef53f7+7,20.0.0-9-g8e1b333,w.2020.28
LSSTDataManagementBasePackage
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. More...
 
std::shared_ptr< afw::geom::SkyWcsgetNewWcs ()
 
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

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 }

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.

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

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

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

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

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

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:
lsst::geom::Box2I::getHeight
int getHeight() const noexcept
Definition: Box.h:188
wcs
table::Key< table::Array< std::uint8_t > > wcs
Definition: SkyWcs.cc:71
std::vector
STL class.
lsst::afw::table._match.matchRaDec
matchRaDec
Definition: _match.py:159
std::vector::size
T size(T... args)
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
crval
table::PointKey< double > crval
Definition: OldWcs.cc:130
lsst::meas::astrom::makeMatchStatisticsInRadians
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.
Definition: makeMatchStatistics.cc:71
lsst::geom::Box2I::isEmpty
bool isEmpty() const noexcept
Return true if the box contains no points.
Definition: Box.h:213
lsst::afw::math::Statistics::getValue
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1056
src
std::shared_ptr< RecordT > src
Definition: Match.cc:48
std::hypot
T hypot(T... args)
lsst::geom::Box2I::getWidth
int getWidth() const noexcept
Definition: Box.h:187
lsst::afw::geom::makeTanSipWcs
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:543
lsst::meas::astrom::sip::CreateWcsWithSip::getNewWcs
std::shared_ptr< afw::geom::SkyWcs > getNewWcs()
Definition: CreateWcsWithSip.h:98
lsst::geom::Box2I::grow
void grow(int buffer)
Increase the size of the box by the given buffer amount in all directions.
Definition: Box.h:249
lsst::geom::Box2I::include
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
Definition: Box.cc:152
crpix
table::PointKey< double > crpix
Definition: OldWcs.cc:131
lsst::meas::astrom::makeMatchStatisticsInPixels
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.
Definition: makeMatchStatistics.cc:48
ptr
uint64_t * ptr
Definition: RangeSet.cc:88
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
LOGL_DEBUG
#define LOGL_DEBUG(logger, message...)
Definition: Log.h:504
std::vector::begin
T begin(T... args)
lsst::geom::Point< int, 2 >
lsst::geom::radians
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
std::vector::empty
T empty(T... args)
std::size_t
lsst::geom::Box2D::getCenter
Point2D const getCenter() const noexcept
Return true if the box contains no points.
Definition: Box.h:549
std::vector::end
T end(T... args)
lsst::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
lsst::geom::Extent< int, 2 >
lsst::meas::astrom::sip::CreateWcsWithSip::CreateWcsWithSip
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.
Definition: CreateWcsWithSip.cc:102
bbox
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
lsst::afw::math::MEDIAN
@ MEDIAN
estimate sample median
Definition: Statistics.h:70
lsst::afw::geom
Definition: frameSetUtils.h:40