Skip to content

Commit

Permalink
Merge pull request #65 from lsst/tickets/DM-9862
Browse files Browse the repository at this point in the history 
DM-9862: functions to manipulate TAN-SIP WCSs
  • Loading branch information
@TallJimbo
TallJimbo committed Apr 20, 2017
2 parents 2013b0d + 7aa0963 commit 0db1a30
Show file tree
Hide file tree
Showing 4 changed files with 303 additions and 1 deletion.
63 changes: 63 additions & 0 deletions include/lsst/meas/astrom/SipTransform.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,9 @@
#define LSST_MEAS_ASTROM_SipTransform_INCLUDED

#include "lsst/afw/geom/LinearTransform.h"
#include "lsst/afw/geom/AffineTransform.h"
#include "lsst/afw/geom/Box.h"
#include "lsst/afw/geom/Angle.h"
#include "lsst/meas/astrom/PolynomialTransform.h"


Expand Down Expand Up @@ -103,6 +106,8 @@ class SipTransformBase {
_poly.swap(other._poly);
}

void transformPixelsInPlace(afw::geom::AffineTransform const & s);

afw::geom::Point2D _pixelOrigin;
afw::geom::LinearTransform _cdMatrix;
PolynomialTransform _poly;
Expand Down Expand Up @@ -185,6 +190,11 @@ class SipForwardTransform : public SipTransformBase {
*/
static SipForwardTransform convert(ScaledPolynomialTransform const & scaled);

/**
* Extract the reverse transform from a TanWcs.
*/
static SipForwardTransform extract(afw::image::TanWcs const & wcs);

/**
* Construct a SipForwardTransform from its components.
*
Expand Down Expand Up @@ -222,6 +232,12 @@ class SipForwardTransform : public SipTransformBase {
*/
afw::geom::Point2D operator()(afw::geom::Point2D const & uv) const;

/**
* Return a new forward SIP transform that includes a transformation of
* the pixel coordinate system by the given affine transform.
*/
SipForwardTransform transformPixels(afw::geom::AffineTransform const & s) const;

};


Expand Down Expand Up @@ -303,6 +319,11 @@ class SipReverseTransform : public SipTransformBase {
*/
static SipReverseTransform convert(ScaledPolynomialTransform const & scaled);

/**
* Extract the reverse transform from a TanWcs.
*/
static SipReverseTransform extract(afw::image::TanWcs const & wcs);

/**
* Construct a SipReverseTransform from its components.
*
Expand Down Expand Up @@ -341,6 +362,12 @@ class SipReverseTransform : public SipTransformBase {
*/
afw::geom::Point2D operator()(afw::geom::Point2D const & xy) const;

/**
* Return a new reverse SIP transform that includes a transformation of
* the pixel coordinate system by the given affine transform.
*/
SipReverseTransform transformPixels(afw::geom::AffineTransform const & s) const;

private:
friend class PolynomialTransform;
friend class ScaledPolynomialTransform;
Expand All @@ -367,6 +394,42 @@ std::shared_ptr<afw::image::TanWcs> makeWcs(
afw::coord::Coord const & skyOrigin
);

/**
* Create a new TanWcs whose pixel coordinate system has been transformed
* via an affine transform.
*
* @param[in] wcs Original TanWcs object.
* @param[in] s AffineTransform to apply to the pixel coordinate
* system.
*
* @return a new Wcs that satisfies the following:
* @code
* newWcs = transformWcsPixels(wcs, s);
* assert(newWcs.skyToPixel(sky), s(wcs.skyToPixel(sky)));
* assert(newWcs.pixelToSky(pixel), wcs.pixelToSky(s.invert()(pixel)));
* @endcode
* for all sky coordinates @c sky and pixel coordinates @c pixel.
*/
std::shared_ptr<afw::image::TanWcs> transformWcsPixels(
afw::image::TanWcs const & wcs,
afw::geom::AffineTransform const & s
);

/**
* Return a new TanWcs that represents a rotation of the image
* it corresponds to about the image's center.
*
* @param[in] wcs Original TanWcs to be rotated.
* @param[in] nQuarter Number of 90 degree rotations (positive is counterclockwise).
* @param[in] dimensions Width and height of the image.
*/
std::shared_ptr<afw::image::TanWcs> rotateWcsPixelsBy90(
afw::image::TanWcs const & wcs,
int nQuarter,
afw::geom::Extent2I const & dimensions
);


}}} // namespace lsst::meas::astrom

