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

DM-17426 #1

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
+173 −48
Open

DM-17426 #1

Changes from 1 commit
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
7ce5656
DM-17098
michitaro Feb 13, 2019
b2ebc0a
Implement MaskObjectsTask
michitaro Apr 1, 2019
832078c
Use MaskObjectsTask in constructSky
michitaro Apr 1, 2019
773e3a2
Use MaskObjectsTask in skyCorrection
michitaro Apr 1, 2019
a078e08
2D interpolation in FocalPlaneBackground
michitaro Apr 1, 2019
876d95d
add an extra focalplane sky-subtraction
michitaro Apr 1, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 58 additions & 0 deletions python/lsst/pipe/drivers/background.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -468,6 +468,9 @@ class FocalPlaneBackgroundConfig(Config):
"NONE": "No background estimation is to be attempted",
},
)
doSmooth = Field(dtype=bool, default=False, doc="Do smoothing?")
smoothScale = Field(dtype=float, doc="Smoothing scale")
Copy link

Choose a reason for hiding this comment

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

No default?

smoothWindowSize = Field(dtype=int, default=15, doc="Window size for smoothing")
binning = Field(dtype=int, default=64, doc="Binning to use for CCD background model (pixels)")


Expand Down Expand Up @@ -721,6 +724,12 @@ def getStatsImage(self):
values /= self._numbers
thresh = self.config.minFrac*self.config.xSize*self.config.ySize
isBad = self._numbers.getArray() < thresh
if self.config.doSmooth:
array = values.getArray()
array[isBad] = numpy.nan
gridSize = min(self.config.xSize, self.config.ySize)
array[:] = NanSafeSmoothing.gaussianSmoothing(array, self.config.smoothWindowSize, self.config.smoothScale / gridSize)
Copy link

Choose a reason for hiding this comment

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

Line length.

isBad = numpy.isnan(values.array)
interpolateBadPixels(values.getArray(), isBad, self.config.interpolation)
return values

Expand Down Expand Up @@ -785,3 +794,52 @@ def run(self, exp):
detected = mask.array & mask.getPlaneBitMask(['DETECTED']) > 0
exp.maskedImage.image.array[detected] = smooth.getImage().getArray()[detected]


class NanSafeSmoothing:
Copy link

Choose a reason for hiding this comment

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

I think this section needs a different implementation.
It looks like you're using class here just as a namespace, which feels dirty. I don't think it adds anything beyond making a stand-alone function with a couple of embedded functions.
A direct convolution implemented in python is going to be crazy slow (especially since you're not making use of the fact that the kernel is separable). Is there a reason you can't use the convolution functionality in afw, as is done in SourceDetectionTask? If you're worried about NaNs, you can always (temporarily?) replace them with zeros when you do the convolution.

Copy link
Author

@michitaro michitaro Apr 2, 2019 •
edited
Loading

Choose a reason for hiding this comment

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

If you're worried about NaNs, you can always (temporarily?) replace them with zeros when you do the convolution.

This is why I didn't use the existing convolution code. Replacing NaNs with 0 is not equivalent to my code.
Sometimes super-pixels on the edge of the field are bit high. In that case it is better that the super-pixels out of the field are as the are extrapolated (this is what I intended to) than assuming that they are zero (replacing NaNs with 0). In other words, replacing NaNs with 0 leads low sky estimation near the bright edge.

Screen Shot 2019-04-02 at 17 17 16

A direct convolution in Python is crazy slow, but the targets of this code are relatively small images such as FocalPlaneBackground._values whose dimensions are about 140x150. Convolutions on such images end in about 0.5 second on my mac.

It looks like you're using class here just as a namespace, which feels dirty. I don't think it adds anything beyond making a stand-alone function with a couple of embedded functions.

I totally agree with you. I will put these functions flat in the module.

Copy link

Choose a reason for hiding this comment

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

I figured out the difference between your convolution and the vanilla convolution. This will allow us to condense the code, and probably run it faster:

def safeConvolve(array, sigma=2):
    bad = np.isnan(array)
    safe = np.where(bad, 0.0, array)
    convolved = gaussian_filter(safe, sigma, mode="constant", cval=0.0)
    corr = np.where(bad, 0.0, 1.0)
    ones = np.ones_like(array)
    numerator = gaussian_filter(ones, sigma, mode="constant", cval=0.0)
    denominator = gaussian_filter(corr, sigma, mode="constant", cval=0.0)
    return convolved*numerator/denominator

This reproduces your result:
image

Copy link
Author

Choose a reason for hiding this comment

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

Excellent solution.

'''
Smooth image dealing with NaN pixels
'''

@classmethod
def gaussianSmoothing(cls, array, windowSize, sigma):
return cls._safeConvolve2d(array, cls._gaussianKernel(windowSize, sigma))

@classmethod
def _gaussianKernel(cls, windowSize, sigma):
''' Returns 2D gaussian kernel '''
s = sigma
r = windowSize
X, Y = numpy.meshgrid(
numpy.linspace(-r, r, 2 * r + 1),
numpy.linspace(-r, r, 2 * r + 1),
)
kernel = cls._normalDist(X, s) * cls._normalDist(Y, s)
# cut off
kernel[X**2 + Y**2 > (r + 0.5)**2] = 0.
return kernel / kernel.sum()

@staticmethod
def _normalDist(x, s=1., m=0.):
''' Normal Distribution '''
return 1. / (s * numpy.sqrt(2. * numpy.pi)) * numpy.exp(-(x-m)**2/(2*s**2)) / (s * numpy.sqrt(2*numpy.pi))
Copy link

Choose a reason for hiding this comment

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

Why two normalisation factors of (s * numpy.sqrt(2. * numpy.pi))?


@staticmethod
def _safeConvolve2d(image, kernel):
''' Convolve 2D safely dealing with NaNs in `image` '''
assert numpy.ndim(image) == 2
assert numpy.ndim(kernel) == 2
assert kernel.shape[0] % 2 == 1 and kernel.shape[1] % 2 == 1
ks = kernel.shape
kl = (ks[0] - 1) // 2, \
(ks[1] - 1) // 2
image2 = numpy.pad(image, ((kl[0], kl[0]), (kl[1], kl[1])), 'constant', constant_values=numpy.nan)
convolved = numpy.empty_like(image)
convolved.fill(numpy.nan)
for yi in range(convolved.shape[0]):
for xi in range(convolved.shape[1]):
patch = image2[yi : yi + ks[0], xi : xi + ks[1]]
c = patch * kernel
ok = numpy.isfinite(c)
if numpy.any(ok):
convolved[yi, xi] = numpy.nansum(c) / kernel[ok].sum()
return convolved