Skip to content

DM-53312: Robustify the calibrateImage background estimation #1187

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Nov 20, 2025
Merged

DM-53312: Robustify the calibrateImage background estimation #1187

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8469eb4
Remove commented out code from testing
laurenam Nov 18, 2025
cd08d84
Move up setting of diffraction spike mask
laurenam Nov 18, 2025
e7279b0
Robustify the calibrateImage background estimation
laurenam Nov 18, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 92 additions & 14 deletions python/lsst/pipe/tasks/calibrateImage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -324,6 +324,24 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali
target=lsst.meas.algorithms.SubtractBackgroundTask,
doc="Task to perform final background subtraction, just before photoCal.",
)
star_background_min_footprints = pexConfig.Field(
dtype=int,
default=3,
doc="Minimum number of footprints in the detection mask for star_background measurement. "
"This number will get adjusted to the fraction config.star_background_peak_fraction "
"of the detected peaks if that number is larger. If the number of footprints is less "
"than the minimum, the detection threshold is iteratively increased until the "
"threshold is met.",
)
star_background_peak_fraction = pexConfig.Field(
dtype=float,
default=0.01,
doc="The minimum number of footprints in the detection mask for star_background measuremen "
"gets set to the maximum of this fraction of the detected peaks and the value set in "
"config.star_background_min_footprints. If the number of footprints is less than the "
"current minimum set, the detection threshold is iteratively increased until the "
"threshold is met.",
)
star_detection = pexConfig.ConfigurableField(
target=lsst.meas.algorithms.SourceDetectionTask,
doc="Task to detect stars to return in the output catalog."
Expand Down Expand Up @@ -997,6 +1015,11 @@ def run(
astrometry_meta)
result.psf_stars = result.psf_stars_footprints.asAstropy()

# Add diffraction spike mask here so it can be used (i.e. avoided)
# in the adaptive background estimation.
if self.config.doMaskDiffractionSpikes:
self.diffractionSpikeMask.run(result.exposure)

if self.config.do_adaptive_threshold_detection:
self._remeasure_star_background(
result,
Expand Down Expand Up @@ -1049,8 +1072,6 @@ def run(
if "photometry_matches" in self.config.optional_outputs:
result.photometry_matches = lsst.meas.astrom.denormalizeMatches(photometry_matches,
photometry_meta)
if self.config.doMaskDiffractionSpikes:
self.diffractionSpikeMask.run(result.exposure)
if "mask" in self.config.optional_outputs:
result.mask = result.exposure.mask.clone()
except pipeBase.AlgorithmError:
Expand Down Expand Up @@ -1766,8 +1787,8 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
# using an adaptive threshold detection iteration to ensure a
# "Goldilocks Zone" for the fraction of detected pixels.
backgroundOrig = result.background.clone()
median_background = np.median(backgroundOrig.getImage().array)
self.log.warning("Original median_background = %.2f", median_background)
median_background_orig = np.median(backgroundOrig.getImage().array)
self.log.info("Original median_background = %.2f", median_background_orig)
Copy link
Contributor

@erykoff erykoff Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thank you for converting to an info on this one! Also ... have we checked at how long this operation takes? In my experience background.getImage() is surprisingly slow. (Not to solve today; but something to consider for the future).

Copy link
Contributor Author

@laurenam laurenam Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok. On the TODO list...

result.exposure.image.array += result.background.getImage().array
result.background = afwMath.BackgroundList()

Expand Down Expand Up @@ -1819,27 +1840,70 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
"highest_detected_fraction_per_amp = %.3f",
n_above_max_per_amp, highest_detected_fraction_per_amp)

bgIgnoreMasksToAdd = ["SAT", "SUSPECT", "SPIKE"]
detected_fraction = 1.0
nFootprintTemp = 1e12
starBackgroundDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()
for maskName in bgIgnoreMasksToAdd:
if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
and maskName not in starBackgroundDetectionConfig.background.ignoredPixelMask):
starBackgroundDetectionConfig.background.ignoredPixelMask += [maskName]
starBackgroundDetectionConfig.doTempLocalBackground = False
starBackgroundDetectionConfig.nSigmaToGrow = 70.0
starBackgroundDetectionConfig.reEstimateBackground = False
starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0
starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background)
starBackgroundDetectionConfig.thresholdType = "pixel_stdev" # "stdev"
starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background_orig)
starBackgroundDetectionConfig.thresholdType = "pixel_stdev"

n_above_max_per_amp = -99
highest_detected_fraction_per_amp = float("nan")
doCheckPerAmpDetFraction = True

