Skip to content

DM-53304: Pass attribute result struct into _compute_psf #1201

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 5 commits into from
Dec 4, 2025
Merged

DM-53304: Pass attribute result struct into _compute_psf #1201

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
2e8c4f1
Fix typos
leeskelvin Dec 2, 2025
88150cc
Pass attribute struct into _compute_psf
leeskelvin Dec 2, 2025
2e76341
Add missing docstring for adaptive_det_res_struct
leeskelvin Dec 2, 2025
ceb55c8
Fix typos and too-long docstrings
leeskelvin Dec 3, 2025
a98e4ea
Update unit tests for new _compute_psf API
leeskelvin Dec 3, 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
110 changes: 61 additions & 49 deletions python/lsst/pipe/tasks/calibrateImage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,7 @@ class AllCentroidsFlaggedError(pipeBase.AlgorithmError):
"""
def __init__(self, n_sources, psf_shape_ixx, psf_shape_iyy, psf_shape_ixy, psf_size):
msg = (f"All source centroids (out of {n_sources}) flagged during PSF fitting. "
"Original image PSF is likely unuseable; best-fit PSF shape parameters: "
"Original image PSF is likely unusable; best-fit PSF shape parameters: "
f"Ixx={psf_shape_ixx}, Iyy={psf_shape_iyy}, Ixy={psf_shape_ixy}, size={psf_size}"
)
super().__init__(msg)
Expand Down Expand Up @@ -284,7 +284,7 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali
)
psf_subtract_background = pexConfig.ConfigurableField(
target=lsst.meas.algorithms.SubtractBackgroundTask,
doc="Task to perform intial background subtraction, before first detection pass.",
doc="Task to perform initial background subtraction, before first detection pass.",
)
psf_detection = pexConfig.ConfigurableField(
target=lsst.meas.algorithms.SourceDetectionTask,
Expand Down Expand Up @@ -337,7 +337,7 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali
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 "
doc="The minimum number of footprints in the detection mask for star_background measurement. "
"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 "
Expand Down Expand Up @@ -494,7 +494,7 @@ class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=Cali
def setDefaults(self):
super().setDefaults()

# Use a very broad PSF here, to throughly reject CRs.
# Use a very broad PSF here, to thoroughly reject CRs.
# TODO investigation: a large initial psf guess may make stars look
# like CRs for very good seeing images.
self.install_simple_psf.fwhm = 4
Expand Down Expand Up @@ -534,7 +534,7 @@ def setDefaults(self):
self.psf_measure_psf.psfDeterminer["psfex"].photometricFluxField = \
"base_CircularApertureFlux_12_0_instFlux"

# No extendeness information available: we need the aperture
# No extendedness information available: we need the aperture
# corrections to determine that.
self.measure_aperture_correction.sourceSelector["science"].doUnresolved = False
self.measure_aperture_correction.sourceSelector["science"].flags.good = ["calib_psf_used"]
Expand Down Expand Up @@ -785,7 +785,7 @@ def __init__(self, initial_stars_schema=None, **kwargs):

# The final catalog will have calibrated flux columns, which we add to
# the init-output schema by calibrating our zero-length catalog with an
# arbitrary dummy PhotoCalib. We also use this schema to initialze
# arbitrary dummy PhotoCalib. We also use this schema to initialize
# the stars catalog in order to ensure it's the same even when we hit
# an error (and write partial outputs) before calibrating the catalog
# - note that calibrateCatalog will happily reuse existing output
Expand Down Expand Up @@ -973,6 +973,7 @@ def run(
result.exposure.detector.getId())

result.background = None
result.background_to_photometric_ratio = None
summary_stat_catalog = None
# Some exposure components are set to initial placeholder objects
# while we try to bootstrap them. If we fail before we fit for them,
Expand All @@ -988,10 +989,9 @@ def run(
illumination_correction,
)

result.psf_stars_footprints, result.background, _, adaptive_det_res_struct = self._compute_psf(
result.exposure,
result.psf_stars_footprints, _, adaptive_det_res_struct = self._compute_psf(
result,
id_generator,
background_to_photometric_ratio=result.background_to_photometric_ratio,
)
have_fit_psf = True

Expand All @@ -1014,7 +1014,7 @@ def run(
self._measure_aperture_correction(result.exposure, result.psf_stars_footprints)
result.psf_stars = result.psf_stars_footprints.asAstropy()
# Run astrometry using PSF candidate stars.
# Update "the psf_stars" source cooordinates with the current wcs.
# Update "the psf_stars" source coordinates with the current wcs.
afwTable.updateSourceCoords(
result.exposure.wcs,
sourceList=result.psf_stars_footprints,
Expand Down Expand Up @@ -1055,7 +1055,7 @@ def run(
self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints,
psfSigma=psfSigma)

# Update the "stars" source cooordinates with the current wcs.
# Update the "stars" source coordinates with the current wcs.
afwTable.updateSourceCoords(
result.exposure.wcs,
sourceList=result.stars_footprints,
Expand Down Expand Up @@ -1203,7 +1203,7 @@ def _apply_illumination_correction(self, exposure, background_flat, illumination

return background_to_photometric_ratio

def _compute_psf(self, exposure, id_generator, background_to_photometric_ratio=None):
def _compute_psf(self, result, id_generator):
"""Find bright sources detected on an exposure and fit a PSF model to
them, repairing likely cosmic rays before detection.

