Coverage for python/lsst/ip/diffim/computeSpatiallySampledMetrics.py: 91%
198 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-21 09:17 +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.science = dataclasses.replace(self.science, storageClass="VisitImage")
90 self.template = dataclasses.replace(self.template, storageClass="DifferenceImage")
91 self.difference = dataclasses.replace(self.difference, storageClass="DifferenceImage")
92 self.psfMatchingKernel = dataclasses.replace(self.psfMatchingKernel,
93 storageClass="ConvolutionKernel")
96class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig,
97 pipelineConnections=SpatiallySampledMetricsConnections):
98 """Config for SpatiallySampledMetricsTask
99 """
100 metricsMaskPlanes = lsst.pex.config.ListField(
101 dtype=str,
102 doc="List of mask planes to include in metrics. "
103 "A mask plane that is not registered reports 0 pixels set.",
104 default=('BAD', 'CLIPPED', 'CR', 'DETECTED', 'DETECTED_NEGATIVE', 'EDGE',
105 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED',
106 'NO_DATA', 'REJECTED', 'SAT', 'SAT_TEMPLATE', 'SENSOR_EDGE', 'STREAK', 'SUSPECT',
107 'UNMASKEDNAN',
108 ),
109 )
110 metricSources = pexConfig.ConfigurableField(
111 target=SkyObjectsTask,
112 doc="Generate QA metric sources",
113 )
114 image_type = pexConfig.ChoiceField[str](
115 doc="Type of the input images and PSF matching kernel. run() takes the"
116 " legacy types either way.",
117 allowed={
118 "legacy": "Read lsst.afw.image.ExposureF and lsst.afw.math.Kernel.",
119 "future": "Read lsst.images.VisitImage for the science image,"
120 " lsst.images.DifferenceImage for the template and difference"
121 " image, and lsst.images.convolution_kernels.ConvolutionKernel.",
122 },
123 optional=False,
124 default="legacy",
125 )
127 def setDefaults(self):
128 self.metricSources.avoidMask = ["NO_DATA", "EDGE"]
131class SpatiallySampledMetricsTask(lsst.pipe.base.PipelineTask):
132 """Detect and measure sources on a difference image.
133 """
134 ConfigClass = SpatiallySampledMetricsConfig
135 _DefaultName = "spatiallySampledMetrics"
137 def __init__(self, **kwargs):
138 super().__init__(**kwargs)
140 self.makeSubtask("metricSources")
141 self.schema = afwTable.SourceTable.makeMinimalSchema()
142 self.schema.addField(
143 "x", "F",
144 "X location of the metric evaluation.",
145 units="pixel")
146 self.schema.addField(
147 "y", "F",
148 "Y location of the metric evaluation.",
149 units="pixel")
150 self.metricSources.skySourceKey = self.schema.addField("sky_source", type="Flag",
151 doc="Metric evaluation objects.")
152 self.schema.addField(
153 "source_density", "F",
154 "Density of diaSources at location.",
155 units="count/degree^2")
156 self.schema.addField(
157 "dipole_density", "F",
158 "Density of dipoles at location.",
159 units="count/degree^2")
160 self.schema.addField(
161 "dipole_direction", "F",
162 "Mean dipole orientation.",
163 units="radian")
164 self.schema.addField(
165 "dipole_separation", "F",
166 "Mean dipole separation.",
167 units="pixel")
168 self.schema.addField(
169 "template_value", "F",
170 "Median of template at location.",
171 units="nJy")
172 self.schema.addField(
173 "template_variance", "F",
174 "Median of template variance at location.",
175 units="nJy^2")
176 self.schema.addField(
177 "science_value", "F",
178 "Median of science at location.",
179 units="nJy")
180 self.schema.addField(
181 "science_variance", "F",
182 "Median of science variance at location.",
183 units="nJy^2")
184 self.schema.addField(
185 "diffim_value", "F",
186 "Median of diffim at location.",
187 units="nJy")
188 self.schema.addField(
189 "diffim_variance", "F",
190 "Median of diffim variance at location.",
191 units="nJy^2")
192 self.schema.addField(
193 "diffim_chi2PerPix", "F",
194 "Robust normalized noise of diffim at location:"
195 " (1.4826*MAD(image))^2 / median(variance), evaluated on background"
196 " pixels (DETECTED, DETECTED_NEGATIVE, BAD, SAT, EDGE, NO_DATA excluded)."
197 " Expected ~1.0 for a well-decorrelated diffim; values >>1 indicate"
198 " residual structure, <<1 indicates over-estimated variance.")
199 self.schema.addField(
200 "science_psfSize", "F",
201 "Width of the science image PSF at location.",
202 units="pixel")
203 self.schema.addField(
204 "template_psfSize", "F",
205 "Width of the template image PSF at location.",
206 units="pixel")
207 for maskPlane in self.config.metricsMaskPlanes:
208 self.schema.addField(
209 "%s_mask_fraction"%maskPlane.lower(), "F",
210 "Fraction of pixels with %s mask"%maskPlane
211 )
212 self.schema.addField(
213 "psfMatchingKernel_sum", "F",
214 "PSF matching kernel sum at location.")
215 self.schema.addField(
216 "psfMatchingKernel_dx", "F",
217 "PSF matching kernel centroid offset in x at location.",
218 units="pixel")
219 self.schema.addField(
220 "psfMatchingKernel_dy", "F",
221 "PSF matching kernel centroid offset in y at location.",
222 units="pixel")
223 self.schema.addField(
224 "psfMatchingKernel_length", "F",
225 "PSF matching kernel centroid offset module.",
226 units="arcsecond")
227 self.schema.addField(
228 "psfMatchingKernel_position_angle", "F",
229 "PSF matching kernel centroid offset position angle.",
230 units="radian")
231 self.schema.addField(
232 "psfMatchingKernel_direction", "F",
233 "PSF matching kernel centroid offset direction in detector plane.",
234 units="radian")
235 self.schema.addField(
236 "psfMatchingKernel_residualNorm", "F",
237 "Shape-only PSF match residual at location:"
238 "L2 norm of (K-convolved template PSF - science PSF),"
239 " relative to the science PSF L2 norm. Larger values indicate worse"
240 " PSF matching. Assumes the kernel was solved to convolve the template.")
242 def runQuantum(self, butlerQC, inputRefs, outputRefs):
243 inputs = butlerQC.get(inputRefs)
244 if self.config.image_type == "future":
245 inputs["science"] = inputs["science"].to_legacy()
246 inputs["template"] = inputs["template"].to_legacy()
247 inputs["difference"] = inputs["difference"].to_legacy()
248 inputs["psfMatchingKernel"] = inputs["psfMatchingKernel"].to_legacy()
249 outputs = self.run(**inputs)
250 butlerQC.put(outputs, outputRefs)
252 @timeMethod
253 def run(self, science, template, difference, diaSources, psfMatchingKernel):
254 """Calculate difference image metrics on specific locations across the images
256 Parameters
257 ----------
258 science : `lsst.afw.image.ExposureF`
259 Science exposure that the template was subtracted from.
260 template : `lsst.afw.image.ExposureF`
261 Warped and non PSF-matched template that was used produce
262 the difference image.
263 difference : `lsst.afw.image.ExposureF`
264 Result of subtracting template from the science image.
265 diaSources : `lsst.afw.table.SourceCatalog`
266 The catalog of detected sources.
267 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
268 The PSF matching kernel of the subtraction to evaluate.
270 Returns
271 -------
272 results : `lsst.pipe.base.Struct`
273 ``spatiallySampledMetrics`` : `astropy.table.Table`
274 Image quality metrics spatially sampled locations.
275 """
277 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory()
279 spatiallySampledMetrics = afwTable.SourceCatalog(self.schema)
280 spatiallySampledMetrics.getTable().setIdFactory(idFactory)
282 self.metricSources.run(mask=science.mask, seed=difference.info.id, catalog=spatiallySampledMetrics)
284 # Mask plane bits are looked up in the global afw dict, not the
285 # image's own. A plane that is not registered in this process has no
286 # pixels set in any image.
287 metricsMaskPlanes = {}
288 for maskPlane in self.config.metricsMaskPlanes:
289 try:
290 metricsMaskPlanes[maskPlane] = afwImage.Mask.getPlaneBitMask(maskPlane)
291 except InvalidParameterError:
292 self.log.info("Mask plane %s is not registered; its metrics will be zero.", maskPlane)
293 metricsMaskPlanes[maskPlane] = 0
295 for src in spatiallySampledMetrics:
296 self._evaluateLocalMetric(src, science, template, difference, diaSources,
297 metricsMaskPlanes=metricsMaskPlanes,
298 psfMatchingKernel=psfMatchingKernel)
299 spatiallySampledMetrics = spatiallySampledMetrics.copy(deep=True).asAstropy()
300 return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics)
302 def _evaluateLocalMetric(self, src, science, template, difference, diaSources,
303 metricsMaskPlanes, psfMatchingKernel):
304 """Calculate image quality metrics at spatially sampled locations.
306 Parameters
307 ----------
308 src : `lsst.afw.table.SourceRecord`
309 The source record to be updated with metric calculations.
310 diaSources : `lsst.afw.table.SourceCatalog`
311 The catalog of detected sources.
312 science : `lsst.afw.image.Exposure`
313 The science image.
314 difference : `lsst.afw.image.Exposure`
315 Result of subtracting template from the science image.
316 metricsMaskPlanes : `dict` [`str`, `int`]
317 Mask planes to calculate metrics from, mapped to their bitmask.
318 The bitmask is 0 for a plane that is not registered.
319 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel`
320 The PSF matching kernel of the subtraction to evaluate.
321 """
322 bbox = src.getFootprint().getBBox()
323 pix = bbox.getCenter()
324 src.set('science_psfSize', getPsfFwhm(science.psf, position=pix))
325 try:
326 src.set('template_psfSize', getPsfFwhm(template.psf, position=pix))
327 except (InvalidParameterError, RangeError):
328 src.set('template_psfSize', np.nan)
330 metricRegionSize = 100
331 bbox.grow(metricRegionSize)
332 bbox = bbox.clippedTo(science.getBBox())
333 nPix = bbox.getArea()
334 pixScale = science.wcs.getPixelScale(bbox.getCenter())
335 area = nPix*pixScale.asDegrees()**2
336 peak = src.getFootprint().getPeaks()[0]
337 src.set('x', peak['i_x'])
338 src.set('y', peak['i_y'])
339 src.setCoord(science.wcs.pixelToSky(peak['i_x'], peak['i_y']))
340 selectSources = diaSources[bbox.contains(diaSources['x'], diaSources['y'])]
341 sourceDensity = len(selectSources)/area
342 dipoleSources = selectSources[selectSources["isDipole"]]
343 dipoleDensity = len(dipoleSources)/area
345 if dipoleSources: 345 ↛ 346line 345 didn't jump to line 346 because the condition on line 345 was never true
346 meanDipoleOrientation = angleMean(dipoleSources["dipoleAngle"])
347 src.set('dipole_direction', meanDipoleOrientation)
348 meanDipoleSeparation = np.mean(dipoleSources["dipoleLength"])
349 src.set('dipole_separation', meanDipoleSeparation)
351 templateVal = np.median(template[bbox].image.array)
352 templateVar = np.median(template[bbox].variance.array)
353 scienceVal = np.median(science[bbox].image.array)
354 scienceVar = np.median(science[bbox].variance.array)
355 diffimVal = np.median(difference[bbox].image.array)
356 diffimVar = np.median(difference[bbox].variance.array)
357 src.set('source_density', sourceDensity)
358 src.set('dipole_density', dipoleDensity)
359 src.set('template_value', templateVal)
360 src.set('template_variance', templateVar)
361 src.set('science_value', scienceVal)
362 src.set('science_variance', scienceVar)
363 src.set('diffim_value', diffimVal)
364 src.set('diffim_variance', diffimVar)
365 src.set('diffim_chi2PerPix', self._diffimChi2PerPix(difference[bbox]))
366 differenceMask = difference.mask[bbox].array
367 for maskPlane, bitMask in metricsMaskPlanes.items():
368 src.set("%s_mask_fraction"%maskPlane.lower(),
369 np.count_nonzero(differenceMask & bitMask)/differenceMask.size
370 )
372 krnlSum, dx, dy, direction, length = getKernelCenterDisplacement(
373 psfMatchingKernel, src.get('x'), src.get('y'))
375 point1 = lsst.geom.SpherePoint(
376 src.get('coord_ra'), src.get('coord_dec'),
377 lsst.geom.radians)
378 point2 = science.wcs.pixelToSky(src.get('x') + dx, src.get('y') + dy)
379 bearing = point1.bearingTo(point2)
380 pa_ref_angle = lsst.geom.Angle(np.pi/2, lsst.geom.radians)
381 pa = pa_ref_angle - bearing
382 # Wrap around to get Delta_RA from -pi to +pi
383 pa = pa.wrapCtr()
384 position_angle = pa.asRadians()
386 src.set('psfMatchingKernel_sum', krnlSum)
387 src.set('psfMatchingKernel_dx', dx)
388 src.set('psfMatchingKernel_dy', dy)
389 src.set('psfMatchingKernel_length', length*pixScale.asArcseconds())
390 src.set('psfMatchingKernel_position_angle', position_angle) # in E of N position angle
391 src.set('psfMatchingKernel_direction', direction) # direction offset in detector
393 src.set('psfMatchingKernel_residualNorm',
394 self._psfMatchResidualNorm(psfMatchingKernel, science.psf, template.psf,
395 src.get('x'), src.get('y')))
397 def _diffimChi2PerPix(self, difference):
398 """Robust normalized noise of the difference image.
400 Computes ``(1.4826 * MAD(image))^2 / median(variance)`` on
401 background-only pixels. The 1.4826 factor rescales the
402 Median Absolute Deviation to a Gaussian-sigma estimate, so the
403 ratio has expectation 1.0 on pure noise.
405 Returns NaN if no usable pixels remain or the median variance
406 is non-positive.
407 """
408 image = difference.image.array
409 variance = difference.variance.array
410 mask = difference.mask
411 excludePlanes = [p for p in ("DETECTED", "DETECTED_NEGATIVE", "BAD",
412 "SAT", "EDGE", "NO_DATA")
413 if p in mask.getMaskPlaneDict()]
414 if excludePlanes: 414 ↛ 418line 414 didn't jump to line 418 because the condition on line 414 was always true
415 badBits = mask.getPlaneBitMask(excludePlanes)
416 good = (mask.array & badBits) == 0
417 else:
418 good = np.ones(image.shape, dtype=bool)
419 good &= np.isfinite(image) & np.isfinite(variance) & (variance > 0)
420 if not np.any(good): 420 ↛ 421line 420 didn't jump to line 421 because the condition on line 420 was never true
421 return np.nan
422 varMed = np.median(variance[good])
423 if not np.isfinite(varMed) or varMed <= 0: 423 ↛ 424line 423 didn't jump to line 424 because the condition on line 423 was never true
424 return np.nan
425 imgGood = image[good]
426 mad = np.median(np.abs(imgGood - np.median(imgGood)))
427 return float((1.4826*mad)**2/varMed)
429 def _psfMatchResidualNorm(self, kernel, sciencePsf, templatePsf, x, y):
430 """Relative L2 norm of the PSF-match residual at (x, y).
432 Convolves the template PSF with the matching kernel evaluated at
433 (x, y) and compares to the science PSF at the same position, both
434 renormalized to unit sum so the result captures shape mismatch only.
435 The kernel is assumed to convolve the template.
437 Returns NaN if either PSF cannot be evaluated at the position or
438 the resulting images cannot be normalized.
439 """
440 point = lsst.geom.Point2D(x, y)
441 try:
442 psfSci = sciencePsf.computeKernelImage(point).array
443 psfTmp = templatePsf.computeKernelImage(point).array
444 except (InvalidParameterError, RangeError):
445 return np.nan
447 kImage = afwImage.ImageD(kernel.getDimensions())
448 kernel.computeImage(kImage, doNormalize=True, x=x, y=y)
449 matched = scipy.signal.fftconvolve(psfTmp, kImage.array, mode='same')
451 matchedSum = matched.sum()
452 psfSciSum = psfSci.sum()
453 if matchedSum <= 0 or psfSciSum <= 0: 453 ↛ 454line 453 didn't jump to line 454 because the condition on line 453 was never true
454 return np.nan
455 matched = matched/matchedSum
456 psfSci = psfSci/psfSciSum
458 # PSF stamps may differ in size between the two exposures; crop
459 # both to a common centered region before differencing.
460 h = min(matched.shape[0], psfSci.shape[0])
461 w = min(matched.shape[1], psfSci.shape[1])
463 def _crop(a):
464 sy = (a.shape[0] - h)//2
465 sx = (a.shape[1] - w)//2
466 return a[sy:sy + h, sx:sx + w]
467 matched = _crop(matched)
468 psfSci = _crop(psfSci)
470 sciNorm = np.sqrt(np.sum(psfSci**2))
471 if sciNorm == 0: 471 ↛ 472line 471 didn't jump to line 472 because the condition on line 471 was never true
472 return np.nan
473 return float(np.sqrt(np.sum((matched - psfSci)**2))/sciNorm)