Review for meas_base DM-2307: transformations for shape measurements

include/lsst/meas/base/SdssShape.h

@@ -225,6 +225,38 @@ class SdssShapeResult : public ShapeResult, public CentroidResult, public FluxRe
 
 };
 
+/**
+ *  Transformation for SdssShape measurements.
+ *
+ *  SdssShape measures not just shape but also flux and centroid. This
+ *  transform operates on the first directly, and delegates to the Flux and
+ *  Centroid transforms for the other two.
+ */
+class SdssShapeTransform : public BaseTransform {
+public:
+    typedef SdssShapeControl Control;
+
+    SdssShapeTransform(Control const & ctrl, std::string const & name, afw::table::SchemaMapper & mapper);
+
+    /*
+     * @brief Perform transformation from inputCatalog to outputCatalog.
+     *
+     * @param[in]     inputCatalog   Source of data to be transformed
+     * @param[in,out] outputCatalog  Container for transformed results
+     * @param[in]     wcs            World coordinate system under which transformation will take place
+     * @param[in]     calib          Photometric calibration under which transformation will take place
+     * @throws        LengthError    Catalog sizes do not match
+     */
+    virtual void operator()(afw::table::SourceCatalog const & inputCatalog,
+                            afw::table::BaseCatalog & outputCatalog,
+                            afw::image::Wcs const & wcs,
+                            afw::image::Calib const & calib) const;
+private:
+    FluxTransform _fluxTransform;
+    CentroidTransform _centroidTransform;
+    ShapeResultKey _outShapeKey;
+};
+
 }}} // namespace lsst::meas::base
 
 #endif // !LSST_MEAS_BASE_SdssShape_h_INCLUDED

include/lsst/meas/base/ShapeUtilities.h

@@ -1,7 +1,7 @@
 // -*- lsst-c++ -*-
 /*
  * LSST Data Management System
- * Copyright 2008-2014 LSST Corporation.
+ * Copyright 2008-2015 AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -130,7 +130,8 @@ class ShapeResultKey : public afw::table::FunctorKey<ShapeResult>  {
         afw::table::Schema & schema,
         std::string const & name,
         std::string const & doc,
-        UncertaintyEnum uncertainty
+        UncertaintyEnum uncertainty,
+        afw::table::CoordinateType coordType=afw::table::CoordinateType::PIXEL
     );
 
     /// Default constructor; instance will not be usuable unless subsequently assigned to.
@@ -192,6 +193,28 @@ class ShapeResultKey : public afw::table::FunctorKey<ShapeResult>  {
     afw::table::CovarianceMatrixKey<ErrElement,3> _shapeErr;
 };
 
+/**
+ *  @brief Construct a matrix suitable for transforming second moments.
+ *
+ *  Given an LinearTransform which maps from positions (x, y) to (a, d),
+ *  returns a 3-by-3 matrix which transforms (xx, yy, xy) to (aa, dd, ad).
+ *
+ *  That is, given an input transform described by the matrix
+ *
+ *  | A11 | A12 |
+ *  | A21 | A22 |
+ *
+ *  we return the matrix
+ *
+ *  | A11*A11 | A12*A12 |         2*A11*A12 |
+ *  | A21*A21 | A22*A22 |         2*A21*A22 |
+ *  | A11*A21 | A12*A22 | A11*A22 + A12*A21 |
+ *
+ *  @param[in]  xform  LinearTransform describing the coordinate mapping
+ *  @returns    A 3-by-3 transformation matrix for the second order moments
+ */
+ShapeTrMatrix makeShapeTransformMatrix(afw::geom::LinearTransform const & xform);
+
 }}} // lsst::meas::base
 
 #endif // !LSST_MEAS_BASE_ShapeUtilities_h_INCLUDED

include/lsst/meas/base/constants.h

@@ -57,6 +57,7 @@ typedef afw::geom::Point<CentroidElement,2> Centroid;
 typedef Eigen::Matrix<ErrElement,2,2,Eigen::DontAlign> CentroidCov;
 typedef afw::geom::ellipses::Quadrupole Shape;
 typedef Eigen::Matrix<ErrElement,3,3,Eigen::DontAlign> ShapeCov;
+typedef Eigen::Matrix<ShapeElement,3,3,Eigen::DontAlign> ShapeTrMatrix;
 //@}
 
 }}} // lsst::meas::base

python/lsst/meas/base/plugins.py