Expand All @@ -1212,22 +1212,35 @@ def _compute_psf(self, exposure, id_generator, background_to_photometric_ratio=N

Parameters
----------
exposure : `lsst.afw.image.Exposure`
Exposure to detect and measure bright stars on.
result : `lsst.pipe.base.Struct`
Result struct that is modified to allow saving of partial outputs
for some failure conditions. Should contain at least the following
attributes:

- exposure : `lsst.afw.image.Exposure`
Exposure to detect and measure bright stars on.
- background : `lsst.afw.math.BackgroundList` | `None`
Background that was fit to the exposure during detection.
- background_to_photometric_ratio : `lsst.afw.image.Image` | `None`
Image to convert photometric-flattened image to
background-flattened image.
id_generator : `lsst.meas.base.IdGenerator`
Object that generates source IDs and provides random seeds.
background_to_photometric_ratio : `lsst.afw.image.Image`, optional
Image to convert photometric-flattened image to
background-flattened image.

Returns
-------
sources : `lsst.afw.table.SourceCatalog`
Catalog of detected bright sources.
background : `lsst.afw.math.BackgroundList`
Background that was fit to the exposure during detection.
cell_set : `lsst.afw.math.SpatialCellSet`
PSF candidates returned by the psf determiner.
adaptive_det_res_struct : `lsst.pipe.base.Struct`
Result struct from the adaptive threshold detection.

Notes
-----
This method modifies the exposure, background and
background_to_photometric_ratio attributes of the result struct
in-place.
"""
def log_psf(msg, addToMetadata=False):
"""Log the parameters of the psf and background, with a prepended
Expand All @@ -1242,28 +1255,28 @@ def log_psf(msg, addToMetadata=False):
Whether to add the final psf sigma value to the task
metadata (the default is False).
"""
position = exposure.psf.getAveragePosition()
sigma = exposure.psf.computeShape(position).getDeterminantRadius()
dimensions = exposure.psf.computeImage(position).getDimensions()
if background is not None:
median_background = np.median(background.getImage().array)
position = result.exposure.psf.getAveragePosition()
sigma = result.exposure.psf.computeShape(position).getDeterminantRadius()
dimensions = result.exposure.psf.computeImage(position).getDimensions()
if result.background is not None:
median_background = np.median(result.background.getImage().array)
else:
median_background = 0.0
self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f",
msg, sigma, dimensions, median_background)
if addToMetadata:
self.metadata["final_psf_sigma"] = sigma

self.log.info("First pass detection with Guassian PSF FWHM=%s pixels",
self.log.info("First pass detection with Gaussian PSF FWHM=%s pixels",
self.config.install_simple_psf.fwhm)
self.install_simple_psf.run(exposure=exposure)
self.install_simple_psf.run(exposure=result.exposure)

background = self.psf_subtract_background.run(
exposure=exposure,
backgroundToPhotometricRatio=background_to_photometric_ratio,
result.background = self.psf_subtract_background.run(
exposure=result.exposure,
backgroundToPhotometricRatio=result.background_to_photometric_ratio,
).background
log_psf("Initial PSF:")
self.psf_repair.run(exposure=exposure, keepCRs=True)
self.psf_repair.run(exposure=result.exposure, keepCRs=True)

table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory())
if not self.config.do_adaptive_threshold_detection:
Expand All @@ -1272,15 +1285,15 @@ def log_psf(msg, addToMetadata=False):
# measurement uses the most accurate background-subtraction.
detections = self.psf_detection.run(
table=table,
exposure=exposure,
background=background,
backgroundToPhotometricRatio=background_to_photometric_ratio,
exposure=result.exposure,
background=result.background,
backgroundToPhotometricRatio=result.background_to_photometric_ratio,
)
else:
initialThreshold = self.config.psf_detection.thresholdValue
initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier
adaptive_det_res_struct = self.psf_adaptive_threshold_detection.run(
table, exposure,
table, result.exposure,
initialThreshold=initialThreshold,
initialThresholdMultiplier=initialThresholdMultiplier,
)
Expand All @@ -1290,8 +1303,8 @@ def log_psf(msg, addToMetadata=False):
self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg
self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks
self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks
self.psf_source_measurement.run(detections.sources, exposure)
psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)
self.psf_source_measurement.run(detections.sources, result.exposure)
psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)

# This 2nd round of PSF fitting has been deemed unnecessary (and
# sometimes even causing harm), so it is being skipped for the
Expand All @@ -1302,7 +1315,7 @@ def log_psf(msg, addToMetadata=False):
# Replace the initial PSF with something simpler for the second
# repair/detect/measure/measure_psf step: this can help it
# converge.
self.install_simple_psf.run(exposure=exposure)
self.install_simple_psf.run(exposure=result.exposure)

log_psf("Rerunning with simple PSF:")
# TODO investigation: Should we only re-run repair here, to use the
Expand All @@ -1312,32 +1325,31 @@ def log_psf(msg, addToMetadata=False):
# for the post-psf_measure_psf step, since we only want to do
# PsfFlux and GaussianFlux *after* we have a PSF? Maybe that's not
# relevant once DM-39203 is merged?
self.psf_repair.run(exposure=exposure, keepCRs=True)
self.psf_repair.run(exposure=result.exposure, keepCRs=True)
# Re-estimate the background during this detection step, so that
# measurement uses the most accurate background-subtraction.
detections = self.psf_detection.run(
table=table,
exposure=exposure,
background=background,
backgroundToPhotometricRatio=background_to_photometric_ratio,
exposure=result.exposure,
background=result.background,
backgroundToPhotometricRatio=result.background_to_photometric_ratio,
)
self.psf_source_measurement.run(detections.sources, exposure)
psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)

self.psf_source_measurement.run(detections.sources, result.exposure)
psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)
self.metadata["simple_psf_positive_footprint_count"] = detections.numPos
self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg
self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks
self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks
log_psf("Final PSF:", addToMetadata=True)

# Final repair with final PSF, removing cosmic rays this time.
self.psf_repair.run(exposure=exposure)
self.psf_repair.run(exposure=result.exposure)
# Final measurement with the CRs removed.
self.psf_source_measurement.run(detections.sources, exposure)
self.psf_source_measurement.run(detections.sources, result.exposure)

# PSF is set on exposure; candidates are returned to use for
# calibration flux normalization and aperture corrections.
return detections.sources, background, psf_result.cellSet, adaptive_det_res_struct
return detections.sources, psf_result.cellSet, adaptive_det_res_struct

def _measure_aperture_correction(self, exposure, bright_sources):
"""Measure and set the ApCorrMap on the Exposure, using
Expand Down Expand Up @@ -1851,7 +1863,7 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
result : `lsst.pipe.base.Struct`
The modified result Struct with the new background subtracted.
"""
# Restore the previously measured backgroud and remeasure it
# Restore the previously measured background and remeasure it
# using an adaptive threshold detection iteration to ensure a
# "Goldilocks Zone" for the fraction of detected pixels.
backgroundOrig = result.background.clone()
Expand Down Expand Up @@ -1992,7 +2004,7 @@ def _remeasure_star_background(self, result, background_to_photometric_ratio=Non
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) nFooprint = %d (current min is %d)",
"(max is %.3f; min is %.3f) nFootprint = %d (current min is %d)",
nIter, starBackgroundDetectionConfig.thresholdValue,
detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,
nFootprintTemp, minFootprints)
Expand Down
Loading