Skip to content

Commit

Permalink
IsrTask: add option to set variance from empirical read noise
Browse files Browse the repository at this point in the history
The read noise is measured from the overscan.

This might be useful when working on detectors whose read noise hasn't
been measured, or drifts due to hardware tests.
  • Loading branch information
PaulPrice committed Apr 24, 2018
1 parent 2c715da commit 2b7fdf1
Show file tree
Hide file tree
Showing 3 changed files with 211 additions and 16 deletions.
2 changes: 1 addition & 1 deletion python/lsst/ip/isr/isrFunctions.py
Original file line number Diff line number Diff line change
Expand Up @@ -437,7 +437,7 @@ def overscanCorrection(ampMaskedImage, overscanImage, fitType='MEDIAN', order=1,

offImage = ampImage.Factory(ampImage.getDimensions())
offArray = offImage.getArray()
overscanFit = overscanImage.Factory(overscanImage.getDimensions())
overscanFit = afwImage.ImageF(overscanImage.getDimensions())
overscanArray = overscanFit.getArray()
if shortInd == 1:
offArray[:, :] = fitBiasArr[:, numpy.newaxis]
Expand Down
76 changes: 61 additions & 15 deletions python/lsst/ip/isr/isrTask.py
Original file line number Diff line number Diff line change
Expand Up @@ -292,6 +292,11 @@ class IsrTaskConfig(pexConfig.Config):
default=True,
doc="Load and use transmission_atmosphere (if doAttachTransmissionCurve is True)?"
)
doEmpiricalReadNoise = pexConfig.Field(
dtype=bool,
default=False,
doc="Calculate empirical read noise instead of value from AmpInfo data?"
)

