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

Tickets/DM-16864: New frequency regularization using relative model #107

Merged
merged 2 commits into from
Feb 1, 2019
Merged

Tickets/DM-16864: New frequency regularization using relative model #107

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8ac52ee
New frequency regularization using relative model
isullivan Dec 17, 2018
e1d653a
Bug fix: supply missing BBox
isullivan Jan 28, 2019
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
79 changes: 50 additions & 29 deletions python/lsst/ip/diffim/dcrModel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,19 +40,9 @@ class DcrModel:
----------
dcrNumSubfilters : `int`
Number of sub-filters used to model chromatic effects within a band.
filterInfo : `lsst.afw.image.Filter`
The filter definition, set in the current instruments' obs package.
modelImages : `list` of `lsst.afw.image.MaskedImage`
A list of masked images, each containing the model for one subfilter

Parameters
----------
modelImages : `list` of `lsst.afw.image.MaskedImage`
A list of masked images, each containing the model for one subfilter.
filterInfo : `lsst.afw.image.Filter`, optional
The filter definition, set in the current instruments' obs package.
Required for any calculation of DCR, including making matched templates.

Notes
-----
The ``DcrModel`` contains an estimate of the true sky, at a higher
Expand Down Expand Up @@ -250,7 +240,7 @@ def mask(self):
return self[0].mask

def getReferenceImage(self, bbox=None):
"""Create a simple template from the DCR model.
"""Calculate a reference image from the average of the subfilter images.

Parameters
----------
Expand All @@ -259,8 +249,8 @@ def getReferenceImage(self, bbox=None):

Returns
-------
templateImage : `numpy.ndarray`
The template with no chromatic effects applied.
refImage : `numpy.ndarray`
The reference image with no chromatic effects applied.
"""
bbox = bbox or self.bbox
return np.mean([model[bbox].image.array for model in self], axis=0)
Expand Down Expand Up @@ -426,11 +416,12 @@ def regularizeModelIter(self, subfilter, newModel, bbox, regularizationFactor,
refImage = self[subfilter][bbox].image.array
highThreshold = np.abs(refImage)*regularizationFactor
lowThreshold = refImage/regularizationFactor
self.applyImageThresholds(newModel, highThreshold=highThreshold, lowThreshold=lowThreshold,
newImage = newModel.image.array
self.applyImageThresholds(newImage, highThreshold=highThreshold, lowThreshold=lowThreshold,
regularizationWidth=regularizationWidth)

def regularizeModelFreq(self, modelImages, bbox, regularizationFactor,
regularizationWidth=2):
def regularizeModelFreq(self, modelImages, bbox, statsCtrl, regularizationFactor,
regularizationWidth=2, mask=None, convergenceMaskPlanes="DETECTED"):
"""Restrict large variations in the model between subfilters.

