fix variance calculation to use masks #26
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -22,8 +22,7 @@ | |
| # | ||
|
|
||
| import numpy as np | ||
| import scipy.fftpack | ||
|
|
||
| import lsst.afw.image as afwImage | ||
| import lsst.meas.algorithms as measAlg | ||
|
|
@@ -39,10 +38,13 @@ class DecorrelateALKernelConfig(pexConfig.Config): | |
| \anchor DecorrelateALKernelConfig_ | ||
|
|
||
| \brief Configuration parameters for the DecorrelateALKernelTask | ||
| """ | ||
|
|
||
| ignoreMaskPlanes = pexConfig.ListField( | ||
| dtype = str, | ||
| doc = """Mask planes to ignore for sigma-clipped statistics""", | ||
| default = ("INTRP", "EDGE", "DETECTED", "SAT", "CR", "BAD", "NO_DATA", "DETECTED_NEGATIVE") | ||
| ) | ||
|
|
||
| ## \addtogroup LSST_task_documentation | ||
| ## \{ | ||
|
|
@@ -125,24 +127,42 @@ def __init__(self, *args, **kwargs): | |
| self.statsControl = afwMath.StatisticsControl() | ||
| self.statsControl.setNumSigmaClip(3.) | ||
| self.statsControl.setNumIter(3) | ||
| self.statsControl.setAndMask(afwImage.MaskU.getPlaneBitMask(self.config.ignoreMaskPlanes)) | ||
|
|
||
| def computeVarianceMean(self, exposure): | ||
| statObj = afwMath.makeStatistics(exposure.getMaskedImage().getVariance(), | ||
| exposure.getMaskedImage().getMask(), | ||
| afwMath.MEANCLIP, self.statsControl) | ||
| var = statObj.getValue(afwMath.MEANCLIP) | ||
| return var | ||
|
|
||
| @pipeBase.timeMethod | ||
| def run(self, exposure, templateExposure, subtractedExposure, psfMatchingKernel, | ||
| xcen=None, ycen=None, svar=None, tvar=None): | ||
| """! Perform decorrelation of an image difference exposure. | ||
|
|
||
| Decorrelates the diffim due to the convolution of the templateExposure with the | ||
| A&L PSF matching kernel. Currently can accept a spatially varying matching kernel but in | ||
| this case it simply uses a static kernel from the center of the exposure. The decorrelation | ||
| is described in [DMTN-021, Equation 1](http://dmtn-021.lsst.io/#equation-1), where | ||
| `exposure` is I_1; templateExposure is I_2; `subtractedExposure` is D(k); | ||
| `psfMatchingKernel` is kappa; and svar and tvar are their respective | ||
| variances (see below). | ||
|
|
||
| @param[in] exposure the science afwImage.Exposure used for PSF matching | ||
| @param[in] templateExposure the template afwImage.Exposure used for PSF matching | ||
| @param[in] subtractedExposure the subtracted exposure produced by | ||
| `ip_diffim.ImagePsfMatchTask.subtractExposures()` | ||
| @param[in] psfMatchingKernel an (optionally spatially-varying) PSF matching kernel produced | ||
| by `ip_diffim.ImagePsfMatchTask.subtractExposures()` | ||
| @param[in] xcen X-pixel coordinate to use for computing constant matching kernel to use | ||
| If `None` (default), then use the center of the image. | ||
| @param[in] ycen Y-pixel coordinate to use for computing constant matching kernel to use | ||
| If `None` (default), then use the center of the image. | ||
| @param[in] svar image variance for science image | ||
| If `None` (default) then compute the variance over the entire input science image. | ||
| @param[in] tvar image variance for template image | ||
| If `None` (default) then compute the variance over the entire input template image. | ||
|
|
||
| @return a `pipeBase.Struct` containing: | ||
| * `correctedExposure`: the decorrelated diffim | ||
|
|
@@ -159,79 +179,89 @@ def run(self, templateExposure, exposure, subtractedExposure, psfMatchingKernel) | |
| """ | ||
| self.log.info("Starting.") | ||
| kimg = None | ||
|
|
||
| spatialKernel = psfMatchingKernel | ||
| kimg = afwImage.ImageD(spatialKernel.getDimensions()) | ||
| bbox = subtractedExposure.getBBox() | ||
| if xcen is None: | ||
| xcen = (bbox.getBeginX() + bbox.getEndX()) / 2. | ||
| if ycen is None: | ||
| ycen = (bbox.getBeginY() + bbox.getEndY()) / 2. | ||
| self.log.info("Using matching kernel computed at (%d, %d)" % (xcen, ycen)) | ||
| spatialKernel.computeImage(kimg, True, xcen, ycen) | ||
|
|
||
| if False: # debug code to save spatially varying kernel for analysis | ||
| import cPickle | ||
| import gzip | ||
| cPickle.dump(spatialKernel, gzip.GzipFile('spatialKernel.p.gz', 'wb')) | ||
|
|
||
| if svar is None: | ||
| svar = self.computeVarianceMean(exposure) | ||
| self.log.info("Variance (science): %f" % svar) | ||
| if tvar is None: | ||
| tvar = self.computeVarianceMean(templateExposure) | ||
| self.log.info("Variance (template): %f" % tvar) | ||
|
|
||
| var = self.computeVarianceMean(subtractedExposure) | ||
| self.log.info("Variance (uncorrected diffim): %f" % var) | ||
|
|
||
| corrKernel = DecorrelateALKernelTask._computeDecorrelationKernel(kimg.getArray(), svar, tvar) | ||
| self.log.info("Convolving.") | ||
| correctedExposure, corrKern = DecorrelateALKernelTask._doConvolve(subtractedExposure, corrKernel) | ||
| self.log.info("Updating correctedExposure and its PSF.") | ||
|
|
||
| # Compute the subtracted exposure's updated psf | ||
| psf = subtractedExposure.getPsf().computeImage().getArray() | ||
| psfc = DecorrelateALKernelTask.computeCorrectedDiffimPsf(corrKernel, psf, svar=svar, tvar=tvar) | ||
| psfcI = afwImage.ImageD(psfc.shape[0], psfc.shape[1]) | ||
| psfcI.getArray()[:, :] = psfc | ||
| psfcK = afwMath.FixedKernel(psfcI) | ||
| psfNew = measAlg.KernelPsf(psfcK) | ||
|
There was a problem hiding this comment. Another if False. Please just comment out the block, with the comment as to why above it. There was a problem hiding this comment. Standards have long been that There was a problem hiding this comment. Whose standards? It's very unpythonic, and can be a cause of confusion (commented out code shows up differently in an editor, whereas "if False" just looks like any other code). I couldn't find anything in the coding standards about not commenting out blocks. I even have a community post about this: https://community.lsst.org/t/use-of-if-true-if-false-in-python-and-if-0-in-c/455/3 There was a problem hiding this comment. It's a pythonification of a C++ standard. Granted that we don't actually have a corresponding python standard, but 'we' have long interpreted this as blessing There was a problem hiding this comment. I understand where it comes from in C, but I very strongly advocate that we stop using it in python, especially since |
||
| correctedExposure.setPsf(psfNew) | ||
| self.log.info("Complete.") | ||
|
|
||
| var = self.computeVarianceMean(correctedExposure) | ||
| self.log.info("Variance (corrected diffim): %f" % var) | ||
|
|
||
| return pipeBase.Struct(correctedExposure=correctedExposure, correctionKernel=corrKern) | ||
|
|
||
| @staticmethod | ||
| def _computeDecorrelationKernel(kappa, svar=0.04, tvar=0.04): | ||
| """! Compute the Lupton/ZOGY post-conv. kernel for decorrelating an | ||
| image difference, based on the PSF-matching kernel. | ||
| @param kappa A matching kernel 2-d numpy.array derived from Alard & Lupton PSF matching | ||
| @param svar Average variance of science image used for PSF matching | ||
| @param tvar Average variance of template image used for PSF matching | ||
| @return a 2-d numpy.array containing the correction kernel | ||
|
|
||
| @note As currently implemented, kappa is a static (single, non-spatially-varying) kernel. | ||
| """ | ||
| kappa = DecorrelateALKernelTask._fixOddKernel(kappa) | ||
| kft = scipy.fftpack.fft2(kappa) | ||
| kft = np.sqrt((svar + tvar) / (svar + tvar * kft**2)) | ||
| pck = scipy.fftpack.ifft2(kft) | ||
| pck = scipy.fftpack.ifftshift(pck.real) | ||
| fkernel = DecorrelateALKernelTask._fixEvenKernel(pck) | ||
|
|
||
| # I think we may need to "reverse" the PSF, as in the ZOGY (and Kaiser) papers... | ||
| # This is the same as taking the complex conjugate in Fourier space before FFT-ing back to real space. | ||
| if False: # TBD: figure this out. For now, we are turning it off. | ||
| fkernel = fkernel[::-1, :] | ||
|
|
||
| return fkernel | ||
|
|
||
| @staticmethod | ||
| def computeCorrectedDiffimPsf(kappa, psf, svar=0.04, tvar=0.04): | ||
| """! Compute the (decorrelated) difference image's new PSF. | ||
| new_psf = psf(k) * sqrt((svar + tvar) / (svar + tvar * kappa_ft(k)**2)) | ||
|
|
||
| @param kappa A matching kernel array derived from Alard & Lupton PSF matching | ||
| @param psf The uncorrected psf array of the science image (and also of the diffim) | ||
| @param svar Average variance of science image used for PSF matching | ||
| @param tvar Average variance of template image used for PSF matching | ||
| @return a 2-d numpy.array containing the new PSF | ||
| """ | ||
| def post_conv_psf_ft2(psf, kernel, svar, tvar): | ||
| # Pad psf or kernel symmetrically to make them the same size! | ||
| # Note this assumes they are both square (width == height) | ||
| if psf.shape[0] < kernel.shape[0]: | ||
|
|
@@ -240,20 +270,41 @@ def post_conv_psf_ft2(psf, kernel, sig1=1., sig2=1.): | |
| elif psf.shape[0] > kernel.shape[0]: | ||
| diff = (psf.shape[0] - kernel.shape[0]) // 2 | ||
| kernel = np.pad(kernel, (diff, diff), mode='constant') | ||
| psf_ft = scipy.fftpack.fft2(psf) | ||
| kft = scipy.fftpack.fft2(kernel) | ||
| out = psf_ft * np.sqrt((svar + tvar) / (svar + tvar * kft**2)) | ||
| return out | ||
|
|
||
| def post_conv_psf(psf, kernel, svar, tvar): | ||
| kft = post_conv_psf_ft2(psf, kernel, svar, tvar) | ||
| out = scipy.fftpack.ifft2(kft) | ||
| return out | ||
|
|
||
| pcf = post_conv_psf(psf=psf, kernel=kappa, svar=svar, tvar=tvar) | ||
| pcf = pcf.real / pcf.real.sum() | ||
| return pcf | ||
|
|
||
| @staticmethod | ||
| def _fixOddKernel(kernel): | ||
|
There was a problem hiding this comment. Not sure this should be a static method: it looks more like a stand-alone, top-level function. There was a problem hiding this comment. I'm going to go with the philosophy of this comment: http://stackoverflow.com/a/16079946 Note that a subclass may want to override some of these methods and I want to make that an explicit choice. |
||
| """! Take a kernel with odd dimensions and make them even for FFT | ||
|
|
||
| @param kernel a numpy.array | ||
| @return a fixed kernel numpy.array. Returns a copy if the dimensions needed to change; | ||
| otherwise just return the input kernel. | ||
| """ | ||
| # Note this works best for the FFT if we left-pad | ||
| out = kernel | ||
|
There was a problem hiding this comment. It looks like the goal is to not modify the kernel in place, so you probably want |
||
| changed = False | ||
| if (out.shape[0] % 2) == 1: | ||
| out = np.pad(out, ((1, 0), (0, 0)), mode='constant') | ||
| changed = True | ||
| if (out.shape[1] % 2) == 1: | ||
| out = np.pad(out, ((0, 0), (1, 0)), mode='constant') | ||
| changed = True | ||
| if changed: | ||
| out *= (np.mean(kernel) / np.mean(out)) # need to re-scale to same mean for FFT | ||
| return out | ||
|
|
||
| @staticmethod | ||
|
There was a problem hiding this comment. Another unlikely static method that should just be a top level function. Actually, I think all your staticmethods in this class should probably just be top level. |
||
| def _fixEvenKernel(kernel): | ||
| """! Take a kernel with even dimensions and make them odd, centered correctly. | ||
|
|
@@ -288,12 +339,10 @@ def _doConvolve(exposure, kernel): | |
| @note We use afwMath.convolve() but keep scipy.convolve for debugging. | ||
| @note We re-center the kernel if necessary and return the possibly re-centered kernel | ||
| """ | ||
| kernelImg = afwImage.ImageD(kernel.shape[0], kernel.shape[1]) | ||
| kernelImg.getArray()[:, :] = kernel | ||
| kern = afwMath.FixedKernel(kernelImg) | ||
| maxloc = np.unravel_index(np.argmax(kernel), kernel.shape) | ||
| kern.setCtrX(maxloc[0]) | ||
| kern.setCtrY(maxloc[1]) | ||
| outExp = exposure.clone() # Do this to keep WCS, PSF, masks, etc. | ||
|
|
||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Might be handy to note in the docstring that exposure, template match equation whatever in some paper.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's not in a paper, yet. But I will reference the (still in progress) DMTN-021 with those details.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I see the reference. Thanks. Maybe I_1 and I_2 should be I_s and I_t in DMTN-021, to make it explicit?