#endif // !LSST_MEAS_ASTROM_SipTransform_INCLUDED
6 changes: 6 additions & 0 deletions python/lsst/meas/astrom/sipTransform.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,8 +64,10 @@ static void declareSipForwardTransform(py::module &mod) {
cls.def_static("convert",
(SipForwardTransform(*)(ScaledPolynomialTransform const &)) & SipForwardTransform::convert,
"scaled"_a);
cls.def_static("extract", &SipForwardTransform::extract, "wcs"_a);

cls.def("__call__", &SipForwardTransform::operator(), "in"_a);
cls.def("transformPixels", &SipForwardTransform::transformPixels, "s"_a);

cls.def("linearize", &SipForwardTransform::linearize);
}
Expand All @@ -91,8 +93,10 @@ static void declareSipReverseTransform(py::module &mod) {
cls.def_static("convert",
(SipReverseTransform(*)(ScaledPolynomialTransform const &)) & SipReverseTransform::convert,
"scaled"_a);
cls.def_static("extract", &SipReverseTransform::extract, "wcs"_a);

cls.def("__call__", &SipReverseTransform::operator(), "in"_a);
cls.def("transformPixels", &SipReverseTransform::transformPixels, "s"_a);

cls.def("linearize", &SipReverseTransform::linearize);
}
Expand All @@ -107,6 +111,8 @@ PYBIND11_PLUGIN(sipTransform) {
declareSipReverseTransform(mod);

mod.def("makeWcs", &makeWcs);
mod.def("transformWcsPixels", &transformWcsPixels, "wcs"_a, "s"_a);
mod.def("rotateWcsPixelsBy90", &rotateWcsPixelsBy90, "wcs"_a, "nQuarter"_a, "dimensions"_a);

return mod.ptr();
}
Expand Down
110 changes: 110 additions & 0 deletions src/SipTransform.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,18 @@

