From 8cec98e0048993687a56840ab06297ff1ae7d496 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Tue, 10 Nov 2026 14:33:11 -0800 Subject: [PATCH] Fix bug when a dipole negative lobe is brighter --- python/lsst/ip/diffim/dipoleFitTask.py | 22 ++++---- tests/test_dipoleFitter.py | 78 ++++++++++++++++++++++++-- 2 files changed, 83 insertions(+), 17 deletions(-) diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index e525c3e47..91f0976cb 100644 --- a/python/lsst/ip/diffim/dipoleFitTask.py +++ b/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, diff --git a/tests/test_dipoleFitter.py b/tests/test_dipoleFitter.py index bc5c5a74d..06b6bb640 100644 --- a/tests/test_dipoleFitter.py +++ b/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,11 +108,6 @@ 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) @@ -104,6 +115,37 @@ def testDipoleAlgorithm(self): 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