Inheritance diagram for lsst.ip.isr.ampOffset.AmpOffsetTask:
Public Member Functions | |
| __init__ (self, *args, **kwargs) | |
| run (self, exposure) | |
| getAmpAssociations (self, amps) | |
| getNeighbors (self, ampIds, ampId) | |
| getAmpOffsets (self, im, amps, associations, sides) | |
| getAmpEdges (self, im, amps, ampSides) | |
| getInterfaceOffset (self, ampIdA, ampIdB, edgeA, edgeB) | |
Public Attributes | |
| shortAmpSide | |
| ampDims | |
| interfaceLengthLookupBySide | |
Static Public Attributes | |
| ConfigClass = AmpOffsetConfig | |
Static Protected Attributes | |
| str | _DefaultName = "isrAmpOffset" |
Detailed Description
Calculate and apply amp offset corrections to an exposure.
Definition at line 119 of file ampOffset.py.
Constructor & Destructor Documentation
◆ __init__()
| lsst.ip.isr.ampOffset.AmpOffsetTask.__init__ | ( | self, | |
| * | args, | ||
| ** | kwargs ) |
Definition at line 125 of file ampOffset.py.
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.
133 self.shortAmpSide = 0
134
Member Function Documentation
◆ getAmpAssociations()
| lsst.ip.isr.ampOffset.AmpOffsetTask.getAmpAssociations | ( | self, | |
| amps ) |
Determine amp geometry and amp associations from a list of
amplifiers.
Parse an input list of amplifiers to determine the layout of amps
within a detector, and identify all amp sides (i.e., the
horizontal and vertical junctions between amps).
Returns a matrix with a shape corresponding to the geometry of the amps
in the detector.
Parameters
----------
amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
List of amplifier objects used to deduce associations.
Returns
-------
ampAssociations : `numpy.ndarray`
An N x N matrix (N = number of amplifiers) that illustrates the
connections between amplifiers within the detector layout. Each row
and column index corresponds to the ampIds of a specific pair of
amplifiers, and the matrix elements indicate their associations as
follows:
* 0: No association
* -1: Association exists (direction specified in the ampSides
matrix)
* n >= 1: Diagonal elements indicate the number of neighboring
amplifiers for the corresponding ampId==row==column number.
ampSides : `numpy.ndarray`
An N x N matrix (N = the number of amplifiers) representing the amp
side information corresponding to the `ampAssociations`
matrix. The elements are integers defined as below:
* -1: No side due to no association or the same amp (diagonals)
* 0: Side on the bottom
* 1: Side on the right
* 2: Side on the top
* 3: Side on the left
Definition at line 246 of file ampOffset.py.
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)
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
◆ getAmpEdges()
| lsst.ip.isr.ampOffset.AmpOffsetTask.getAmpEdges | ( | self, | |
| im, | |||
| amps, | |||
| ampSides ) |
Calculate the amp edges for all amplifiers.
Parameters
----------
im : `lsst.afw.image._image.ImageF`
Amplifier image to extract data from.
amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
List of amplifier objects.
ampSides : `numpy.ndarray`
An N x N matrix containing amp side information, where N is the
number of amplifiers.
Returns
-------
ampEdges : `dict` [`int`, `dict` [`int`, `numpy.ndarray`]]
A dictionary containing amp edge(s) for each amplifier,
corresponding to one or more potential sides, where each edge is
associated with a side. The outer dictionary has integer keys
representing amplifier IDs, and the inner dictionary has integer
keys representing side IDs for each amplifier and values that are
1D arrays of floats representing the 1D medianified strips from the
amp image, referred to as "amp edge":
{ampID: {sideID: numpy.ndarray}, ...}
Definition at line 402 of file ampOffset.py.
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
◆ getAmpOffsets()
| lsst.ip.isr.ampOffset.AmpOffsetTask.getAmpOffsets | ( | self, | |
| im, | |||
| amps, | |||
| associations, | |||
| sides ) |
Calculate the amp offsets for all amplifiers.
Parameters
----------
im : `lsst.afw.image._image.ImageF`
Amplifier image to extract data from.
amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
List of amplifier objects.
associations : numpy.ndarray
An N x N matrix containing amp association information, where N is
the number of amplifiers.
sides : numpy.ndarray
An N x N matrix containing amp side information, where N is the
number of amplifiers.
Returns
-------
ampsOffsets : `numpy.ndarray`
1D float array containing the calculated amp offsets for all
amplifiers.
Definition at line 359 of file ampOffset.py.
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
◆ getInterfaceOffset()
| lsst.ip.isr.ampOffset.AmpOffsetTask.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.
Definition at line 449 of file ampOffset.py.
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
◆ getNeighbors()
| lsst.ip.isr.ampOffset.AmpOffsetTask.getNeighbors | ( | self, | |
| ampIds, | |||
| ampId ) |
Get the neighbor amplifiers and their sides for a given
amplifier.
Parameters
----------
ampIds : `numpy.ndarray`
Matrix with amp side association information.
ampId : `int`
The amplifier ID for which neighbor amplifiers and side IDs
are to be found.
Returns
-------
neighbors : `list` [`int`]
List of neighbor amplifier IDs.
sides : `list` [`int`]
List of side IDs, with each ID corresponding to its respective
neighbor amplifier.
Definition at line 324 of file ampOffset.py.
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
◆ run()
| lsst.ip.isr.ampOffset.AmpOffsetTask.run | ( | self, | |
| exposure ) |
Calculate amp offset values, determine corrective pedestals for each
amp, and update the input exposure in-place.
Parameters
----------
exposure: `lsst.afw.image.Exposure`
Exposure to be corrected for amp offsets.
Definition at line 135 of file ampOffset.py.
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
Member Data Documentation
◆ _DefaultName
|
staticprotected |
Definition at line 123 of file ampOffset.py.
◆ ampDims
| lsst.ip.isr.ampOffset.AmpOffsetTask.ampDims |
Definition at line 156 of file ampOffset.py.
◆ ConfigClass
|
static |
Definition at line 122 of file ampOffset.py.
◆ interfaceLengthLookupBySide
| lsst.ip.isr.ampOffset.AmpOffsetTask.interfaceLengthLookupBySide |
Definition at line 159 of file ampOffset.py.
◆ shortAmpSide
| lsst.ip.isr.ampOffset.AmpOffsetTask.shortAmpSide |
Definition at line 133 of file ampOffset.py.
The documentation for this class was generated from the following file:
- /j/snowflake/release/lsstsw/stack/lsst-scipipe-8.0.0/Linux64/ip_isr/g02d81e74bb+86cf3d8bc9/python/lsst/ip/isr/ampOffset.py
