LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
imageDecorrelation.py
Go to the documentation of this file.
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#
22
23import numpy as np
24
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
33
34from .imageMapReduce import (ImageMapReduceConfig, ImageMapReduceTask,
35 ImageMapper)
36from .utils import computeAveragePsf
37
38__all__ = ("DecorrelateALKernelTask", "DecorrelateALKernelConfig",
39 "DecorrelateALKernelMapper", "DecorrelateALKernelMapReduceConfig",
40 "DecorrelateALKernelSpatialConfig", "DecorrelateALKernelSpatialTask")
41
42
43class DecorrelateALKernelConfig(pexConfig.Config):
44 """Configuration parameters for the DecorrelateALKernelTask
45 """
46
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 )
58
59
60class DecorrelateALKernelTask(pipeBase.Task):
61 """Decorrelate the effect of convolution by Alard-Lupton matching kernel in image difference
62
63 """
64 ConfigClass = DecorrelateALKernelConfig
65 _DefaultName = "ip_diffim_decorrelateALKernel"
66
67 def __init__(self, *args, **kwargs):
68 """Create the image decorrelation Task
69
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)
78
80 self.statsControl.setNumSigmaClip(3.)
81 self.statsControl.setNumIter(3)
82 self.statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.ignoreMaskPlanes))
83
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
90
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.
96
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.
101
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.
140
141 Returns
142 -------
143 result : `lsst.pipe.base.Struct`
144 - ``correctedExposure`` : the decorrelated diffim
145
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.
153
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.
157
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.
163
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.
170
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).
173
174 We are also using a constant accross-the-image measure of sigma (sqrt(variance)) to compute
175 the decorrelation kernel.
176
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.")
183
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)
193
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)
201
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)
208
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, )
219
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)
250
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)
258
259 oldVarMean = self.computeVarianceMean(subtractedExposure)
260 self.log.info("Variance plane mean of uncorrected diffim: %f", oldVarMean)
261
262 kArr = kimg.array
263 diffimShape = subtractedExposure.image.array.shape
264 psfShape = psfImg.array.shape
265
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)
275
276 correctedImage = self.computeCorrectedImage(corr.corrft, subtractedExposure.image.array)
277 correctedPsf = self.computeCorrectedDiffimPsf(corr.corrft, psfImg.array)
278
279 # The subtracted exposure variance plane is already correlated, we cannot propagate
280 # it through another convolution; instead we need to use the uncorrelated originals
281 # The whitening should scale it to varianceMean + targetVarianceMean on average
282 if self.config.completeVarPlanePropagation:
283 self.log.debug("Using full variance plane calculation in decorrelation")
284 correctedVariance = self.calculateVariancePlane(
285 variance, targetVariance,
286 varianceMean, targetVarianceMean, corr.cnft, corr.crft)
287 else:
288 self.log.debug("Using estimated variance plane calculation in decorrelation")
289 correctedVariance = self.estimateVariancePlane(
290 variance, targetVariance,
291 corr.cnft, corr.crft)
292
293 # Determine the common shape
294 kSum = np.sum(kArr)
295 kSumSq = kSum*kSum
296 self.log.debug("Matching kernel sum: %.3e", kSum)
297 if not templateMatched:
298 # ImagePsfMatch.subtractExposures re-scales the difference in
299 # the science image convolution mode
300 correctedVariance /= kSumSq
301 subtractedExposure.image.array[...] = correctedImage # Allow for numpy type casting
302 subtractedExposure.variance.array[...] = correctedVariance
303 subtractedExposure.setPsf(correctedPsf)
304
305 newVarMean = self.computeVarianceMean(subtractedExposure)
306 self.log.info("Variance plane mean of corrected diffim: %.5e", newVarMean)
307
308 # TODO DM-23857 As part of the spatially varying correction implementation
309 # consider whether returning a Struct is still necessary.
310 return pipeBase.Struct(correctedExposure=subtractedExposure, )
311
312 def computeCommonShape(self, *shapes):
313 """Calculate the common shape for FFT operations. Set `self.freqSpaceShape`
314 internally.
315
316 Parameters
317 ----------
318 shapes : one or more `tuple` of `int`
319 Shapes of the arrays. All must have the same dimensionality.
320 At least one shape must be provided.
321
322 Returns
323 -------
324 None.
325
326 Notes
327 -----
328 For each dimension, gets the smallest even number greater than or equal to
329 `N1+N2-1` where `N1` and `N2` are the two largest values.
330 In case of only one shape given, rounds up to even each dimension value.
331 """
332 S = np.array(shapes, dtype=int)
333 if len(shapes) > 2:
334 S.sort(axis=0)
335 S = S[-2:]
336 if len(shapes) > 1:
337 commonShape = np.sum(S, axis=0) - 1
338 else:
339 commonShape = S[0]
340 commonShape[commonShape % 2 != 0] += 1
341 self.freqSpaceShape = tuple(commonShape)
342 self.log.info("Common frequency space shape %s", self.freqSpaceShape)
343
344 @staticmethod
345 def padCenterOriginArray(A, newShape: tuple, useInverse=False):
346 """Zero pad an image where the origin is at the center and replace the
347 origin to the corner as required by the periodic input of FFT. Implement also
348 the inverse operation, crop the padding and re-center data.
349
350 Parameters
351 ----------
352 A : `numpy.ndarray`
353 An array to copy from.
354 newShape : `tuple` of `int`
355 The dimensions of the resulting array. For padding, the resulting array
356 must be larger than A in each dimension. For the inverse operation this
357 must be the original, before padding size of the array.
358 useInverse : bool, optional
359 Selector of forward, add padding, operation (False)
360 or its inverse, crop padding, operation (True).
361
362 Returns
363 -------
364 R : `numpy.ndarray`
365 The padded or unpadded array with shape of `newShape` and the same dtype as A.
366
367 Notes
368 -----
369 For odd dimensions, the splitting is rounded to
370 put the center pixel into the new corner origin (0,0). This is to be consistent
371 e.g. for a dirac delta kernel that is originally located at the center pixel.
372 """
373
374 # The forward and inverse operations should round odd dimension halves at the opposite
375 # sides to get the pixels back to their original positions.
376 if not useInverse:
377 # Forward operation: First and second halves with respect to the axes of A.
378 firstHalves = [x//2 for x in A.shape]
379 secondHalves = [x-y for x, y in zip(A.shape, firstHalves)]
380 else:
381 # Inverse operation: Opposite rounding
382 secondHalves = [x//2 for x in newShape]
383 firstHalves = [x-y for x, y in zip(newShape, secondHalves)]
384
385 R = np.zeros_like(A, shape=newShape)
386 R[-firstHalves[0]:, -firstHalves[1]:] = A[:firstHalves[0], :firstHalves[1]]
387 R[:secondHalves[0], -firstHalves[1]:] = A[-secondHalves[0]:, :firstHalves[1]]
388 R[:secondHalves[0], :secondHalves[1]] = A[-secondHalves[0]:, -secondHalves[1]:]
389 R[-firstHalves[0]:, :secondHalves[1]] = A[:firstHalves[0], -secondHalves[1]:]
390 return R
391
392 def computeDiffimCorrection(self, kappa, svar, tvar):
393 """Compute the Lupton decorrelation post-convolution kernel for decorrelating an
394 image difference, based on the PSF-matching kernel.
395
396 Parameters
397 ----------
398 kappa : `numpy.ndarray` of `float`
399 A matching kernel 2-d numpy.array derived from Alard & Lupton PSF matching.
400 svar : `float` > 0.
401 Average variance of science image used for PSF matching.
402 tvar : `float` > 0.
403 Average variance of the template (matched) image used for PSF matching.
404
405 Returns
406 -------
407 corrft : `numpy.ndarray` of `float`
408 The frequency space representation of the correction. The array is real (dtype float).
409 Shape is `self.freqSpaceShape`.
410
411 cnft, crft : `numpy.ndarray` of `complex`
412 The overall convolution (pre-conv, PSF matching, noise correction) kernel
413 for the science and template images, respectively for the variance plane
414 calculations. These are intermediate results in frequency space.
415
416 Notes
417 -----
418 The maximum correction factor converges to `sqrt(tvar/svar)` towards high frequencies.
419 This should be a plausible value.
420 """
421 kSum = np.sum(kappa) # We scale the decorrelation to preserve fluxes
422 kappa = self.padCenterOriginArray(kappa, self.freqSpaceShape)
423 kft = np.fft.fft2(kappa)
424 kftAbsSq = np.real(np.conj(kft) * kft)
425
426 denom = svar + tvar * kftAbsSq
427 corrft = np.sqrt((svar + tvar * kSum*kSum) / denom)
428 cnft = corrft
429 crft = kft*corrft
430 return pipeBase.Struct(corrft=corrft, cnft=cnft, crft=crft)
431
432 def computeScoreCorrection(self, kappa, svar, tvar, preConvArr):
433 """Compute the correction kernel for a score image.
434
435 Parameters
436 ----------
437 kappa : `numpy.ndarray`
438 A matching kernel 2-d numpy.array derived from Alard & Lupton PSF matching.
439 svar : `float`
440 Average variance of science image used for PSF matching (before pre-convolution).
441 tvar : `float`
442 Average variance of the template (matched) image used for PSF matching.
443 preConvArr : `numpy.ndarray`
444 The pre-convolution kernel of the science image. It should be the PSF
445 of the science image or an approximation of it. It must be normed to sum 1.
446
447 Returns
448 -------
449 corrft : `numpy.ndarray` of `float`
450 The frequency space representation of the correction. The array is real (dtype float).
451 Shape is `self.freqSpaceShape`.
452 cnft, crft : `numpy.ndarray` of `complex`
453 The overall convolution (pre-conv, PSF matching, noise correction) kernel
454 for the science and template images, respectively for the variance plane
455 calculations. These are intermediate results in frequency space.
456
457 Notes
458 -----
459 To be precise, the science image should be _correlated_ by ``preConvArray`` but this
460 does not matter for this calculation.
461
462 ``cnft``, ``crft`` contain the scaling factor as well.
463
464 """
465 kSum = np.sum(kappa)
466 kappa = self.padCenterOriginArray(kappa, self.freqSpaceShape)
467 kft = np.fft.fft2(kappa)
468 preConvArr = self.padCenterOriginArray(preConvArr, self.freqSpaceShape)
469 preFt = np.fft.fft2(preConvArr)
470 preFtAbsSq = np.real(np.conj(preFt) * preFt)
471 kftAbsSq = np.real(np.conj(kft) * kft)
472 # Avoid zero division, though we don't normally approach `tiny`.
473 # We have numerical noise instead.
474 tiny = np.finfo(preFtAbsSq.dtype).tiny * 1000.
475 flt = preFtAbsSq < tiny
476 # If we pre-convolve to avoid deconvolution in AL, then kftAbsSq / preFtAbsSq
477 # theoretically expected to diverge to +inf. But we don't care about the convergence
478 # properties here, S goes to 0 at these frequencies anyway.
479 preFtAbsSq[flt] = tiny
480 denom = svar + tvar * kftAbsSq / preFtAbsSq
481 corrft = (svar + tvar * kSum*kSum) / denom
482 cnft = np.conj(preFt)*corrft
483 crft = kft*corrft
484 return pipeBase.Struct(corrft=corrft, cnft=cnft, crft=crft)
485
486 @staticmethod
487 def estimateVariancePlane(vplane1, vplane2, c1ft, c2ft):
488 """Estimate the variance planes.
489
490 The estimation assumes that around each pixel the surrounding
491 pixels' sigmas within the convolution kernel are the same.
492
493 Parameters
494 ----------
495 vplane1, vplane2 : `numpy.ndarray` of `float`
496 Variance planes of the original (before pre-convolution or matching)
497 exposures.
498 c1ft, c2ft : `numpy.ndarray` of `complex`
499 The overall convolution that includes the matching and the
500 afterburner in frequency space. The result of either
501 ``computeScoreCorrection`` or ``computeDiffimCorrection``.
502
503 Returns
504 -------
505 vplaneD : `numpy.ndarray` of `float`
506 The estimated variance plane of the difference/score image
507 as a weighted sum of the input variances.
508
509 Notes
510 -----
511 See DMTN-179 Section 5 about the variance plane calculations.
512 """
513 w1 = np.sum(np.real(np.conj(c1ft)*c1ft)) / c1ft.size
514 w2 = np.sum(np.real(np.conj(c2ft)*c2ft)) / c2ft.size
515 # w1, w2: the frequency space sum of abs(c1)^2 is the same as in image
516 # space.
517 return vplane1*w1 + vplane2*w2
518
519 def calculateVariancePlane(self, vplane1, vplane2, varMean1, varMean2, c1ft, c2ft):
520 """Full propagation of the variance planes of the original exposures.
521
522 The original variance planes of independent pixels are convolved with the
523 image space square of the overall kernels.
524
525 Parameters
526 ----------
527 vplane1, vplane2 : `numpy.ndarray` of `float`
528 Variance planes of the original (before pre-convolution or matching)
529 exposures.
530 varMean1, varMean2 : `float`
531 Replacement average values for non-finite ``vplane1`` and ``vplane2`` values respectively.
532
533 c1ft, c2ft : `numpy.ndarray` of `complex`
534 The overall convolution that includes the matching and the
535 afterburner in frequency space. The result of either
536 ``computeScoreCorrection`` or ``computeDiffimCorrection``.
537
538 Returns
539 -------
540 vplaneD : `numpy.ndarray` of `float`
541 The variance plane of the difference/score images.
542
543 Notes
544 -----
545 See DMTN-179 Section 5 about the variance plane calculations.
546
547 Infs and NaNs are allowed and kept in the returned array.
548 """
549 D = np.real(np.fft.ifft2(c1ft))
550 c1SqFt = np.fft.fft2(D*D)
551
552 v1shape = vplane1.shape
553 filtInf = np.isinf(vplane1)
554 filtNan = np.isnan(vplane1)
555 # This copy could be eliminated if inf/nan handling were go into padCenterOriginArray
556 vplane1 = np.copy(vplane1)
557 vplane1[filtInf | filtNan] = varMean1
558 D = self.padCenterOriginArray(vplane1, self.freqSpaceShape)
559 v1 = np.real(np.fft.ifft2(np.fft.fft2(D) * c1SqFt))
560 v1 = self.padCenterOriginArray(v1, v1shape, useInverse=True)
561 v1[filtNan] = np.nan
562 v1[filtInf] = np.inf
563
564 D = np.real(np.fft.ifft2(c2ft))
565 c2SqFt = np.fft.fft2(D*D)
566
567 v2shape = vplane2.shape
568 filtInf = np.isinf(vplane2)
569 filtNan = np.isnan(vplane2)
570 vplane2 = np.copy(vplane2)
571 vplane2[filtInf | filtNan] = varMean2
572 D = self.padCenterOriginArray(vplane2, self.freqSpaceShape)
573 v2 = np.real(np.fft.ifft2(np.fft.fft2(D) * c2SqFt))
574 v2 = self.padCenterOriginArray(v2, v2shape, useInverse=True)
575 v2[filtNan] = np.nan
576 v2[filtInf] = np.inf
577
578 return v1 + v2
579
580 def computeCorrectedDiffimPsf(self, corrft, psfOld):
581 """Compute the (decorrelated) difference image's new PSF.
582
583 Parameters
584 ----------
585 corrft : `numpy.ndarray`
586 The frequency space representation of the correction calculated by
587 `computeCorrection`. Shape must be `self.freqSpaceShape`.
588 psfOld : `numpy.ndarray`
589 The psf of the difference image to be corrected.
590
591 Returns
592 -------
593 correctedPsf : `lsst.meas.algorithms.KernelPsf`
594 The corrected psf, same shape as `psfOld`, sum normed to 1.
595
596 Notes
597 -----
598 There is no algorithmic guarantee that the corrected psf can
599 meaningfully fit to the same size as the original one.
600 """
601 psfShape = psfOld.shape
602 psfNew = self.padCenterOriginArray(psfOld, self.freqSpaceShape)
603 psfNew = np.fft.fft2(psfNew)
604 psfNew *= corrft
605 psfNew = np.fft.ifft2(psfNew)
606 psfNew = psfNew.real
607 psfNew = self.padCenterOriginArray(psfNew, psfShape, useInverse=True)
608 psfNew = psfNew/psfNew.sum()
609
610 psfcI = afwImage.ImageD(geom.Extent2I(psfShape[1], psfShape[0]))
611 psfcI.array = psfNew
612 psfcK = afwMath.FixedKernel(psfcI)
613 correctedPsf = measAlg.KernelPsf(psfcK)
614 return correctedPsf
615
616 def computeCorrectedImage(self, corrft, imgOld):
617 """Compute the decorrelated difference image.
618
619 Parameters
620 ----------
621 corrft : `numpy.ndarray`
622 The frequency space representation of the correction calculated by
623 `computeCorrection`. Shape must be `self.freqSpaceShape`.
624 imgOld : `numpy.ndarray`
625 The difference image to be corrected.
626
627 Returns
628 -------
629 imgNew : `numpy.ndarray`
630 The corrected image, same size as the input.
631 """
632 expShape = imgOld.shape
633 imgNew = np.copy(imgOld)
634 filtInf = np.isinf(imgNew)
635 filtNan = np.isnan(imgNew)
636 imgNew[filtInf] = np.nan
637 imgNew[filtInf | filtNan] = np.nanmean(imgNew)
638 imgNew = self.padCenterOriginArray(imgNew, self.freqSpaceShape)
639 imgNew = np.fft.fft2(imgNew)
640 imgNew *= corrft
641 imgNew = np.fft.ifft2(imgNew)
642 imgNew = imgNew.real
643 imgNew = self.padCenterOriginArray(imgNew, expShape, useInverse=True)
644 imgNew[filtNan] = np.nan
645 imgNew[filtInf] = np.inf
646 return imgNew
647
648
650 """Task to be used as an ImageMapper for performing
651 A&L decorrelation on subimages on a grid across a A&L difference image.
652
653 This task subclasses DecorrelateALKernelTask in order to implement
654 all of that task's configuration parameters, as well as its `run` method.
655 """
656
657 ConfigClass = DecorrelateALKernelConfig
658 _DefaultName = 'ip_diffim_decorrelateALKernelMapper'
659
660 def __init__(self, *args, **kwargs):
661 DecorrelateALKernelTask.__init__(self, *args, **kwargs)
662
663 def run(self, subExposure, expandedSubExposure, fullBBox,
664 template, science, alTaskResult=None, psfMatchingKernel=None,
665 preConvKernel=None, **kwargs):
666 """Perform decorrelation operation on `subExposure`, using
667 `expandedSubExposure` to allow for invalid edge pixels arising from
668 convolutions.
669
670 This method performs A&L decorrelation on `subExposure` using
671 local measures for image variances and PSF. `subExposure` is a
672 sub-exposure of the non-decorrelated A&L diffim. It also
673 requires the corresponding sub-exposures of the template
674 (`template`) and science (`science`) exposures.
675
676 Parameters
677 ----------
678 subExposure : `lsst.afw.image.Exposure`
679 the sub-exposure of the diffim
680 expandedSubExposure : `lsst.afw.image.Exposure`
681 the expanded sub-exposure upon which to operate
682 fullBBox : `lsst.geom.Box2I`
683 the bounding box of the original exposure
684 template : `lsst.afw.image.Exposure`
685 the corresponding sub-exposure of the template exposure
686 science : `lsst.afw.image.Exposure`
687 the corresponding sub-exposure of the science exposure
688 alTaskResult : `lsst.pipe.base.Struct`
689 the result of A&L image differencing on `science` and
690 `template`, importantly containing the resulting
691 `psfMatchingKernel`. Can be `None`, only if
692 `psfMatchingKernel` is not `None`.
693 psfMatchingKernel : Alternative parameter for passing the
694 A&L `psfMatchingKernel` directly.
695 preConvKernel : If not None, then pre-filtering was applied
696 to science exposure, and this is the pre-convolution
697 kernel.
698 kwargs :
699 additional keyword arguments propagated from
700 `ImageMapReduceTask.run`.
701
702 Returns
703 -------
704 A `pipeBase.Struct` containing:
705
706 - ``subExposure`` : the result of the `subExposure` processing.
707 - ``decorrelationKernel`` : the decorrelation kernel, currently
708 not used.
709
710 Notes
711 -----
712 This `run` method accepts parameters identical to those of
713 `ImageMapper.run`, since it is called from the
714 `ImageMapperTask`. See that class for more information.
715 """
716 templateExposure = template # input template
717 scienceExposure = science # input science image
718 if alTaskResult is None and psfMatchingKernel is None:
719 raise RuntimeError('Both alTaskResult and psfMatchingKernel cannot be None')
720 psfMatchingKernel = alTaskResult.psfMatchingKernel if alTaskResult is not None else psfMatchingKernel
721
722 # subExp and expandedSubExp are subimages of the (un-decorrelated) diffim!
723 # So here we compute corresponding subimages of templateExposure and scienceExposure
724 subExp2 = scienceExposure.Factory(scienceExposure, expandedSubExposure.getBBox())
725 subExp1 = templateExposure.Factory(templateExposure, expandedSubExposure.getBBox())
726
727 # Prevent too much log INFO verbosity from DecorrelateALKernelTask.run
728 logLevel = self.log.level
729 self.log.setLevel(self.log.WARNING)
730 res = DecorrelateALKernelTask.run(self, subExp2, subExp1, expandedSubExposure,
731 psfMatchingKernel, preConvKernel, **kwargs)
732 self.log.setLevel(logLevel) # reset the log level
733
734 diffim = res.correctedExposure.Factory(res.correctedExposure, subExposure.getBBox())
735 out = pipeBase.Struct(subExposure=diffim, )
736 return out
737
738
740 """Configuration parameters for the ImageMapReduceTask to direct it to use
741 DecorrelateALKernelMapper as its mapper for A&L decorrelation.
742 """
743 mapper = pexConfig.ConfigurableField(
744 doc='A&L decorrelation task to run on each sub-image',
745 target=DecorrelateALKernelMapper
746 )
747
748
749class DecorrelateALKernelSpatialConfig(pexConfig.Config):
750 """Configuration parameters for the DecorrelateALKernelSpatialTask.
751 """
752 decorrelateConfig = pexConfig.ConfigField(
753 dtype=DecorrelateALKernelConfig,
754 doc='DecorrelateALKernel config to use when running on complete exposure (non spatially-varying)',
755 )
756
757 decorrelateMapReduceConfig = pexConfig.ConfigField(
758 dtype=DecorrelateALKernelMapReduceConfig,
759 doc='DecorrelateALKernelMapReduce config to use when running on each sub-image (spatially-varying)',
760 )
761
762 ignoreMaskPlanes = pexConfig.ListField(
763 dtype=str,
764 doc="""Mask planes to ignore for sigma-clipped statistics""",
765 default=("INTRP", "EDGE", "DETECTED", "SAT", "CR", "BAD", "NO_DATA", "DETECTED_NEGATIVE")
766 )
767
768 def setDefaults(self):
769 self.decorrelateMapReduceConfig.gridStepX = self.decorrelateMapReduceConfig.gridStepY = 40
770 self.decorrelateMapReduceConfig.cellSizeX = self.decorrelateMapReduceConfig.cellSizeY = 41
771 self.decorrelateMapReduceConfig.borderSizeX = self.decorrelateMapReduceConfig.borderSizeY = 8
772 self.decorrelateMapReduceConfig.reducer.reduceOperation = 'average'
773
774
776 """Decorrelate the effect of convolution by Alard-Lupton matching kernel in image difference
777
778 """
779 ConfigClass = DecorrelateALKernelSpatialConfig
780 _DefaultName = "ip_diffim_decorrelateALKernelSpatial"
781
782 def __init__(self, *args, **kwargs):
783 """Create the image decorrelation Task
784
785 Parameters
786 ----------
787 args :
788 arguments to be passed to
789 `lsst.pipe.base.task.Task.__init__`
790 kwargs :
791 additional keyword arguments to be passed to
792 `lsst.pipe.base.task.Task.__init__`
793 """
794 pipeBase.Task.__init__(self, *args, **kwargs)
795
797 self.statsControl.setNumSigmaClip(3.)
798 self.statsControl.setNumIter(3)
799 self.statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.ignoreMaskPlanes))
800
801 def computeVarianceMean(self, exposure):
802 """Compute the mean of the variance plane of `exposure`.
803 """
804 statObj = afwMath.makeStatistics(exposure.variance,
805 exposure.mask,
806 afwMath.MEANCLIP, self.statsControl)
807 var = statObj.getValue(afwMath.MEANCLIP)
808 return var
809
810 def run(self, scienceExposure, templateExposure, subtractedExposure, psfMatchingKernel,
811 spatiallyVarying=True, preConvKernel=None, templateMatched=True, preConvMode=False):
812 """Perform decorrelation of an image difference exposure.
813
814 Decorrelates the diffim due to the convolution of the
815 templateExposure with the A&L psfMatchingKernel. If
816 `spatiallyVarying` is True, it utilizes the spatially varying
817 matching kernel via the `imageMapReduce` framework to perform
818 spatially-varying decorrelation on a grid of subExposures.
819
820 Parameters
821 ----------
822 scienceExposure : `lsst.afw.image.Exposure`
823 the science Exposure used for PSF matching
824 templateExposure : `lsst.afw.image.Exposure`
825 the template Exposure used for PSF matching
826 subtractedExposure : `lsst.afw.image.Exposure`
827 the subtracted Exposure produced by `ip_diffim.ImagePsfMatchTask.subtractExposures()`
828 psfMatchingKernel : an (optionally spatially-varying) PSF matching kernel produced
829 by `ip_diffim.ImagePsfMatchTask.subtractExposures()`
830 spatiallyVarying : `bool`
831 if True, perform the spatially-varying operation
832 preConvKernel : `lsst.meas.algorithms.Psf`
833 if not none, the scienceExposure has been pre-filtered with this kernel. (Currently
834 this option is experimental.)
835 templateMatched : `bool`, optional
836 If True, the template exposure was matched (convolved) to the science exposure.
837 preConvMode : `bool`, optional
838 If True, ``subtractedExposure`` is assumed to be a likelihood difference image
839 and will be noise corrected as a likelihood image.
840
841 Returns
842 -------
843 results : `lsst.pipe.base.Struct`
844 a structure containing:
845 - ``correctedExposure`` : the decorrelated diffim
846 """
847 self.log.info('Running A&L decorrelation: spatiallyVarying=%r', spatiallyVarying)
848
849 svar = self.computeVarianceMean(scienceExposure)
850 tvar = self.computeVarianceMean(templateExposure)
851 if np.isnan(svar) or np.isnan(tvar): # Should not happen unless entire image has been masked.
852 # Double check that one of the exposures is all NaNs
853 if (np.all(np.isnan(scienceExposure.image.array))
854 or np.all(np.isnan(templateExposure.image.array))):
855 self.log.warning('Template or science image is entirely NaNs: skipping decorrelation.')
856 if np.isnan(svar):
857 svar = 1e-9
858 if np.isnan(tvar):
859 tvar = 1e-9
860
861 var = self.computeVarianceMean(subtractedExposure)
862
863 if spatiallyVarying:
864 self.log.info("Variance (science, template): (%f, %f)", svar, tvar)
865 self.log.info("Variance (uncorrected diffim): %f", var)
866 config = self.config.decorrelateMapReduceConfig
867 task = ImageMapReduceTask(config=config)
868 results = task.run(subtractedExposure, science=scienceExposure,
869 template=templateExposure, psfMatchingKernel=psfMatchingKernel,
870 preConvKernel=preConvKernel, forceEvenSized=True,
871 templateMatched=templateMatched, preConvMode=preConvMode)
872 results.correctedExposure = results.exposure
873
874 # Make sure masks of input image are propagated to diffim
875 def gm(exp):
876 return exp.mask
877 gm(results.correctedExposure)[:, :] = gm(subtractedExposure)
878
879 var = self.computeVarianceMean(results.correctedExposure)
880 self.log.info("Variance (corrected diffim): %f", var)
881
882 else:
883 config = self.config.decorrelateConfig
884 task = DecorrelateALKernelTask(config=config)
885 results = task.run(scienceExposure, templateExposure,
886 subtractedExposure, psfMatchingKernel, preConvKernel=preConvKernel,
887 templateMatched=templateMatched, preConvMode=preConvMode)
888
889 return results
A kernel created from an Image.
Definition Kernel.h:471
Pass parameters to a Statistics object.
Definition Statistics.h:83
run(self, subExposure, expandedSubExposure, fullBBox, template, science, alTaskResult=None, psfMatchingKernel=None, preConvKernel=None, **kwargs)
run(self, scienceExposure, templateExposure, subtractedExposure, psfMatchingKernel, spatiallyVarying=True, preConvKernel=None, templateMatched=True, preConvMode=False)
run(self, scienceExposure, templateExposure, subtractedExposure, psfMatchingKernel, preConvKernel=None, xcen=None, ycen=None, svar=None, tvar=None, templateMatched=True, preConvMode=False, **kwargs)
calculateVariancePlane(self, vplane1, vplane2, varMean1, varMean2, c1ft, c2ft)
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361