LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
ampOffset.py
Go to the documentation of this file.
1# This file is part of ip_isr.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22__all__ = ["AmpOffsetConfig", "AmpOffsetTask"]
23
24import warnings
25
26import numpy as np
27from lsst.afw.math import MEANCLIP, StatisticsControl, makeStatistics
28from lsst.afw.table import SourceTable
29from lsst.meas.algorithms import SourceDetectionTask, SubtractBackgroundTask
30from lsst.pex.config import Config, ConfigurableField, Field
31from lsst.pipe.base import Struct, Task
32
33
35 """Configuration parameters for AmpOffsetTask."""
36
37 def setDefaults(self):
38 self.background.algorithm = "AKIMA_SPLINE"
39 self.background.useApprox = False
40 self.background.ignoredPixelMask = [
41 "BAD",
42 "SAT",
43 "INTRP",
44 "CR",
45 "EDGE",
46 "DETECTED",
47 "DETECTED_NEGATIVE",
48 "SUSPECT",
49 "NO_DATA",
50 ]
51 self.detection.reEstimateBackground = False
52
53 # This maintains existing behavior and test values after DM-39796.
54 self.detection.thresholdType = "stdev"
55
56 ampEdgeInset = Field(
57 doc="Number of pixels the amp edge strip is inset from the amp edge. A thin strip of pixels running "
58 "parallel to the edge of the amp is used to characterize the average flux level at the amp edge.",
59 dtype=int,
60 default=5,
61 )
62 ampEdgeWidth = Field(
63 doc="Pixel width of the amp edge strip, starting at ampEdgeInset and extending inwards.",
64 dtype=int,
65 default=64,
66 )
67 ampEdgeMinFrac = Field(
68 doc="Minimum allowed fraction of viable pixel rows along an amp edge. No amp offset estimate will be "
69 "generated for amp edges that do not have at least this fraction of unmasked pixel rows.",
70 dtype=float,
71 default=0.5,
72 )
73 ampEdgeMaxOffset = Field(
74 doc="Maximum allowed amp offset ADU value. If a measured amp offset value is larger than this, the "
75 "result will be discarded and therefore not used to determine amp pedestal corrections.",
76 dtype=float,
77 default=5.0,
78 )
79 ampEdgeWindowFrac = Field(
80 doc="Fraction of the amp edge lengths utilized as the sliding window for generating rolling average "
81 "amp offset values. It should be reconfigured for every instrument (HSC, LSSTCam, etc.) and should "
82 "not exceed 1. If not provided, it defaults to the fraction that recovers the pixel size of the "
83 "sliding window used in obs_subaru for compatibility with existing HSC data.",
84 dtype=float,
85 default=512 / 4176,
86 )
87 doBackground = Field(
88 doc="Estimate and subtract background prior to amp offset estimation?",
89 dtype=bool,
90 default=True,
91 )
92 background = ConfigurableField(
93 doc="An initial background estimation step run prior to amp offset calculation.",
94 target=SubtractBackgroundTask,
95 )
96 backgroundFractionSample = Field(
97 doc="The fraction of the shorter side of the amplifier used for background binning.",
98 dtype=float,
99 default=1.0,
100 )
101 doDetection = Field(
102 doc="Detect sources and update cloned exposure prior to amp offset estimation?",
103 dtype=bool,
104 default=True,
105 )
106 detection = ConfigurableField(
107 doc="Source detection to add temporary detection footprints prior to amp offset calculation.",
108 target=SourceDetectionTask,
109 )
110 applyWeights = Field(
111 doc="Weights the amp offset calculation by the length of the interface between amplifiers. Applying "
112 "weights does not affect outcomes for amplifiers in a 2D grid with square-shaped amplifiers or in "
113 "any 1D layout on a detector, regardless of whether the amplifiers are square.",
114 dtype=bool,
115 default=True,
116 )
117
118
119class AmpOffsetTask(Task):
120 """Calculate and apply amp offset corrections to an exposure."""
121
122 ConfigClass = AmpOffsetConfig
123 _DefaultName = "isrAmpOffset"
124
125 def __init__(self, *args, **kwargs):
126 super().__init__(*args, **kwargs)
127 # Always load background subtask, even if doBackground=False;
128 # this allows for default plane bit masks to be defined.
129 self.makeSubtask("background")
130 if self.config.doDetection:
131 self.makeSubtask("detection")
132 # Initialize all of the instance variables here.
134
135 def run(self, exposure):
136 """Calculate amp offset values, determine corrective pedestals for each
137 amp, and update the input exposure in-place.
138
139 Parameters
140 ----------
141 exposure: `lsst.afw.image.Exposure`
142 Exposure to be corrected for amp offsets.
143 """
144
145 # Generate an exposure clone to work on and establish the bit mask.
146 exp = exposure.clone()
147 bitMask = exp.mask.getPlaneBitMask(self.background.config.ignoredPixelMask)
148 amps = exp.getDetector().getAmplifiers()
149
150 # Check that all amps have the same gemotry.
151 ampDims = [amp.getBBox().getDimensions() for amp in amps]
152 if not all(dim == ampDims[0] for dim in ampDims):
153 raise RuntimeError("All amps should have the same geometry.")
154 else:
155 # The zeroth amp is representative of all amps in the detector.
156 self.ampDims = ampDims[0]
157 # Dictionary mapping side numbers to interface lengths.
158 # See `getAmpAssociations()` for details about sides.
159 self.interfaceLengthLookupBySide = {i: self.ampDims[i % 2] for i in range(4)}
160
161 # Determine amplifier geometry.
162 ampWidths = {amp.getBBox().getWidth() for amp in amps}
163 ampHeights = {amp.getBBox().getHeight() for amp in amps}
164 if len(ampWidths) > 1 or len(ampHeights) > 1:
165 raise NotImplementedError(
166 "Amp offset correction is not yet implemented for detectors with differing amp sizes."
167 )
168
169 # Assuming all the amps have the same geometry.
170 self.shortAmpSide = np.min(ampDims[0])
171
172 # Check that the edge width and inset are not too large.
173 if self.config.ampEdgeWidth >= self.shortAmpSide - 2 * self.config.ampEdgeInset:
174 raise RuntimeError(
175 f"The edge width ({self.config.ampEdgeWidth}) plus insets ({self.config.ampEdgeInset}) "
176 f"exceed the amp's short side ({self.shortAmpSide}). This setup leads to incorrect results."
177 )
178
179 # Fit and subtract background.
180 if self.config.doBackground:
181 maskedImage = exp.getMaskedImage()
182 # Assuming all the detectors are the same.
183 nX = exp.getWidth() // (self.shortAmpSide * self.config.backgroundFractionSample) + 1
184 nY = exp.getHeight() // (self.shortAmpSide * self.config.backgroundFractionSample) + 1
185 # This ensures that the `binSize` is as large as possible,
186 # preventing background subtraction from inadvertently removing the
187 # amp offset signature. Here it's set to the shorter dimension of
188 # the amplifier by default (`backgroundFractionSample` = 1), which
189 # seems reasonable.
190 bg = self.background.fitBackground(maskedImage, nx=int(nX), ny=int(nY))
191 bgImage = bg.getImageF(self.background.config.algorithm, self.background.config.undersampleStyle)
192 maskedImage -= bgImage
193
194 # Detect sources and update cloned exposure mask planes in-place.
195 if self.config.doDetection:
196 schema = SourceTable.makeMinimalSchema()
197 table = SourceTable.make(schema)
198 # Detection sigma, used for smoothing and to grow detections, is
199 # normally measured from the PSF of the exposure. As the PSF hasn't
200 # been measured at this stage of processing, sigma is instead
201 # set to an approximate value here (which should be sufficient).
202 _ = self.detection.run(table=table, exposure=exp, sigma=2)
203
204 # Safety check: do any pixels remain for amp offset estimation?
205 if (exp.mask.array & bitMask).all():
206 self.log.warning(
207 "All pixels masked: cannot calculate any amp offset corrections. All pedestals are being set "
208 "to zero."
209 )
210 pedestals = np.zeros(len(amps))
211 else:
212 # Set up amp offset inputs.
213 im = exp.image
214 im.array[(exp.mask.array & bitMask) > 0] = np.nan
215
216 if self.config.ampEdgeWindowFrac > 1:
217 raise RuntimeError(
218 f"The specified fraction (`ampEdgeWindowFrac`={self.config.ampEdgeWindowFrac}) of the "
219 "edge length exceeds 1. This leads to complications downstream, after convolution in "
220 "the `getInterfaceOffset()` method. Please modify the `ampEdgeWindowFrac` value in the "
221 "config to be 1 or less and rerun."
222 )
223
224 # Obtain association and offset matrices.
225 A, sides = self.getAmpAssociations(amps)
226 B = self.getAmpOffsets(im, amps, A, sides)
227
228 # If least-squares minimization fails, convert NaNs to zeroes,
229 # ensuring that no values are erroneously added/subtracted.
230 pedestals = np.nan_to_num(np.linalg.lstsq(A, B, rcond=None)[0])
231
232 metadata = exposure.getMetadata()
233 for amp, pedestal in zip(amps, pedestals):
234 ampIm = exposure.image[amp.getBBox()].array
235 ampIm -= pedestal
236 ampName = amp.getName()
237 metadata.set(
238 f"LSST ISR AMPOFFSET PEDESTAL {ampName}",
239 float(pedestal),
240 f"Pedestal level subtracted from amp {ampName}",
241 )
242 self.log.info(f"amp pedestal values: {', '.join([f'{x:.4f}' for x in pedestals])}")
243
244 return Struct(pedestals=pedestals)
245
246 def getAmpAssociations(self, amps):
247 """Determine amp geometry and amp associations from a list of
248 amplifiers.
249
250 Parse an input list of amplifiers to determine the layout of amps
251 within a detector, and identify all amp sides (i.e., the
252 horizontal and vertical junctions between amps).
253
254 Returns a matrix with a shape corresponding to the geometry of the amps
255 in the detector.
256
257 Parameters
258 ----------
259 amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
260 List of amplifier objects used to deduce associations.
261
262 Returns
263 -------
264 ampAssociations : `numpy.ndarray`
265 An N x N matrix (N = number of amplifiers) that illustrates the
266 connections between amplifiers within the detector layout. Each row
267 and column index corresponds to the ampIds of a specific pair of
268 amplifiers, and the matrix elements indicate their associations as
269 follows:
270
271 * 0: No association
272 * -1: Association exists (direction specified in the ampSides
273 matrix)
274 * n >= 1: Diagonal elements indicate the number of neighboring
275 amplifiers for the corresponding ampId==row==column number.
276
277 ampSides : `numpy.ndarray`
278 An N x N matrix (N = the number of amplifiers) representing the amp
279 side information corresponding to the `ampAssociations`
280 matrix. The elements are integers defined as below:
281
282 * -1: No side due to no association or the same amp (diagonals)
283 * 0: Side on the bottom
284 * 1: Side on the right
285 * 2: Side on the top
286 * 3: Side on the left
287 """
288 xCenters = [amp.getBBox().getCenterX() for amp in amps]
289 yCenters = [amp.getBBox().getCenterY() for amp in amps]
290 xIndices = np.ceil(xCenters / np.min(xCenters) / 2).astype(int) - 1
291 yIndices = np.ceil(yCenters / np.min(yCenters) / 2).astype(int) - 1
292
293 nAmps = len(amps)
294 ampIds = np.zeros((len(set(yIndices)), len(set(xIndices))), dtype=int)
295
296 for ampId, xIndex, yIndex in zip(np.arange(nAmps), xIndices, yIndices):
297 ampIds[yIndex, xIndex] = ampId
298
299 ampAssociations = np.zeros((nAmps, nAmps), dtype=int)
300 ampSides = np.full_like(ampAssociations, -1)
301
302 for ampId in ampIds.ravel():
303 neighbors, sides = self.getNeighbors(ampIds, ampId)
304 interfaceWeights = (
305 1
306 if not self.config.applyWeights
307 else np.array([self.interfaceLengthLookupBySide[side] for side in sides])
308 )
309 ampAssociations[ampId, neighbors] = -1 * interfaceWeights
310 ampSides[ampId, neighbors] = sides
311 ampAssociations[ampId, ampId] = -ampAssociations[ampId].sum()
312
313 if ampAssociations.sum() != 0:
314 raise RuntimeError("The `ampAssociations` array does not sum to zero.")
315
316 if not np.all(ampAssociations == ampAssociations.T):
317 raise RuntimeError("The `ampAssociations` is not symmetric about the diagonal.")
318
319 self.log.debug("amp associations:\n%s", ampAssociations)
320 self.log.debug("amp sides:\n%s", ampSides)
321
322 return ampAssociations, ampSides
323
324 def getNeighbors(self, ampIds, ampId):
325 """Get the neighbor amplifiers and their sides for a given
326 amplifier.
327
328 Parameters
329 ----------
330 ampIds : `numpy.ndarray`
331 Matrix with amp side association information.
332 ampId : `int`
333 The amplifier ID for which neighbor amplifiers and side IDs
334 are to be found.
335
336 Returns
337 -------
338 neighbors : `list` [`int`]
339 List of neighbor amplifier IDs.
340 sides : `list` [`int`]
341 List of side IDs, with each ID corresponding to its respective
342 neighbor amplifier.
343 """
344 m, n = ampIds.shape
345 r, c = np.ravel(np.where(ampIds == ampId))
346 neighbors, sides = [], []
347 sideLookup = {
348 0: (r + 1, c),
349 1: (r, c + 1),
350 2: (r - 1, c),
351 3: (r, c - 1),
352 }
353 for side, (row, column) in sideLookup.items():
354 if 0 <= row < m and 0 <= column < n:
355 neighbors.append(ampIds[row][column])
356 sides.append(side)
357 return neighbors, sides
358
359 def getAmpOffsets(self, im, amps, associations, sides):
360 """Calculate the amp offsets for all amplifiers.
361
362 Parameters
363 ----------
364 im : `lsst.afw.image._image.ImageF`
365 Amplifier image to extract data from.
366 amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
367 List of amplifier objects.
368 associations : numpy.ndarray
369 An N x N matrix containing amp association information, where N is
370 the number of amplifiers.
371 sides : numpy.ndarray
372 An N x N matrix containing amp side information, where N is the
373 number of amplifiers.
374
375 Returns
376 -------
377 ampsOffsets : `numpy.ndarray`
378 1D float array containing the calculated amp offsets for all
379 amplifiers.
380 """
381 ampsOffsets = np.zeros(len(amps))
382 ampsEdges = self.getAmpEdges(im, amps, sides)
383 interfaceOffsetLookup = {}
384
385 for ampId, ampAssociations in enumerate(associations):
386 ampNeighbors = np.ravel(np.where(ampAssociations < 0))
387 for ampNeighbor in ampNeighbors:
388 ampSide = sides[ampId][ampNeighbor]
389 interfaceWeight = (
390 1 if not self.config.applyWeights else self.interfaceLengthLookupBySide[ampSide]
391 )
392 edgeA = ampsEdges[ampId][ampSide]
393 edgeB = ampsEdges[ampNeighbor][(ampSide + 2) % 4]
394 if ampId < ampNeighbor:
395 interfaceOffset = self.getInterfaceOffset(ampId, ampNeighbor, edgeA, edgeB)
396 interfaceOffsetLookup[f"{ampId}{ampNeighbor}"] = interfaceOffset
397 else:
398 interfaceOffset = -interfaceOffsetLookup[f"{ampNeighbor}{ampId}"]
399 ampsOffsets[ampId] += interfaceWeight * interfaceOffset
400 return ampsOffsets
401
402 def getAmpEdges(self, im, amps, ampSides):
403 """Calculate the amp edges for all amplifiers.
404
405 Parameters
406 ----------
407 im : `lsst.afw.image._image.ImageF`
408 Amplifier image to extract data from.
409 amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
410 List of amplifier objects.
411 ampSides : `numpy.ndarray`
412 An N x N matrix containing amp side information, where N is the
413 number of amplifiers.
414
415 Returns
416 -------
417 ampEdges : `dict` [`int`, `dict` [`int`, `numpy.ndarray`]]
418 A dictionary containing amp edge(s) for each amplifier,
419 corresponding to one or more potential sides, where each edge is
420 associated with a side. The outer dictionary has integer keys
421 representing amplifier IDs, and the inner dictionary has integer
422 keys representing side IDs for each amplifier and values that are
423 1D arrays of floats representing the 1D medianified strips from the
424 amp image, referred to as "amp edge":
425 {ampID: {sideID: numpy.ndarray}, ...}
426 """
427 ampEdgeOuter = self.config.ampEdgeInset + self.config.ampEdgeWidth
428 ampEdges = {}
429 slice_map = {
430 0: (slice(-ampEdgeOuter, -self.config.ampEdgeInset), slice(None)),
431 1: (slice(None), slice(-ampEdgeOuter, -self.config.ampEdgeInset)),
432 2: (slice(self.config.ampEdgeInset, ampEdgeOuter), slice(None)),
433 3: (slice(None), slice(self.config.ampEdgeInset, ampEdgeOuter)),
434 }
435 for ampId, (amp, ampSides) in enumerate(zip(amps, ampSides)):
436 ampEdges[ampId] = {}
437 ampIm = im[amp.getBBox()].array
438 # Loop over identified sides.
439 for ampSide in ampSides:
440 if ampSide < 0:
441 continue
442 strip = ampIm[slice_map[ampSide]]
443 # Catch warnings to prevent all-NaN slice RuntimeWarning.
444 with warnings.catch_warnings():
445 warnings.filterwarnings("ignore", r"All-NaN (slice|axis) encountered")
446 ampEdges[ampId][ampSide] = np.nanmedian(strip, axis=ampSide % 2) # 1D medianified strip
447 return ampEdges
448
449 def getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB):
450 """Calculate the amp offset for a given interface between two
451 amplifiers.
452
453 Parameters
454 ----------
455 ampIdA : int
456 ID of the first amplifier.
457 ampIdB : int
458 ID of the second amplifier.
459 edgeA : numpy.ndarray
460 Amp edge for the first amplifier.
461 edgeB : numpy.ndarray
462 Amp edge for the second amplifier.
463
464 Returns
465 -------
466 interfaceOffset : float
467 The calculated amp offset value for the given interface between
468 amps A and B.
469 """
470 interfaceId = f"{ampIdA}{ampIdB}"
471 sctrl = StatisticsControl()
472 # NOTE: Taking the difference with the order below fixes the sign flip
473 # in the B matrix.
474 edgeDiff = edgeA - edgeB
475 window = int(self.config.ampEdgeWindowFrac * len(edgeDiff))
476 # Compute rolling averages.
477 edgeDiffSum = np.convolve(np.nan_to_num(edgeDiff), np.ones(window), "same")
478 edgeDiffNum = np.convolve(~np.isnan(edgeDiff), np.ones(window), "same")
479 edgeDiffAvg = edgeDiffSum / np.clip(edgeDiffNum, 1, None)
480 edgeDiffAvg[np.isnan(edgeDiff)] = np.nan
481 # Take clipped mean of rolling average data as amp offset value.
482 interfaceOffset = makeStatistics(edgeDiffAvg, MEANCLIP, sctrl).getValue()
483 # Perform a couple of do-no-harm safety checks:
484 # a) The fraction of unmasked pixel rows is > ampEdgeMinFrac,
485 # b) The absolute offset ADU value is < ampEdgeMaxOffset.
486 ampEdgeGoodFrac = 1 - (np.sum(np.isnan(edgeDiffAvg)) / len(edgeDiffAvg))
487 minFracFail = ampEdgeGoodFrac < self.config.ampEdgeMinFrac
488 maxOffsetFail = np.abs(interfaceOffset) > self.config.ampEdgeMaxOffset
489 if minFracFail or maxOffsetFail:
490 interfaceOffset = 0
491 if minFracFail:
492 self.log.warning(
493 f"The fraction of unmasked pixels for amp interface {interfaceId} is below the threshold "
494 f"({ampEdgeGoodFrac:.2f} < {self.config.ampEdgeMinFrac}). Setting the interface offset "
495 f"to {interfaceOffset}."
496 )
497 if maxOffsetFail:
498 self.log.warning(
499 "The absolute offset value exceeds the limit "
500 f"({np.abs(interfaceOffset):.2f} > {self.config.ampEdgeMaxOffset} ADU). Setting the "
501 f"interface offset to {interfaceOffset}."
502 )
503 self.log.debug(
504 f"amp interface {interfaceId} : "
505 f"viable edge difference frac = {ampEdgeGoodFrac}, "
506 f"interface offset = {interfaceOffset:.3f}"
507 )
508 return interfaceOffset
Pass parameters to a Statistics object.
Definition Statistics.h:83
getAmpOffsets(self, im, amps, associations, sides)
Definition ampOffset.py:359
__init__(self, *args, **kwargs)
Definition ampOffset.py:125
getAmpEdges(self, im, amps, ampSides)
Definition ampOffset.py:402
getNeighbors(self, ampIds, ampId)
Definition ampOffset.py:324
getInterfaceOffset(self, ampIdA, ampIdB, edgeA, edgeB)
Definition ampOffset.py:449
daf::base::PropertySet * set
Definition fits.cc:931