minFootprints = self.config.star_background_min_footprints
maxIter = 40
for nIter in range(maxIter):
currentThresh = starBackgroundDetectionConfig.thresholdValue
if detected_fraction > maxDetFracForFinalBg:
nZeroEncountered = 0
if nFootprintTemp == 0:
zeroFactor = min(0.98, 0.9 + 0.01*nZeroEncountered)
starBackgroundDetectionConfig.thresholdValue = zeroFactor*currentThresh
self.log.warning("No footprints detected. Decreasing threshold to %.2f and rerunning.",
starBackgroundDetectionConfig.thresholdValue)
starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
config=starBackgroundDetectionConfig)
table = afwTable.SourceTable.make(self.initial_stars_schema.schema)
tempDetections = starBackgroundDetectionTask.run(
table=table, exposure=result.exposure, clearMask=True)
nFootprintTemp = len(tempDetections.sources)
minFootprints = max(self.config.star_background_min_footprints,
int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))
minFootprints = min(200, minFootprints)
nZeroEncountered += 1
if nFootprintTemp >= minFootprints:
detected_fraction = self._compute_mask_fraction(result.exposure.mask,
detected_mask_planes,
bad_mask_planes)
self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "
"DETECTED_NEGATIVE in star_background_detection = %.3f "
"(max is %.3f; min is %.3f) nFootprint = %d (current min is %d)",
nIter, starBackgroundDetectionConfig.thresholdValue,
detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,
nFootprintTemp, minFootprints)
break
else:
# Still not enough footprints, so make sure this loop is
# entered again.
if nFootprintTemp > 0 and nFootprintTemp < minFootprints:
nFootprintTemp = 0
continue
if detected_fraction > maxDetFracForFinalBg or nFootprintTemp <= minFootprints:
starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh
if nFootprintTemp < 3 and detected_fraction > 0.9*maxDetFracForFinalBg:
starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh
if nFootprintTemp < minFootprints and detected_fraction > 0.9*maxDetFracForFinalBg:
if nFootprintTemp == 1:
starBackgroundDetectionConfig.thresholdValue = 1.4*currentThresh
Copy link
Contributor

@erykoff erykoff Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where does 1.4 (and 1.2) come from?

Copy link
Contributor Author

@laurenam laurenam Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Trial and error and a trade-off between "get there faster" and "don't overshoot". The detection plane can be disturbingly sensitive to small tweaks in the threshold, but when there's only a single footprint, a bigger increase is usually needed.

Copy link
Contributor

@erykoff erykoff Nov 19, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't love that but ... sure at the moment...

else:
starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh

if n_above_max_per_amp > 1:
starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh
if detected_fraction < minDetFracForFinalBg:
Expand All @@ -1851,13 +1915,17 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
table=table, exposure=result.exposure, clearMask=True)
result.exposure.mask |= dilatedMask
nFootprintTemp = len(tempDetections.sources)
minFootprints = max(self.config.star_background_min_footprints,
int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))
minFootprints = min(200, minFootprints)
detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,
bad_mask_planes)
self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "
"DETECTED_NEGATIVE in star_background_detection = %.3f "
"(max is %.3f; min is %.3f)",
"(max is %.3f; min is %.3f) nFooprint = %d (current min is %d)",
Copy link
Contributor

@erykoff erykoff Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nFooprint -> nFootprint.

nIter, starBackgroundDetectionConfig.thresholdValue,
detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg)
detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,
nFootprintTemp, minFootprints)

n_amp = len(result.exposure.detector.getAmplifiers())
if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg:
Expand All @@ -1870,7 +1938,8 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non

if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg
and n_above_max_per_amp < int(0.75*n_amp)
and no_zero_det_amps):
and no_zero_det_amps
and nFootprintTemp >= minFootprints):
if (n_above_max_per_amp < max(1, int(0.15*n_amp))
or detected_fraction < 0.85*maxDetFracForFinalBg):
break
Expand All @@ -1879,6 +1948,8 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh
self.log.warning("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp))

detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,
bad_mask_planes)
self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f "
"(highest per amp section = %.5f)",
detected_fraction, highest_detected_fraction_per_amp)
Expand All @@ -1895,6 +1966,9 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
).background
result.background.append(star_background[0])

median_background = np.median(result.background.getImage().array)
self.log.info("New initial median_background = %.2f", median_background)

# Perform one more round of background subtraction that is just an
# overall pedestal (order = 0). This is intended to account for
# any potential gross oversubtraction imposed by the higher-order
Expand Down Expand Up @@ -1923,6 +1997,10 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
detected_fraction_orig, detected_fraction_dilated)

pedestalBackgroundConfig = lsst.meas.algorithms.SubtractBackgroundConfig()
for maskName in bgIgnoreMasksToAdd:
if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
and maskName not in pedestalBackgroundConfig.ignoredPixelMask):
pedestalBackgroundConfig.ignoredPixelMask += [maskName]
pedestalBackgroundConfig.statisticsProperty = "MEDIAN"
pedestalBackgroundConfig.approxOrderX = 0
pedestalBackgroundConfig.binSize = 64
Expand All @@ -1932,11 +2010,11 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
background=result.background,
backgroundToPhotometricRatio=background_to_photometric_ratio,
).background
# Isolate the final pedestal background to log the computed value
# Isolate the final pedestal background to log the computed value.
pedestalBackground = afwMath.BackgroundList()
pedestalBackground.append(pedestalBackgroundList[1])
pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0]
self.log.warning("Subtracted pedestalBackgroundLevel = %.4f", pedestalBackgroundLevel)
self.log.info("Subtracted pedestalBackgroundLevel = %.4f", pedestalBackgroundLevel)
Copy link
Contributor

@erykoff erykoff Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So ... I'm confused how this even works now. Does this change the value inside the background list?

        pedestalBackground = afwMath.BackgroundList()
        pedestalBackground.append(pedestalBackgroundList[1])
        pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0]

And then ... it's already been subtracted by the task ... so this won't work.

Copy link
Contributor Author

@laurenam laurenam Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're so right (and thanks for catching this!). Going with dropping the failsafes for now (will follow-up after the next DRP run).


# Clear detected mask plane before final round of detection
mask = result.exposure.mask
Expand Down