Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DM-33822 Fix WCS bug in DCR templates #208

Merged
merged 5 commits into from
Mar 1, 2022
Merged

DM-33822 Fix WCS bug in DCR templates #208

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
47264f2
Remove Gen 2 support
isullivan Feb 18, 2022
4860255
Support loading a single subfilter of the DCR model at a time
isullivan Feb 18, 2022
5fba208
DCR calculation should use model WCS not science.
isullivan Feb 26, 2022
04f06ee
Prefiltering is not allowed for order 1 spline
isullivan Feb 27, 2022
ede0eca
Fix filterLabel typo
isullivan Feb 28, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
135 changes: 53 additions & 82 deletions python/lsst/ip/diffim/dcrModel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,20 +50,22 @@ class DcrModel:
"""

def __init__(self, modelImages, effectiveWavelength, bandwidth, filterLabel=None, psf=None,
mask=None, variance=None, photoCalib=None):
bbox=None, wcs=None, mask=None, variance=None, photoCalib=None):
self.dcrNumSubfilters = len(modelImages)
self.modelImages = modelImages
self._filterLabel = filterLabel
self._effectiveWavelength = effectiveWavelength
self._bandwidth = bandwidth
self._psf = psf
self._bbox = bbox
self._wcs = wcs
self._mask = mask
self._variance = variance
self.photoCalib = photoCalib

@classmethod
def fromImage(cls, maskedImage, dcrNumSubfilters, effectiveWavelength, bandwidth,
filterLabel=None, psf=None, photoCalib=None):
wcs=None, filterLabel=None, psf=None, photoCalib=None):
"""Initialize a DcrModel by dividing a coadd between the subfilters.

