Initial version of the near edge fallback on SdssCentroid.

lsst · Nov 16, 2022 · 930f269 · 930f269
1 parent e79e97c
commit 930f269
Show file tree

Hide file tree

Showing 2 changed files with 25 additions and 4 deletions.
diff --git a/include/lsst/meas/base/SdssCentroid.h b/include/lsst/meas/base/SdssCentroid.h
@@ -75,6 +75,7 @@ class SdssCentroidAlgorithm : public SimpleAlgorithm {
     static FlagDefinition const NO_SECOND_DERIVATIVE;
     static FlagDefinition const ALMOST_NO_SECOND_DERIVATIVE;
     static FlagDefinition const NOT_AT_MAXIMUM;
+    static FlagDefinition const NEAR_EDGE;
 
     /// A typedef to the Control object for this algorithm, defined above.
     /// The control object contains the configuration parameters for this algorithm.

diff --git a/src/SdssCentroid.cc b/src/SdssCentroid.cc
@@ -26,6 +26,7 @@
 #include "ndarray/eigen.h"
 #include "lsst/geom/Angle.h"
 #include "lsst/afw/detection/Psf.h"
+#include "lsst/afw/detection/GaussianPsf.h"
 #include "lsst/pex/exceptions.h"
 #include "lsst/afw/math/ConvolveImage.h"
 #include "lsst/afw/math/offsetImage.h"
@@ -41,13 +42,15 @@ FlagDefinitionList flagDefinitions;
 
 FlagDefinition const SdssCentroidAlgorithm::FAILURE = flagDefinitions.addFailureFlag();
 FlagDefinition const SdssCentroidAlgorithm::EDGE =
-        flagDefinitions.add("flag_edge", "Object too close to edge");
+        flagDefinitions.add("flag_edge", "Object too close to edge; peak used.");
 FlagDefinition const SdssCentroidAlgorithm::NO_SECOND_DERIVATIVE =
         flagDefinitions.add("flag_noSecondDerivative", "Vanishing second derivative");
 FlagDefinition const SdssCentroidAlgorithm::ALMOST_NO_SECOND_DERIVATIVE =
         flagDefinitions.add("flag_almostNoSecondDerivative", "Almost vanishing second derivative");
 FlagDefinition const SdssCentroidAlgorithm::NOT_AT_MAXIMUM =
         flagDefinitions.add("flag_notAtMaximum", "Object is not at a maximum");
+FlagDefinition const SdssCentroidAlgorithm::NEAR_EDGE =
+        flagDefinitions.add("flag_near_edge", "Object close to edge; fallback kernel used.");
 
 FlagDefinitionList const &SdssCentroidAlgorithm::getFlagDefinitions() { return flagDefinitions; }
 
@@ -434,14 +437,27 @@ void SdssCentroidAlgorithm::measure(afw::table::SourceRecord &measRecord,
     int binX = 1;
     int binY = 1;
     double xc = 0., yc = 0., dxc = 0., dyc = 0.;  // estimated centre and error therein
+    bool stopBinning = false;
     for (int binsize = 1; binsize <= _ctrl.binmax; binsize *= 2) {
         std::tuple<MaskedImageT, double, int> smoothResult =
             smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
         int errorFlag = std::get<2>(smoothResult);
-        if (errorFlag > 0) {
+        if (errorFlag == EDGE.number) {
+            psf = std::make_shared<afw::detection::GaussianPsf>(5, 5, 0.5);
+            smoothResult = smoothAndBinImage(psf, x, y, mimage, binX, binY, _flagHandler);
+            stopBinning = true;
+            errorFlag = std::get<2>(smoothResult);
+            if (errorFlag == 0) {
+                errorFlag = NEAR_EDGE.number;
+            }
+        }
+        if (errorFlag > 0 && errorFlag != NEAR_EDGE.number) {
             _flagHandler.setValue(measRecord, errorFlag, true);
-            _flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
-            return;
+            // if NEAR_EDGE is not a fatal error we continue otherwise return
+            if (errorFlag != NEAR_EDGE.number) {
+                _flagHandler.setValue(measRecord, SdssCentroidAlgorithm::FAILURE.number, true);
+                return;
+            }
         }
         MaskedImageT const smoothedImage = std::get<0>(smoothResult);
         double const smoothingSigma = std::get<1>(smoothResult);
@@ -474,6 +490,10 @@ void SdssCentroidAlgorithm::measure(afw::table::SourceRecord &measRecord,
         xc += x;  // xc, yc are measured relative to pixel (x, y)
         yc += y;
 
+        if (stopBinning) {
+            break;
+        }
+
         double const fac = _ctrl.wfac * (1 + smoothingSigma * smoothingSigma);
         double const facX2 = fac * binX * binX;
         double const facY2 = fac * binY * binY;