@@ -52,7 +52,8 @@
 wrapSimpleAlgorithm(SdssCentroidAlgorithm, Control=SdssCentroidControl,
                     TransformClass=SdssCentroidTransform, executionOrder=0.0)
 wrapSimpleAlgorithm(PixelFlagsAlgorithm, Control=PixelFlagsControl, executionOrder=2.0)
-wrapSimpleAlgorithm(SdssShapeAlgorithm, Control=SdssShapeControl, executionOrder=1.0)
+wrapSimpleAlgorithm(SdssShapeAlgorithm, Control=SdssShapeControl,
+                    TransformClass=SdssShapeTransform, executionOrder=1.0)
 
 wrapSimpleAlgorithm(CircularApertureFluxAlgorithm, needsMetadata=True, Control=ApertureFluxControl,
                     TransformClass=ApertureFluxTransform, executionOrder=2.0)

python/lsst/meas/base/pluginsLib.i

@@ -44,6 +44,8 @@
 %convertConfig(lsst::meas::base, GaussianCentroidTransform)
 %convertConfig(lsst::meas::base, SdssCentroidTransform)
 
+%convertConfig(lsst::meas::base, SdssShapeTransform)
+
 %include "lsst/meas/base/Transform.h"
 
 // flux algorithms

python/lsst/meas/base/tests.py

@@ -157,7 +157,7 @@ def makeMinimalSchema(cls):
                 schema, "truth", "true simulated centroid", "pixels"
             )
             cls.keys["shape"] = lsst.afw.table.QuadrupoleKey.addFields(
-                schema, "truth", "true shape after PSF convolution", "pixels^2"
+                schema, "truth", "true shape after PSF convolution", lsst.afw.table.PIXEL
             )
             cls.keys["isStar"] = schema.addField("truth_isStar", type="Flag",
                                                  doc="set if the object is a star")

python/lsst/meas/base/utilities.i

