Coverage for python/lsst/ip/diffim/computeSpatiallySampledMetrics.py: 91%
192 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-14 19:10 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import dataclasses
24import numpy as np
25import scipy.signal
27import lsst.geom
29import lsst.afw.image as afwImage
30import lsst.afw.table as afwTable
31import lsst.pipe.base as pipeBase
32import lsst.pex.config as pexConfig
34from lsst.ip.diffim.utils import getPsfFwhm, angleMean, getKernelCenterDisplacement
35from lsst.meas.algorithms import SkyObjectsTask
36from lsst.pex.exceptions import InvalidParameterError, RangeError
37from lsst.utils.timer import timeMethod
39import lsst.utils
41__all__ = ["SpatiallySampledMetricsConfig", "SpatiallySampledMetricsTask"]
44class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections,
45 dimensions=("instrument", "visit", "detector"),
46 defaultTemplates={"coaddName": "deep",
47 "warpTypeSuffix": "",
48 "fakesType": ""}):
49 science = pipeBase.connectionTypes.Input(
50 doc="Input science exposure.",
51 dimensions=("instrument", "visit", "detector"),
52 storageClass="ExposureF",
53 name="{fakesType}calexp"
54 )
55 template = pipeBase.connectionTypes.Input(
56 doc="Warped and not PSF-matched template used to create the difference image.",
57 dimensions=("instrument", "visit", "detector"),
58 storageClass="ExposureF",
59 name="{fakesType}{coaddName}Diff_templateExp",
60 )
61 difference = pipeBase.connectionTypes.Input(
62 doc="Difference image with detection mask plane filled in.",
63 dimensions=("instrument", "visit", "detector"),
64 storageClass="ExposureF",
65 name="{fakesType}{coaddName}Diff_differenceExp",
66 )
67 diaSources = pipeBase.connectionTypes.Input(
68 doc="Filtered diaSources on the difference image.",
69 dimensions=("instrument", "visit", "detector"),
70 storageClass="ArrowAstropy",
71 name="{fakesType}dia_source_detector",
72 )
73 psfMatchingKernel = pipeBase.connectionTypes.Input(
74 doc="Kernel used to PSF match the science and template images.",
75 dimensions=("instrument", "visit", "detector"),
76 storageClass="MatchingKernel",
77 name="{fakesType}{coaddName}Diff_psfMatchKernel",
78 )
79 spatiallySampledMetrics = pipeBase.connectionTypes.Output(
80 doc="Summary metrics computed at randomized locations.",
81 dimensions=("instrument", "visit", "detector"),
82 storageClass="ArrowAstropy",
83 name="{fakesType}{coaddName}Diff_spatiallySampledMetrics",
84 )
86 def __init__(self, *, config=None):
87 super().__init__(config=config)
88 if config.image_type == "future":
89 self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel,
90 storageClass="ConvolutionKernel")
93class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig,
94 pipelineConnections=SpatiallySampledMetricsConnections):
95 """Config for SpatiallySampledMetricsTask
96 """
97 metricsMaskPlanes = lsst.pex.config.ListField(
98 dtype=str,
99 doc="List of mask planes to include in metrics. "
100 "A mask plane that is not registered reports 0 pixels set.",
101 default=('BAD', 'CLIPPED', 'CR', 'DETECTED', 'DETECTED_NEGATIVE', 'EDGE',
102 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED',
103 'NO_DATA', 'REJECTED', 'SAT', 'SAT_TEMPLATE', 'SENSOR_EDGE', 'STREAK', 'SUSPECT',
104 'UNMASKEDNAN',
105 ),
106 )
107 metricSources = pexConfig.ConfigurableField(
108 target=SkyObjectsTask,
109 doc="Generate QA metric sources",
110 )
111 image_type = pexConfig.ChoiceField[str](
112 doc="Type of the input PSF matching kernel. The input images are read"
113 " as lsst.afw.image.ExposureF either way.",
114 allowed={
115 "legacy": "Read lsst.afw.math.Kernel.",
116 "future": "Read lsst.images.convolution_kernels.ConvolutionKernel.",
117 },
118 optional=False,
119 default="legacy",
120 )
122 def setDefaults(self):
123 self.metricSources.avoidMask = ["NO_DATA", "EDGE"]
126class SpatiallySampledMetricsTask(lsst.pipe.base.PipelineTask):
127 """Detect and measure sources on a difference image.
128 """
129 ConfigClass = SpatiallySampledMetricsConfig
130 _DefaultName = "spatiallySampledMetrics"
132 def __init__(self, **kwargs):
133 super().__init__(**kwargs)
135 self.makeSubtask("metricSources")
136 self.schema = afwTable.SourceTable.makeMinimalSchema()
137 self.schema.addField(
138 "x", "F",
139 "X location of the metric evaluation.",
140 units="pixel")
141 self.schema.addField(
142 "y", "F",
143 "Y location of the metric evaluation.",
144 units="pixel")
145 self.metricSources.skySourceKey = self.schema.addField("sky_source", type="Flag",
146 doc="Metric evaluation objects.")
147 self.schema.addField(
148 "source_density", "F",
149 "Density of diaSources at location.",
150 units="count/degree^2")
151 self.schema.addField(
152 "dipole_density", "F",
153 "Density of dipoles at location.",
154 units="count/degree^2")
155 self.schema.addField(
156 "dipole_direction", "F",
157 "Mean dipole orientation.",
158 units="radian")
159 self.schema.addField(
160 "dipole_separation", "F",
161 "Mean dipole separation.",
162 units="pixel")
163 self.schema.addField(
164 "template_value", "F",
165 "Median of template at location.",
166 units="nJy")
167 self.schema.addField(
168 "template_variance", "F",
169 "Median of template variance at location.",
170 units="nJy^2")
171 self.schema.addField(
172 "science_value", "F",
173 "Median of science at location.",
174 units="nJy")
175 self.schema.addField(
176 "science_variance", "F",
177 "Median of science variance at location.",
178 units="nJy^2")
179 self.schema.addField(
180 "diffim_value", "F",
181 "Median of diffim at location.",
182 units="nJy")
183 self.schema.addField(
184 "diffim_variance", "F",
185 "Median of diffim variance at location.",
186 units="nJy^2")
187 self.schema.addField(
188 "diffim_chi2PerPix", "F",
189 "Robust normalized noise of diffim at location:"
190 " (1.4826*MAD(image))^2 / median(variance), evaluated on background"
191 " pixels (DETECTED, DETECTED_NEGATIVE, BAD, SAT, EDGE, NO_DATA excluded)."
192 " Expected ~1.0 for a well-decorrelated diffim; values >>1 indicate"
193 " residual structure, <<1 indicates over-estimated variance.")
194 self.schema.addField(
195 "science_psfSize", "F",
196 "Width of the science image PSF at location.",
197 units="pixel")
198 self.schema.addField(
199 "template_psfSize", "F",
200 "Width of the template image PSF at location.",
201 units="pixel")
202 for maskPlane in self.config.metricsMaskPlanes:
203 self.schema.addField(
204 "%s_mask_fraction"%maskPlane.lower(), "F",
205 "Fraction of pixels with %s mask"%maskPlane
206 )
207 self.schema.addField(
208 "psfMatchingKernel_sum", "F",
209 "PSF matching kernel sum at location.")
210 self.schema.addField(
211 "psfMatchingKernel_dx", "F",
212 "PSF matching kernel centroid offset in x at location.",
213 units="pixel")
214 self.schema.addField(
215 "psfMatchingKernel_dy", "F",
216 "PSF matching kernel centroid offset in y at location.",
217 units="pixel")
218 self.schema.addField(
219 "psfMatchingKernel_length", "F",
220 "PSF matching kernel centroid offset module.",
221 units="arcsecond")
222 self.schema.addField(
223 "psfMatchingKernel_position_angle", "F",
224 "PSF matching kernel centroid offset position angle.",
225 units="radian")
226 self.schema.addField(
227 "psfMatchingKernel_direction", "F",
228 "PSF matching kernel centroid offset direction in detector plane.",
229 units="radian")
230 self.schema.addField(
231 "psfMatchingKernel_residualNorm", "F",
232 "Shape-only PSF match residual at location:"
233 "L2 norm of (K-convolved template PSF - science PSF),"
234 " relative to the science PSF L2 norm. Larger values indicate worse"
235 " PSF matching. Assumes the kernel was solved to convolve the template.")
237 def runQuantum(self, butlerQC, inputRefs, outputRefs):
238 inputs = butlerQC.get(inputRefs)
239 if self.config.image_type == "future":
240 inputs["psfMatchingKernel"] = inputs["psfMatchingKernel"].to_legacy()
241 outputs = self.run(**inputs)
242 butlerQC.put(outputs, outputRefs)
244 @timeMethod
245 def run(self, science, template, difference, diaSources, psfMatchingKernel):
246 """Calculate difference image metrics on specific locations across the images
248 Parameters
249 ----------
250 science : `lsst.afw.image.ExposureF`
251 Science exposure that the template was subtracted from.
252 template : `lsst.afw.image.ExposureF`
253 Warped and non PSF-matched template that was used produce
254 the difference image.
255 difference : `lsst.afw.image.ExposureF`
256 Result of subtracting template from the science image.
257 diaSources : `lsst.afw.table.SourceCatalog`
258 The catalog of detected sources.
259 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
260 The PSF matching kernel of the subtraction to evaluate.
262 Returns
263 -------
264 results : `lsst.pipe.base.Struct`
265 ``spatiallySampledMetrics`` : `astropy.table.Table`
266 Image quality metrics spatially sampled locations.
267 """
269 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
271 spatiallySampledMetrics = afwTable.SourceCatalog(self.schema)
272 spatiallySampledMetrics.getTable().setIdFactory(idFactory)
274 self.metricSources.run(mask=science.mask, seed=difference.info.id, catalog=spatiallySampledMetrics)
276 # Mask plane bits are looked up in the global afw dict, not the
277 # image's own. A plane that is not registered in this process has no
278 # pixels set in any image.
279 metricsMaskPlanes = {}
280 for maskPlane in self.config.metricsMaskPlanes:
281 try:
282 metricsMaskPlanes[maskPlane] = afwImage.Mask.getPlaneBitMask(maskPlane)
283 except InvalidParameterError:
284 self.log.info("Mask plane %s is not registered; its metrics will be zero.", maskPlane)
285 metricsMaskPlanes[maskPlane] = 0
287 for src in spatiallySampledMetrics:
288 self._evaluateLocalMetric(src, science, template, difference, diaSources,
289 metricsMaskPlanes=metricsMaskPlanes,
290 psfMatchingKernel=psfMatchingKernel)
291 spatiallySampledMetrics = spatiallySampledMetrics.copy(deep=True).asAstropy()
292 return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics)
294 def _evaluateLocalMetric(self, src, science, template, difference, diaSources,
295 metricsMaskPlanes, psfMatchingKernel):
296 """Calculate image quality metrics at spatially sampled locations.
298 Parameters
299 ----------
300 src : `lsst.afw.table.SourceRecord`
301 The source record to be updated with metric calculations.
302 diaSources : `lsst.afw.table.SourceCatalog`
303 The catalog of detected sources.
304 science : `lsst.afw.image.Exposure`
305 The science image.
306 difference : `lsst.afw.image.Exposure`
307 Result of subtracting template from the science image.
308 metricsMaskPlanes : `dict` [`str`, `int`]
309 Mask planes to calculate metrics from, mapped to their bitmask.
310 The bitmask is 0 for a plane that is not registered.
311 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
312 The PSF matching kernel of the subtraction to evaluate.
313 """
314 bbox = src.getFootprint().getBBox()
315 pix = bbox.getCenter()
316 src.set('science_psfSize', getPsfFwhm(science.psf, position=pix))
317 try:
318 src.set('template_psfSize', getPsfFwhm(template.psf, position=pix))
319 except (InvalidParameterError, RangeError):
320 src.set('template_psfSize', np.nan)
322 metricRegionSize = 100
323 bbox.grow(metricRegionSize)
324 bbox = bbox.clippedTo(science.getBBox())
325 nPix = bbox.getArea()
326 pixScale = science.wcs.getPixelScale(bbox.getCenter())
327 area = nPix*pixScale.asDegrees()**2
328 peak = src.getFootprint().getPeaks()[0]
329 src.set('x', peak['i_x'])
330 src.set('y', peak['i_y'])
331 src.setCoord(science.wcs.pixelToSky(peak['i_x'], peak['i_y']))
332 selectSources = diaSources[bbox.contains(diaSources['x'], diaSources['y'])]
333 sourceDensity = len(selectSources)/area
334 dipoleSources = selectSources[selectSources["isDipole"]]
335 dipoleDensity = len(dipoleSources)/area
337 if dipoleSources: 337 ↛ 338line 337 didn't jump to line 338 because the condition on line 337 was never true
338 meanDipoleOrientation = angleMean(dipoleSources["dipoleAngle"])
339 src.set('dipole_direction', meanDipoleOrientation)
340 meanDipoleSeparation = np.mean(dipoleSources["dipoleLength"])
341 src.set('dipole_separation', meanDipoleSeparation)
343 templateVal = np.median(template[bbox].image.array)
344 templateVar = np.median(template[bbox].variance.array)
345 scienceVal = np.median(science[bbox].image.array)
346 scienceVar = np.median(science[bbox].variance.array)
347 diffimVal = np.median(difference[bbox].image.array)
348 diffimVar = np.median(difference[bbox].variance.array)
349 src.set('source_density', sourceDensity)
350 src.set('dipole_density', dipoleDensity)
351 src.set('template_value', templateVal)
352 src.set('template_variance', templateVar)
353 src.set('science_value', scienceVal)
354 src.set('science_variance', scienceVar)
355 src.set('diffim_value', diffimVal)
356 src.set('diffim_variance', diffimVar)
357 src.set('diffim_chi2PerPix', self._diffimChi2PerPix(difference[bbox]))
358 differenceMask = difference.mask[bbox].array
359 for maskPlane, bitMask in metricsMaskPlanes.items():
360 src.set("%s_mask_fraction"%maskPlane.lower(),
361 np.count_nonzero(differenceMask & bitMask)/differenceMask.size
362 )
364 krnlSum, dx, dy, direction, length = getKernelCenterDisplacement(
365 psfMatchingKernel, src.get('x'), src.get('y'))
367 point1 = lsst.geom.SpherePoint(
368 src.get('coord_ra'), src.get('coord_dec'),
369 lsst.geom.radians)
370 point2 = science.wcs.pixelToSky(src.get('x') + dx, src.get('y') + dy)
371 bearing = point1.bearingTo(point2)
372 pa_ref_angle = lsst.geom.Angle(np.pi/2, lsst.geom.radians)
373 pa = pa_ref_angle - bearing
374 # Wrap around to get Delta_RA from -pi to +pi
375 pa = pa.wrapCtr()
376 position_angle = pa.asRadians()
378 src.set('psfMatchingKernel_sum', krnlSum)
379 src.set('psfMatchingKernel_dx', dx)
380 src.set('psfMatchingKernel_dy', dy)
381 src.set('psfMatchingKernel_length', length*pixScale.asArcseconds())
382 src.set('psfMatchingKernel_position_angle', position_angle) # in E of N position angle
383 src.set('psfMatchingKernel_direction', direction) # direction offset in detector
385 src.set('psfMatchingKernel_residualNorm',
386 self._psfMatchResidualNorm(psfMatchingKernel, science.psf, template.psf,
387 src.get('x'), src.get('y')))
389 def _diffimChi2PerPix(self, difference):
390 """Robust normalized noise of the difference image.
392 Computes ``(1.4826 * MAD(image))^2 / median(variance)`` on
393 background-only pixels. The 1.4826 factor rescales the
394 Median Absolute Deviation to a Gaussian-sigma estimate, so the
395 ratio has expectation 1.0 on pure noise.
397 Returns NaN if no usable pixels remain or the median variance
398 is non-positive.
399 """
400 image = difference.image.array
401 variance = difference.variance.array
402 mask = difference.mask
403 excludePlanes = [p for p in ("DETECTED", "DETECTED_NEGATIVE", "BAD",
404 "SAT", "EDGE", "NO_DATA")
405 if p in mask.getMaskPlaneDict()]
406 if excludePlanes: 406 ↛ 410line 406 didn't jump to line 410 because the condition on line 406 was always true
407 badBits = mask.getPlaneBitMask(excludePlanes)
408 good = (mask.array & badBits) == 0
409 else:
410 good = np.ones(image.shape, dtype=bool)
411 good &= np.isfinite(image) & np.isfinite(variance) & (variance > 0)
412 if not np.any(good): 412 ↛ 413line 412 didn't jump to line 413 because the condition on line 412 was never true
413 return np.nan
414 varMed = np.median(variance[good])
415 if not np.isfinite(varMed) or varMed <= 0: 415 ↛ 416line 415 didn't jump to line 416 because the condition on line 415 was never true
416 return np.nan
417 imgGood = image[good]
418 mad = np.median(np.abs(imgGood - np.median(imgGood)))
419 return float((1.4826*mad)**2/varMed)
421 def _psfMatchResidualNorm(self, kernel, sciencePsf, templatePsf, x, y):
422 """Relative L2 norm of the PSF-match residual at (x, y).
424 Convolves the template PSF with the matching kernel evaluated at
425 (x, y) and compares to the science PSF at the same position, both
426 renormalized to unit sum so the result captures shape mismatch only.
427 The kernel is assumed to convolve the template.
429 Returns NaN if either PSF cannot be evaluated at the position or
430 the resulting images cannot be normalized.
431 """
432 point = lsst.geom.Point2D(x, y)
433 try:
434 psfSci = sciencePsf.computeKernelImage(point).array
435 psfTmp = templatePsf.computeKernelImage(point).array
436 except (InvalidParameterError, RangeError):
437 return np.nan
439 kImage = afwImage.ImageD(kernel.getDimensions())
440 kernel.computeImage(kImage, doNormalize=True, x=x, y=y)
441 matched = scipy.signal.fftconvolve(psfTmp, kImage.array, mode='same')
443 matchedSum = matched.sum()
444 psfSciSum = psfSci.sum()
445 if matchedSum <= 0 or psfSciSum <= 0: 445 ↛ 446line 445 didn't jump to line 446 because the condition on line 445 was never true
446 return np.nan
447 matched = matched/matchedSum
448 psfSci = psfSci/psfSciSum
450 # PSF stamps may differ in size between the two exposures; crop
451 # both to a common centered region before differencing.
452 h = min(matched.shape[0], psfSci.shape[0])
453 w = min(matched.shape[1], psfSci.shape[1])
455 def _crop(a):
456 sy = (a.shape[0] - h)//2
457 sx = (a.shape[1] - w)//2
458 return a[sy:sy + h, sx:sx + w]
459 matched = _crop(matched)
460 psfSci = _crop(psfSci)
462 sciNorm = np.sqrt(np.sum(psfSci**2))
463 if sciNorm == 0: 463 ↛ 464line 463 didn't jump to line 464 because the condition on line 463 was never true
464 return np.nan
465 return float(np.sqrt(np.sum((matched - psfSci)**2))/sciNorm)