matchOptimisticB: correct for distortion DM-3492

include/lsst/meas/astrom/matchOptimisticB.h

@@ -52,9 +52,14 @@ namespace astrom {
 
     typedef std::vector<RecordProxy> ProxyVector;
 
-    ProxyVector makeProxies(lsst::afw::table::SourceCatalog const & sourceCat);
+    ProxyVector makeProxies(lsst::afw::table::SourceCatalog const & sourceCat,
+                            afw::image::Wcs const& distortedWcs,
+                            afw::image::Wcs const& tanWcs
+        );
 
-    ProxyVector makeProxies(lsst::afw::table::SimpleCatalog const & posRefCat);
+    ProxyVector makeProxies(lsst::afw::table::SimpleCatalog const & posRefCat,
+                            afw::image::Wcs const& tanWcs
+        );
 
     struct ProxyPair {
         RecordProxy first;
@@ -111,13 +116,13 @@ namespace astrom {
     Optimistic Pattern Matching is described in V. Tabur 2007, PASA, 24, 189-198
     "Fast Algorithms for Matching CCD Images to a Stellar Catalogue"
 
-    @param[in] posRefCat  catalog of position reference stars
-        fields that are used:
-        coord
-        centroid
-        hasCentroid
-        control.refFluxField  flux for desired filter
-    @param[in] sourceCat  catalog of detected sources
+    @param[in] posRefCat  catalog of position reference stars; fields read:
+        - "coord"
+        - control.refFluxField
+    @param[in] sourceCat  catalog of detected sources; fields read:
+        - "Centroid_x"
+        - "Centroid_y"
+        - control.refFluxField
     @param[in] wcs  estimated WCS
     @param[in] control  control object
     @param[in] posRefBegInd  index of first start to use in posRefCat
@@ -128,6 +133,7 @@ namespace astrom {
         lsst::afw::table::SimpleCatalog const &posRefCat,
         lsst::afw::table::SourceCatalog const &sourceCat,
         MatchOptimisticBControl const &control,
+        afw::image::Wcs const& wcs,
         int posRefBegInd=0,
         bool verbose = false
     );

python/lsst/meas/astrom/fitTanSipWcs.py

@@ -1,5 +1,8 @@
 from __future__ import absolute_import, division, print_function
 
+import numpy
+import lsst.afw.image as afwImage
+import lsst.afw.coord as afwCoord
 import lsst.afw.geom as afwGeom
 import lsst.afw.table as afwTable
 import lsst.pex.config as pexConfig
@@ -23,6 +26,18 @@ class FitTanSipWcsConfig(pexConfig.Config):
         default = 3,
         min = 1,
     )
+    numRejIter = pexConfig.RangeField(
+        doc = "number of rejection iterations",
+        dtype = int,
+        default = 1,
+        min = 0,
+    )
+    rejSigma = pexConfig.RangeField(
+        doc = "Number of standard deviations for clipping level",
+        dtype = float,
+        default = 3.0,
+        min = 0.0,
+    )
     maxScatterArcsec = pexConfig.RangeField(
         doc = "maximum median scatter of a WCS fit beyond which the fit fails (arcsec); " +
             "be generous, as this is only intended to catch catastrophic failures",
@@ -115,10 +130,27 @@ def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None):
         """
         if bbox is None:
             bbox = afwGeom.Box2I()
-        wcs = initWcs
-        for i in range(self.config.numIter):
-            sipObject = makeCreateWcsWithSip(matches, wcs, self.config.order, bbox)
+
+        import lsstDebug
+        debug = lsstDebug.Info(__name__)
+
+        wcs = self.initialWcs(matches, initWcs)
+        rejected = numpy.zeros(len(matches), dtype=bool)
+        for rej in range(self.config.numRejIter):
+            sipObject = self._fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
             wcs = sipObject.getNewWcs()
+            rejected = self.rejectMatches(matches, wcs, rejected)
+            if rejected.sum() == len(rejected):
+                raise RuntimeError("All matches rejected in iteration %d" % (rej + 1,))
+            if debug.plot:
+                print("Plotting fit after rejection iteration %d/%d" % (rej + 1, self.config.numRejIter))
+                self.plotFit(matches, wcs, rejected)
+        # Final fit after rejection
+        sipObject = self._fitWcs([mm for i, mm in enumerate(matches) if not rejected[i]], wcs)
+        wcs = sipObject.getNewWcs()
+        if debug.plot:
+            print("Plotting final fit")
+            self.plotFit(matches, wcs, rejected)
 
         if refCat is not None:
             self.log.info("Updating centroids in refCat")
@@ -149,6 +181,46 @@ def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None):
             scatterOnSky = scatterOnSky,
         )
 
+    def initialWcs(self, matches, wcs):
+        """Generate a guess Wcs from the astrometric matches
+
+        We create a Wcs anchored at the center of the matches, with the scale
+        of the input Wcs.  This is necessary because matching returns only
+        matches with no estimated Wcs, and the input Wcs is a wild guess.
+        We're using the best of each: positions from the matches, and scale
+        from the input Wcs.
+        """
+        crpix = afwGeom.Extent2D(0, 0)
+        crval = afwGeom.Extent3D(0, 0, 0)
+        for mm in matches:
+            crpix += afwGeom.Extent2D(mm.second.getCentroid())
+            crval += afwGeom.Extent3D(mm.first.getCoord().toIcrs().getVector())
+        crpix /= len(matches)
+        crval /= len(matches)
+        newWcs = afwImage.Wcs(afwCoord.IcrsCoord(afwGeom.Point3D(crval)).getPosition(),
+                              afwGeom.Point2D(crpix), wcs.getCDMatrix())
+        return newWcs
+
+    def _fitWcs(self, matches, wcs):
+        """Fit a Wcs based on the matches and a guess Wcs"""
+        for i in range(self.config.numIter):
+            sipObject = makeCreateWcsWithSip(matches, wcs, self.config.order)
+            wcs = sipObject.getNewWcs()
+        return sipObject
+
+    def rejectMatches(self, matches, wcs, rejected):
+        """Flag deviant matches
+
+        We return a boolean numpy array indicating whether the corresponding
+        match should be rejected.  The previous list of rejections is used
+        so we can calculate uncontaminated statistics.
+        """
+        fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
+        dx = numpy.array([ff.getX() - mm.second.getCentroid().getX() for ff, mm in zip(fit, matches)])
+        dy = numpy.array([ff.getY() - mm.second.getCentroid().getY() for ff, mm in zip(fit, matches)])
+        good = numpy.logical_not(rejected)
+        return (dx > self.config.rejSigma*dx[good].std()) | (dy > self.config.rejSigma*dy[good].std())
+
     @staticmethod
     def updateRefCentroids(wcs, refList):
         """Update centroids in a collection of reference objects, given a WCS
@@ -165,10 +237,55 @@ def updateRefCentroids(wcs, refList):
     def updateSourceCoords(wcs, sourceList):
         """Update coords in a collection of sources, given a WCS
         """
-
         if len(sourceList) < 1:
             return
         schema = sourceList[1].schema
         srcCoordKey = afwTable.CoordKey(schema["coord"])
         for src in sourceList:
             src.set(srcCoordKey, wcs.pixelToSky(src.getCentroid()))
+
+    def plotFit(self, matches, wcs, rejected):
+        """Plot the fit
+
+        We create four plots, for all combinations of (dx, dy) against
+        (x, y).  Good points are black, while rejected points are red.
+        """
+        import matplotlib.pyplot as plt
+
+        fit = [wcs.skyToPixel(m.first.getCoord()) for m in matches]
+        x1 = numpy.array([ff.getX() for ff in fit])
+        y1 = numpy.array([ff.getY() for ff in fit])
+        x2 = numpy.array([m.second.getCentroid().getX() for m in matches])
+        y2 = numpy.array([m.second.getCentroid().getY() for m in matches])
+
+        dx = x1 - x2
+        dy = y1 - y2
+
+        good = numpy.logical_not(rejected)
+
+        figure = plt.figure()
+        axes = figure.add_subplot(2, 2, 1)
+        axes.plot(x2[good], dx[good], 'ko')
+        axes.plot(x2[rejected], dx[rejected], 'ro')
+        axes.set_xlabel("x")
+        axes.set_ylabel("dx")
+
+        axes = figure.add_subplot(2, 2, 2)
+        axes.plot(x2[good], dy[good], 'ko')
+        axes.plot(x2[rejected], dy[rejected], 'ro')
+        axes.set_xlabel("x")
+        axes.set_ylabel("dy")
+
+        axes = figure.add_subplot(2, 2, 3)
+        axes.plot(y2[good], dx[good], 'ko')
+        axes.plot(y2[rejected], dx[rejected], 'ro')
+        axes.set_xlabel("y")
+        axes.set_ylabel("dx")
+
+        axes = figure.add_subplot(2, 2, 4)
+        axes.plot(y2[good], dy[good], 'ko')
+        axes.plot(y2[rejected], dy[rejected], 'ro')
+        axes.set_xlabel("y")
+        axes.set_ylabel("dy")
+
+        plt.show()

python/lsst/meas/astrom/matchOptimisticB.py

@@ -269,7 +269,7 @@ def matchObjectsToSources(self, refCat, sourceCat, wcs, refFluxField, maxMatchDi
         usableSourceCat = afwTable.SourceCatalog(sourceCat.table)
         usableSourceCat.extend(s for s in sourceCat if sourceInfo.isUsable(s))
         numUsableSources = len(usableSourceCat)
-        self.log.info("Purged %d unsuable sources, leaving %d usable sources" % \
+        self.log.info("Purged %d unusable sources, leaving %d usable sources" % \
             (numSources - numUsableSources, numUsableSources))
 
         if len(usableSourceCat) == 0:
@@ -364,6 +364,7 @@ def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources, minMa
                         refCat,
                         sourceCat,
                         matchControl,
+                        wcs,
                         posRefBegInd,
                         verbose,
                     )

src/matchOptimisticB.cc

@@ -17,6 +17,7 @@
 #include "lsst/utils/ieee.h"
 #include "lsst/pex/exceptions.h"
 #include "lsst/afw/image/Wcs.h"
+#include "lsst/afw/image/DistortedTanWcs.h"
 #include "lsst/afw/geom/Angle.h"
 #include "lsst/meas/astrom/matchOptimisticB.h"
 
@@ -488,27 +489,31 @@ namespace astrom {
         }
     }
 
-    ProxyVector makeProxies(afwTable::SourceCatalog const & sourceCat) {
+    ProxyVector makeProxies(afwTable::SourceCatalog const & sourceCat,
+                            afw::image::Wcs const& distortedWcs,
+                            afw::image::Wcs const& tanWcs
+                            )
+    {
         ProxyVector r;
         r.reserve(sourceCat.size());
         for (afwTable::SourceCatalog::const_iterator sourcePtr = sourceCat.begin();
             sourcePtr != sourceCat.end(); ++sourcePtr) {
-            r.push_back(RecordProxy(sourcePtr, sourcePtr->getCentroid()));
+            r.push_back(RecordProxy(sourcePtr,
+                                    tanWcs.skyToPixel(*distortedWcs.pixelToSky(sourcePtr->getCentroid()))));
         }
         return r;
     }
 
-    ProxyVector makeProxies(afwTable::SimpleCatalog const & posRefCat) {
-        auto centroidKey = afwTable::PointKey<double>(posRefCat.getSchema()["centroid"]);
-        auto hasCentroidKey = posRefCat.getSchema().find<afwTable::Flag>("hasCentroid").key;
+    ProxyVector makeProxies(afwTable::SimpleCatalog const & posRefCat,
+                            afw::image::Wcs const& tanWcs
+                            )
+    {
+        auto coordKey = afwTable::CoordKey(posRefCat.getSchema()["coord"]);
         ProxyVector r;
         r.reserve(posRefCat.size());
         for (afwTable::SimpleCatalog::const_iterator posRefPtr = posRefCat.begin();
             posRefPtr != posRefCat.end(); ++posRefPtr) {
-            if (!posRefPtr->get(hasCentroidKey)) {
-                throw LSST_EXCEPT(pexExcept::InvalidParameterError, "unknown centroid");
-            }
-            r.push_back(RecordProxy(posRefPtr, posRefPtr->get(centroidKey)));
+            r.push_back(RecordProxy(posRefPtr, tanWcs.skyToPixel(posRefPtr->get(coordKey))));
         }
         return r;
     }
@@ -517,6 +522,7 @@ namespace astrom {
         afwTable::SimpleCatalog const &posRefCat,
         afwTable::SourceCatalog const &sourceCat,
         MatchOptimisticBControl const &control,
+        afw::image::Wcs const& wcs,
         int posRefBegInd,
         bool verbose
     ) {
@@ -528,10 +534,32 @@ namespace astrom {
             throw LSST_EXCEPT(pexExcept::InvalidParameterError, "posRefBegInd too big");
         }
         double const maxRotationRad = afw::geom::degToRad(control.maxRotationDeg);
-        double const allowedNonperpRad = afw::geom::degToRad(control.allowedNonperpDeg);
 
-        ProxyVector posRefProxyCat = makeProxies(posRefCat);
-        ProxyVector sourceProxyCat = makeProxies(sourceCat);
+        CONST_PTR(afw::image::Wcs) tanWcs; // Undistorted Wcs, providing tangent plane
+        if (wcs.hasDistortion()) {
+            // Create an undistorted Wcs to project everything with
+            // We'll anchor it at the center of the area.
+            afw::geom::Extent2D srcCenter(0, 0);
+            for (auto iter = sourceCat.begin(); iter != sourceCat.end(); ++iter) {
+                srcCenter += afw::geom::Extent2D(iter->getCentroid());
+            }
+            srcCenter /= sourceCat.size();
+
+            afw::geom::Extent3D refCenter(0, 0, 0);
+            for (auto iter = posRefCat.begin(); iter != posRefCat.end(); ++iter) {
+                refCenter += afw::geom::Extent3D(iter->getCoord().toIcrs().getVector());
+            }
+            refCenter /= posRefCat.size();
+
+            tanWcs = boost::make_shared<afw::image::Wcs>(
+                afw::coord::IcrsCoord(afw::geom::Point3D(refCenter)).getPosition(),
+                afw::geom::Point2D(srcCenter),
+                wcs.getCDMatrix()
+                );
+        }
+
+        ProxyVector posRefProxyCat = makeProxies(posRefCat, tanWcs ? *tanWcs : wcs);
+        ProxyVector sourceProxyCat = makeProxies(sourceCat, wcs, tanWcs ? *tanWcs : wcs);
 
         // sourceSubCat contains at most the numBrightStars brightest sources, sorted by decreasing flux
         ProxyVector sourceSubCat = selectPoint(
@@ -668,18 +696,22 @@ namespace astrom {
                         double b = coeff[2];
                         double c = coeff[4];
                         double d = coeff[5];
-                        double theta = std::acos((a*b+c*d)/(std::sqrt(a*a+c*c)*std::sqrt(b*b+d*d))) /
-                            M_PI * 180.0;
+                        afw::geom::Angle const theta(std::acos((a*b+c*d)/
+                                                               (std::sqrt(a*a+c*c)*std::sqrt(b*b+d*d))),
+                                                     afw::geom::radians);
                         if (verbose) {
                             std::cout << "Linear fit from match:" << std::endl;
                             std::cout << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
                             std::cout << coeff[3] << " " << coeff[4] << " " << coeff[5] << std::endl;
+                            std::cout << "Determinant (max " << control.maxDeterminant << "): ";
                             std::cout << coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1. << std::endl;
-                            std::cout << theta << std::endl;
+                            std::cout << "Angle between axes (deg; allowed 90 +/- ";
+                            std::cout << control.allowedNonperpDeg << "): ";
+                            std::cout << theta.asDegrees() << std::endl;
                         }
                         if (std::fabs(coeff[1] * coeff[5] - coeff[2] * coeff[4] - 1.) 
                                 > control.maxDeterminant ||
-                            std::fabs(theta - 90.) > allowedNonperpRad ||
+                            std::fabs(theta.asDegrees() - 90) > control.allowedNonperpDeg ||
                             std::fabs(coeff[0]) > control.maxOffsetPix ||
                             std::fabs(coeff[3]) > control.maxOffsetPix) {
                             if (verbose) {

tests/testAstrometryTask.py

@@ -91,6 +91,7 @@ def doTest(self, pixelsToTanPixels, order=3):
         sourceCat = self.makeSourceCat(distortedWcs)
         config = measAstrom.AstrometryTask.ConfigClass()
         config.wcsFitter.order = order
+        config.wcsFitter.numRejIter = 0
         solver = measAstrom.AstrometryTask(config=config)
         results = solver.run(
             sourceCat = sourceCat,