@@ -1,7 +1,7 @@
 // -*- lsst-c++ -*-
 /*
  * LSST Data Management System
- * Copyright 2008-2014 LSST Corporation.
+ * Copyright 2008-2015 AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -37,6 +37,7 @@
 
 %declareNumPyConverters(lsst::meas::base::CentroidCov);
 %declareNumPyConverters(lsst::meas::base::ShapeCov);
+%declareNumPyConverters(lsst::meas::base::ShapeTrMatrix);
 
 %declareFunctorKey(FluxResult, lsst::meas::base::FluxResult)
 %shared_ptr(lsst::meas::base::FluxResultKey)

src/SdssShape.cc

@@ -882,5 +882,55 @@ INSTANTIATE_PIXEL(int);
 INSTANTIATE_PIXEL(float);
 INSTANTIATE_PIXEL(double);
 
-}}} // end namespace lsst::meas::base
+SdssShapeTransform::SdssShapeTransform(
+    Control const & ctrl,
+    std::string const & name,
+    afw::table::SchemaMapper & mapper
+) :
+    BaseTransform{name},
+    _fluxTransform{name, mapper},
+    _centroidTransform{name, mapper}
+{
+    for (auto flag = flagDefs.begin() + 1; flag < flagDefs.end(); flag++) {
+        mapper.addMapping(mapper.getInputSchema().find<afw::table::Flag>(
+                          mapper.getInputSchema().join(name, flag->name)).key);
+    }
+
+    _outShapeKey = ShapeResultKey::addFields(mapper.editOutputSchema(), name, "Shape in celestial moments",
+                                             FULL_COVARIANCE, afw::table::CoordinateType::CELESTIAL);
+}
+
+void SdssShapeTransform::operator()(
+    afw::table::SourceCatalog const & inputCatalog,
+    afw::table::BaseCatalog & outputCatalog,
+    afw::image::Wcs const & wcs,
+    afw::image::Calib const & calib
+) const {
+    // The flux and cetroid transforms will check that the catalog lengths
+    // match and throw if not, so we don't repeat the test here.
+    _fluxTransform(inputCatalog, outputCatalog, wcs, calib);
+    _centroidTransform(inputCatalog, outputCatalog, wcs, calib);
+
+    CentroidResultKey centroidKey(inputCatalog.getSchema()[_name]);
+    ShapeResultKey inShapeKey(inputCatalog.getSchema()[_name]);
+
+    afw::table::SourceCatalog::const_iterator inSrc = inputCatalog.begin();
+    afw::table::BaseCatalog::iterator outSrc = outputCatalog.begin();
+    for (; inSrc < inputCatalog.end(); ++inSrc, ++outSrc) {
+        ShapeResult inShape = inShapeKey.get(*inSrc);
+        ShapeResult outShape;
+
+        // The transformation from the (x, y) to the (Ra, Dec) basis.
+        afw::geom::AffineTransform crdTr = wcs.linearizePixelToSky(centroidKey.get(*inSrc).getCentroid(),
+                                                                   afw::geom::radians);
+        outShape.setShape(inShape.getShape().transform(crdTr.getLinear()));
+
+        // Transformation matrix from pixel to celestial basis.
+        ShapeTrMatrix m = makeShapeTransformMatrix(crdTr.getLinear());
+        outShape.setShapeErr((m*inShape.getShapeErr().cast<double>()*m.transpose()).cast<ErrElement>());
+
+        _outShapeKey.set(*outSrc, outShape);
+    }
+}
 
+}}} // end namespace lsst::meas::base

src/ShapeUtilities.cc

@@ -1,6 +1,6 @@
 /*
  * LSST Data Management System
- * Copyright 2008-2014 LSST Corporation.
+ * Copyright 2008-2015 AURA/LSST.
  *
  * This product includes software developed by the
  * LSST Project (http://www.lsst.org/).
@@ -76,41 +76,48 @@ ShapeResultKey ShapeResultKey::addFields(
     afw::table::Schema & schema,
     std::string const & name,
     std::string const & doc,
-    UncertaintyEnum uncertainty
+    UncertaintyEnum uncertainty,
+    afw::table::CoordinateType coordType
 ) {
     ShapeResultKey r;
     r._shape = afw::table::QuadrupoleKey::addFields(
         schema,
         name,
         doc,
-        "pixels^2"
+        coordType
     );
     if (uncertainty != NO_UNCERTAINTY) {
         std::vector< afw::table::Key<ErrElement> > sigma(3);
         std::vector< afw::table::Key<ErrElement> > cov;
         sigma[0] = schema.addField<ErrElement>(
-            schema.join(name, "xxSigma"), "1-sigma uncertainty on xx moment", "pixels^2"
+            schema.join(name, "xxSigma"), "1-sigma uncertainty on xx moment",
+            coordType == afw::table::CoordinateType::PIXEL ? "pixels^2" : "radians^2"
         );
         sigma[1] = schema.addField<ErrElement>(
-            schema.join(name, "yySigma"), "1-sigma uncertainty on yy moment", "pixels^2"
+            schema.join(name, "yySigma"), "1-sigma uncertainty on yy moment",
+            coordType == afw::table::CoordinateType::PIXEL ? "pixels^2" : "radians^2"
         );
         sigma[2] = schema.addField<ErrElement>(
-            schema.join(name, "xySigma"), "1-sigma uncertainty on xy moment", "pixels^2"
+            schema.join(name, "xySigma"), "1-sigma uncertainty on xy moment",
+            coordType == afw::table::CoordinateType::PIXEL ? "pixels^2" : "radians^2"
         );
         if (uncertainty == FULL_COVARIANCE) {
             cov.push_back(
                 schema.addField<ErrElement>(
-                    schema.join(name, "xx_yy_Cov"), "uncertainty covariance in xx and yy", "pixels^4"
+                    schema.join(name, "xx_yy_Cov"), "uncertainty covariance in xx and yy",
+                    coordType == afw::table::CoordinateType::PIXEL ? "pixels^4" : "radians^4"
                 )
             );
             cov.push_back(
                 schema.addField<ErrElement>(
-                    schema.join(name, "xx_xy_Cov"), "uncertainty covariance in xx and xy", "pixels^4"
+                    schema.join(name, "xx_xy_Cov"), "uncertainty covariance in xx and xy",
+                    coordType == afw::table::CoordinateType::PIXEL ? "pixels^4" : "radians^4"
                 )
             );
             cov.push_back(
                 schema.addField<ErrElement>(
-                    schema.join(name, "yy_xy_Cov"), "uncertainty covariance in yy and xy", "pixels^4"
+                    schema.join(name, "yy_xy_Cov"), "uncertainty covariance in yy and xy",
+                    coordType == afw::table::CoordinateType::PIXEL ? "pixels^4" : "radians^4"
                 )
             );
         }
@@ -156,4 +163,14 @@ void ShapeResultKey::set(afw::table::BaseRecord & record, ShapeResult const & va
     }
 }
 
+ShapeTrMatrix makeShapeTransformMatrix(afw::geom::LinearTransform const & xform) {
+    typedef afw::geom::LinearTransform LT;
+    Eigen::Matrix<ShapeElement,3,3,Eigen::DontAlign> m;
+    m << xform[LT::XX]*xform[LT::XX], xform[LT::XY]*xform[LT::XY], 2*xform[LT::XX]*xform[LT::XY],
+         xform[LT::YX]*xform[LT::YX], xform[LT::YY]*xform[LT::YY], 2*xform[LT::YX]*xform[LT::YY],
+         xform[LT::XX]*xform[LT::YX], xform[LT::XY]*xform[LT::YY],
+         xform[LT::XX]*xform[LT::YY] + xform[LT::XY]*xform[LT::YX];
+    return m;
+}
+
 }}} // lsst::meas::base

tests/testSdssShape.py

@@ -27,8 +27,10 @@
 
 import lsst.utils.tests
 import lsst.meas.base.tests
+from lsst.meas.base.tests import (AlgorithmTestCase, FluxTransformTestCase,
+                                  CentroidTransformTestCase, SingleFramePluginTransformSetupHelper)
 
-class SdssShapeTestCase(lsst.meas.base.tests.AlgorithmTestCase):
+class SdssShapeTestCase(AlgorithmTestCase):
 
     def setUp(self):
         self.bbox = lsst.afw.geom.Box2I(lsst.afw.geom.Point2I(-20, -30),
@@ -73,13 +75,51 @@ def testGaussians(self):
             self.assertFinite(result.yy_xy_Cov)
 
 
+class SdssShapeTransformTestCase(FluxTransformTestCase, CentroidTransformTestCase,
+                                 SingleFramePluginTransformSetupHelper):
+    name = "sdssShape"
+    controlClass = lsst.meas.base.SdssShapeControl
+    algorithmClass = lsst.meas.base.SdssShapeAlgorithm
+    transformClass = lsst.meas.base.SdssShapeTransform
+    flagNames = ("flag", "flag_unweighted", "flag_unweightedBad", "flag_shift", "flag_maxIter")
+    singleFramePlugins = ('base_SdssShape',)
+    forcedPlugins = ('base_SdssShape',)
+
+    def _setFieldsInRecord(self, record, name):
+        FluxTransformTestCase._setFieldsInRecord(self, record, name)
+        CentroidTransformTestCase._setFieldsInRecord(self, record, name)
+        for field in ('xx', 'yy', 'xy', 'xxSigma', 'yySigma', 'xySigma'):
+            record[record.schema.join(name, field)] = numpy.random.random()
+
+    def _compareFieldsInRecords(self, inSrc, outSrc, name):
+        FluxTransformTestCase._compareFieldsInRecords(self, inSrc, outSrc, name)
+        CentroidTransformTestCase._compareFieldsInRecords(self, inSrc, outSrc, name)
+
+        inShape = lsst.meas.base.ShapeResultKey(inSrc.schema[name]).get(inSrc)
+        outShape = lsst.meas.base.ShapeResultKey(outSrc.schema[name]).get(outSrc)
+
+        centroid = lsst.meas.base.CentroidResultKey(inSrc.schema[name]).get(inSrc).getCentroid()
+        xform = self.calexp.getWcs().linearizePixelToSky(centroid, lsst.afw.geom.radians)
+
+        trInShape = inShape.getShape().transform(xform.getLinear())
+        self.assertEqual(trInShape.getIxx(), outShape.getShape().getIxx())
+        self.assertEqual(trInShape.getIyy(), outShape.getShape().getIyy())
+        self.assertEqual(trInShape.getIxy(), outShape.getShape().getIxy())
+
+        m = lsst.meas.base.makeShapeTransformMatrix(xform.getLinear())
+        numpy.testing.assert_array_almost_equal(
+            numpy.dot(numpy.dot(m, inShape.getShapeErr()), m.transpose()), outShape.getShapeErr()
+        )
+
+
 def suite():
     """Returns a suite containing all the test cases in this module."""
 
     lsst.utils.tests.init()
 
     suites = []
     suites += unittest.makeSuite(SdssShapeTestCase)
+    suites += unittest.makeSuite(SdssShapeTransformTestCase)
     suites += unittest.makeSuite(lsst.utils.tests.MemoryTestCase)
     return unittest.TestSuite(suites)