namespace lsst { namespace meas { namespace astrom {

void SipTransformBase::transformPixelsInPlace(afw::geom::AffineTransform const & s) {
// The implementation for transformPixels is identical for
// SipForwardTransform and SipReverseTransform. That's pretty obvious for
// the pixel origin and CD matrix, which are the same in both cases, but
// it wasn't obvious to me until I did the math that the polynomial
// transforms are composed with the affine transform the same way.
auto sInv = s.invert();
_pixelOrigin = s.getLinear()(_pixelOrigin - sInv.getTranslation());
_cdMatrix = _cdMatrix * sInv.getLinear();
_poly = compose(s.getLinear(), compose(getPoly(), sInv.getLinear()));
}

SipForwardTransform SipForwardTransform::convert(
PolynomialTransform const & poly,
afw::geom::Point2D const & pixelOrigin,
Expand Down Expand Up @@ -75,6 +87,24 @@ SipForwardTransform SipForwardTransform::convert(ScaledPolynomialTransform const
return convert(scaled, pixelOrigin, cdMatrix);
}

SipForwardTransform SipForwardTransform::extract(afw::image::TanWcs const & wcs) {
if (!wcs.hasDistortion()) {
throw LSST_EXCEPT(
pex::exceptions::LogicError,
"Constructing a SipForwardTransform from a TanWcs with no distortions is not implemented."
);
}
ndarray::Array<double,2,2> xCoeffs = ndarray::allocate(wcs.getSipA().rows(), wcs.getSipA().cols());
ndarray::Array<double,2,2> yCoeffs = ndarray::allocate(wcs.getSipB().rows(), wcs.getSipB().cols());
xCoeffs.asEigen() = wcs.getSipA();
yCoeffs.asEigen() = wcs.getSipB();
return SipForwardTransform(
wcs.getPixelOrigin(),
afw::geom::LinearTransform(wcs.getCDMatrix()),
PolynomialTransform(xCoeffs, yCoeffs)
);
}

afw::geom::AffineTransform SipForwardTransform::linearize(afw::geom::Point2D const & in) const {
afw::geom::AffineTransform tail(-afw::geom::Extent2D(getPixelOrigin()));
return afw::geom::AffineTransform(_cdMatrix)
Expand All @@ -87,6 +117,12 @@ afw::geom::Point2D SipForwardTransform::operator()(afw::geom::Point2D const & uv
return getCDMatrix()(afw::geom::Extent2D(duv) + getPoly()(duv));
}

SipForwardTransform SipForwardTransform::transformPixels(afw::geom::AffineTransform const & s) const {
SipForwardTransform result(*this);
result.transformPixelsInPlace(s);
return result;
}

SipReverseTransform SipReverseTransform::convert(
PolynomialTransform const & poly,
afw::geom::Point2D const & pixelOrigin,
Expand Down Expand Up @@ -134,6 +170,31 @@ SipReverseTransform SipReverseTransform::convert(ScaledPolynomialTransform const
);
}

SipReverseTransform SipReverseTransform::extract(afw::image::TanWcs const & wcs) {
if (!wcs.hasDistortion()) {
throw LSST_EXCEPT(
pex::exceptions::LogicError,
"Constructing a SipReverseTransform from a TanWcs with no distortions is not implemented."
);
}
ndarray::Array<double,2,2> xCoeffs = ndarray::allocate(wcs.getSipAp().rows(), wcs.getSipAp().cols());
ndarray::Array<double,2,2> yCoeffs = ndarray::allocate(wcs.getSipBp().rows(), wcs.getSipBp().cols());
xCoeffs.asEigen() = wcs.getSipAp();
yCoeffs.asEigen() = wcs.getSipBp();
return SipReverseTransform(
wcs.getPixelOrigin(),
afw::geom::LinearTransform(wcs.getCDMatrix()),
PolynomialTransform(xCoeffs, yCoeffs)
);
}

SipReverseTransform SipReverseTransform::transformPixels(afw::geom::AffineTransform const & s) const {
SipReverseTransform result(*this);
result.transformPixelsInPlace(s);
result._cdInverse = result._cdMatrix.invert();
return result;
}

afw::geom::AffineTransform SipReverseTransform::linearize(afw::geom::Point2D const & in) const {
return afw::geom::AffineTransform(afw::geom::Extent2D(getPixelOrigin()))
* (afw::geom::AffineTransform() + _poly.linearize(_cdInverse(in)))
Expand Down Expand Up @@ -199,4 +260,53 @@ PTR(afw::image::TanWcs) makeWcs(
);
}

std::shared_ptr<afw::image::TanWcs> transformWcsPixels(
afw::image::TanWcs const & wcs,
afw::geom::AffineTransform const & s
) {
if (wcs.hasDistortion()) {
auto fwd = SipForwardTransform::extract(wcs).transformPixels(s);
auto rev = SipReverseTransform::extract(wcs).transformPixels(s);
return makeWcs(fwd, rev, *wcs.getSkyOrigin());
} else {
auto sInv = s.invert();
auto pixelOrigin = s.getLinear()(wcs.getPixelOrigin() - sInv.getTranslation());
Eigen::Matrix2d cdMatrix = wcs.getCDMatrix() * sInv.getLinear().getMatrix();
return std::make_shared<afw::image::TanWcs>(
wcs.getSkyOrigin()->toIcrs().getPosition(afw::geom::degrees),
pixelOrigin,
cdMatrix,
wcs.getEquinox(),
"ICRS"
);
}
}

std::shared_ptr<afw::image::TanWcs> rotateWcsPixelsBy90(
afw::image::TanWcs const & wcs,
int nQuarter,
afw::geom::Extent2I const & dimensions
) {
afw::geom::Extent2D offset;
switch(nQuarter % 4) {
case 0:
offset = afw::geom::Extent2D(0, 0);
break;
case 1:
offset = afw::geom::Extent2D(dimensions.getY() - 1, 0);
break;
case 2:
offset = afw::geom::Extent2D(dimensions - afw::geom::Extent2I(1, 1));
break;
case 3:
offset = afw::geom::Extent2D(0, dimensions.getX() - 1);
break;
}
auto rot = afw::geom::LinearTransform::makeRotation(nQuarter*90.0*afw::geom::degrees);
return transformWcsPixels(
wcs,
afw::geom::AffineTransform(rot, offset)
);
}

}}} // namespace lsst::meas::astrom

0 comments on commit 0db1a30

Please sign in to comment.