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