Skip to content

Commit

Permalink
Merge pull request #107 from lsst/tickets/DM-16864
Browse files Browse the repository at this point in the history 
Tickets/DM-16864: New frequency regularization using relative model
  • Loading branch information
@isullivan
isullivan committed Feb 1, 2019
2 parents e8afd39 + e1d653a commit 42ca1a7
Show file tree
Hide file tree
Showing 2 changed files with 74 additions and 58 deletions.
79 changes: 50 additions & 29 deletions python/lsst/ip/diffim/dcrModel.py
View file
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
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

0 comments on commit 42ca1a7

Please sign in to comment.