## @addtogroup LSST_task_documentation
## @{
Expand Down Expand Up @@ -521,12 +526,16 @@ def run(self, ccdExposure, bias=None, linearizer=None, dark=None, flat=None, def
assert not self.config.doAssembleCcd, "You need a Detector to run assembleCcd"
ccd = [FakeAmp(ccdExposure, self.config)]

overscans = []
for amp in ccd:
# if ccdExposure is one amp, check for coverage to prevent performing ops multiple times
if ccdExposure.getBBox().contains(amp.getBBox()):
self.saturationDetection(ccdExposure, amp)
self.suspectDetection(ccdExposure, amp)
self.overscanCorrection(ccdExposure, amp)
overscanResults = self.overscanCorrection(ccdExposure, amp)
overscans.append(overscanResults.overscanImage if overscanResults is not None else None)
else:
overscans.append(None)

if self.config.doCrosstalk:
self.crosstalk.run(ccdExposure, crosstalkSources)
Expand All @@ -542,11 +551,12 @@ def run(self, ccdExposure, bias=None, linearizer=None, dark=None, flat=None, def
if self.doLinearize(ccd):
linearizer(image=ccdExposure.getMaskedImage().getImage(), detector=ccd, log=self.log)

for amp in ccd:
assert len(ccd) == len(overscans)
for amp, overscanImage in zip(ccd, overscans):
# if ccdExposure is one amp, check for coverage to prevent performing ops multiple times
if ccdExposure.getBBox().contains(amp.getBBox()):
ampExposure = ccdExposure.Factory(ccdExposure, amp.getBBox())
self.updateVariance(ampExposure, amp)
self.updateVariance(ampExposure, amp, overscanImage)

interpolationDone = False

Expand Down Expand Up @@ -699,12 +709,23 @@ def doLinearize(self, detector):
return self.config.doLinearize and \
detector.getAmpInfoCatalog()[0].getLinearityType() != NullLinearityType

def updateVariance(self, ampExposure, amp):
"""!Set the variance plane based on the image plane, plus amplifier gain and read noise
def updateVariance(self, ampExposure, amp, overscanImage=None):
"""Set the variance plane using the amplifier gain and read noise
@param[in,out] ampExposure exposure to process
@param[in] amp amplifier detector information
The read noise is calculated from the ``overscanImage`` if the
``doEmpiricalReadNoise`` option is set in the configuration; otherwise
the value from the amplifier data is used.
Parameters
----------
ampExposure : `lsst.afw.image.Exposure`
Exposure to process.
amp : `lsst.afw.table.AmpInfoRecord` or `FakeAmp`
Amplifier detector data.
overscanImage : `lsst.afw.image.MaskedImage`, optional.
Image of overscan, required only for empirical read noise.
"""
maskPlanes = [self.config.saturatedMaskName, self.config.suspectMaskName]
gain = amp.getGain()
if not math.isnan(gain):
if gain <= 0:
Expand All @@ -713,10 +734,18 @@ def updateVariance(self, ampExposure, amp):
(amp.getName(), gain, patchedGain))
gain = patchedGain

if self.config.doEmpiricalReadNoise and overscanImage is not None:
stats = afwMath.StatisticsControl()
stats.setAndMask(overscanImage.mask.getPlaneBitMask(maskPlanes))
readNoise = afwMath.makeStatistics(overscanImage, afwMath.STDEVCLIP, stats).getValue()
self.log.info("Calculated empirical read noise for amp %s: %f", amp.getName(), readNoise)
else:
readNoise = amp.getReadNoise()

isrFunctions.updateVariance(
maskedImage=ampExposure.getMaskedImage(),
gain=gain,
readNoise=amp.getReadNoise(),
readNoise=readNoise,
)

def flatCorrection(self, exposure, flatExposure, invert=False):
Expand Down Expand Up @@ -870,10 +899,27 @@ def maskAndInterpNan(self, exposure):
)

def overscanCorrection(self, exposure, amp):
"""!Apply overscan correction, in place
"""Apply overscan correction, in-place
@param[in,out] exposure exposure to process; must include both DataSec and BiasSec pixels
@param[in] amp amplifier device data
Parameters
----------
exposure : `lsst.afw.image.Exposure`
Exposure to process; must include both data and bias regions.
amp : `lsst.afw.table.AmpInfoRecord`
Amplifier device data.
Results
-------
result : `lsst.pipe.base.Struct` or `NoneType`
`None` if there is no overscan; otherwise, this is a
result struct with components:
- ``imageFit``: Value(s) removed from image (scalar or
`lsst.afw.image.Image`).
- ``overscanFit``: Value(s) removed from overscan (scalar or
`lsst.afw.image.Image`).
- ``overscanImage``: Image of the overscan, post-subtraction
(`lsst.afw.image.Image`).
"""
if not amp.getHasRawInfo():
raise RuntimeError("This method must be executed on an amp with raw information.")
Expand All @@ -884,17 +930,17 @@ def overscanCorrection(self, exposure, amp):

maskedImage = exposure.getMaskedImage()
dataView = maskedImage.Factory(maskedImage, amp.getRawDataBBox())
overscanImage = maskedImage.Factory(maskedImage, amp.getRawHorizontalOverscanBBox())

expImage = exposure.getMaskedImage().getImage()
overscanImage = expImage.Factory(expImage, amp.getRawHorizontalOverscanBBox())

return isrFunctions.overscanCorrection(
results = isrFunctions.overscanCorrection(
ampMaskedImage=dataView,
overscanImage=overscanImage,
fitType=self.config.overscanFitType,
order=self.config.overscanOrder,
collapseRej=self.config.overscanRej,
)
results.overscanImage = overscanImage
return results

def addDistortionModel(self, exposure, camera):
"""!Update the WCS in exposure with a distortion model based on camera geometry
Expand Down
149 changes: 149 additions & 0 deletions tests/test_empiricalVariance.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,149 @@
#
# LSST Data Management System
# Copyright 2018 LSST Corporation.
#
# This product includes software developed by the
# LSST Project (http://www.lsst.org/).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the LSST License Statement and
# the GNU General Public License along with this program. If not,
# see <http://www.lsstcorp.org/LegalNotices/>.
#
import numpy as np
import unittest

import lsst.utils.tests

from lsst.daf.base import PropertyList
from lsst.afw.cameraGeom import Detector, SCIENCE, Orientation, TransformMap, CameraSys
from lsst.afw.table import AmpInfoCatalog, AmpInfoTable
from lsst.afw.geom import Point2I, Extent2I, Box2I, Extent2D
from lsst.afw.image import ExposureF, VisitInfo
from lsst.pipe.base import Struct

from lsst.ip.isr.isrTask import IsrTask


def makeAmplifier(catalog, name, bbox, rawImageBox, overscanBox, gain, readNoise, saturation):
amp = catalog.addNew()
amp.setName(name)
amp.setBBox(bbox)
amp.setRawDataBBox(rawImageBox)
amp.setRawHorizontalOverscanBBox(overscanBox)
amp.setHasRawInfo(True)
amp.setGain(gain)
amp.setReadNoise(readNoise)
amp.setSaturation(saturation)
amp.setSuspectLevel(np.nan)
return amp

class EmpiricalVarianceTestCast(lsst.utils.tests.TestCase):
def setUp(self):
"""Constructs a CCD with two amplifiers and prepares for ISR"""
np.random.seed(12345)
baseValue = 100.0
gain = 1.0
readNoise = 123456789.0
saturation = 987654321.0
height = 234
imageSize = Extent2I(123, height)
overscanSize = Extent2I(16, height)
self.sigma = 1.234

# Set up the various regions
overscan1 = Box2I(Point2I(0, 0), overscanSize)
image1 = Box2I(Point2I(overscanSize[0], 0), imageSize)
image2 = Box2I(Point2I(overscanSize[0] + imageSize[0], 0), imageSize)
overscan2 = Box2I(Point2I(overscanSize[0] + 2*imageSize[0], 0), overscanSize)

leftBox = Box2I(overscan1.getMin(), Extent2I(overscan1.getWidth() + image1.getWidth(), height))
rightBox = Box2I(image2.getMin(), Extent2I(image2.getWidth() + overscan2.getWidth(), height))

target1 = Box2I(Point2I(0, 0), imageSize)
target2 = Box2I(Point2I(image1.getWidth(), 0), imageSize)

# Set the pixels
exposure = ExposureF(Box2I(Point2I(0, 0), Extent2I(imageSize[0]*2 + overscanSize[0]*2, height)))
yy = np.arange(0, height, 1, dtype=np.float32)
leftImage = ExposureF(exposure, leftBox)
leftImage.image.array[:] = baseValue + yy[:, np.newaxis]
rightImage = ExposureF(exposure, rightBox)
rightImage.image.array[:] = baseValue - yy[:, np.newaxis]

leftOverscan = ExposureF(exposure, overscan1)
leftOverscan.image.array += np.random.normal(0.0, self.sigma, leftOverscan.image.array.shape)
rightOverscan = ExposureF(exposure, overscan2)
rightOverscan.image.array += np.random.normal(0.0, self.sigma, leftOverscan.image.array.shape)
exposure.mask.array[:] = 0.0
exposure.variance.array[:] = np.nan

# Construct the detectors
amps = AmpInfoCatalog(AmpInfoTable.makeMinimalSchema())
makeAmplifier(amps, "left", target1, image1, overscan1, gain, readNoise, saturation)
makeAmplifier(amps, "right", target2, image2, overscan2, gain, readNoise, saturation)
ccdBox = Box2I(Point2I(0, 0), Extent2I(image1.getWidth() + image2.getWidth(), height))
ccd = Detector("detector", 1, SCIENCE, "det1", ccdBox, amps, Orientation(), Extent2D(1.0, 1.0), {})
exposure.setDetector(ccd)
header = PropertyList()
header.add("EXPTIME", 0.0)
exposure.getInfo().setVisitInfo(VisitInfo(header))

self.exposure = exposure
self.config = IsrTask.ConfigClass()

# Disable everything we don't care about
self.config.doBias = False
self.config.doDark = False
self.config.doFlat = False
self.config.doFringe = False
self.config.doDefect = False
self.config.doAddDistortionModel = False
self.config.doWrite = False
self.config.expectWcs = False
self.config.doLinearize = False
self.config.doCrosstalk = False
self.config.doBrighterFatter = False
self.config.doAttachTransmissionCurve = False

# Set the things that match our test setup
self.config.overscanFitType = "CHEB"
self.config.overscanOrder = 1
self.config.doEmpiricalReadNoise = True

self.task = IsrTask(config=self.config)

def tearDown(self):
del self.exposure

def testEmpiricalVariance(self):
results = self.task.run(self.exposure)
exposure = results.exposure
self.assertFloatsEqual(exposure.mask.array, 0)
# Image is not exactly zero because the noise in the overscan (required to be able to set
# the empirical variance) leads to a slight misestimate in the polynomial fit.
self.assertFloatsAlmostEqual(exposure.image.array, 0.0, atol=5.0e-2)
self.assertFloatsAlmostEqual(exposure.variance.array, self.sigma**2, rtol=5.0e-2)


class MemoryTester(lsst.utils.tests.MemoryTestCase):
pass


def setup_module(module):
lsst.utils.tests.init()


if __name__ == "__main__":
import sys
setup_module(sys.modules[__name__])
unittest.main()

0 comments on commit 2b7fdf1

Please sign in to comment.