Parameters
Expand All @@ -77,6 +79,8 @@ def fromImage(cls, maskedImage, dcrNumSubfilters, effectiveWavelength, bandwidth
The effective wavelengths of the current filter, in nanometers.
bandwidth : `float`
The bandwidth of the current filter, in nanometers.
wcs : `lsst.afw.geom.SkyWcs`
Coordinate system definition (wcs) for the exposure.
filterLabel : `lsst.afw.image.FilterLabel`, optional
The filter label, set in the current instruments' obs package.
Required for any calculation of DCR, including making matched
Expand All @@ -97,6 +101,7 @@ def fromImage(cls, maskedImage, dcrNumSubfilters, effectiveWavelength, bandwidth
# depending on the shift or convolution type used.
model = maskedImage.image.clone()
mask = maskedImage.mask.clone()
bbox = maskedImage.getBBox()
# We divide the variance by N and not N**2 because we will assume each
# subfilter is independent. That means that the significance of
# detected sources will be lower by a factor of sqrt(N) in the
Expand All @@ -109,61 +114,8 @@ def fromImage(cls, maskedImage, dcrNumSubfilters, effectiveWavelength, bandwidth
for subfilter in range(1, dcrNumSubfilters):
modelImages.append(model.clone())
return cls(modelImages, effectiveWavelength, bandwidth,
filterLabel=filterLabel, psf=psf, mask=mask, variance=variance, photoCalib=photoCalib)

@classmethod
def fromDataRef(cls, dataRef, effectiveWavelength, bandwidth, datasetType="dcrCoadd", numSubfilters=None,
**kwargs):
"""Load an existing DcrModel from a Gen 2 repository.

Parameters
----------
dataRef : `lsst.daf.persistence.ButlerDataRef`
Data reference defining the patch for coaddition and the
reference Warp
effectiveWavelength : `float`
The effective wavelengths of the current filter, in nanometers.
bandwidth : `float`
The bandwidth of the current filter, in nanometers.
datasetType : `str`, optional
Name of the DcrModel in the registry {"dcrCoadd", "dcrCoadd_sub"}
numSubfilters : `int`
Number of sub-filters used to model chromatic effects within a
band.
**kwargs
Additional keyword arguments to pass to look up the model in the
data registry.
Common keywords and their types include: ``tract``:`str`,
``patch``:`str`, ``bbox``:`lsst.afw.geom.Box2I`

Returns
-------
dcrModel : `lsst.pipe.tasks.DcrModel`
Best fit model of the true sky after correcting chromatic effects.
"""
modelImages = []
filterLabel = None
psf = None
mask = None
variance = None
photoCalib = None
if "subfilter" in kwargs:
kwargs.pop("subfilter")
for subfilter in range(numSubfilters):
dcrCoadd = dataRef.get(datasetType, subfilter=subfilter,
numSubfilters=numSubfilters, **kwargs)
if filterLabel is None:
filterLabel = dcrCoadd.getFilterLabel()
if psf is None:
psf = dcrCoadd.getPsf()
if mask is None:
mask = dcrCoadd.mask
if variance is None:
variance = dcrCoadd.variance
if photoCalib is None:
photoCalib = dcrCoadd.getPhotoCalib()
modelImages.append(dcrCoadd.image)
return cls(modelImages, effectiveWavelength, bandwidth, filterLabel, psf, mask, variance, photoCalib)
filterLabel=filterLabel, psf=psf, bbox=bbox, wcs=wcs,
mask=mask, variance=variance, photoCalib=photoCalib)

@classmethod
def fromQuantum(cls, availableCoaddRefs, effectiveWavelength, bandwidth):
Expand All @@ -186,6 +138,8 @@ def fromQuantum(cls, availableCoaddRefs, effectiveWavelength, bandwidth):
"""
filterLabel = None
psf = None
bbox = None
wcs = None
mask = None
variance = None
photoCalib = None
Expand All @@ -198,14 +152,19 @@ def fromQuantum(cls, availableCoaddRefs, effectiveWavelength, bandwidth):
filterLabel = dcrCoadd.getFilterLabel()
if psf is None:
psf = dcrCoadd.getPsf()
if bbox is None:
bbox = dcrCoadd.getBBox()
if wcs is None:
wcs = dcrCoadd.wcs
if mask is None:
mask = dcrCoadd.mask
if variance is None:
variance = dcrCoadd.variance
if photoCalib is None:
photoCalib = dcrCoadd.getPhotoCalib()
modelImages[subfilter] = dcrCoadd.image
return cls(modelImages, effectiveWavelength, bandwidth, filterLabel, psf, mask, variance, photoCalib)
return cls(modelImages, effectiveWavelength, bandwidth, filterLabel,
psf, bbox, wcs, mask, variance, photoCalib)

def __len__(self):
"""Return the number of subfilters.
Expand Down Expand Up @@ -319,7 +278,18 @@ def bbox(self):
bbox : `lsst.afw.geom.Box2I`
Bounding box of the DCR model.
"""
return self[0].getBBox()
return self._bbox

@property
def wcs(self):
"""Return the WCS of each subfilter image.

Returns
-------
bbox : `lsst.afw.geom.SkyWcs`
Coordinate system definition (wcs) for the exposure.
"""
return self._wcs

@property
def mask(self):
Expand Down Expand Up @@ -384,7 +354,7 @@ def assign(self, dcrSubModel, bbox=None):
model.assign(subModel[bbox], bbox)

def buildMatchedTemplate(self, exposure=None, order=3,
visitInfo=None, bbox=None, wcs=None, mask=None,
visitInfo=None, bbox=None, mask=None,
splitSubfilters=True, splitThreshold=0., amplifyModel=1.):
"""Create a DCR-matched template image for an exposure.

Expand All @@ -398,10 +368,7 @@ def buildMatchedTemplate(self, exposure=None, order=3,
visitInfo : `lsst.afw.image.VisitInfo`, optional
Metadata for the exposure. Ignored if ``exposure`` is set.
bbox : `lsst.afw.geom.Box2I`, optional
Sub-region of the coadd. Ignored if ``exposure`` is set.
wcs : `lsst.afw.geom.SkyWcs`, optional
Coordinate system definition (wcs) for the exposure.
Ignored if ``exposure`` is set.
Sub-region of the coadd, or use the entire coadd if not supplied.
mask : `lsst.afw.image.Mask`, optional
reference mask to use for the template image.
splitSubfilters : `bool`, optional
Expand All @@ -422,24 +389,31 @@ def buildMatchedTemplate(self, exposure=None, order=3,
Raises
------
ValueError
If neither ``exposure`` or all of ``visitInfo``, ``bbox``, and
``wcs`` are set.
If neither ``exposure`` or ``visitInfo`` are set.
"""
if self.effectiveWavelength is None or self.bandwidth is None:
raise ValueError("'effectiveWavelength' and 'bandwidth' must be set for the DcrModel in order "
"to calculate DCR.")
if exposure is not None:
visitInfo = exposure.getInfo().getVisitInfo()
bbox = exposure.getBBox()
wcs = exposure.getInfo().getWcs()
elif visitInfo is None or bbox is None or wcs is None:
raise ValueError("Either exposure or visitInfo, bbox, and wcs must be set.")
dcrShift = calculateDcr(visitInfo, wcs, self.effectiveWavelength, self.bandwidth, len(self),
elif visitInfo is None:
raise ValueError("Either exposure or visitInfo must be set.")
if bbox is None:
bbox = self.bbox
dcrShift = calculateDcr(visitInfo, self.wcs, self.effectiveWavelength, self.bandwidth, len(self),
splitSubfilters=splitSubfilters)
templateImage = afwImage.ImageF(bbox)
refModel = self.getReferenceImage(bbox)
refModel = None
for subfilter, dcr in enumerate(dcrShift):
if self[subfilter] is None:
# It is possible to load only a single DCR subfilter at a time.
self.log.debug("Skipping missing DCR model subfilter %d", subfilter)
continue
if amplifyModel > 1:
if refModel is None:
# amplifyModel is only an option while constructing the DcrModel,
# and we don't want to calculate a reference image during image differencing.
refModel = self.getReferenceImage(bbox)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not quite following the logic here. Before one would run self.getReferenceImage(bbox) once before the for loop. Now there is code to run it only the first time through the loop. What's the advantage? The only one I see is in case the first time through has a None subfilter, you would skip it. Am I missing something else?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The advantage is that getReferenceImage() does not have to be run during image differencing. amplifyModel is an option that is only used while building the DcrCoadds. I'll add a comment to make it more clear.

model = (self[subfilter][bbox].array - refModel)*amplifyModel + refModel
else:
model = self[subfilter][bbox].array
Expand All @@ -448,7 +422,7 @@ def buildMatchedTemplate(self, exposure=None, order=3,
return templateImage

def buildMatchedExposure(self, exposure=None,
visitInfo=None, bbox=None, wcs=None, mask=None):
visitInfo=None, bbox=None, mask=None):
"""Wrapper to create an exposure from a template image.

Parameters
Expand All @@ -459,10 +433,7 @@ def buildMatchedExposure(self, exposure=None,
visitInfo : `lsst.afw.image.VisitInfo`, optional
Metadata for the exposure. Ignored if ``exposure`` is set.
bbox : `lsst.afw.geom.Box2I`, optional
Sub-region of the coadd. Ignored if ``exposure`` is set.
wcs : `lsst.afw.geom.SkyWcs`, optional
Coordinate system definition (wcs) for the exposure.
Ignored if ``exposure`` is set.
Sub-region of the coadd, or use the entire coadd if not supplied.
mask : `lsst.afw.image.Mask`, optional
reference mask to use for the template image.

Expand All @@ -477,20 +448,20 @@ def buildMatchedExposure(self, exposure=None,
If no `photcCalib` is set.
"""
if bbox is None:
bbox = exposure.getBBox()
bbox = self.bbox
templateImage = self.buildMatchedTemplate(exposure=exposure, visitInfo=visitInfo,
bbox=bbox, wcs=wcs, mask=mask)
bbox=bbox, mask=mask)
maskedImage = afwImage.MaskedImageF(bbox)
maskedImage.image = templateImage[bbox]
maskedImage.mask = self.mask[bbox]
maskedImage.variance = self.variance[bbox]
# The variance of the stacked image will be `dcrNumSubfilters`
# times the variance of the individual subfilters.
maskedImage.variance *= self.dcrNumSubfilters
templateExposure = afwImage.ExposureF(bbox, wcs)
templateExposure = afwImage.ExposureF(bbox, self.wcs)
templateExposure.setMaskedImage(maskedImage[bbox])
templateExposure.setPsf(self.psf)
templateExposure.setFilterLabel(self.filterLabel)
templateExposure.setFilterLabel(self.filter)
if self.photoCalib is None:
raise RuntimeError("No PhotoCalib set for the DcrModel. "
"If the DcrModel was created from a masked image"
Expand Down Expand Up @@ -727,7 +698,7 @@ def applyDcr(image, dcr, useInverse=False, splitSubfilters=False, splitThreshold
shiftedImage : `numpy.ndarray`
A copy of the input image with the specified shift applied.
"""
if doPrefilter:
if doPrefilter and order > 1:
prefilteredImage = ndimage.spline_filter(image, order=order)
else:
prefilteredImage = image
Expand Down
1 change: 0 additions & 1 deletion python/lsst/ip/diffim/getTemplate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -350,7 +350,6 @@ def run(self, tractInfo, patchList, skyCorners, availableCoaddRefs,
dcrBBox.grow(-self.config.templateBorderSize)
dcrBBox.include(patchInnerBBox)
coaddPatch = dcrModel.buildMatchedExposure(bbox=dcrBBox,
wcs=coaddWcs,
visitInfo=visitInfo)
else:
if sensorRef is None:
Expand Down