Coverage for python / lsst / ip / diffim / imageDecorrelation.py: 15%
288 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-18 09:01 +0000
1#
2# LSST Data Management System
3# Copyright 2016 AURA/LSST.
4#
5# This product includes software developed by the
6# LSST Project (http://www.lsst.org/).
7#
8# This program is free software: you can redistribute it and/or modify
9# it under the terms of the GNU General Public License as published by
10# the Free Software Foundation, either version 3 of the License, or
11# (at your option) any later version.
12#
13# This program is distributed in the hope that it will be useful,
14# but WITHOUT ANY WARRANTY; without even the implied warranty of
15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16# GNU General Public License for more details.
17#
18# You should have received a copy of the LSST License Statement and
19# the GNU General Public License along with this program. If not,
20# see <https://www.lsstcorp.org/LegalNotices/>.
21#
23import numpy as np
25import lsst.afw.image as afwImage
26import lsst.afw.math as afwMath
27import lsst.geom as geom
28import lsst.meas.algorithms as measAlg
29import lsst.pex.config as pexConfig
30import lsst.pipe.base as pipeBase
31from lsst.pex.exceptions import InvalidParameterError
32from lsst.utils.timer import timeMethod
34from .imageMapReduce import (ImageMapReduceConfig, ImageMapReduceTask,
35 ImageMapper)
36from .utils import computeAveragePsf
38__all__ = ("DecorrelateALKernelTask", "DecorrelateALKernelConfig",
39 "DecorrelateALKernelMapper", "DecorrelateALKernelMapReduceConfig",
40 "DecorrelateALKernelSpatialConfig", "DecorrelateALKernelSpatialTask")
43class DecorrelateALKernelConfig(pexConfig.Config):
44 """Configuration parameters for the DecorrelateALKernelTask
45 """
47 ignoreMaskPlanes = pexConfig.ListField(
48 dtype=str,
49 doc="""Mask planes to ignore for sigma-clipped statistics""",
50 default=("INTRP", "EDGE", "DETECTED", "SAT", "CR", "BAD", "NO_DATA", "DETECTED_NEGATIVE")
51 )
52 completeVarPlanePropagation = pexConfig.Field(
53 dtype=bool,
54 default=False,
55 doc="Compute the full effect of the decorrelated matching kernel on the variance plane."
56 " Otherwise use a model weighed sum of the input variances."
57 )
60class DecorrelateALKernelTask(pipeBase.Task):
61 """Decorrelate the effect of convolution by Alard-Lupton matching kernel in image difference
63 """
64 ConfigClass = DecorrelateALKernelConfig
65 _DefaultName = "ip_diffim_decorrelateALKernel"
67 def __init__(self, *args, **kwargs):
68 """Create the image decorrelation Task
70 Parameters
71 ----------
72 args :
73 arguments to be passed to ``lsst.pipe.base.task.Task.__init__``
74 kwargs :
75 keyword arguments to be passed to ``lsst.pipe.base.task.Task.__init__``
76 """
77 pipeBase.Task.__init__(self, *args, **kwargs)
79 self.statsControl = afwMath.StatisticsControl()
80 self.statsControl.setNumSigmaClip(3.)
81 self.statsControl.setNumIter(3)
82 self.statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.ignoreMaskPlanes))
84 def computeVarianceMean(self, exposure):
85 statObj = afwMath.makeStatistics(exposure.variance,
86 exposure.mask,
87 afwMath.MEANCLIP, self.statsControl)
88 var = statObj.getValue(afwMath.MEANCLIP)
89 return var
91 @timeMethod
92 def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatchingKernel,
93 preConvKernel=None, xcen=None, ycen=None, svar=None, tvar=None,
94 templateMatched=True, preConvMode=False, **kwargs):
95 """Perform decorrelation of an image difference or of a score difference exposure.
97 Corrects the difference or score image due to the convolution of the
98 templateExposure with the A&L PSF matching kernel.
99 See [DMTN-021, Equation 1](http://dmtn-021.lsst.io/#equation-1) and
100 [DMTN-179](http://dmtn-179.lsst.io/) for details.
102 Parameters
103 ----------
104 scienceExposure : `lsst.afw.image.Exposure`
105 The original science exposure (before pre-convolution, if ``preConvMode==True``).
106 templateExposure : `lsst.afw.image.Exposure`
107 The original template exposure warped, but not psf-matched, to the science exposure.
108 subtractedExposure : `lsst.afw.image.Exposure`
109 the subtracted exposure produced by
110 `ip_diffim.ImagePsfMatchTask.subtractExposures()`. The `subtractedExposure` must
111 inherit its PSF from `exposure`, see notes below.
112 psfMatchingKernel : `lsst.afw.detection.Psf`
113 An (optionally spatially-varying) PSF matching kernel produced
114 by `ip_diffim.ImagePsfMatchTask.subtractExposures()`.
115 preConvKernel : `lsst.afw.math.Kernel`, optional
116 If not `None`, then the `scienceExposure` was pre-convolved with (the reflection of)
117 this kernel. Must be normalized to sum to 1.
118 Allowed only if ``templateMatched==True`` and ``preConvMode==True``.
119 Defaults to the PSF of the science exposure at the image center.
120 xcen : `float`, optional
121 X-pixel coordinate to use for computing constant matching kernel to use
122 If `None` (default), then use the center of the image.
123 ycen : `float`, optional
124 Y-pixel coordinate to use for computing constant matching kernel to use
125 If `None` (default), then use the center of the image.
126 svar : `float`, optional
127 Image variance for science image
128 If `None` (default) then compute the variance over the entire input science image.
129 tvar : `float`, optional
130 Image variance for template image
131 If `None` (default) then compute the variance over the entire input template image.
132 templateMatched : `bool`, optional
133 If True, the template exposure was matched (convolved) to the science exposure.
134 See also notes below.
135 preConvMode : `bool`, optional
136 If True, ``subtractedExposure`` is assumed to be a likelihood difference image
137 and will be noise corrected as a likelihood image.
138 **kwargs
139 Additional keyword arguments propagated from DecorrelateALKernelSpatialTask.
141 Returns
142 -------
143 result : `lsst.pipe.base.Struct`
144 - ``correctedExposure`` : the decorrelated diffim
146 Notes
147 -----
148 If ``preConvMode==True``, ``subtractedExposure`` is assumed to be a
149 score image and the noise correction for likelihood images
150 is applied. The resulting image is an optimal detection likelihood image
151 when the templateExposure has noise. (See DMTN-179) If ``preConvKernel`` is
152 not specified, the PSF of ``scienceExposure`` is assumed as pre-convolution kernel.
154 The ``subtractedExposure`` is NOT updated. The returned ``correctedExposure``
155 has an updated but spatially fixed PSF. It is calculated as the center of
156 image PSF corrected by the center of image matching kernel.
158 If ``templateMatched==True``, the templateExposure was matched (convolved)
159 to the ``scienceExposure`` by ``psfMatchingKernel`` during image differencing.
160 Otherwise the ``scienceExposure`` was matched (convolved) by ``psfMatchingKernel``.
161 In either case, note that the original template and science images are required,
162 not the psf-matched version.
164 This task discards the variance plane of ``subtractedExposure`` and re-computes
165 it from the variance planes of ``scienceExposure`` and ``templateExposure``.
166 The image plane of ``subtractedExposure`` must be at the photometric level
167 set by the AL PSF matching in `ImagePsfMatchTask.subtractExposures`.
168 The assumptions about the photometric level are controlled by the
169 `templateMatched` option in this task.
171 Here we currently convert a spatially-varying matching kernel into a constant kernel,
172 just by computing it at the center of the image (tickets DM-6243, DM-6244).
174 We are also using a constant accross-the-image measure of sigma (sqrt(variance)) to compute
175 the decorrelation kernel.
177 TODO DM-23857 As part of the spatially varying correction implementation
178 consider whether returning a Struct is still necessary.
179 """
180 if preConvKernel is not None and not (templateMatched and preConvMode):
181 raise ValueError("Pre-convolution kernel is allowed only if "
182 "preConvMode==True and templateMatched==True.")
184 spatialKernel = psfMatchingKernel
185 kimg = afwImage.ImageD(spatialKernel.getDimensions())
186 bbox = subtractedExposure.getBBox()
187 if xcen is None:
188 xcen = (bbox.getBeginX() + bbox.getEndX()) / 2.
189 if ycen is None:
190 ycen = (bbox.getBeginY() + bbox.getEndY()) / 2.
191 self.log.info("Using matching kernel computed at (%d, %d)", xcen, ycen)
192 spatialKernel.computeImage(kimg, False, xcen, ycen)
194 preConvImg = None
195 if preConvMode:
196 if preConvKernel is None:
197 pos = scienceExposure.getPsf().getAveragePosition()
198 preConvKernel = scienceExposure.getPsf().getLocalKernel(pos)
199 preConvImg = afwImage.ImageD(preConvKernel.getDimensions())
200 preConvKernel.computeImage(preConvImg, True)
202 if svar is None:
203 svar = self.computeVarianceMean(scienceExposure)
204 if tvar is None:
205 tvar = self.computeVarianceMean(templateExposure)
206 self.log.info("Original variance plane means. Science:%.5e, warped template:%.5e)",
207 svar, tvar)
209 # Should not happen unless entire image has been masked, which could happen
210 # if this is a small subimage of the main exposure. In this case, just return a full NaN
211 # exposure
212 if np.isnan(svar) or np.isnan(tvar):
213 # Double check that one of the exposures is all NaNs
214 if (np.all(np.isnan(scienceExposure.image.array))
215 or np.all(np.isnan(templateExposure.image.array))):
216 self.log.warning('Template or science image is entirely NaNs: skipping decorrelation.')
217 outExposure = subtractedExposure.clone()
218 return pipeBase.Struct(correctedExposure=outExposure, )
220 if templateMatched:
221 # Regular subtraction, we convolved the template
222 self.log.info("Decorrelation after template image convolution")
223 varianceMean = svar
224 targetVarianceMean = tvar
225 # variance plane of the image that is not convolved
226 variance = scienceExposure.variance.array
227 # Variance plane of the convolved image, before convolution.
228 targetVariance = templateExposure.variance.array
229 # If the template is convolved, the PSF of the difference image is
230 # that of the science image
231 psfImg = scienceExposure.getPsf().computeKernelImage(geom.Point2D(xcen, ycen))
232 else:
233 # We convolved the science image
234 self.log.info("Decorrelation after science image convolution")
235 varianceMean = tvar
236 targetVarianceMean = svar
237 # variance plane of the image that is not convolved
238 variance = templateExposure.variance.array
239 # Variance plane of the convolved image, before convolution.
240 targetVariance = scienceExposure.variance.array
241 # If the science image is convolved, the PSF of the difference image
242 # is that of the template image, and will be a CoaddPsf which might
243 # not be defined for the entire image.
244 # Try the simple calculation first, and fall back on a slower method
245 # if the PSF of the template is not defined in the center.
246 try:
247 psfImg = templateExposure.getPsf().computeKernelImage(geom.Point2D(xcen, ycen))
248 except InvalidParameterError:
249 psfImg = computeAveragePsf(templateExposure, psfExposureBuffer=0.05, psfExposureGrid=100)
251 # The maximal correction value converges to sqrt(targetVarianceMean/varianceMean).
252 # Correction divergence warning if the correction exceeds 4 orders of magnitude.
253 mOverExpVar = targetVarianceMean/varianceMean
254 if mOverExpVar > 1e8:
255 self.log.warning("Diverging correction: matched image variance is "
256 " much larger than the unconvolved one's"
257 ", targetVarianceMean/varianceMean:%.2e", mOverExpVar)
259 oldVarMean = self.computeVarianceMean(subtractedExposure)
260 self.log.info("Variance plane mean of uncorrected diffim: %f", oldVarMean)
262 kArr = kimg.array
263 diffimShape = subtractedExposure.image.array.shape
264 psfShape = psfImg.array.shape
266 if preConvMode:
267 self.log.info("Decorrelation of likelihood image")
268 self.computeCommonShape(preConvImg.array.shape, kArr.shape,
269 psfShape, diffimShape)
270 corr = self.computeScoreCorrection(kArr, varianceMean, targetVarianceMean, preConvImg.array)
271 else:
272 self.log.info("Decorrelation of difference image")
273 self.computeCommonShape(kArr.shape, psfShape, diffimShape)
274 corr = self.computeDiffimCorrection(kArr, varianceMean, targetVarianceMean)
276 correctedImage = self.computeCorrectedImage(corr.corrft, subtractedExposure.image.array)
277 # TODO DM-47461: only decorrelate the PSF if it is calculated for the
278 # difference image. If it is taken directly from the science (or template)
279 # image the decorrelation calculation is incorrect.
280 # correctedPsf = self.computeCorrectedDiffimPsf(corr.corrft, psfImg.array)
282 # The subtracted exposure variance plane is already correlated, we cannot propagate
283 # it through another convolution; instead we need to use the uncorrelated originals
284 # The whitening should scale it to varianceMean + targetVarianceMean on average
285 if self.config.completeVarPlanePropagation:
286 self.log.debug("Using full variance plane calculation in decorrelation")
287 correctedVariance = self.calculateVariancePlane(
288 variance, targetVariance,
289 varianceMean, targetVarianceMean, corr.cnft, corr.crft)
290 else:
291 self.log.debug("Using estimated variance plane calculation in decorrelation")
292 correctedVariance = self.estimateVariancePlane(
293 variance, targetVariance,
294 corr.cnft, corr.crft)
296 # Determine the common shape
297 kSum = np.sum(kArr)
298 kSumSq = kSum*kSum
299 self.log.debug("Matching kernel sum: %.3e", kSum)
300 if not templateMatched:
301 # ImagePsfMatch.subtractExposures re-scales the difference in
302 # the science image convolution mode
303 correctedVariance /= kSumSq
304 subtractedExposure.image.array[...] = correctedImage # Allow for numpy type casting
305 subtractedExposure.variance.array[...] = correctedVariance
306 # TODO DM-47461
307 # subtractedExposure.setPsf(correctedPsf)
309 newVarMean = self.computeVarianceMean(subtractedExposure)
310 self.log.info("Variance plane mean of corrected diffim: %.5e", newVarMean)
312 # TODO DM-23857 As part of the spatially varying correction implementation
313 # consider whether returning a Struct is still necessary.
314 return pipeBase.Struct(correctedExposure=subtractedExposure, )
316 def computeCommonShape(self, *shapes):
317 """Calculate the common shape for FFT operations. Set `self.freqSpaceShape`
318 internally.
320 Parameters
321 ----------
322 shapes : one or more `tuple` of `int`
323 Shapes of the arrays. All must have the same dimensionality.
324 At least one shape must be provided.
326 Returns
327 -------
328 None.
330 Notes
331 -----
332 For each dimension, gets the smallest even number greater than or equal to
333 `N1+N2-1` where `N1` and `N2` are the two largest values.
334 In case of only one shape given, rounds up to even each dimension value.
335 """
336 S = np.array(shapes, dtype=int)
337 if len(shapes) > 2:
338 S.sort(axis=0)
339 S = S[-2:]
340 if len(shapes) > 1:
341 commonShape = np.sum(S, axis=0) - 1
342 else:
343 commonShape = S[0]
344 commonShape[commonShape % 2 != 0] += 1
345 self.freqSpaceShape = tuple(commonShape)
346 self.log.info("Common frequency space shape %s", self.freqSpaceShape)
348 @staticmethod
349 def padCenterOriginArray(A, newShape: tuple, useInverse=False):
350 """Zero pad an image where the origin is at the center and replace the
351 origin to the corner as required by the periodic input of FFT. Implement also
352 the inverse operation, crop the padding and re-center data.
354 Parameters
355 ----------
356 A : `numpy.ndarray`
357 An array to copy from.
358 newShape : `tuple` of `int`
359 The dimensions of the resulting array. For padding, the resulting array
360 must be larger than A in each dimension. For the inverse operation this
361 must be the original, before padding size of the array.
362 useInverse : bool, optional
363 Selector of forward, add padding, operation (False)
364 or its inverse, crop padding, operation (True).
366 Returns
367 -------
368 R : `numpy.ndarray`
369 The padded or unpadded array with shape of `newShape` and the same dtype as A.
371 Notes
372 -----
373 For odd dimensions, the splitting is rounded to
374 put the center pixel into the new corner origin (0,0). This is to be consistent
375 e.g. for a dirac delta kernel that is originally located at the center pixel.
376 """
378 # The forward and inverse operations should round odd dimension halves at the opposite
379 # sides to get the pixels back to their original positions.
380 if not useInverse:
381 # Forward operation: First and second halves with respect to the axes of A.
382 firstHalves = [x//2 for x in A.shape]
383 secondHalves = [x-y for x, y in zip(A.shape, firstHalves)]
384 else:
385 # Inverse operation: Opposite rounding
386 secondHalves = [x//2 for x in newShape]
387 firstHalves = [x-y for x, y in zip(newShape, secondHalves)]
389 R = np.zeros_like(A, shape=newShape)
390 R[-firstHalves[0]:, -firstHalves[1]:] = A[:firstHalves[0], :firstHalves[1]]
391 R[:secondHalves[0], -firstHalves[1]:] = A[-secondHalves[0]:, :firstHalves[1]]
392 R[:secondHalves[0], :secondHalves[1]] = A[-secondHalves[0]:, -secondHalves[1]:]
393 R[-firstHalves[0]:, :secondHalves[1]] = A[:firstHalves[0], -secondHalves[1]:]
394 return R
396 def computeDiffimCorrection(self, kappa, svar, tvar):
397 """Compute the Lupton decorrelation post-convolution kernel for decorrelating an
398 image difference, based on the PSF-matching kernel.
400 Parameters
401 ----------
402 kappa : `numpy.ndarray` of `float`
403 A matching kernel 2-d numpy.array derived from Alard & Lupton PSF matching.
404 svar : `float` > 0.
405 Average variance of science image used for PSF matching.
406 tvar : `float` > 0.
407 Average variance of the template (matched) image used for PSF matching.
409 Returns
410 -------
411 corrft : `numpy.ndarray` of `float`
412 The frequency space representation of the correction. The array is real (dtype float).
413 Shape is `self.freqSpaceShape`.
415 cnft, crft : `numpy.ndarray` of `complex`
416 The overall convolution (pre-conv, PSF matching, noise correction) kernel
417 for the science and template images, respectively for the variance plane
418 calculations. These are intermediate results in frequency space.
420 Notes
421 -----
422 The maximum correction factor converges to `sqrt(tvar/svar)` towards high frequencies.
423 This should be a plausible value.
424 """
425 kSum = np.sum(kappa) # We scale the decorrelation to preserve fluxes
426 kappa = self.padCenterOriginArray(kappa, self.freqSpaceShape)
427 kft = np.fft.fft2(kappa)
428 kftAbsSq = np.real(np.conj(kft) * kft)
430 denom = svar + tvar * kftAbsSq
431 corrft = np.sqrt((svar + tvar * kSum*kSum) / denom)
432 cnft = corrft
433 crft = kft*corrft
434 return pipeBase.Struct(corrft=corrft, cnft=cnft, crft=crft)
436 def computeScoreCorrection(self, kappa, svar, tvar, preConvArr):
437 """Compute the correction kernel for a score image.
439 Parameters
440 ----------
441 kappa : `numpy.ndarray`
442 A matching kernel 2-d numpy.array derived from Alard & Lupton PSF matching.
443 svar : `float`
444 Average variance of science image used for PSF matching (before pre-convolution).
445 tvar : `float`
446 Average variance of the template (matched) image used for PSF matching.
447 preConvArr : `numpy.ndarray`
448 The pre-convolution kernel of the science image. It should be the PSF
449 of the science image or an approximation of it. It must be normed to sum 1.
451 Returns
452 -------
453 corrft : `numpy.ndarray` of `float`
454 The frequency space representation of the correction. The array is real (dtype float).
455 Shape is `self.freqSpaceShape`.
456 cnft, crft : `numpy.ndarray` of `complex`
457 The overall convolution (pre-conv, PSF matching, noise correction) kernel
458 for the science and template images, respectively for the variance plane
459 calculations. These are intermediate results in frequency space.
461 Notes
462 -----
463 To be precise, the science image should be _correlated_ by ``preConvArray`` but this
464 does not matter for this calculation.
466 ``cnft``, ``crft`` contain the scaling factor as well.
468 """
469 kSum = np.sum(kappa)
470 kappa = self.padCenterOriginArray(kappa, self.freqSpaceShape)
471 kft = np.fft.fft2(kappa)
472 preConvArr = self.padCenterOriginArray(preConvArr, self.freqSpaceShape)
473 preFt = np.fft.fft2(preConvArr)
474 preFtAbsSq = np.real(np.conj(preFt) * preFt)
475 kftAbsSq = np.real(np.conj(kft) * kft)
476 # Avoid zero division, though we don't normally approach `tiny`.
477 # We have numerical noise instead.
478 tiny = np.finfo(preFtAbsSq.dtype).tiny * 1000.
479 flt = preFtAbsSq < tiny
480 # If we pre-convolve to avoid deconvolution in AL, then kftAbsSq / preFtAbsSq
481 # theoretically expected to diverge to +inf. But we don't care about the convergence
482 # properties here, S goes to 0 at these frequencies anyway.
483 preFtAbsSq[flt] = tiny
484 denom = svar + tvar * kftAbsSq / preFtAbsSq
485 corrft = (svar + tvar * kSum*kSum) / denom
486 cnft = np.conj(preFt)*corrft
487 crft = kft*corrft
488 return pipeBase.Struct(corrft=corrft, cnft=cnft, crft=crft)
490 @staticmethod
491 def estimateVariancePlane(vplane1, vplane2, c1ft, c2ft):
492 """Estimate the variance planes.
494 The estimation assumes that around each pixel the surrounding
495 pixels' sigmas within the convolution kernel are the same.
497 Parameters
498 ----------
499 vplane1, vplane2 : `numpy.ndarray` of `float`
500 Variance planes of the original (before pre-convolution or matching)
501 exposures.
502 c1ft, c2ft : `numpy.ndarray` of `complex`
503 The overall convolution that includes the matching and the
504 afterburner in frequency space. The result of either
505 ``computeScoreCorrection`` or ``computeDiffimCorrection``.
507 Returns
508 -------
509 vplaneD : `numpy.ndarray` of `float`
510 The estimated variance plane of the difference/score image
511 as a weighted sum of the input variances.
513 Notes
514 -----
515 See DMTN-179 Section 5 about the variance plane calculations.
516 """
517 w1 = np.sum(np.real(np.conj(c1ft)*c1ft)) / c1ft.size
518 w2 = np.sum(np.real(np.conj(c2ft)*c2ft)) / c2ft.size
519 # w1, w2: the frequency space sum of abs(c1)^2 is the same as in image
520 # space.
521 return vplane1*w1 + vplane2*w2
523 def calculateVariancePlane(self, vplane1, vplane2, varMean1, varMean2, c1ft, c2ft):
524 """Full propagation of the variance planes of the original exposures.
526 The original variance planes of independent pixels are convolved with the
527 image space square of the overall kernels.
529 Parameters
530 ----------
531 vplane1, vplane2 : `numpy.ndarray` of `float`
532 Variance planes of the original (before pre-convolution or matching)
533 exposures.
534 varMean1, varMean2 : `float`
535 Replacement average values for non-finite ``vplane1`` and ``vplane2`` values respectively.
537 c1ft, c2ft : `numpy.ndarray` of `complex`
538 The overall convolution that includes the matching and the
539 afterburner in frequency space. The result of either
540 ``computeScoreCorrection`` or ``computeDiffimCorrection``.
542 Returns
543 -------
544 vplaneD : `numpy.ndarray` of `float`
545 The variance plane of the difference/score images.
547 Notes
548 -----
549 See DMTN-179 Section 5 about the variance plane calculations.
551 Infs and NaNs are allowed and kept in the returned array.
552 """
553 D = np.real(np.fft.ifft2(c1ft))
554 c1SqFt = np.fft.fft2(D*D)
556 v1shape = vplane1.shape
557 filtInf = np.isinf(vplane1)
558 filtNan = np.isnan(vplane1)
559 # This copy could be eliminated if inf/nan handling were go into padCenterOriginArray
560 vplane1 = np.copy(vplane1)
561 vplane1[filtInf | filtNan] = varMean1
562 D = self.padCenterOriginArray(vplane1, self.freqSpaceShape)
563 v1 = np.real(np.fft.ifft2(np.fft.fft2(D) * c1SqFt))
564 v1 = self.padCenterOriginArray(v1, v1shape, useInverse=True)
565 v1[filtNan] = np.nan
566 v1[filtInf] = np.inf
568 D = np.real(np.fft.ifft2(c2ft))
569 c2SqFt = np.fft.fft2(D*D)
571 v2shape = vplane2.shape
572 filtInf = np.isinf(vplane2)
573 filtNan = np.isnan(vplane2)
574 vplane2 = np.copy(vplane2)
575 vplane2[filtInf | filtNan] = varMean2
576 D = self.padCenterOriginArray(vplane2, self.freqSpaceShape)
577 v2 = np.real(np.fft.ifft2(np.fft.fft2(D) * c2SqFt))
578 v2 = self.padCenterOriginArray(v2, v2shape, useInverse=True)
579 v2[filtNan] = np.nan
580 v2[filtInf] = np.inf
582 return v1 + v2
584 def computeCorrectedDiffimPsf(self, corrft, psfOld):
585 """Compute the (decorrelated) difference image's new PSF.
587 Parameters
588 ----------
589 corrft : `numpy.ndarray`
590 The frequency space representation of the correction calculated by
591 `computeCorrection`. Shape must be `self.freqSpaceShape`.
592 psfOld : `numpy.ndarray`
593 The psf of the difference image to be corrected.
595 Returns
596 -------
597 correctedPsf : `lsst.meas.algorithms.KernelPsf`
598 The corrected psf, same shape as `psfOld`, sum normed to 1.
600 Notes
601 -----
602 There is no algorithmic guarantee that the corrected psf can
603 meaningfully fit to the same size as the original one.
604 """
605 psfShape = psfOld.shape
606 psfNew = self.padCenterOriginArray(psfOld, self.freqSpaceShape)
607 psfNew = np.fft.fft2(psfNew)
608 psfNew *= corrft
609 psfNew = np.fft.ifft2(psfNew)
610 psfNew = psfNew.real
611 psfNew = self.padCenterOriginArray(psfNew, psfShape, useInverse=True)
612 psfNew = psfNew/psfNew.sum()
614 psfcI = afwImage.ImageD(geom.Extent2I(psfShape[1], psfShape[0]))
615 psfcI.array = psfNew
616 psfcK = afwMath.FixedKernel(psfcI)
617 correctedPsf = measAlg.KernelPsf(psfcK)
618 return correctedPsf
620 def computeCorrectedImage(self, corrft, imgOld):
621 """Compute the decorrelated difference image.
623 Parameters
624 ----------
625 corrft : `numpy.ndarray`
626 The frequency space representation of the correction calculated by
627 `computeCorrection`. Shape must be `self.freqSpaceShape`.
628 imgOld : `numpy.ndarray`
629 The difference image to be corrected.
631 Returns
632 -------
633 imgNew : `numpy.ndarray`
634 The corrected image, same size as the input.
635 """
636 expShape = imgOld.shape
637 imgNew = np.copy(imgOld)
638 filtInf = np.isinf(imgNew)
639 filtNan = np.isnan(imgNew)
640 imgNew[filtInf] = np.nan
641 imgNew[filtInf | filtNan] = np.nanmean(imgNew)
642 imgNew = self.padCenterOriginArray(imgNew, self.freqSpaceShape)
643 imgNew = np.fft.fft2(imgNew)
644 imgNew *= corrft
645 imgNew = np.fft.ifft2(imgNew)
646 imgNew = imgNew.real
647 imgNew = self.padCenterOriginArray(imgNew, expShape, useInverse=True)
648 imgNew[filtNan] = np.nan
649 imgNew[filtInf] = np.inf
650 return imgNew
653class DecorrelateALKernelMapper(DecorrelateALKernelTask, ImageMapper):
654 """Task to be used as an ImageMapper for performing
655 A&L decorrelation on subimages on a grid across a A&L difference image.
657 This task subclasses DecorrelateALKernelTask in order to implement
658 all of that task's configuration parameters, as well as its `run` method.
659 """
661 ConfigClass = DecorrelateALKernelConfig
662 _DefaultName = 'ip_diffim_decorrelateALKernelMapper'
664 def __init__(self, *args, **kwargs):
665 DecorrelateALKernelTask.__init__(self, *args, **kwargs)
667 def run(self, subExposure, expandedSubExposure, fullBBox,
668 template, science, alTaskResult=None, psfMatchingKernel=None,
669 preConvKernel=None, **kwargs):
670 """Perform decorrelation operation on `subExposure`, using
671 `expandedSubExposure` to allow for invalid edge pixels arising from
672 convolutions.
674 This method performs A&L decorrelation on `subExposure` using
675 local measures for image variances and PSF. `subExposure` is a
676 sub-exposure of the non-decorrelated A&L diffim. It also
677 requires the corresponding sub-exposures of the template
678 (`template`) and science (`science`) exposures.
680 Parameters
681 ----------
682 subExposure : `lsst.afw.image.Exposure`
683 the sub-exposure of the diffim
684 expandedSubExposure : `lsst.afw.image.Exposure`
685 the expanded sub-exposure upon which to operate
686 fullBBox : `lsst.geom.Box2I`
687 the bounding box of the original exposure
688 template : `lsst.afw.image.Exposure`
689 the corresponding sub-exposure of the template exposure
690 science : `lsst.afw.image.Exposure`
691 the corresponding sub-exposure of the science exposure
692 alTaskResult : `lsst.pipe.base.Struct`
693 the result of A&L image differencing on `science` and
694 `template`, importantly containing the resulting
695 `psfMatchingKernel`. Can be `None`, only if
696 `psfMatchingKernel` is not `None`.
697 psfMatchingKernel : Alternative parameter for passing the
698 A&L `psfMatchingKernel` directly.
699 preConvKernel : If not None, then pre-filtering was applied
700 to science exposure, and this is the pre-convolution
701 kernel.
702 kwargs :
703 additional keyword arguments propagated from
704 `ImageMapReduceTask.run`.
706 Returns
707 -------
708 A `pipeBase.Struct` containing:
710 - ``subExposure`` : the result of the `subExposure` processing.
711 - ``decorrelationKernel`` : the decorrelation kernel, currently
712 not used.
714 Notes
715 -----
716 This `run` method accepts parameters identical to those of
717 `ImageMapper.run`, since it is called from the
718 `ImageMapperTask`. See that class for more information.
719 """
720 templateExposure = template # input template
721 scienceExposure = science # input science image
722 if alTaskResult is None and psfMatchingKernel is None:
723 raise RuntimeError('Both alTaskResult and psfMatchingKernel cannot be None')
724 psfMatchingKernel = alTaskResult.psfMatchingKernel if alTaskResult is not None else psfMatchingKernel
726 # subExp and expandedSubExp are subimages of the (un-decorrelated) diffim!
727 # So here we compute corresponding subimages of templateExposure and scienceExposure
728 subExp2 = scienceExposure.Factory(scienceExposure, expandedSubExposure.getBBox())
729 subExp1 = templateExposure.Factory(templateExposure, expandedSubExposure.getBBox())
731 # Prevent too much log INFO verbosity from DecorrelateALKernelTask.run
732 logLevel = self.log.level
733 self.log.setLevel(self.log.WARNING)
734 res = DecorrelateALKernelTask.run(self, subExp2, subExp1, expandedSubExposure,
735 psfMatchingKernel, preConvKernel, **kwargs)
736 self.log.setLevel(logLevel) # reset the log level
738 diffim = res.correctedExposure.Factory(res.correctedExposure, subExposure.getBBox())
739 out = pipeBase.Struct(subExposure=diffim, )
740 return out
743class DecorrelateALKernelMapReduceConfig(ImageMapReduceConfig):
744 """Configuration parameters for the ImageMapReduceTask to direct it to use
745 DecorrelateALKernelMapper as its mapper for A&L decorrelation.
746 """
747 mapper = pexConfig.ConfigurableField(
748 doc='A&L decorrelation task to run on each sub-image',
749 target=DecorrelateALKernelMapper
750 )
753class DecorrelateALKernelSpatialConfig(pexConfig.Config):
754 """Configuration parameters for the DecorrelateALKernelSpatialTask.
755 """
756 decorrelateConfig = pexConfig.ConfigField(
757 dtype=DecorrelateALKernelConfig,
758 doc='DecorrelateALKernel config to use when running on complete exposure (non spatially-varying)',
759 )
761 decorrelateMapReduceConfig = pexConfig.ConfigField(
762 dtype=DecorrelateALKernelMapReduceConfig,
763 doc='DecorrelateALKernelMapReduce config to use when running on each sub-image (spatially-varying)',
764 )
766 ignoreMaskPlanes = pexConfig.ListField(
767 dtype=str,
768 doc="""Mask planes to ignore for sigma-clipped statistics""",
769 default=("INTRP", "EDGE", "DETECTED", "SAT", "CR", "BAD", "NO_DATA", "DETECTED_NEGATIVE")
770 )
772 def setDefaults(self):
773 self.decorrelateMapReduceConfig.gridStepX = self.decorrelateMapReduceConfig.gridStepY = 40
774 self.decorrelateMapReduceConfig.cellSizeX = self.decorrelateMapReduceConfig.cellSizeY = 41
775 self.decorrelateMapReduceConfig.borderSizeX = self.decorrelateMapReduceConfig.borderSizeY = 8
776 self.decorrelateMapReduceConfig.reducer.reduceOperation = 'average'
779class DecorrelateALKernelSpatialTask(pipeBase.Task):
780 """Decorrelate the effect of convolution by Alard-Lupton matching kernel in image difference
782 """
783 ConfigClass = DecorrelateALKernelSpatialConfig
784 _DefaultName = "ip_diffim_decorrelateALKernelSpatial"
786 def __init__(self, *args, **kwargs):
787 """Create the image decorrelation Task
789 Parameters
790 ----------
791 args :
792 arguments to be passed to
793 `lsst.pipe.base.task.Task.__init__`
794 kwargs :
795 additional keyword arguments to be passed to
796 `lsst.pipe.base.task.Task.__init__`
797 """
798 pipeBase.Task.__init__(self, *args, **kwargs)
800 self.statsControl = afwMath.StatisticsControl()
801 self.statsControl.setNumSigmaClip(3.)
802 self.statsControl.setNumIter(3)
803 self.statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.ignoreMaskPlanes))
805 def computeVarianceMean(self, exposure):
806 """Compute the mean of the variance plane of `exposure`.
807 """
808 statObj = afwMath.makeStatistics(exposure.variance,
809 exposure.mask,
810 afwMath.MEANCLIP, self.statsControl)
811 var = statObj.getValue(afwMath.MEANCLIP)
812 return var
814 def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatchingKernel,
815 spatiallyVarying=True, preConvKernel=None, templateMatched=True, preConvMode=False):
816 """Perform decorrelation of an image difference exposure.
818 Decorrelates the diffim due to the convolution of the
819 templateExposure with the A&L psfMatchingKernel. If
820 `spatiallyVarying` is True, it utilizes the spatially varying
821 matching kernel via the `imageMapReduce` framework to perform
822 spatially-varying decorrelation on a grid of subExposures.
824 Parameters
825 ----------
826 scienceExposure : `lsst.afw.image.Exposure`
827 the science Exposure used for PSF matching
828 templateExposure : `lsst.afw.image.Exposure`
829 the template Exposure used for PSF matching
830 subtractedExposure : `lsst.afw.image.Exposure`
831 the subtracted Exposure produced by `ip_diffim.ImagePsfMatchTask.subtractExposures()`
832 psfMatchingKernel : an (optionally spatially-varying) PSF matching kernel produced
833 by `ip_diffim.ImagePsfMatchTask.subtractExposures()`
834 spatiallyVarying : `bool`
835 if True, perform the spatially-varying operation
836 preConvKernel : `lsst.meas.algorithms.Psf`
837 if not none, the scienceExposure has been pre-filtered with this kernel. (Currently
838 this option is experimental.)
839 templateMatched : `bool`, optional
840 If True, the template exposure was matched (convolved) to the science exposure.
841 preConvMode : `bool`, optional
842 If True, ``subtractedExposure`` is assumed to be a likelihood difference image
843 and will be noise corrected as a likelihood image.
845 Returns
846 -------
847 results : `lsst.pipe.base.Struct`
848 a structure containing:
849 - ``correctedExposure`` : the decorrelated diffim
850 """
851 self.log.info('Running A&L decorrelation: spatiallyVarying=%r', spatiallyVarying)
853 svar = self.computeVarianceMean(scienceExposure)
854 tvar = self.computeVarianceMean(templateExposure)
855 if np.isnan(svar) or np.isnan(tvar): # Should not happen unless entire image has been masked.
856 # Double check that one of the exposures is all NaNs
857 if (np.all(np.isnan(scienceExposure.image.array))
858 or np.all(np.isnan(templateExposure.image.array))):
859 self.log.warning('Template or science image is entirely NaNs: skipping decorrelation.')
860 if np.isnan(svar):
861 svar = 1e-9
862 if np.isnan(tvar):
863 tvar = 1e-9
865 var = self.computeVarianceMean(subtractedExposure)
867 if spatiallyVarying:
868 self.log.info("Variance (science, template): (%f, %f)", svar, tvar)
869 self.log.info("Variance (uncorrected diffim): %f", var)
870 config = self.config.decorrelateMapReduceConfig
871 task = ImageMapReduceTask(config=config)
872 results = task.run(subtractedExposure, science=scienceExposure,
873 template=templateExposure, psfMatchingKernel=psfMatchingKernel,
874 preConvKernel=preConvKernel, forceEvenSized=True,
875 templateMatched=templateMatched, preConvMode=preConvMode)
876 results.correctedExposure = results.exposure
878 # Make sure masks of input image are propagated to diffim
879 def gm(exp):
880 return exp.mask
881 gm(results.correctedExposure)[:, :] = gm(subtractedExposure)
883 var = self.computeVarianceMean(results.correctedExposure)
884 self.log.info("Variance (corrected diffim): %f", var)
886 else:
887 config = self.config.decorrelateConfig
888 task = DecorrelateALKernelTask(config=config)
889 results = task.run(scienceExposure, templateExposure,
890 subtractedExposure, psfMatchingKernel, preConvKernel=preConvKernel,
891 templateMatched=templateMatched, preConvMode=preConvMode)
893 return results