Parameters
Expand All @@ -440,21 +431,53 @@ def regularizeModelFreq(self, modelImages, bbox, regularizationFactor,
The values will be modified in place.
bbox : `lsst.afw.geom.Box2I`
Sub-region to coadd
statsCtrl : `lsst.afw.math.StatisticsControl`
Statistics control object for coaddition.
regularizationFactor : `float`
Maximum relative change of the model allowed between subfilters.
regularizationWidth : `int`, optional
Minimum radius of a region to include in regularization, in pixels.
mask : `lsst.afw.image.Mask`, optional
Optional alternate mask
convergenceMaskPlanes : `list` of `str`, or `str`, optional
Mask planes to use to calculate convergence.

Notes
-----
This implementation of frequency regularization restricts each subfilter
image to be a smoothly-varying function times a reference image.
"""
# ``regularizationFactor`` is the maximum change between subfilter images, so the maximum difference
# between one subfilter image and the average will be the square root of that.
maxDiff = np.sqrt(regularizationFactor)
refImage = self.getReferenceImage(bbox)

for model in modelImages:
highThreshold = np.abs(refImage)*maxDiff
lowThreshold = refImage/maxDiff
self.applyImageThresholds(model[bbox], highThreshold=highThreshold, lowThreshold=lowThreshold,
noiseLevel = self.calculateNoiseCutoff(modelImages[0], statsCtrl, bufferSize=5, mask=mask, bbox=bbox)
referenceImage = self.getReferenceImage(bbox)
badPixels = np.isnan(referenceImage) | (referenceImage <= 0.)
if np.sum(~badPixels) == 0:
# Skip regularization if there are no valid pixels
return
referenceImage[badPixels] = 0.
filterWidth = regularizationWidth
fwhm = 2.*filterWidth
# The noise should be lower in the smoothed image by sqrt(Nsmooth) ~ fwhm pixels
noiseLevel /= fwhm
smoothRef = ndimage.filters.gaussian_filter(referenceImage, filterWidth) + noiseLevel

baseThresh = np.ones_like(referenceImage)
highThreshold = baseThresh*maxDiff
lowThreshold = baseThresh/maxDiff
for subfilter, model in enumerate(modelImages):
smoothModel = ndimage.filters.gaussian_filter(model.image.array, filterWidth) + noiseLevel
relativeModel = smoothModel/smoothRef
# Now sharpen the smoothed relativeModel using an alpha of 3.
relativeModel2 = ndimage.filters.gaussian_filter(relativeModel, filterWidth/3.)
relativeModel = relativeModel + 3.*(relativeModel - relativeModel2)
self.applyImageThresholds(relativeModel,
highThreshold=highThreshold,
lowThreshold=lowThreshold,
regularizationWidth=regularizationWidth)
relativeModel *= referenceImage
modelImages[subfilter].image.array = relativeModel

def calculateNoiseCutoff(self, maskedImage, statsCtrl, bufferSize,
convergenceMaskPlanes="DETECTED", mask=None, bbox=None):
Expand Down Expand Up @@ -493,8 +516,7 @@ def calculateNoiseCutoff(self, maskedImage, statsCtrl, bufferSize,
noiseCutoff = np.std(maskedImage[bboxShrink].image.array[backgroundPixels])
return noiseCutoff

def applyImageThresholds(self, maskedImage, highThreshold=None, lowThreshold=None,
regularizationWidth=2):
def applyImageThresholds(self, image, highThreshold=None, lowThreshold=None, regularizationWidth=2):
"""Restrict image values to be between upper and lower limits.

This method flags all pixels in an image that are outside of the given
Expand All @@ -507,20 +529,19 @@ def applyImageThresholds(self, maskedImage, highThreshold=None, lowThreshold=Non

Parameters
----------
maskedImage : `lsst.afw.image.MaskedImage`
image : `numpy.ndarray`
The image to apply the thresholds to.
The image plane values will be modified in place.
The values will be modified in place.
highThreshold : `numpy.ndarray`, optional
Array of upper limit values for each pixel of ``maskedImage``.
Array of upper limit values for each pixel of ``image``.
lowThreshold : `numpy.ndarray`, optional
Array of lower limit values for each pixel of ``maskedImage``.
Array of lower limit values for each pixel of ``image``.
regularizationWidth : `int`, optional
Minimum radius of a region to include in regularization, in pixels.
"""
# Generate the structure for binary erosion and dilation, which is used to remove noise-like pixels.
# Groups of pixels with a radius smaller than ``regularizationWidth``
# will be excluded from regularization.
image = maskedImage.image.array
filterStructure = ndimage.iterate_structure(ndimage.generate_binary_structure(2, 1),
regularizationWidth)
if highThreshold is not None:
Expand Down
53 changes: 24 additions & 29 deletions tests/test_dcrModel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,22 @@ def makeTestImages(self, seed=5, nSrc=5, psfSize=2., noiseLevel=5.,
self.mask = modelImages[0].mask
return modelImages

def prepareStats(self):
"""Make a simple statistics object for testing.

Returns
-------
statsCtrl : `lsst.afw.math.StatisticsControl`
Statistics control object for coaddition.
"""
statsCtrl = afwMath.StatisticsControl()
statsCtrl.setNumSigmaClip(5)
statsCtrl.setNumIter(3)
statsCtrl.setNanSafe(True)
statsCtrl.setWeighted(True)
statsCtrl.setCalcErrorFromInputVariance(False)
return statsCtrl

def makeDummyWcs(self, rotAngle, pixelScale, crval, flipX=True):
"""Make a World Coordinate System object for testing.

Expand Down Expand Up @@ -377,30 +393,20 @@ def testConditionDcrModelWithChange(self):
refModel.image.array[:] *= 2.
self.assertMaskedImagesAlmostEqual(refModel, newModel)

def testRegularizationLargeClamp(self):
"""Frequency regularization should leave the models unchanged if the clamp factor is large.
"""
clampFrequency = 3.
regularizationWidth = 2
dcrModels = DcrModel(modelImages=self.makeTestImages())
newModels = [model.clone() for model in dcrModels]
dcrModels.regularizeModelFreq(newModels, self.bbox, clampFrequency, regularizationWidth)
for model, refModel in zip(newModels, dcrModels):
self.assertMaskedImagesEqual(model, refModel)

def testRegularizationSmallClamp(self):
"""Test that large variations between model planes are reduced.

This also tests that noise-like pixels are not regularized.
"""
clampFrequency = 1.1
clampFrequency = 2
regularizationWidth = 2
fluxRange = 10.
dcrModels = DcrModel(modelImages=self.makeTestImages(fluxRange=fluxRange))
newModels = [model.clone() for model in dcrModels]
templateImage = dcrModels.getReferenceImage(self.bbox)

dcrModels.regularizeModelFreq(newModels, self.bbox, clampFrequency, regularizationWidth)
statsCtrl = self.prepareStats()
dcrModels.regularizeModelFreq(newModels, self.bbox, statsCtrl, clampFrequency, regularizationWidth)
for model, refModel in zip(newModels, dcrModels):
# The mask and variance planes should be unchanged
self.assertFloatsEqual(model.mask.array, refModel.mask.array)
Expand Down Expand Up @@ -431,6 +437,7 @@ def testRegularizationSidelobes(self):
templateImage = np.mean([model.image.array for model in modelImages], axis=0)
sidelobeImages = self.makeTestImages(seed=5, nSrc=5, psfSize=1.5, noiseLevel=noiseLevel/10.,
detectionSigma=5., sourceSigma=sourceAmplitude*5., fluxRange=2.)
statsCtrl = self.prepareStats()
signList = [-1., 0., 1.]
sidelobeShift = afwGeom.Extent2D(4., 0.)
for model, sidelobe, sign in zip(modelImages, sidelobeImages, signList):
Expand All @@ -441,27 +448,15 @@ def testRegularizationSidelobes(self):
dcrModels = DcrModel(modelImages=modelImages)
refModels = [dcrModels[subfilter].clone() for subfilter in range(self.dcrNumSubfilters)]

dcrModels.regularizeModelFreq(modelImages, self.bbox, clampFrequency,
dcrModels.regularizeModelFreq(modelImages, self.bbox, statsCtrl, clampFrequency,
regularizationWidth=regularizationWidth)
for model, refModel, sign in zip(modelImages, refModels, signList):
# The mask and variance planes should be unchanged
self.assertFloatsEqual(model.mask.array, refModel.mask.array)
self.assertFloatsEqual(model.variance.array, refModel.variance.array)
if sign == 0:
# The center subfilter does not have sidelobes, and should be unaffected.
self.assertFloatsEqual(model.image.array, refModel.image.array)
else:
# Make sure the test parameters do reduce the outliers
self.assertGreater(np.sum(np.abs(refModel.image.array - templateImage)),
np.sum(np.abs(model.image.array - templateImage)))
highThreshold = templateImage*clampFrequency
highPix = model.image.array > highThreshold
highPix = ndimage.morphology.binary_opening(highPix, iterations=regularizationWidth)
self.assertFalse(np.all(highPix))
lowThreshold = templateImage/clampFrequency
lowPix = model.image.array < lowThreshold
lowPix = ndimage.morphology.binary_opening(lowPix, iterations=regularizationWidth)
self.assertFalse(np.all(lowPix))
# Make sure the test parameters do reduce the outliers
self.assertGreater(np.sum(np.abs(refModel.image.array - templateImage)),
np.sum(np.abs(model.image.array - templateImage)))

def testRegularizeModelIter(self):
"""Test that large amplitude changes between iterations are restricted.
Expand Down