LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
CreateWcsWithSip.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 #include <cmath>
25 
26 #include "Eigen/SVD"
27 #include "Eigen/Cholesky"
28 #include "Eigen/LU"
29 
32 #include "lsst/geom/Angle.h"
33 #include "lsst/geom/Box.h"
35 #include "lsst/afw/geom/SkyWcs.h"
36 #include "lsst/log/Log.h"
38 
39 namespace {
40 LOG_LOGGER _log = LOG_GET("lsst.meas.astrom.sip");
41 }
42 
43 namespace lsst {
44 namespace meas {
45 namespace astrom {
46 namespace sip {
47 
48 namespace {
49 
50 double const MAX_DISTANCE_CRPIX_TO_BBOXCTR = 1000;
51 
52 /*
53  * Given an index and a SIP order, calculate p and q for the index'th term u^p v^q
54  * (Cf. Eqn 2 in http://fits.gsfc.nasa.gov/registry/sip/SIP_distortion_v1_0.pdf)
55  */
56 std::pair<int, int> indexToPQ(int const index, int const order) {
57  int p = 0, q = index;
58 
59  for (int decrement = order; q >= decrement && decrement > 0; --decrement) {
60  q -= decrement;
61  p++;
62  }
63 
64  return std::make_pair(p, q);
65 }
66 
67 Eigen::MatrixXd calculateCMatrix(Eigen::VectorXd const& axis1, Eigen::VectorXd const& axis2,
68  int const order) {
69  int nTerms = 0.5 * order * (order + 1);
70 
71  int const n = axis1.size();
72  Eigen::MatrixXd C = Eigen::MatrixXd::Zero(n, nTerms);
73  for (int i = 0; i < n; ++i) {
74  for (int j = 0; j < nTerms; ++j) {
75  std::pair<int, int> pq = indexToPQ(j, order);
76  int p = pq.first, q = pq.second;
77 
78  assert(p + q < order);
79  C(i, j) = ::pow(axis1[i], p) * ::pow(axis2[i], q);
80  }
81  }
82 
83  return C;
84 }
85 
92 Eigen::VectorXd leastSquaresSolve(Eigen::VectorXd b, Eigen::MatrixXd A) {
93  assert(A.rows() == b.rows());
94  Eigen::VectorXd par = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
95  return par;
96 }
97 
98 } // anonymous namespace
99 
101 template <class MatchT>
103  afw::geom::SkyWcs const& linearWcs, int const order,
104  geom::Box2I const& bbox, int const ngrid)
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(
122  str(boost::format("SIP forward order %d exceeds the convention limit of 9") % _sipOrder));
123  }
124  if (_reverseSipOrder > 9) {
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 }
184 
185 template <class MatchT>
187  // Assumes FITS (1-indexed) coordinates.
188  geom::Point2D crpix = _linearWcs->getPixelOrigin();
189 
190  // Calculate u, v and intermediate world coordinates
191  int const nPoints = _matches.size();
192  Eigen::VectorXd u(nPoints), v(nPoints), iwc1(nPoints), iwc2(nPoints);
193 
194  int i = 0;
195  auto linearIwcToSky = getIntermediateWorldCoordsToSky(*_linearWcs);
196  for (typename std::vector<MatchT>::const_iterator ptr = _matches.begin(); ptr != _matches.end();
197  ++ptr, ++i) {
198  afw::table::ReferenceMatch const& match = *ptr;
199 
200  // iwc: intermediate world coordinate positions of catalogue objects
201  auto c = match.first->getCoord();
202  geom::Point2D p = linearIwcToSky->applyInverse(c);
203  iwc1[i] = p[0];
204  iwc2[i] = p[1];
205  // u and v are intermediate pixel coordinates of observed (distorted) positions
206  u[i] = match.second->getX() - crpix[0];
207  v[i] = match.second->getY() - crpix[1];
208  }
209  // Scale u and v down to [-1,,+1] in order to avoid too large numbers in the polynomials
210  double uMax = u.cwiseAbs().maxCoeff();
211  double vMax = v.cwiseAbs().maxCoeff();
212  double norm = std::max(uMax, vMax);
213  u = u / norm;
214  v = v / norm;
215 
216  // Forward transform
217  int ord = _sipOrder;
218  Eigen::MatrixXd forwardC = calculateCMatrix(u, v, ord);
219  Eigen::VectorXd mu = leastSquaresSolve(iwc1, forwardC);
220  Eigen::VectorXd nu = leastSquaresSolve(iwc2, forwardC);
221 
222  // Use mu and nu to refine CD
223 
224  // Given the implementation of indexToPQ(), the refined values
225  // of the elements of the CD matrices are in elements 1 and "_sipOrder" of mu and nu
226  // If the implementation of indexToPQ() changes, these assertions
227  // will catch that change.
228  assert((indexToPQ(0, ord) == std::pair<int, int>(0, 0)));
229  assert((indexToPQ(1, ord) == std::pair<int, int>(0, 1)));
230  assert((indexToPQ(ord, ord) == std::pair<int, int>(1, 0)));
231 
232  // Scale back CD matrix
233  Eigen::Matrix2d CD;
234  CD(1, 0) = nu[ord] / norm;
235  CD(1, 1) = nu[1] / norm;
236  CD(0, 0) = mu[ord] / norm;
237  CD(0, 1) = mu[1] / norm;
238 
239  Eigen::Matrix2d CDinv = CD.inverse(); // Direct inverse OK for 2x2 matrix in Eigen
240 
241  // The zeroth elements correspond to a shift in crpix
242  crpix[0] -= mu[0] * CDinv(0, 0) + nu[0] * CDinv(0, 1);
243  crpix[1] -= mu[0] * CDinv(1, 0) + nu[0] * CDinv(1, 1);
244 
245  auto const crval = _linearWcs->getSkyOrigin();
246 
247  _linearWcs = afw::geom::makeSkyWcs(crpix, crval, CD);
248 
249  // Get Sip terms
250 
251  // The rest of the elements correspond to
252  // mu[i] == CD11*Apq + CD12*Bpq and
253  // nu[i] == CD21*Apq + CD22*Bpq and
254  //
255  // We solve for Apq and Bpq with the equation
256  // (Apq) = (CD11 CD12)-1 * (mu[i])
257  // (Bpq) (CD21 CD22) (nu[i])
258 
259  for (int i = 1; i < mu.rows(); ++i) {
260  std::pair<int, int> pq = indexToPQ(i, ord);
261  int p = pq.first, q = pq.second;
262 
263  if (p + q > 1 && p + q < ord) {
264  Eigen::Vector2d munu(2, 1);
265  munu(0) = mu(i);
266  munu(1) = nu(i);
267  Eigen::Vector2d AB = CDinv * munu;
268  // Scale back sip coefficients
269  _sipA(p, q) = AB[0] / ::pow(norm, p + q);
270  _sipB(p, q) = AB[1] / ::pow(norm, p + q);
271  }
272  }
273 }
274 
275 template <class MatchT>
276 void CreateWcsWithSip<MatchT>::_calculateReverseMatrices() {
277  int const ngrid2 = _ngrid * _ngrid;
278 
279  Eigen::VectorXd U(ngrid2), V(ngrid2);
280  Eigen::VectorXd delta1(ngrid2), delta2(ngrid2);
281 
282  int const x0 = _bbox.getMinX();
283  double const dx = _bbox.getWidth() / (double)(_ngrid - 1);
284  int const y0 = _bbox.getMinY();
285  double const dy = _bbox.getHeight() / (double)(_ngrid - 1);
286 
287  // wcs->getPixelOrigin() returns LSST-style (0-indexed) pixel coords.
288  geom::Point2D crpix = _newWcs->getPixelOrigin();
289 
290  LOGL_DEBUG(_log, "_calcReverseMatrices: x0,y0 %i,%i, W,H %i,%i, ngrid %i, dx,dy %g,%g, CRPIX %g,%g", x0,
291  y0, _bbox.getWidth(), _bbox.getHeight(), _ngrid, dx, dy, crpix[0], crpix[1]);
292 
293  auto tanWcs = _newWcs->getTanWcs(_newWcs->getPixelOrigin());
294  auto applySipAB = afw::geom::makeWcsPairTransform(*_newWcs, *tanWcs);
295  int k = 0;
296  for (int i = 0; i < _ngrid; ++i) {
297  double const y = y0 + i * dy;
298  std::vector<geom::Point2D> beforeSipABPoints;
299  beforeSipABPoints.reserve(_ngrid);
300  for (int j = 0; j < _ngrid; ++j) {
301  double const x = x0 + j * dx;
302  beforeSipABPoints.emplace_back(geom::Point2D(x, y));
303  }
304  auto const afterSipABPoints = applySipAB->applyForward(beforeSipABPoints);
305  for (int j = 0; j < _ngrid; ++j, ++k) {
306  double const x = x0 + j * dx;
307  double u, v;
308  // u and v are intermediate pixel coordinates on a grid of positions
309  u = x - crpix[0];
310  v = y - crpix[1];
311 
312  // U and V are the result of applying the "forward" (A,B) SIP coefficients
313  geom::Point2D const xy = afterSipABPoints[j];
314  U[k] = xy[0] - crpix[0];
315  V[k] = xy[1] - crpix[1];
316 
317  if ((i == 0 || i == (_ngrid - 1) || i == (_ngrid / 2)) &&
318  (j == 0 || j == (_ngrid - 1) || j == (_ngrid / 2))) {
319  LOGL_DEBUG(_log, " x,y (%.1f, %.1f), u,v (%.1f, %.1f), U,V (%.1f, %.1f)", x, y, u, v, U[k],
320  V[k]);
321  }
322 
323  delta1[k] = u - U[k];
324  delta2[k] = v - V[k];
325  }
326  }
327 
328  // Scale down U and V in order to avoid too large numbers in the polynomials
329  double UMax = U.cwiseAbs().maxCoeff();
330  double VMax = V.cwiseAbs().maxCoeff();
331  double norm = (UMax > VMax) ? UMax : VMax;
332  U = U / norm;
333  V = V / norm;
334 
335  // Reverse transform
336  int const ord = _reverseSipOrder;
337  Eigen::MatrixXd reverseC = calculateCMatrix(U, V, ord);
338  Eigen::VectorXd tmpA = leastSquaresSolve(delta1, reverseC);
339  Eigen::VectorXd tmpB = leastSquaresSolve(delta2, reverseC);
340 
341  assert(tmpA.rows() == tmpB.rows());
342  for (int j = 0; j < tmpA.rows(); ++j) {
343  std::pair<int, int> pq = indexToPQ(j, ord);
344  int p = pq.first, q = pq.second;
345  // Scale back sip coefficients
346  _sipAp(p, q) = tmpA[j] / ::pow(norm, p + q);
347  _sipBp(p, q) = tmpB[j] / ::pow(norm, p + q);
348  }
349 }
350 
351 template <class MatchT>
353  assert(_newWcs.get());
354  return makeMatchStatisticsInPixels(*_newWcs, _matches, afw::math::MEDIAN).getValue();
355 }
356 
357 template <class MatchT>
359  assert(_linearWcs.get());
360  return makeMatchStatisticsInPixels(*_linearWcs, _matches, afw::math::MEDIAN).getValue();
361 }
362 
363 template <class MatchT>
365  assert(_newWcs.get());
367 }
368 
369 template <class MatchT>
371  assert(_linearWcs.get());
372  return makeMatchStatisticsInRadians(*_linearWcs, _matches, afw::math::MEDIAN).getValue() * geom::radians;
373 }
374 
375 #define INSTANTIATE(MATCH) template class CreateWcsWithSip<MATCH>;
376 
379 
380 } // namespace sip
381 } // namespace astrom
382 } // namespace meas
383 } // namespace lsst
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
LSST DM logging module built on log4cxx.
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
#define LOG_LOGGER
Definition: Log.h:714
#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:88
int y
Definition: SpanSet.cc:48
table::Key< int > b
T begin(T... args)
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1047
Record class that contains measurements made on a single exposure.
Definition: Source.h:78
A class representing an angle.
Definition: Angle.h:127
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
An integer coordinate rectangle.
Definition: Box.h:55
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
Measure the distortions in an image plane and express them a SIP polynomials.
double getLinearScatterInPixels() const
Compute the median radial separation between items in this object's match list.
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.
geom::Angle getLinearScatterOnSky() const
Compute the median on-sky separation between items in this object's match list,.
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.
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
Reports attempts to access elements outside a valid range of indices.
Definition: Runtime.h:89
T emplace_back(T... args)
T empty(T... args)
T end(T... args)
T hypot(T... args)
T make_pair(T... args)
T max(T... args)
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
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
std::shared_ptr< TransformPoint2ToSpherePoint > getIntermediateWorldCoordsToSky(SkyWcs const &wcs, bool simplify=true)
Return a transform from intermediate world coordinates to sky.
Definition: SkyWcs.cc:553
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
Definition: SkyWcs.cc:146
T norm(const T &x)
Definition: Integrate.h:160
@ MEDIAN
estimate sample median
Definition: Statistics.h:69
Match< SimpleRecord, SourceRecord > ReferenceMatch
Definition: fwd.h:104
Match< SourceRecord, SourceRecord > SourceMatch
Definition: fwd.h:105
constexpr AngleUnit radians
constant with units of radians
Definition: Angle.h:108
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.
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.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
STL namespace.
T reserve(T... args)
T size(T... args)
#define INSTANTIATE(MATCH)
Lightweight representation of a geometric match between two records.
Definition: Match.h:67
std::shared_ptr< Record2 > second
Definition: Match.h:69
std::shared_ptr< Record1 > first
Definition: Match.h:68
table::Key< int > order