27 #include "Eigen/Cholesky"
43 namespace except = lsst::pex::exceptions;
44 namespace afwCoord = lsst::afw::coord;
45 namespace afwGeom = lsst::afw::geom;
47 namespace afwDet = lsst::afw::detection;
48 namespace afwMath = lsst::afw::math;
49 namespace afwTable = lsst::afw::table;
50 namespace pexLog = lsst::pex::logging;
58 indexToPQ(
int const index,
int const order)
62 for (
int decrement = order; q >= decrement && decrement > 0; --decrement) {
67 return std::make_pair(p, q);
71 calculateCMatrix(Eigen::VectorXd
const& axis1, Eigen::VectorXd
const& axis2,
int const order)
73 int nTerms = 0.5*order*(order+1);
75 int const n = axis1.size();
76 Eigen::MatrixXd C = Eigen::MatrixXd::Zero(n, nTerms);
77 for (
int i = 0; i < n; ++i) {
78 for (
int j = 0; j < nTerms; ++j) {
79 std::pair<int, int> pq = indexToPQ(j, order);
80 int p = pq.first, q = pq.second;
82 assert(p + q < order);
83 C(i, j) = ::pow(axis1[i], p)*::pow(axis2[i], q);
97 leastSquaresSolve(Eigen::VectorXd
b, Eigen::MatrixXd A) {
98 assert(A.rows() == b.rows());
99 Eigen::VectorXd par = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
106 template<
class MatchT>
108 std::vector<MatchT>
const & matches,
114 _log(pexLog::Log(pexLog::Log::getDefaultLog(),
"meas.astrom.sip")),
118 _linearWcs(linearWcs.clone()),
120 _reverseSipOrder(order+2),
121 _sipA(Eigen::MatrixXd::Zero(_sipOrder, _sipOrder)),
122 _sipB(Eigen::MatrixXd::Zero(_sipOrder, _sipOrder)),
123 _sipAp(Eigen::MatrixXd::Zero(_reverseSipOrder, _reverseSipOrder)),
124 _sipBp(Eigen::MatrixXd::Zero(_reverseSipOrder, _reverseSipOrder)),
128 throw LSST_EXCEPT(except::OutOfRangeError,
"SIP must be at least 2nd order");
132 str(
boost::format(
"SIP forward order %d exceeds the convention limit of 9") %
137 str(
boost::format(
"SIP reverse order %d exceeds the convention limit of 9") %
142 throw LSST_EXCEPT(except::LengthError,
"Number of matches less than requested sip order");
158 typename std::vector<MatchT>::const_iterator ptr =
_matches.begin();
165 float const borderFrac = 1/::sqrt(
_matches.size());
176 Eigen::MatrixXd CD =
_linearWcs->getCDMatrix();
187 template<
class MatchT>
195 int const nPoints = _matches.size();
196 Eigen::VectorXd u(nPoints), v(nPoints), iwc1(nPoints), iwc2(nPoints);
200 typename std::vector<MatchT>::const_iterator ptr = _matches.begin();
201 ptr != _matches.end();
212 u[i] = match.
second->getX() - crpix[0];
213 v[i] = match.
second->getY() - crpix[1];
216 double uMax = u.cwiseAbs().maxCoeff();
217 double vMax = v.cwiseAbs().maxCoeff();
218 double norm = std::max(uMax, vMax);
224 Eigen::MatrixXd forwardC = calculateCMatrix(u, v, ord);
225 Eigen::VectorXd mu = leastSquaresSolve(iwc1, forwardC);
226 Eigen::VectorXd nu = leastSquaresSolve(iwc2, forwardC);
234 assert ((indexToPQ(0, ord) == std::pair<int, int>(0, 0)));
235 assert ((indexToPQ(1, ord) == std::pair<int, int>(0, 1)));
236 assert ((indexToPQ(ord, ord) == std::pair<int, int>(1, 0)));
240 CD(1,0) = nu[ord]/
norm;
241 CD(1,1) = nu[1]/
norm;
242 CD(0,0) = mu[ord]/
norm;
243 CD(0,1) = mu[1]/
norm;
245 Eigen::Matrix2d CDinv = CD.inverse();
248 crpix[0] -= mu[0]*CDinv(0,0) + nu[0]*CDinv(0,1);
249 crpix[1] -= mu[0]*CDinv(1,0) + nu[0]*CDinv(1,1);
265 for(
int i=1; i< mu.rows(); ++i) {
266 std::pair<int, int> pq = indexToPQ(i, ord);
267 int p = pq.first, q = pq.second;
269 if(p + q > 1 && p + q < ord) {
270 Eigen::Vector2d munu(2,1);
273 Eigen::Vector2d AB = CDinv*munu;
275 _sipA(p,q) = AB[0]/::pow(norm,p+q);
276 _sipB(p,q) = AB[1]/::pow(norm,p+q);
281 template<
class MatchT>
283 int const ngrid2 = _ngrid*_ngrid;
285 Eigen::VectorXd U(ngrid2), V(ngrid2);
286 Eigen::VectorXd delta1(ngrid2), delta2(ngrid2);
296 _log.debugf(
"_calcReverseMatrices: x0,y0 %i,%i, W,H %i,%i, ngrid %i, dx,dy %g,%g, CRPIX %g,%g",
300 for (
int i = 0; i < _ngrid; ++i) {
301 double const y = y0 + i*dy;
302 for (
int j = 0; j < _ngrid; ++j, ++k) {
303 double const x = x0 + j*dx;
315 U[k] = xy[0] - 1 - crpix[0];
316 V[k] = xy[1] - 1 - crpix[1];
318 if ((i == 0 || i == (_ngrid-1) || i == (_ngrid/2)) &&
319 (j == 0 || j == (_ngrid-1) || j == (_ngrid/2))) {
320 _log.debugf(
" x,y (%.1f, %.1f), u,v (%.1f, %.1f), U,V (%.1f, %.1f)", x, y, u, v, U[k], V[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) {
343 std::pair<int, int> pq = indexToPQ(j, ord);
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());
370 template<
class MatchT>
372 assert(_linearWcs.get());
377 template<
class MatchT>
384 #define INSTANTIATE(MATCH) \
385 template class CreateWcsWithSip<MATCH>;
lsst::afw::geom::Point2D getPosition(lsst::afw::geom::AngleUnit unit=lsst::afw::geom::degrees) const
Return our contents in a Point2D object.
void _calculateForwardMatrices()
boost::shared_ptr< afw::image::TanWcs > _newWcs
boost::shared_ptr< Record1 > first
void _calculateReverseMatrices()
boost::shared_ptr< Record2 > second
void include(Point2I const &point)
Expand this to ensure that this->contains(point).
double getY() const
Return the centroid slot y coordinate.
afw::geom::Angle getScatterOnSky() const
table::PointKey< double > crval
Implementation of the WCS standard for a any projection.
double getX() const
Return the centroid slot x coordinate.
An integer coordinate rectangle.
table::Key< table::Array< Kernel::Pixel > > image
virtual Fk5Coord toFk5(double const epoch) const
Convert ourself to Fk5 (ie. a no-op): RA, Dec (precess to new epoch)
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
afw::geom::Point2D _getCrvalAsGeomPoint() const
int const _reverseSipOrder
boost::shared_ptr< afw::image::Wcs const > _linearWcs
Implementation of the WCS standard for the special case of the Gnomonic (tangent plane) projection...
CreateWcsWithSip(std::vector< MatchT > const &matches, afw::image::Wcs const &linearWcs, int const order, afw::geom::Box2I const &bbox=afw::geom::Box2I(), int const ngrid=0)
Constructor.
AngleUnit const radians
constant with units of radians
Lightweight representation of a geometric match between two records.
geom::Box2I const & _bbox
afw::math::Statistics makeMatchStatisticsInPixels(afw::image::Wcs const &wcs, std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
bool isEmpty() const
Return true if the box contains no points.
double getScatterInPixels() const
boost::shared_ptr< Wcs > Ptr
A class to handle Fk5 coordinates (inherits from Coord)
#define LSST_EXCEPT(type,...)
void grow(int buffer)
Increase the size of the box by the given buffer amount in all directions.
afw::geom::Angle getLinearScatterOnSky() const
double getLinearScatterInPixels() const
afw::math::Statistics makeMatchStatisticsInRadians(afw::image::Wcs const &wcs, std::vector< MatchT > const &matchList, int const flags, afw::math::StatisticsControl const &sctrl=afw::math::StatisticsControl())
Compute Image Statistics.
afw::table::Key< double > b
Record class that contains measurements made on a single exposure.
#define INSTANTIATE(MATCH)
A class to handle Icrs coordinates (inherits from Coord)
std::vector< MatchT > const _matches
table::PointKey< double > crpix