27 #include "Eigen/Cholesky"
50 double const MAX_DISTANCE_CRPIX_TO_BBOXCTR = 1000;
59 for (
int decrement =
order; q >= decrement && decrement > 0; --decrement) {
67 Eigen::MatrixXd calculateCMatrix(Eigen::VectorXd
const& axis1, Eigen::VectorXd
const& axis2,
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) {
76 int p = pq.first, q = pq.second;
78 assert(p + q <
order);
79 C(i, j) = ::pow(axis1[i], p) * ::pow(axis2[i], q);
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);
101 template <
class MatchT>
109 _sipOrder(
order + 1),
110 _reverseSipOrder(
order + 2),
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)),
122 str(
boost::format(
"SIP forward order %d exceeds the convention limit of 9") % _sipOrder));
124 if (_reverseSipOrder > 9) {
126 str(
boost::format(
"SIP reverse order %d exceeds the convention limit of 9") %
135 _ngrid = 5 * _sipOrder;
151 float const borderFrac = 1 / ::sqrt(_matches.
size());
158 auto const initialCrpix = _linearWcs->getPixelOrigin();
160 if (
std::hypot(initialCrpix[0] - bboxCenter[0], initialCrpix[1] - bboxCenter[1]) >
161 MAX_DISTANCE_CRPIX_TO_BBOXCTR) {
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],
166 _linearWcs = _linearWcs->getTanWcs(bboxCenter);
170 _calculateForwardMatrices();
173 auto const crval = _linearWcs->getSkyOrigin();
174 auto const crpix = _linearWcs->getPixelOrigin();
175 Eigen::MatrixXd cdMatrix = _linearWcs->getCdMatrix();
179 _calculateReverseMatrices();
185 template <
class MatchT>
191 int const nPoints = _matches.size();
192 Eigen::VectorXd u(nPoints), v(nPoints), iwc1(nPoints), iwc2(nPoints);
201 auto c = match.
first->getCoord();
210 double uMax = u.cwiseAbs().maxCoeff();
211 double vMax = v.cwiseAbs().maxCoeff();
218 Eigen::MatrixXd forwardC = calculateCMatrix(u, v, ord);
219 Eigen::VectorXd mu = leastSquaresSolve(iwc1, forwardC);
220 Eigen::VectorXd nu = leastSquaresSolve(iwc2, forwardC);
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;
239 Eigen::Matrix2d CDinv = CD.inverse();
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);
245 auto const crval = _linearWcs->getSkyOrigin();
259 for (
int i = 1; i < mu.rows(); ++i) {
261 int p = pq.first, q = pq.second;
263 if (p + q > 1 && p + q < ord) {
264 Eigen::Vector2d munu(2, 1);
267 Eigen::Vector2d AB = CDinv * munu;
269 _sipA(p, q) = AB[0] / ::pow(
norm, p + q);
270 _sipB(p, q) = AB[1] / ::pow(
norm, p + q);
275 template <
class MatchT>
276 void CreateWcsWithSip<MatchT>::_calculateReverseMatrices() {
277 int const ngrid2 = _ngrid * _ngrid;
279 Eigen::VectorXd U(ngrid2), V(ngrid2);
280 Eigen::VectorXd delta1(ngrid2), delta2(ngrid2);
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);
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]);
293 auto tanWcs = _newWcs->getTanWcs(_newWcs->getPixelOrigin());
296 for (
int i = 0; i < _ngrid; ++i) {
297 double const y = y0 + i * dy;
299 beforeSipABPoints.
reserve(_ngrid);
300 for (
int j = 0; j < _ngrid; ++j) {
301 double const x = x0 + j * dx;
304 auto const afterSipABPoints = applySipAB->applyForward(beforeSipABPoints);
305 for (
int j = 0; j < _ngrid; ++j, ++k) {
306 double const x = x0 + j * dx;
314 U[k] = xy[0] -
crpix[0];
315 V[k] = xy[1] -
crpix[1];
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],
323 delta1[k] = u - U[k];
324 delta2[k] = v - V[k];
329 double UMax = U.cwiseAbs().maxCoeff();
330 double VMax = V.cwiseAbs().maxCoeff();
331 double norm = (UMax > VMax) ? UMax : VMax;
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);
341 assert(tmpA.rows() == tmpB.rows());
342 for (
int j = 0; j < tmpA.rows(); ++j) {
344 int p = pq.first, q = pq.second;
346 _sipAp(p, q) = tmpA[j] / ::pow(
norm, p + q);
347 _sipBp(p, q) = tmpB[j] / ::pow(
norm, p + q);
351 template <
class MatchT>
353 assert(_newWcs.get());
357 template <
class MatchT>
359 assert(_linearWcs.get());
363 template <
class MatchT>
365 assert(_newWcs.get());
369 template <
class MatchT>
371 assert(_linearWcs.get());
375 #define INSTANTIATE(MATCH) template class CreateWcsWithSip<MATCH>;
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
LSST DM logging module built on log4cxx.
#define LOG_GET(logger)
Returns a Log object associated with logger.
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
std::shared_ptr< RecordT > src
table::PointKey< double > crpix
table::PointKey< double > crval
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Record class that contains measurements made on a single exposure.
A class representing an angle.
A floating-point coordinate rectangle geometry.
Point2D const getCenter() const noexcept
Return true if the box contains no points.
An integer coordinate rectangle.
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
int getHeight() const noexcept
bool isEmpty() const noexcept
Return true if the box contains no points.
void grow(int buffer)
Increase the size of the box by the given buffer amount in all directions.
int getWidth() const noexcept
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.
Reports attempts to access elements outside a valid range of indices.
T emplace_back(T... args)
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
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.
std::shared_ptr< TransformPoint2ToSpherePoint > getIntermediateWorldCoordsToSky(SkyWcs const &wcs, bool simplify=true)
Return a transform from intermediate world coordinates to sky.
std::shared_ptr< TransformPoint2ToPoint2 > makeWcsPairTransform(SkyWcs const &src, SkyWcs const &dst)
A Transform obtained by putting two SkyWcs objects "back to back".
@ MEDIAN
estimate sample median
Match< SimpleRecord, SourceRecord > ReferenceMatch
Match< SourceRecord, SourceRecord > SourceMatch
constexpr AngleUnit radians
constant with units of radians
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)
A base class for image defects.
#define INSTANTIATE(MATCH)
Lightweight representation of a geometric match between two records.
std::shared_ptr< Record2 > second
std::shared_ptr< Record1 > first