DM-57045: Fix bug when a dipole negative lobe is brighter

python/lsst/ip/diffim/dipoleFitTask.py

@@ -643,10 +643,12 @@ def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None
 
         pks = fp.getPeaks()
 
+        # If individual (merged) peaks were detected, start the lobes at the
+        # highest and the lowest peak. The peaks may be in any order.
         if len(pks) >= 1:
-            cenPos = pks[0].getF()    # if individual (merged) peaks were detected, use those
-        if len(pks) >= 2:    # peaks are already sorted by centroid flux so take the most negative one
-            cenNeg = pks[-1].getF()
+            cenPos = max(pks, key=lambda peak: peak.getPeakValue()).getF()
+        if len(pks) >= 2:
+            cenNeg = min(pks, key=lambda peak: peak.getPeakValue()).getF()
 
         # For close/faint dipoles the starting locs (min/max) might be way off, let's help them a bit.
         # First assume dipole is not separated by more than 5*psfSigma.
@@ -656,7 +658,7 @@ def dipoleModelFunctor(x, flux, xcenPos, ycenPos, xcenNeg, ycenNeg, fluxNeg=None
         if np.sum(np.sqrt((np.array(cenPos) - fpCentroid)**2.)) > maxSep:
             cenPos = fpCentroid
         if np.sum(np.sqrt((np.array(cenNeg) - fpCentroid)**2.)) > maxSep:
-            cenPos = fpCentroid
+            cenNeg = fpCentroid
 
         # parameter hints/constraints: https://lmfit.github.io/lmfit-py/model.html#model-param-hints-section
         # might make sense to not use bounds -- see http://lmfit.github.io/lmfit-py/bounds.html
@@ -1130,16 +1132,14 @@ def measureDipoles(self, measRecord, exposure, posExp=None, negExp=None):
         testing (@see `tests/testDipoleFitter.py`)
         """
         result = None
-        pks = measRecord.getFootprint().getPeaks()
+        peakValues = [peak.getPeakValue() for peak in measRecord.getFootprint().getPeaks()]
 
         # Check if the footprint consists of a putative dipole - else don't fit it.
         if (
-                # One peak in the footprint (not a dipole)
-                ((nPeaks := len(pks)) <= 1)
-                # Peaks are the same sign (not a dipole); peaks are ordered
-                # from highest to lowest.
-                or (nPeaks > 1 and (np.sign(pks[0].getPeakValue())
-                    == np.sign(pks[-1].getPeakValue())))
+                # One peak in the footprint (not a dipole).
+                len(peakValues) <= 1
+                # Peaks are the same sign (not a dipole).
+                or np.sign(max(peakValues)) == np.sign(min(peakValues))
         ):
             if not self.config.fitAllDiaSources:
                 # Non-dipoles fall back on the centroid slot for positions,

tests/test_dipoleFitter.py

@@ -28,11 +28,27 @@
 import numpy as np
 
 import lsst.utils.tests
+import lsst.afw.detection as afwDetection
 import lsst.afw.table as afwTable
 from lsst.ip.diffim.dipoleFitTask import (DipoleFitAlgorithm, DipoleFitTask)
 import utils
 
 
+def reversePeaks(footprint):
+    """Reverse the order of the peaks of a footprint.
+
+    Parameters
+    ----------
+    footprint : `lsst.afw.detection.Footprint`
+        Footprint to modify in place.
+    """
+    peaks = footprint.getPeaks()
+    reordered = afwDetection.PeakCatalog(peaks.table)
+    reordered.extend(list(peaks)[::-1])
+    peaks.clear()
+    footprint.setPeakCatalog(reordered)
+
+
 class DipoleTestImage:
     """Create a test dipole image and store the parameters used to make it,
     for comparison with the fitted results.
@@ -92,18 +108,44 @@ def testDipoleAlgorithm(self):
         Test that the resulting fluxes/centroids are very close to the
         input values for both dipoles in the image.
         """
-        # Display (plot) the output dipole thumbnails with matplotlib.
-        display = False
-        # Be verbose during fitting, including the lmfit internal details.
-        verbose = False
-
         dipoleTestImage = DipoleTestImage()
         catalog = dipoleTestImage.testImage.detectDipoleSources(minBinSize=32)
 
         for s in catalog:
             fp = s.getFootprint()
             self.assertTrue(len(fp.getPeaks()) == 2)
 
+        self._checkAlgorithmOutput(dipoleTestImage, catalog)
+
+    def testDipoleAlgorithmNegativePeakFirst(self):
+        """Test that the fit is the same when the negative peak of each dipole
+        comes before the positive peak.
+        """
+        dipoleTestImage = DipoleTestImage()
+        catalog = dipoleTestImage.testImage.detectDipoleSources(minBinSize=32)
+
+        for s in catalog:
+            reversePeaks(s.getFootprint())
+            self.assertLess(s.getFootprint().getPeaks()[0].getPeakValue(), 0)
+
+        self._checkAlgorithmOutput(dipoleTestImage, catalog)
+
+    def _checkAlgorithmOutput(self, dipoleTestImage, catalog):
+        """Fit each dipole in ``catalog`` and compare the resulting
+        fluxes/centroids with the input values.
+
+        Parameters
+        ----------
+        dipoleTestImage : `DipoleTestImage`
+            Image that was generated for the test.
+        catalog : `lsst.afw.table.SourceCatalog`
+            Sources detected on the test image.
+        """
+        # Display (plot) the output dipole thumbnails with matplotlib.
+        display = False
+        # Be verbose during fitting, including the lmfit internal details.
+        verbose = False
+
         rtol = dipoleTestImage.rtol
         offsets = dipoleTestImage.offsets
         testImage = dipoleTestImage.testImage
@@ -120,11 +162,13 @@ def testDipoleAlgorithm(self):
             self.assertFloatsAlmostEqual(result.negCentroid.x, dipoleTestImage.xc[i] - offsets[i], rtol=rtol)
             self.assertFloatsAlmostEqual(result.negCentroid.y, dipoleTestImage.yc[i] - offsets[i], rtol=rtol)
 
-    def _runDetection(self, dipoleTestImage, maxFootprintArea=None):
+    def _runDetection(self, dipoleTestImage, maxFootprintArea=None, editFootprint=None):
         """Run 'diaSource' detection on the diffim, including merging of
         positive and negative sources.
 
         Then run DipoleFitTask on the image and return the resulting catalog.
+        If ``editFootprint`` is given, it is called with the footprint of
+        each source before the sources are measured.
         """
         # Create the various tasks and schema -- avoid code reuse.
         testImage = dipoleTestImage.testImage
@@ -143,6 +187,9 @@ def _runDetection(self, dipoleTestImage, maxFootprintArea=None):
         fpSet.merge(detectResult.negative, 2, 2, False)
         sources = afwTable.SourceCatalog(table)
         fpSet.makeSources(sources)
+        if editFootprint is not None:
+            for source in sources:
+                editFootprint(source.getFootprint())
 
         measureTask.run(sources, testImage.diffim, testImage.posImage, testImage.negImage)
         return sources
@@ -253,6 +300,25 @@ def testDipoleTask(self):
         sources = self._runDetection(dipoleTestImage)
         self._checkTaskOutput(dipoleTestImage, sources)
 
+    def testDipoleTaskExtraPeak(self):
+        """Test that a dipole is fit when a faint positive peak comes after
+        the peaks of its two lobes.
+        """
+        def addFaintPeak(footprint):
+            peak = footprint.getPeaks()[0]
+            footprint.addPeak(peak.getFx() + 1, peak.getFy(), 0.01*peak.getPeakValue())
+
+        dipoleTestImage = DipoleTestImage()
+        sources = self._runDetection(dipoleTestImage, editFootprint=addFaintPeak)
+
+        for i, record in enumerate(sources):
+            peakValues = [peak.getPeakValue() for peak in record.getFootprint().getPeaks()]
+            self.assertEqual(list(np.sign(peakValues)), [1, -1, 1])
+            self.assertFloatsAlmostEqual((record['ip_diffim_DipoleFit_pos_instFlux']
+                                          + abs(record['ip_diffim_DipoleFit_neg_instFlux']))/2.,
+                                         dipoleTestImage.flux[i], rtol=dipoleTestImage.rtol)
+            self.assertTrue(record['ip_diffim_DipoleFit_classification'])
+
     def testDipoleTaskNoPosImage(self):
         """Test the dipole fitting singleFramePlugin in the case where no
         `posImage` is provided. It should be the same as above because