Skip to content

Commit

Permalink
Address review comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@enourbakhsh
enourbakhsh committed Aug 21, 2023
1 parent 299615a commit baed1a8
Show file tree
Hide file tree
Showing 2 changed files with 182 additions and 128 deletions.
177 changes: 91 additions & 86 deletions python/lsst/ip/isr/ampOffset.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@
from lsst.afw.table import SourceTable
from lsst.meas.algorithms import SourceDetectionTask, SubtractBackgroundTask
from lsst.pex.config import Config, ConfigurableField, Field
from lsst.pipe.base import Task
from lsst.pipe.base import Task, Struct


class AmpOffsetConfig(Config):
Expand Down Expand Up @@ -73,12 +73,13 @@ def setDefaults(self):
dtype=float,
default=5.0,
)
ampEdgeWindow = Field(
doc="Pixel size of the sliding window used to generate rolling average amp offset values. It should "
"should be reconfigured for every instrument (HSC, LSSTCam, etc.) and should not exceed the shorter "
"dimension of amps. If not provided, it defaults to the shorter dimension of amps.",
dtype=int,
default=0,
ampEdgeWindowFrac = Field(
doc="Fraction of the amp edge lengths utilized as the sliding window for generating rolling average "
"amp offset values. It should be reconfigured for every instrument (HSC, LSSTCam, etc.) and should "
"not exceed 1. If not provided, it defaults to the fraction that recovers the pixel size of the "
"sliding window used in obs_subaru for compatibility with existing HSC data.",
dtype=float,
default=512 / 4176,
)
doBackground = Field(
doc="Estimate and subtract background prior to amp offset estimation?",
Expand Down Expand Up @@ -118,6 +119,8 @@ def __init__(self, *args, **kwargs):
self.makeSubtask("background")
if self.config.doDetection:
self.makeSubtask("detection")
# Initialize all of the instance variables here.
self.shortAmpSide = 0

def run(self, exposure):
"""Calculate amp offset values, determine corrective pedestals for each
Expand All @@ -133,9 +136,14 @@ def run(self, exposure):
exp = exposure.clone()
bitMask = exp.mask.getPlaneBitMask(self.background.config.ignoredPixelMask)
amps = exp.getDetector().getAmplifiers()

# Check that all amps have the same gemotry.
ampDims = [amp.getBBox().getDimensions() for amp in amps]
if not all(dim == ampDims[0] for dim in ampDims):
raise RuntimeError("All amps should have the same geometry.")

# Assuming all the amps have the same geometry.
self.shortAmpSide = np.min(amps[0].getBBox().getDimensions())
self.longAmpSide = np.max(amps[0].getBBox().getDimensions())
self.shortAmpSide = np.min(ampDims[0])

# Fit and subtract background.
if self.config.doBackground:
Expand Down Expand Up @@ -165,18 +173,18 @@ def run(self, exposure):
# Safety check: do any pixels remain for amp offset estimation?
if (exp.mask.array & bitMask).all():
self.log.warning("All pixels masked: cannot calculate any amp offset corrections.")
pedestals = np.zeros(len(amps))
else:
# Set up amp offset inputs.
im = exp.image
im.array[(exp.mask.array & bitMask) > 0] = np.nan

if self.config.ampEdgeWindow > self.shortAmpSide:
if self.config.ampEdgeWindowFrac > 1:
raise RuntimeError(
f"The sliding window ({self.config.ampEdgeWindow} pixels) is larger than the smaller "
f"dimension in amplifiers ({self.shortAmpSide} pixels). This can lead to downstream "
"issues after convolution in the `getSideAmpOffset()` method. Please adjust "
f"`ampEdgeWindow` in the config to a value equal to or less than {self.shortAmpSide} "
"pixels and rerun."
f"The specified fraction (`ampEdgeWindowFrac`={self.config.ampEdgeWindowFrac}) of the "
"edge length exceeds 1. This leads to complications downstream, after convolution in "
"the `getSideAmpOffset()` method. Please modify the `ampEdgeWindowFrac` value in the "
"config to be 1 or less and rerun."
)

# Determine amplifier geometry.
Expand All @@ -193,17 +201,20 @@ def run(self, exposure):
# If least-squares minimization fails, convert NaNs to zeroes,
# ensuring that no values are erroneously added/subtracted.
pedestals = np.nan_to_num(np.linalg.lstsq(A, B, rcond=None)[0])
metadata = exposure.getMetadata()
for amp, pedestal in zip(amps, pedestals):
ampIm = exposure.image[amp.getBBox()].array
ampIm -= pedestal
ampName = amp.getName()
metadata.set(
f"LSST ISR AMPOFFSET PEDESTAL {ampName}",
float(pedestal),
f"Pedestal level subtracted from amp {ampName}"
)
self.log.info(f"amp pedestal values: {', '.join([f'{x:.4f}' for x in pedestals])}")

metadata = exposure.getMetadata()
for amp, pedestal in zip(amps, pedestals):
ampIm = exposure.image[amp.getBBox()].array
ampIm -= pedestal
ampName = amp.getName()
metadata.set(
f"LSST ISR AMPOFFSET PEDESTAL {ampName}",
float(pedestal),
f"Pedestal level subtracted from amp {ampName}",
)
self.log.info(f"amp pedestal values: {', '.join([f'{x:.4f}' for x in pedestals])}")

return Struct(pedestals=pedestals)

def getAmpAssociations(self, amps):
"""Determine amp geometry and amp associations from a list of
Expand Down Expand Up @@ -303,66 +314,6 @@ def getNeighbors(self, ampIds, ampId):
sides.append(side)
return neighbors, sides

def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB):
"""Calculate the amp offset for a given interface between two
amplifiers.
Parameters
----------
ampIdA : int
ID of the first amplifier.
ampIdB : int
ID of the second amplifier.
edgeA : numpy.ndarray
Amp edge for the first amplifier.
edgeB : numpy.ndarray
Amp edge for the second amplifier.
Returns
-------
interfaceOffset : float
The calculated amp offset value for the given interface between
amps A and B.
"""
sctrl = StatisticsControl()
# NOTE: Taking the difference with the order below fixes the sign flip
# in the B matrix.
edgeDiff = edgeA - edgeB
if self.config.ampEdgeWindow == 0:
if len(edgeDiff) > self.shortAmpSide:
window = self.shortAmpSide
else:
window = int(self.shortAmpSide**2 / self.longAmpSide)
self.log.info(
f"ampEdgeWindow automatically set to {window} pixels, for the "
f"amp {ampIdA} interfacing amp {ampIdB}, as value of 0 or no "
"input is provided for ampEdgeWindow."
)
else:
window = self.config.ampEdgeWindow
# Compute rolling averages.
edgeDiffSum = np.convolve(np.nan_to_num(edgeDiff), np.ones(window), "same")
edgeDiffNum = np.convolve(~np.isnan(edgeDiff), np.ones(window), "same")
edgeDiffAvg = edgeDiffSum / np.clip(edgeDiffNum, 1, None)
edgeDiffAvg[np.isnan(edgeDiff)] = np.nan
# Take clipped mean of rolling average data as amp offset value.
interfaceOffset = makeStatistics(edgeDiffAvg, MEANCLIP, sctrl).getValue()
# Perform a couple of do-no-harm safety checks:
# a) The fraction of unmasked pixel rows is > ampEdgeMinFrac,
# b) The absolute offset ADU value is < ampEdgeMaxOffset.
ampEdgeGoodFrac = 1 - (np.sum(np.isnan(edgeDiffAvg)) / len(edgeDiffAvg))
minFracFail = ampEdgeGoodFrac < self.config.ampEdgeMinFrac
maxOffsetFail = np.abs(interfaceOffset) > self.config.ampEdgeMaxOffset
if minFracFail or maxOffsetFail:
interfaceOffset = 0
interfaceId = f"{ampIdA}{ampIdB}"
self.log.debug(
f"amp interface {interfaceId} : "
f"viable edge difference frac = {ampEdgeGoodFrac}, "
f"interface offset = {interfaceOffset:.3f}"
)
return interfaceOffset

def getAmpOffsets(self, im, amps, associations, sides):
"""Calculate the amp offsets for all amplifiers.
Expand Down Expand Up @@ -448,3 +399,57 @@ def getAmpEdges(self, im, amps, ampSides):
warnings.filterwarnings("ignore", r"All-NaN (slice|axis) encountered")
ampEdges[ampId][ampSide] = np.nanmedian(strip, axis=ampSide % 2) # 1D medianified strip
return ampEdges

def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB):
"""Calculate the amp offset for a given interface between two
amplifiers.
Parameters
----------
ampIdA : int
ID of the first amplifier.
ampIdB : int
ID of the second amplifier.
edgeA : numpy.ndarray
Amp edge for the first amplifier.
edgeB : numpy.ndarray
Amp edge for the second amplifier.
Returns
-------
interfaceOffset : float
The calculated amp offset value for the given interface between
amps A and B.
"""
interfaceId = f"{ampIdA}{ampIdB}"
sctrl = StatisticsControl()
# NOTE: Taking the difference with the order below fixes the sign flip
# in the B matrix.
edgeDiff = edgeA - edgeB
window = int(self.config.ampEdgeWindowFrac * len(edgeDiff))
# Compute rolling averages.
edgeDiffSum = np.convolve(np.nan_to_num(edgeDiff), np.ones(window), "same")
edgeDiffNum = np.convolve(~np.isnan(edgeDiff), np.ones(window), "same")
edgeDiffAvg = edgeDiffSum / np.clip(edgeDiffNum, 1, None)
edgeDiffAvg[np.isnan(edgeDiff)] = np.nan
# Take clipped mean of rolling average data as amp offset value.
interfaceOffset = makeStatistics(edgeDiffAvg, MEANCLIP, sctrl).getValue()
# Perform a couple of do-no-harm safety checks:
# a) The fraction of unmasked pixel rows is > ampEdgeMinFrac,
# b) The absolute offset ADU value is < ampEdgeMaxOffset.
ampEdgeGoodFrac = 1 - (np.sum(np.isnan(edgeDiffAvg)) / len(edgeDiffAvg))
minFracFail = ampEdgeGoodFrac < self.config.ampEdgeMinFrac
maxOffsetFail = np.abs(interfaceOffset) > self.config.ampEdgeMaxOffset
if minFracFail or maxOffsetFail:
interfaceOffset = 0
self.log.warning(
f"The fraction of unmasked pixels for amp interface {interfaceId} is below the threshold "
f"({self.config.ampEdgeMinFrac}) or the absolute offset value exceeds the limit "
f"({self.config.ampEdgeMaxOffset} ADU). Setting the interface offset to 0."
)
self.log.debug(
f"amp interface {interfaceId} : "
f"viable edge difference frac = {ampEdgeGoodFrac}, "
f"interface offset = {interfaceOffset:.3f}"
)
return interfaceOffset
Loading

0 comments on commit baed1a8

Please sign in to comment.