DM-38659: Relax dipole fit tolerance to avoid segfaults, and add cut on footprint size.

python/lsst/ip/diffim/dipoleFitTask.py

@@ -92,6 +92,11 @@ class DipoleFitPluginConfig(measBase.SingleFramePluginConfig):
         Default value means \"Choose a chi2DoF corresponding to a significance level of at most 0.05\"
         (note this is actually a significance, not a chi2 value).""")
 
+    maxFootprintArea = pexConfig.Field(
+        dtype=int, default=5_000,
+        doc=("Maximum area for footprints before they are ignored as large; "
+             "non-positive means no threshold applied"))
+
 
 class DipoleFitTaskConfig(measBase.SingleFrameMeasurementConfig):
     """Measurement of detected diaSources as dipoles
@@ -755,7 +760,7 @@ def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None
                               # see scipy docs for the meaning of these keywords
                               fit_kws={'ftol': tol, 'xtol': tol, 'gtol': tol,
                                        # Our model is float32 internally, so we need a larger epsfcn.
-                                       'epsfcn': 1e-10},
+                                       'epsfcn': 1e-8},
                               psf=self.diffim.getPsf(),  # hereon: kwargs that get passed to makeModel()
                               rel_weight=rel_weight,
                               footprint=fp,
@@ -1089,9 +1094,13 @@ def measure(self, measRecord, exposure, posExp=None, negExp=None):
 
         # Check if the footprint consists of a putative dipole - else don't fit it.
         if (
-                (len(pks) <= 1)  # one peak in the footprint - not a dipole
+                # One peak in the footprint (not a dipole)
+                (len(pks) <= 1)
+                # Peaks are the same sign (not a dipole)
                 or (len(pks) > 1 and (np.sign(pks[0].getPeakValue())
-                    == np.sign(pks[-1].getPeakValue())))  # peaks are same sign - not a dipole
+                    == np.sign(pks[-1].getPeakValue())))
+                # Footprint is too large (not a dipole)
+                or (measRecord.getFootprint().getArea() > self.config.maxFootprintArea)
         ):
             measRecord.set(self.classificationFlagKey, False)
             measRecord.set(self.classificationAttemptedFlagKey, False)

tests/test_dipoleFitter.py

@@ -119,7 +119,7 @@ def testDipoleAlgorithm(self):
             self.assertFloatsAlmostEqual(result.negCentroidX, dipoleTestImage.xc[i] - offsets[i], rtol=rtol)
             self.assertFloatsAlmostEqual(result.negCentroidY, dipoleTestImage.yc[i] - offsets[i], rtol=rtol)
 
-    def _runDetection(self, dipoleTestImage):
+    def _runDetection(self, dipoleTestImage, maxFootprintArea=None):
         """Run 'diaSource' detection on the diffim, including merging of
         positive and negative sources.
 
@@ -150,6 +150,9 @@ def _runDetection(self, dipoleTestImage):
         # Here is where we make the dipole fitting task. It can run the other measurements as well.
         measureTask = DipoleFitTask(config=measureConfig, schema=schema)
 
+        if maxFootprintArea:
+            measureTask.config.plugins["ip_diffim_DipoleFit"].maxFootprintArea = maxFootprintArea
+
         table = afwTable.SourceTable.make(schema)
         detectResult = detectTask.run(table, testImage.diffim)
         # catalog = detectResult.sources
@@ -306,6 +309,15 @@ def testDipoleEdge(self):
             result = s.extract("ip_diffim_DipoleFit*")
             self.assertTrue(result.get("ip_diffim_DipoleFit_flag"))
 
+    def testDipoleFootprintTooLarge(self):
+        """Test that the footprint area cut flags sources."""
+
+        dipoleTestImage = DipoleTestImage()
+        # This area is smaller than the area of the test sources (~750).
+        sources = self._runDetection(dipoleTestImage, maxFootprintArea=500)
+
+        self.assertTrue(np.all(sources["ip_diffim_DipoleFit_flag"]))
+
 
 class TestMemory(lsst.utils.tests.MemoryTestCase):
     pass