Coverage for python / lsst / ip / diffim / dcrModel.py: 13%
267 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-22 08:51 +0000
1# This file is part of ip_diffim.
2#
3# LSST Data Management System
4# This product includes software developed by the
5# LSST Project (http://www.lsst.org/).
6# See COPYRIGHT file at the top of the source tree.
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
24from scipy import ndimage
25from lsst.afw.coord import differentialRefraction
26import lsst.afw.image as afwImage
27import lsst.geom as geom
28import lsst.pipe.base as pipeBase
30__all__ = ["DcrModel", "applyDcr", "calculateDcr", "calculateImageParallacticAngle"]
33class DcrModel:
34 """A model of the true sky after correcting chromatic effects.
36 Attributes
37 ----------
38 dcrNumSubfilters : `int`
39 Number of sub-filters used to model chromatic effects within a band.
40 modelImages : `list` of `lsst.afw.image.Image`
41 A list of masked images, each containing the model for one subfilter
43 Notes
44 -----
45 The ``DcrModel`` contains an estimate of the true sky, at a higher
46 wavelength resolution than the input observations. It can be forward-
47 modeled to produce Differential Chromatic Refraction (DCR) matched
48 templates for a given ``Exposure``, and provides utilities for conditioning
49 the model in ``dcrAssembleCoadd`` to avoid oscillating solutions between
50 iterations of forward modeling or between the subfilters of the model.
51 """
53 def __init__(self, modelImages, effectiveWavelength, bandwidth, filterLabel=None, psf=None,
54 bbox=None, wcs=None, mask=None, variance=None, photoCalib=None):
55 self.dcrNumSubfilters = len(modelImages)
56 self.modelImages = modelImages
57 self._filterLabel = filterLabel
58 self._effectiveWavelength = effectiveWavelength
59 self._bandwidth = bandwidth
60 self._psf = psf
61 self._bbox = bbox
62 self._wcs = wcs
63 self._mask = mask
64 self._variance = variance
65 self.photoCalib = photoCalib
67 @classmethod
68 def fromImage(cls, maskedImage, dcrNumSubfilters, effectiveWavelength, bandwidth,
69 wcs=None, filterLabel=None, psf=None, photoCalib=None):
70 """Initialize a DcrModel by dividing a coadd between the subfilters.
72 Parameters
73 ----------
74 maskedImage : `lsst.afw.image.MaskedImage`
75 Input coadded image to divide equally between the subfilters.
76 dcrNumSubfilters : `int`
77 Number of sub-filters used to model chromatic effects within a
78 band.
79 effectiveWavelength : `float`
80 The effective wavelengths of the current filter, in nanometers.
81 bandwidth : `float`
82 The bandwidth of the current filter, in nanometers.
83 wcs : `lsst.afw.geom.SkyWcs`
84 Coordinate system definition (wcs) for the exposure.
85 filterLabel : `lsst.afw.image.FilterLabel`, optional
86 The filter label, set in the current instruments' obs package.
87 Required for any calculation of DCR, including making matched
88 templates.
89 psf : `lsst.afw.detection.Psf`, optional
90 Point spread function (PSF) of the model.
91 Required if the ``DcrModel`` will be persisted.
92 photoCalib : `lsst.afw.image.PhotoCalib`, optional
93 Calibration to convert instrumental flux and
94 flux error to nanoJansky.
96 Returns
97 -------
98 dcrModel : `lsst.pipe.tasks.DcrModel`
99 Best fit model of the true sky after correcting chromatic effects.
100 """
101 # NANs will potentially contaminate the entire image,
102 # depending on the shift or convolution type used.
103 model = maskedImage.image.clone()
104 mask = maskedImage.mask.clone()
105 bbox = maskedImage.getBBox()
106 # We divide the variance by N and not N**2 because we will assume each
107 # subfilter is independent. That means that the significance of
108 # detected sources will be lower by a factor of sqrt(N) in the
109 # subfilter images, but we will recover it when we combine the
110 # subfilter images to construct matched templates.
111 variance = maskedImage.variance.clone()
112 variance /= dcrNumSubfilters
113 model /= dcrNumSubfilters
114 modelImages = [model, ]
115 for subfilter in range(1, dcrNumSubfilters):
116 modelImages.append(model.clone())
117 return cls(modelImages, effectiveWavelength, bandwidth,
118 filterLabel=filterLabel, psf=psf, bbox=bbox, wcs=wcs,
119 mask=mask, variance=variance, photoCalib=photoCalib)
121 @classmethod
122 def fromQuantum(cls, availableCoaddRefs, effectiveWavelength, bandwidth, numSubfilters):
123 """Load an existing DcrModel from a Gen 3 repository.
125 Parameters
126 ----------
127 availableCoaddRefs : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
128 Dictionary of spatially relevant retrieved coadd patches,
129 indexed by their sequential patch number.
130 effectiveWavelength : `float`
131 The effective wavelengths of the current filter, in nanometers.
132 bandwidth : `float`
133 The bandwidth of the current filter, in nanometers.
134 numSubfilters : `int`
135 Number of subfilters in the DcrCoadd.
137 Returns
138 -------
139 dcrModel : `lsst.pipe.tasks.DcrModel`
140 Best fit model of the true sky after correcting chromatic effects.
141 """
142 filterLabel = None
143 psf = None
144 bbox = None
145 wcs = None
146 mask = None
147 variance = None
148 photoCalib = None
149 modelImages = [None]*numSubfilters
151 for coaddRef in availableCoaddRefs:
152 subfilter = coaddRef.dataId["subfilter"]
153 dcrCoadd = coaddRef.get()
154 if filterLabel is None:
155 filterLabel = dcrCoadd.getFilter()
156 if psf is None:
157 psf = dcrCoadd.getPsf()
158 if bbox is None:
159 bbox = dcrCoadd.getBBox()
160 if wcs is None:
161 wcs = dcrCoadd.wcs
162 if mask is None:
163 mask = dcrCoadd.mask
164 if variance is None:
165 variance = dcrCoadd.variance
166 if photoCalib is None:
167 photoCalib = dcrCoadd.getPhotoCalib()
168 modelImages[subfilter] = dcrCoadd.image
169 return cls(modelImages, effectiveWavelength, bandwidth, filterLabel,
170 psf, bbox, wcs, mask, variance, photoCalib)
172 def __len__(self):
173 """Return the number of subfilters.
175 Returns
176 -------
177 dcrNumSubfilters : `int`
178 The number of DCR subfilters in the model.
179 """
180 return self.dcrNumSubfilters
182 def __getitem__(self, subfilter):
183 """Iterate over the subfilters of the DCR model.
185 Parameters
186 ----------
187 subfilter : `int`
188 Index of the current ``subfilter`` within the full band.
189 Negative indices are allowed, and count in reverse order
190 from the highest ``subfilter``.
192 Returns
193 -------
194 modelImage : `lsst.afw.image.Image`
195 The DCR model for the given ``subfilter``.
197 Raises
198 ------
199 IndexError
200 If the requested ``subfilter`` is greater or equal to the number
201 of subfilters in the model.
202 """
203 if np.abs(subfilter) >= len(self):
204 raise IndexError("subfilter out of bounds.")
205 return self.modelImages[subfilter]
207 def __setitem__(self, subfilter, maskedImage):
208 """Update the model image for one subfilter.
210 Parameters
211 ----------
212 subfilter : `int`
213 Index of the current subfilter within the full band.
214 maskedImage : `lsst.afw.image.Image`
215 The DCR model to set for the given ``subfilter``.
217 Raises
218 ------
219 IndexError
220 If the requested ``subfilter`` is greater or equal to the number
221 of subfilters in the model.
222 ValueError
223 If the bounding box of the new image does not match.
224 """
225 if np.abs(subfilter) >= len(self):
226 raise IndexError("subfilter out of bounds.")
227 if maskedImage.getBBox() != self.bbox:
228 raise ValueError("The bounding box of a subfilter must not change.")
229 self.modelImages[subfilter] = maskedImage
231 @property
232 def effectiveWavelength(self):
233 """Return the effective wavelength of the model.
235 Returns
236 -------
237 effectiveWavelength : `float`
238 The effective wavelength of the current filter, in nanometers.
239 """
240 return self._effectiveWavelength
242 @property
243 def filter(self):
244 """Return the filter label for the model.
246 Returns
247 -------
248 filterLabel : `lsst.afw.image.FilterLabel`
249 The filter used for the input observations.
250 """
251 return self._filterLabel
253 @property
254 def bandwidth(self):
255 """Return the bandwidth of the model.
257 Returns
258 -------
259 bandwidth : `float`
260 The bandwidth of the current filter, in nanometers.
261 """
262 return self._bandwidth
264 @property
265 def psf(self):
266 """Return the psf of the model.
268 Returns
269 -------
270 psf : `lsst.afw.detection.Psf`
271 Point spread function (PSF) of the model.
272 """
273 return self._psf
275 @property
276 def bbox(self):
277 """Return the common bounding box of each subfilter image.
279 Returns
280 -------
281 bbox : `lsst.afw.geom.Box2I`
282 Bounding box of the DCR model.
283 """
284 return self._bbox
286 @property
287 def wcs(self):
288 """Return the WCS of each subfilter image.
290 Returns
291 -------
292 wcs : `lsst.afw.geom.SkyWcs`
293 Coordinate system definition (wcs) for the exposure.
294 """
295 return self._wcs
297 @property
298 def mask(self):
299 """Return the common mask of each subfilter image.
301 Returns
302 -------
303 mask : `lsst.afw.image.Mask`
304 Mask plane of the DCR model.
305 """
306 return self._mask
308 @property
309 def variance(self):
310 """Return the common variance of each subfilter image.
312 Returns
313 -------
314 variance : `lsst.afw.image.Image`
315 Variance plane of the DCR model.
316 """
317 return self._variance
319 def getReferenceImage(self, bbox=None):
320 """Calculate a reference image from the average of the subfilter
321 images.
323 Parameters
324 ----------
325 bbox : `lsst.afw.geom.Box2I`, optional
326 Sub-region of the coadd. Returns the entire image if `None`.
328 Returns
329 -------
330 refImage : `numpy.ndarray`
331 The reference image with no chromatic effects applied.
332 """
333 bbox = bbox or self.bbox
334 return np.mean([model[bbox].array for model in self], axis=0)
336 def assign(self, dcrSubModel, bbox=None):
337 """Update a sub-region of the ``DcrModel`` with new values.
339 Parameters
340 ----------
341 dcrSubModel : `lsst.pipe.tasks.DcrModel`
342 New model of the true scene after correcting chromatic effects.
343 bbox : `lsst.afw.geom.Box2I`, optional
344 Sub-region of the coadd.
345 Defaults to the bounding box of ``dcrSubModel``.
347 Raises
348 ------
349 ValueError
350 If the new model has a different number of subfilters.
351 """
352 if len(dcrSubModel) != len(self):
353 raise ValueError("The number of DCR subfilters must be the same "
354 "between the old and new models.")
355 bbox = bbox or self.bbox
356 for model, subModel in zip(self, dcrSubModel):
357 model.assign(subModel[bbox], bbox)
359 def buildMatchedTemplate(self, exposure=None, order=3,
360 visitInfo=None, bbox=None, mask=None,
361 splitSubfilters=True, splitThreshold=0., amplifyModel=1.):
362 """Create a DCR-matched template image for an exposure.
364 Parameters
365 ----------
366 exposure : `lsst.afw.image.Exposure`, optional
367 The input exposure to build a matched template for.
368 May be omitted if all of the metadata is supplied separately
369 order : `int`, optional
370 Interpolation order of the DCR shift.
371 visitInfo : `lsst.afw.image.VisitInfo`, optional
372 Metadata for the exposure. Ignored if ``exposure`` is set.
373 bbox : `lsst.afw.geom.Box2I`, optional
374 Sub-region of the coadd, or use the entire coadd if not supplied.
375 mask : `lsst.afw.image.Mask`, optional
376 reference mask to use for the template image.
377 splitSubfilters : `bool`, optional
378 Calculate DCR for two evenly-spaced wavelengths in each subfilter,
379 instead of at the midpoint. Default: True
380 splitThreshold : `float`, optional
381 Minimum DCR difference within a subfilter required to use
382 ``splitSubfilters``
383 amplifyModel : `float`, optional
384 Multiplication factor to amplify differences between model planes.
385 Used to speed convergence of iterative forward modeling.
387 Returns
388 -------
389 templateImage : `lsst.afw.image.ImageF`
390 The DCR-matched template
392 Raises
393 ------
394 ValueError
395 If neither ``exposure`` or ``visitInfo`` are set.
396 """
397 if self.effectiveWavelength is None or self.bandwidth is None:
398 raise ValueError("'effectiveWavelength' and 'bandwidth' must be set for the DcrModel in order "
399 "to calculate DCR.")
400 if exposure is not None:
401 visitInfo = exposure.getInfo().getVisitInfo()
402 elif visitInfo is None:
403 raise ValueError("Either exposure or visitInfo must be set.")
404 if bbox is None:
405 bbox = self.bbox
406 else:
407 bbox = bbox.clippedTo(self.bbox)
408 dcrShift = calculateDcr(visitInfo, self.wcs, self.effectiveWavelength, self.bandwidth, len(self),
409 splitSubfilters=splitSubfilters, bbox=bbox)
410 templateImage = afwImage.ImageF(bbox)
411 refModel = None
412 for subfilter, dcr in enumerate(dcrShift):
413 if self[subfilter] is None:
414 # It is possible to load only a single DCR subfilter at a time.
415 self.log.debug("Skipping missing DCR model subfilter %d", subfilter)
416 continue
417 if amplifyModel > 1:
418 if refModel is None:
419 # amplifyModel is only an option while constructing the DcrModel,
420 # and we don't want to calculate a reference image during image differencing.
421 refModel = self.getReferenceImage(bbox)
422 model = (self[subfilter][bbox].array - refModel)*amplifyModel + refModel
423 else:
424 model = self[subfilter][bbox].array
425 templateImage.array += applyDcr(model, dcr, splitSubfilters=splitSubfilters,
426 splitThreshold=splitThreshold, order=order)
427 return templateImage
429 def buildMatchedExposure(self, exposure=None,
430 visitInfo=None, bbox=None, mask=None):
431 """Wrapper to create an exposure from a template image.
433 Parameters
434 ----------
435 exposure : `lsst.afw.image.Exposure`, optional
436 The input exposure to build a matched template for.
437 May be omitted if all of the metadata is supplied separately
438 visitInfo : `lsst.afw.image.VisitInfo`, optional
439 Metadata for the exposure. Ignored if ``exposure`` is set.
440 bbox : `lsst.afw.geom.Box2I`, optional
441 Sub-region of the coadd, or use the entire coadd if not supplied.
442 mask : `lsst.afw.image.Mask`, optional
443 reference mask to use for the template image.
445 Returns
446 -------
447 templateExposure : `lsst.afw.image.exposureF`
448 The DCR-matched template
450 Raises
451 ------
452 RuntimeError
453 If no `photcCalib` is set.
454 """
455 if bbox is None:
456 bbox = self.bbox
457 templateImage = self.buildMatchedTemplate(exposure=exposure, visitInfo=visitInfo,
458 bbox=bbox, mask=mask)
459 maskedImage = afwImage.MaskedImageF(bbox)
460 maskedImage.image = templateImage[bbox]
461 maskedImage.mask = self.mask[bbox]
462 maskedImage.variance = self.variance[bbox]
463 # The variance of the stacked image will be `dcrNumSubfilters`
464 # times the variance of the individual subfilters.
465 maskedImage.variance *= self.dcrNumSubfilters
466 templateExposure = afwImage.ExposureF(bbox, self.wcs)
467 templateExposure.setMaskedImage(maskedImage[bbox])
468 templateExposure.setPsf(self.psf)
469 templateExposure.setFilter(self.filter)
470 if self.photoCalib is None:
471 raise RuntimeError("No PhotoCalib set for the DcrModel. "
472 "If the DcrModel was created from a masked image"
473 " you must also specify the photoCalib.")
474 templateExposure.setPhotoCalib(self.photoCalib)
475 return templateExposure
477 def buildMatchedExposureHandle(self, exposure=None, visitInfo=None, bbox=None, mask=None):
478 """Create in-memory butler dataset reference containing the DCR-matched
479 template.
481 Parameters
482 ----------
483 exposure : `lsst.afw.image.Exposure`, optional
484 The input exposure to build a matched template for.
485 May be omitted if all of the metadata is supplied separately
486 visitInfo : `lsst.afw.image.VisitInfo`, optional
487 Metadata for the exposure. Ignored if ``exposure`` is set.
488 bbox : `lsst.afw.geom.Box2I`, optional
489 Sub-region of the coadd, or use the entire coadd if not supplied.
490 mask : `lsst.afw.image.Mask`, optional
491 reference mask to use for the template image.
493 Returns
494 -------
495 templateExposureHandle: `lsst.pipe.base.InMemoryDatasetHandle`
496 In-memory butler dataset reference containing the DCR-matched
497 template.
499 Raises
500 ------
501 ValueError
502 If no `exposure` or `visitInfo` is set.
503 """
504 if exposure is not None:
505 visitInfo = exposure.visitInfo
506 elif visitInfo is None:
507 raise ValueError("Either exposure or visitInfo must be set.")
508 templateExposure = self.buildMatchedExposure(
509 exposure=exposure, visitInfo=visitInfo, bbox=bbox, mask=mask
510 )
511 templateExposureHandle = pipeBase.InMemoryDatasetHandle(
512 templateExposure,
513 storageClass="ExposureF",
514 copy=False,
515 photoCalib=self.photoCalib
516 )
517 return templateExposureHandle
519 def conditionDcrModel(self, modelImages, bbox, gain=1.):
520 """Average two iterations' solutions to reduce oscillations.
522 Parameters
523 ----------
524 modelImages : `list` of `lsst.afw.image.Image`
525 The new DCR model images from the current iteration.
526 The values will be modified in place.
527 bbox : `lsst.afw.geom.Box2I`
528 Sub-region of the coadd
529 gain : `float`, optional
530 Relative weight to give the new solution when updating the model.
531 Defaults to 1.0, which gives equal weight to both solutions.
532 """
533 # Calculate weighted averages of the images.
534 for model, newModel in zip(self, modelImages):
535 newModel *= gain
536 newModel += model[bbox]
537 newModel /= 1. + gain
539 def regularizeModelIter(self, subfilter, newModel, bbox, regularizationFactor,
540 regularizationWidth=2):
541 """Restrict large variations in the model between iterations.
543 Parameters
544 ----------
545 subfilter : `int`
546 Index of the current subfilter within the full band.
547 newModel : `lsst.afw.image.Image`
548 The new DCR model for one subfilter from the current iteration.
549 Values in ``newModel`` that are extreme compared with the last
550 iteration are modified in place.
551 bbox : `lsst.afw.geom.Box2I`
552 Sub-region to coadd
553 regularizationFactor : `float`
554 Maximum relative change of the model allowed between iterations.
555 regularizationWidth : int, optional
556 Minimum radius of a region to include in regularization, in pixels.
557 """
558 refImage = self[subfilter][bbox].array
559 highThreshold = np.abs(refImage)*regularizationFactor
560 lowThreshold = refImage/regularizationFactor
561 newImage = newModel.array
562 self.applyImageThresholds(newImage, highThreshold=highThreshold, lowThreshold=lowThreshold,
563 regularizationWidth=regularizationWidth)
565 def regularizeModelFreq(self, modelImages, bbox, statsCtrl, regularizationFactor,
566 regularizationWidth=2, mask=None, convergenceMaskPlanes="DETECTED"):
567 """Restrict large variations in the model between subfilters.
569 Parameters
570 ----------
571 modelImages : `list` of `lsst.afw.image.Image`
572 The new DCR model images from the current iteration.
573 The values will be modified in place.
574 bbox : `lsst.afw.geom.Box2I`
575 Sub-region to coadd
576 statsCtrl : `lsst.afw.math.StatisticsControl`
577 Statistics control object for coaddition.
578 regularizationFactor : `float`
579 Maximum relative change of the model allowed between subfilters.
580 regularizationWidth : `int`, optional
581 Minimum radius of a region to include in regularization, in pixels.
582 mask : `lsst.afw.image.Mask`, optional
583 Optional alternate mask
584 convergenceMaskPlanes : `list` of `str`, or `str`, optional
585 Mask planes to use to calculate convergence.
587 Notes
588 -----
589 This implementation of frequency regularization restricts each
590 subfilter image to be a smoothly-varying function times a reference
591 image.
592 """
593 # ``regularizationFactor`` is the maximum change between subfilter
594 # images, so the maximum difference between one subfilter image and the
595 # average will be the square root of that.
596 maxDiff = np.sqrt(regularizationFactor)
597 noiseLevel = self.calculateNoiseCutoff(modelImages[0], statsCtrl, bufferSize=5, mask=mask, bbox=bbox)
598 referenceImage = self.getReferenceImage(bbox)
599 badPixels = np.isnan(referenceImage) | (referenceImage <= 0.)
600 if np.sum(~badPixels) == 0:
601 # Skip regularization if there are no valid pixels
602 return
603 referenceImage[badPixels] = 0.
604 filterWidth = regularizationWidth
605 fwhm = 2.*filterWidth
606 # The noise should be lower in the smoothed image by
607 # sqrt(Nsmooth) ~ fwhm pixels
608 noiseLevel /= fwhm
609 smoothRef = ndimage.gaussian_filter(referenceImage, filterWidth, mode='constant')
610 # Add a three sigma offset to both the reference and model to prevent
611 # dividing by zero. Note that this will also slightly suppress faint
612 # variations in color.
613 smoothRef += 3.*noiseLevel
615 lowThreshold = smoothRef/maxDiff
616 highThreshold = smoothRef*maxDiff
617 for model in modelImages:
618 self.applyImageThresholds(model.array,
619 highThreshold=highThreshold,
620 lowThreshold=lowThreshold,
621 regularizationWidth=regularizationWidth)
622 smoothModel = ndimage.gaussian_filter(model.array, filterWidth, mode='constant')
623 smoothModel += 3.*noiseLevel
624 relativeModel = smoothModel/smoothRef
625 # Now sharpen the smoothed relativeModel using an alpha of 3.
626 alpha = 3.
627 relativeModel2 = ndimage.gaussian_filter(relativeModel, filterWidth/alpha)
628 relativeModel += alpha*(relativeModel - relativeModel2)
629 model.array = relativeModel*referenceImage
631 def calculateNoiseCutoff(self, image, statsCtrl, bufferSize,
632 convergenceMaskPlanes="DETECTED", mask=None, bbox=None):
633 """Helper function to calculate the background noise level of an image.
635 Parameters
636 ----------
637 image : `lsst.afw.image.Image`
638 The input image to evaluate the background noise properties.
639 statsCtrl : `lsst.afw.math.StatisticsControl`
640 Statistics control object for coaddition.
641 bufferSize : `int`
642 Number of additional pixels to exclude
643 from the edges of the bounding box.
644 convergenceMaskPlanes : `list` of `str`, or `str`
645 Mask planes to use to calculate convergence.
646 mask : `lsst.afw.image.Mask`, Optional
647 Optional alternate mask
648 bbox : `lsst.afw.geom.Box2I`, optional
649 Sub-region of the masked image to calculate the noise level over.
651 Returns
652 -------
653 noiseCutoff : `float`
654 The threshold value to treat pixels as noise in an image..
655 """
656 if bbox is None:
657 bbox = self.bbox
658 if mask is None:
659 mask = self.mask[bbox]
660 bboxShrink = geom.Box2I(bbox)
661 bboxShrink.grow(-bufferSize)
662 convergeMask = mask.getPlaneBitMask(convergenceMaskPlanes)
664 backgroundPixels = mask[bboxShrink].array & (statsCtrl.getAndMask() | convergeMask) == 0
665 noiseCutoff = np.std(image[bboxShrink].array[backgroundPixels])
666 return noiseCutoff
668 def applyImageThresholds(self, image, highThreshold=None, lowThreshold=None, regularizationWidth=2):
669 """Restrict image values to be between upper and lower limits.
671 This method flags all pixels in an image that are outside of the given
672 threshold values. The threshold values are taken from a reference
673 image, so noisy pixels are likely to get flagged. In order to exclude
674 those noisy pixels, the array of flags is eroded and dilated, which
675 removes isolated pixels outside of the thresholds from the list of
676 pixels to be modified. Pixels that remain flagged after this operation
677 have their values set to the appropriate upper or lower threshold
678 value.
680 Parameters
681 ----------
682 image : `numpy.ndarray`
683 The image to apply the thresholds to.
684 The values will be modified in place.
685 highThreshold : `numpy.ndarray`, optional
686 Array of upper limit values for each pixel of ``image``.
687 lowThreshold : `numpy.ndarray`, optional
688 Array of lower limit values for each pixel of ``image``.
689 regularizationWidth : `int`, optional
690 Minimum radius of a region to include in regularization, in pixels.
691 """
692 # Generate the structure for binary erosion and dilation, which is used
693 # to remove noise-like pixels. Groups of pixels with a radius smaller
694 # than ``regularizationWidth`` will be excluded from regularization.
695 filterStructure = ndimage.iterate_structure(ndimage.generate_binary_structure(2, 1),
696 regularizationWidth)
697 if highThreshold is not None:
698 highPixels = image > highThreshold
699 if regularizationWidth > 0:
700 # Erode and dilate ``highPixels`` to exclude noisy pixels.
701 highPixels = ndimage.binary_opening(highPixels, structure=filterStructure)
702 image[highPixels] = highThreshold[highPixels]
703 if lowThreshold is not None:
704 lowPixels = image < lowThreshold
705 if regularizationWidth > 0:
706 # Erode and dilate ``lowPixels`` to exclude noisy pixels.
707 lowPixels = ndimage.binary_opening(lowPixels, structure=filterStructure)
708 image[lowPixels] = lowThreshold[lowPixels]
711def applyDcr(image, dcr, useInverse=False, splitSubfilters=False, splitThreshold=0.,
712 doPrefilter=True, order=3):
713 """Shift an image along the X and Y directions.
715 Parameters
716 ----------
717 image : `numpy.ndarray`
718 The input image to shift.
719 dcr : `tuple`
720 Shift calculated with ``calculateDcr``.
721 Uses numpy axes ordering (Y, X).
722 If ``splitSubfilters`` is set, each element is itself a `tuple`
723 of two `float`, corresponding to the DCR shift at the two wavelengths.
724 Otherwise, each element is a `float` corresponding to the DCR shift at
725 the effective wavelength of the subfilter.
726 useInverse : `bool`, optional
727 Apply the shift in the opposite direction. Default: False
728 splitSubfilters : `bool`, optional
729 Calculate DCR for two evenly-spaced wavelengths in each subfilter,
730 instead of at the midpoint. Default: False
731 splitThreshold : `float`, optional
732 Minimum DCR difference within a subfilter required to use
733 ``splitSubfilters``
734 doPrefilter : `bool`, optional
735 Spline filter the image before shifting, if set. Filtering is required,
736 so only set to False if the image is already filtered.
737 Filtering takes ~20% of the time of shifting, so if `applyDcr` will be
738 called repeatedly on the same image it is more efficient to
739 precalculate the filter.
740 order : `int`, optional
741 The order of the spline interpolation, default is 3.
743 Returns
744 -------
745 shiftedImage : `numpy.ndarray`
746 A copy of the input image with the specified shift applied.
747 """
748 if doPrefilter and order > 1:
749 prefilteredImage = ndimage.spline_filter(image, order=order)
750 else:
751 prefilteredImage = image
752 if splitSubfilters:
753 shiftAmp = np.max(np.abs([_dcr0 - _dcr1 for _dcr0, _dcr1 in zip(dcr[0], dcr[1])]))
754 if shiftAmp >= splitThreshold:
755 if useInverse:
756 shift = [-1.*s for s in dcr[0]]
757 shift1 = [-1.*s for s in dcr[1]]
758 else:
759 shift = dcr[0]
760 shift1 = dcr[1]
761 shiftedImage = ndimage.shift(prefilteredImage, shift, prefilter=False, order=order)
762 shiftedImage += ndimage.shift(prefilteredImage, shift1, prefilter=False, order=order)
763 shiftedImage /= 2.
764 return shiftedImage
765 else:
766 # If the difference in the DCR shifts is less than the threshold,
767 # then just use the average shift for efficiency.
768 dcr = (np.mean(dcr[0]), np.mean(dcr[1]))
769 if useInverse:
770 shift = [-1.*s for s in dcr]
771 else:
772 shift = dcr
773 shiftedImage = ndimage.shift(prefilteredImage, shift, prefilter=False, order=order)
774 return shiftedImage
777def calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, dcrNumSubfilters, splitSubfilters=False,
778 bbox=None):
779 """Calculate the shift in pixels of an exposure due to DCR.
781 Parameters
782 ----------
783 visitInfo : `lsst.afw.image.VisitInfo`
784 Metadata for the exposure.
785 wcs : `lsst.afw.geom.SkyWcs`
786 Coordinate system definition (wcs) for the exposure.
787 effectiveWavelength : `float`
788 The effective wavelengths of the current filter, in nanometers.
789 bandwidth : `float`
790 The bandwidth of the current filter, in nanometers.
791 dcrNumSubfilters : `int`
792 Number of sub-filters used to model chromatic effects within a band.
793 splitSubfilters : `bool`, optional
794 Calculate DCR for two evenly-spaced wavelengths in each subfilter,
795 instead of at the midpoint. Default: False
796 bbox : `lsst.afw.geom.Box2I`, optional
797 Bounding box for the region of interest for evaluating the local
798 pixelScale (defaults to the Sky Origin of the ``wcs`` provided if
799 ``bbox`` is None).
801 Returns
802 -------
803 dcrShift : `tuple` of two `float`
804 The 2D shift due to DCR, in pixels.
805 Uses numpy axes ordering (Y, X).
806 """
807 rotation = calculateImageParallacticAngle(visitInfo)
808 dcrShift = []
809 weight = [0.75, 0.25]
810 for wl0, wl1 in wavelengthGenerator(effectiveWavelength, bandwidth, dcrNumSubfilters):
811 # Note that diffRefractAmp can be negative, since it's relative to the
812 # midpoint of the full band
813 diffRefractAmp0 = differentialRefraction(wavelength=wl0, wavelengthRef=effectiveWavelength,
814 elevation=visitInfo.getBoresightAzAlt().getLatitude(),
815 observatory=visitInfo.getObservatory(),
816 weather=visitInfo.getWeather())
817 diffRefractAmp1 = differentialRefraction(wavelength=wl1, wavelengthRef=effectiveWavelength,
818 elevation=visitInfo.getBoresightAzAlt().getLatitude(),
819 observatory=visitInfo.getObservatory(),
820 weather=visitInfo.getWeather())
821 if bbox is not None:
822 pixelScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
823 else:
824 pixelScale = wcs.getPixelScale().asArcseconds()
826 if splitSubfilters:
827 diffRefractPix0 = diffRefractAmp0.asArcseconds()/pixelScale
828 diffRefractPix1 = diffRefractAmp1.asArcseconds()/pixelScale
829 diffRefractArr = [diffRefractPix0*weight[0] + diffRefractPix1*weight[1],
830 diffRefractPix0*weight[1] + diffRefractPix1*weight[0]]
831 shiftX = [diffRefractPix*np.sin(rotation.asRadians()) for diffRefractPix in diffRefractArr]
832 shiftY = [diffRefractPix*np.cos(rotation.asRadians()) for diffRefractPix in diffRefractArr]
833 dcrShift.append(((shiftY[0], shiftX[0]), (shiftY[1], shiftX[1])))
834 else:
835 diffRefractAmp = (diffRefractAmp0 + diffRefractAmp1)/2.
836 diffRefractPix = diffRefractAmp.asArcseconds()/pixelScale
837 shiftX = diffRefractPix*np.sin(rotation.asRadians())
838 shiftY = diffRefractPix*np.cos(rotation.asRadians())
839 dcrShift.append((shiftY, shiftX))
840 return dcrShift
843def calculateImageParallacticAngle(visitInfo, wcs=None):
844 """Calculate the total sky rotation angle of an exposure.
846 Parameters
847 ----------
848 visitInfo : `lsst.afw.image.VisitInfo`
849 Metadata for the exposure.
850 wcs : `lsst.afw.geom.SkyWcs`, optional
851 Coordinate system definition (wcs) for the exposure.
852 Deprecated, will be removed following v30.
854 Returns
855 -------
856 `lsst.geom.Angle`
857 The rotation of the image axis, East from North.
858 Equal to the parallactic angle plus any additional rotation of the
859 coordinate system.
860 A rotation angle of 0 degrees is defined with
861 North along the +y axis and East along the +x axis.
862 A rotation angle of 90 degrees is defined with
863 North along the +x axis and East along the -y axis.
864 """
865 rotAngle = visitInfo.boresightParAngle - visitInfo.boresightRotAngle
866 return rotAngle
869def wavelengthGenerator(effectiveWavelength, bandwidth, dcrNumSubfilters):
870 """Iterate over the wavelength endpoints of subfilters.
872 Parameters
873 ----------
874 effectiveWavelength : `float`
875 The effective wavelength of the current filter, in nanometers.
876 bandwidth : `float`
877 The bandwidth of the current filter, in nanometers.
878 dcrNumSubfilters : `int`
879 Number of sub-filters used to model chromatic effects within a band.
881 Yields
882 ------
883 `tuple` of two `float`
884 The next set of wavelength endpoints for a subfilter, in nanometers.
885 """
886 lambdaMin = effectiveWavelength - bandwidth/2
887 lambdaMax = effectiveWavelength + bandwidth/2
888 wlStep = bandwidth/dcrNumSubfilters
889 for wl in np.linspace(lambdaMin, lambdaMax, dcrNumSubfilters, endpoint=False):
890 yield (wl, wl + wlStep)