LSST Applications g063fba187b+cac8b7c890,g0f08755f38+6aee506743,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+b4475c5878,g1dcb35cd9c+8f9bc1652e,g20f6ffc8e0+6aee506743,g217e2c1bcf+73dee94bd0,g28da252d5a+1f19c529b9,g2bbee38e9b+3f2625acfc,g2bc492864f+3f2625acfc,g3156d2b45e+6e55a43351,g32e5bea42b+1bb94961c2,g347aa1857d+3f2625acfc,g35bb328faa+a8ce1bb630,g3a166c0a6a+3f2625acfc,g3e281a1b8c+c5dd892a6c,g3e8969e208+a8ce1bb630,g414038480c+5927e1bc1e,g41af890bb2+8a9e676b2a,g7af13505b9+809c143d88,g80478fca09+6ef8b1810f,g82479be7b0+f568feb641,g858d7b2824+6aee506743,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+2903d499ea,gb58c049af0+d64f4d3760,gc28159a63d+3f2625acfc,gcab2d0539d+b12535109e,gcf0d15dbbd+46a3f46ba9,gda6a2b7d83+46a3f46ba9,gdaeeff99f8+1711a396fd,ge79ae78c31+3f2625acfc,gef2f8181fd+0a71e47438,gf0baf85859+c1f95f4921,gfa517265be+6aee506743,gfa999e8aa5+17cd334064,w.2024.51
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 doWindowSmoothing = Field(
80 doc="Smooth amp edge differences by taking a rolling average.",
81 dtype=bool,
82 default=True,
83 )
84 ampEdgeWindowFrac = Field(
85 doc="Fraction of the amp edge lengths utilized as the sliding window for generating rolling average "
86 "amp offset values. It should be reconfigured for every instrument (HSC, LSSTCam, etc.) and should "
87 "not exceed 1. If not provided, it defaults to the fraction that recovers the pixel size of the "
88 "sliding window used in obs_subaru for compatibility with existing HSC data. Only relevant if "
89 "`doWindowSmoothing` is set to True.",
90 dtype=float,
91 default=512 / 4176,
92 )
93 doBackground = Field(
94 doc="Estimate and subtract background prior to amp offset estimation?",
95 dtype=bool,
96 default=True,
97 )
98 background = ConfigurableField(
99 doc="An initial background estimation step run prior to amp offset calculation.",
100 target=SubtractBackgroundTask,
101 )
102 backgroundFractionSample = Field(
103 doc="The fraction of the shorter side of the amplifier used for background binning.",
104 dtype=float,
105 default=1.0,
106 )
107 doDetection = Field(
108 doc="Detect sources and update cloned exposure prior to amp offset estimation?",
109 dtype=bool,
110 default=True,
111 )
112 detection = ConfigurableField(
113 doc="Source detection to add temporary detection footprints prior to amp offset calculation.",
114 target=SourceDetectionTask,
115 )
116 applyWeights = Field(
117 doc="Weights the amp offset calculation by the length of the interface between amplifiers. Applying "
118 "weights does not affect outcomes for amplifiers in a 2D grid with square-shaped amplifiers or in "
119 "any 1D layout on a detector, regardless of whether the amplifiers are square.",
120 dtype=bool,
121 default=True,
122 )
123 doApplyAmpOffset = Field(
124 doc="Apply amp offset corrections to the input exposure?",
125 dtype=bool,
126 default=False,
127 )
128
129
130class AmpOffsetTask(Task):
131 """Calculate and apply amp offset corrections to an exposure."""
132
133 ConfigClass = AmpOffsetConfig
134 _DefaultName = "isrAmpOffset"
135
136 def __init__(self, *args, **kwargs):
137 super().__init__(*args, **kwargs)
138 # Always load background subtask, even if doBackground=False;
139 # this allows for default plane bit masks to be defined.
140 self.makeSubtask("background")
141 if self.config.doDetection:
142 self.makeSubtask("detection")
143 # Initialize all of the instance variables here.
145
146 def run(self, exposure):
147 """Calculate amp offset values, determine corrective pedestals for each
148 amp, and update the input exposure in-place.
149
150 Parameters
151 ----------
152 exposure: `lsst.afw.image.Exposure`
153 Exposure to be corrected for amp offsets.
154 """
155
156 # Generate an exposure clone to work on and establish the bit mask.
157 exp = exposure.clone()
158 bitMask = exp.mask.getPlaneBitMask(self.background.config.ignoredPixelMask)
159 amps = exp.getDetector().getAmplifiers()
160
161 # Check that all amps have the same gemotry.
162 ampDims = [amp.getBBox().getDimensions() for amp in amps]
163 if not all(dim == ampDims[0] for dim in ampDims):
164 raise RuntimeError("All amps should have the same geometry.")
165 else:
166 # The zeroth amp is representative of all amps in the detector.
167 self.ampDims = ampDims[0]
168 # Dictionary mapping side numbers to interface lengths.
169 # See `getAmpAssociations()` for details about sides.
170 self.interfaceLengthLookupBySide = {i: self.ampDims[i % 2] for i in range(4)}
171
172 # Determine amplifier geometry.
173 ampWidths = {amp.getBBox().getWidth() for amp in amps}
174 ampHeights = {amp.getBBox().getHeight() for amp in amps}
175 if len(ampWidths) > 1 or len(ampHeights) > 1:
176 raise NotImplementedError(
177 "Amp offset correction is not yet implemented for detectors with differing amp sizes."
178 )
179
180 # Assuming all the amps have the same geometry.
181 self.shortAmpSide = np.min(ampDims[0])
182
183 # Check that the edge width and inset are not too large.
184 if self.config.ampEdgeWidth >= self.shortAmpSide - 2 * self.config.ampEdgeInset:
185 raise RuntimeError(
186 f"The edge width ({self.config.ampEdgeWidth}) plus insets ({self.config.ampEdgeInset}) "
187 f"exceed the amp's short side ({self.shortAmpSide}). This setup leads to incorrect results."
188 )
189
190 # Fit and subtract background.
191 if self.config.doBackground:
192 maskedImage = exp.getMaskedImage()
193 # Assuming all the detectors are the same.
194 nX = exp.getWidth() // (self.shortAmpSide * self.config.backgroundFractionSample) + 1
195 nY = exp.getHeight() // (self.shortAmpSide * self.config.backgroundFractionSample) + 1
196 # This ensures that the `binSize` is as large as possible,
197 # preventing background subtraction from inadvertently removing the
198 # amp offset signature. Here it's set to the shorter dimension of
199 # the amplifier by default (`backgroundFractionSample` = 1), which
200 # seems reasonable.
201 bg = self.background.fitBackground(maskedImage, nx=int(nX), ny=int(nY))
202 bgImage = bg.getImageF(self.background.config.algorithm, self.background.config.undersampleStyle)
203 maskedImage -= bgImage
204
205 # Detect sources and update cloned exposure mask planes in-place.
206 if self.config.doDetection:
207 schema = SourceTable.makeMinimalSchema()
208 table = SourceTable.make(schema)
209 # Detection sigma, used for smoothing and to grow detections, is
210 # normally measured from the PSF of the exposure. As the PSF hasn't
211 # been measured at this stage of processing, sigma is instead
212 # set to an approximate value here (which should be sufficient).
213 _ = self.detection.run(table=table, exposure=exp, sigma=2)
214
215 # Safety check: do any pixels remain for amp offset estimation?
216 if (exp.mask.array & bitMask).all():
217 log_fn = self.log.warning if self.config.doApplyAmpOffset else self.log.info
218 log_fn(
219 "All pixels masked: cannot calculate any amp offset corrections. "
220 "All pedestals are being set to zero."
221 )
222 pedestals = np.zeros(len(amps))
223 else:
224 # Set up amp offset inputs.
225 im = exp.image
226 im.array[(exp.mask.array & bitMask) > 0] = np.nan
227
228 if self.config.ampEdgeWindowFrac > 1:
229 raise RuntimeError(
230 f"The specified fraction (`ampEdgeWindowFrac`={self.config.ampEdgeWindowFrac}) of the "
231 "edge length exceeds 1. This leads to complications downstream, after convolution in "
232 "the `getInterfaceOffset()` method. Please modify the `ampEdgeWindowFrac` value in the "
233 "config to be 1 or less and rerun."
234 )
235
236 # Obtain association and offset matrices.
237 A, sides = self.getAmpAssociations(amps)
238 B, interfaceOffsetDict = self.getAmpOffsets(im, amps, A, sides)
239
240 # If least-squares minimization fails, convert NaNs to zeroes,
241 # ensuring that no values are erroneously added/subtracted.
242 pedestals = np.nan_to_num(np.linalg.lstsq(A, B, rcond=None)[0])
243
244 metadata = exposure.getMetadata() # Exposure metadata.
245 self.metadata["AMPOFFSET_PEDESTALS"] = {} # Task metadata.
246 ampNames = [amp.getName() for amp in amps]
247
248 # Add the amp interface offsets to the exposure metadata.
249 for interfaceId, interfaceOffset in interfaceOffsetDict.items():
250 metadata.set(
251 f"LSST ISR AMPOFFSET INTERFACEOFFSET {interfaceId}",
252 float(interfaceOffset),
253 f"Raw amp interface offset calculated for {interfaceId}",
254 )
255
256 for ampName, amp, pedestal in zip(ampNames, amps, pedestals):
257 # Add the amp pedestal to the exposure metadata.
258 metadata.set(
259 f"LSST ISR AMPOFFSET PEDESTAL {ampName}",
260 float(pedestal),
261 f"Pedestal level calculated for amp {ampName}",
262 )
263 if self.config.doApplyAmpOffset:
264 ampIm = exposure.image[amp.getBBox()].array
265 ampIm -= pedestal
266 # Add the amp pedestal to the "Task" metadata as well.
267 # Needed for Sasquatch/Chronograf!
268 self.metadata["AMPOFFSET_PEDESTALS"][ampName] = float(pedestal)
269 if self.config.doApplyAmpOffset:
270 status = "subtracted from exposure"
271 metadata.set("LSST ISR AMPOFFSET PEDESTAL SUBTRACTED", True, "Amp pedestals have been subtracted")
272 else:
273 status = "not subtracted from exposure"
274 metadata.set(
275 "LSST ISR AMPOFFSET PEDESTAL SUBTRACTED", False, "Amp pedestals have not been subtracted"
276 )
277 ampPedestalReport = ", ".join(
278 [f"{ampName}: {ampPedestal:.4f}" for (ampName, ampPedestal) in zip(ampNames, pedestals)]
279 )
280 self.log.info(f"Amp pedestal values ({status}): {ampPedestalReport}")
281
282 return Struct(pedestals=pedestals)
283
284 def getAmpAssociations(self, amps):
285 """Determine amp geometry and amp associations from a list of
286 amplifiers.
287
288 Parse an input list of amplifiers to determine the layout of amps
289 within a detector, and identify all amp sides (i.e., the
290 horizontal and vertical junctions between amps).
291
292 Returns a matrix with a shape corresponding to the geometry of the amps
293 in the detector.
294
295 Parameters
296 ----------
297 amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
298 List of amplifier objects used to deduce associations.
299
300 Returns
301 -------
302 ampAssociations : `numpy.ndarray`
303 An N x N matrix (N = number of amplifiers) that illustrates the
304 connections between amplifiers within the detector layout. Each row
305 and column index corresponds to the ampIds of a specific pair of
306 amplifiers, and the matrix elements indicate their associations as
307 follows:
308
309 * 0: No association
310 * -1: Association exists (direction specified in the ampSides
311 matrix)
312 * n >= 1: Diagonal elements indicate the number of neighboring
313 amplifiers for the corresponding ampId==row==column number.
314
315 ampSides : `numpy.ndarray`
316 An N x N matrix (N = the number of amplifiers) representing the amp
317 side information corresponding to the `ampAssociations`
318 matrix. The elements are integers defined as below:
319
320 * -1: No side due to no association or the same amp (diagonals)
321 * 0: Side on the bottom
322 * 1: Side on the right
323 * 2: Side on the top
324 * 3: Side on the left
325 """
326 xCenters = [amp.getBBox().getCenterX() for amp in amps]
327 yCenters = [amp.getBBox().getCenterY() for amp in amps]
328 xIndices = np.ceil(xCenters / np.min(xCenters) / 2).astype(int) - 1
329 yIndices = np.ceil(yCenters / np.min(yCenters) / 2).astype(int) - 1
330
331 nAmps = len(amps)
332 ampIds = np.zeros((len(set(yIndices)), len(set(xIndices))), dtype=int)
333
334 for ampId, xIndex, yIndex in zip(np.arange(nAmps), xIndices, yIndices):
335 ampIds[yIndex, xIndex] = ampId
336
337 ampAssociations = np.zeros((nAmps, nAmps), dtype=int)
338 ampSides = np.full_like(ampAssociations, -1)
339
340 for ampId in ampIds.ravel():
341 neighbors, sides = self.getNeighbors(ampIds, ampId)
342 interfaceWeights = (
343 1
344 if not self.config.applyWeights
345 else np.array([self.interfaceLengthLookupBySide[side] for side in sides])
346 )
347 ampAssociations[ampId, neighbors] = -1 * interfaceWeights
348 ampSides[ampId, neighbors] = sides
349 ampAssociations[ampId, ampId] = -ampAssociations[ampId].sum()
350
351 if ampAssociations.sum() != 0:
352 raise RuntimeError("The `ampAssociations` array does not sum to zero.")
353
354 if not np.all(ampAssociations == ampAssociations.T):
355 raise RuntimeError("The `ampAssociations` is not symmetric about the diagonal.")
356
357 with np.printoptions(linewidth=200):
358 self.log.debug("amp associations:\n%s", ampAssociations)
359 self.log.debug("amp sides:\n%s", ampSides)
360
361 return ampAssociations, ampSides
362
363 def getNeighbors(self, ampIds, ampId):
364 """Get the neighbor amplifiers and their sides for a given
365 amplifier.
366
367 Parameters
368 ----------
369 ampIds : `numpy.ndarray`
370 Matrix with amp side association information.
371 ampId : `int`
372 The amplifier ID for which neighbor amplifiers and side IDs
373 are to be found.
374
375 Returns
376 -------
377 neighbors : `list` [`int`]
378 List of neighbor amplifier IDs.
379 sides : `list` [`int`]
380 List of side IDs, with each ID corresponding to its respective
381 neighbor amplifier.
382 """
383 m, n = ampIds.shape
384 r, c = np.ravel(np.where(ampIds == ampId))
385 neighbors, sides = [], []
386 sideLookup = {
387 0: (r + 1, c),
388 1: (r, c + 1),
389 2: (r - 1, c),
390 3: (r, c - 1),
391 }
392 for side, (row, column) in sideLookup.items():
393 if 0 <= row < m and 0 <= column < n:
394 neighbors.append(ampIds[row][column])
395 sides.append(side)
396 return neighbors, sides
397
398 def getAmpOffsets(self, im, amps, associations, sides):
399 """Calculate the amp offsets for all amplifiers.
400
401 Parameters
402 ----------
403 im : `lsst.afw.image._image.ImageF`
404 Amplifier image to extract data from.
405 amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
406 List of amplifier objects.
407 associations : numpy.ndarray
408 An N x N matrix containing amp association information, where N is
409 the number of amplifiers.
410 sides : numpy.ndarray
411 An N x N matrix containing amp side information, where N is the
412 number of amplifiers.
413
414 Returns
415 -------
416 ampsOffsets : `numpy.ndarray`
417 1D float array containing the calculated amp offsets for all
418 amplifiers.
419 interfaceOffsetDict : `dict` [`str`, `float`]
420 Dictionary mapping interface IDs to their corresponding raw
421 (uncapped) offset values.
422 """
423 ampsOffsets = np.zeros(len(amps))
424 ampsEdges = self.getAmpEdges(im, amps, sides)
425 ampsNames = [amp.getName() for amp in amps]
426 interfaceOffsetLookup = {}
427
428 interfaceIds = []
429 interfaceOffsetOriginals = []
430 ampEdgeGoodFracs = []
431 minFracFails = []
432 maxOffsetFails = []
433 for ampId, ampAssociations in enumerate(associations):
434 ampNeighbors = np.ravel(np.where(ampAssociations < 0))
435 for ampNeighbor in ampNeighbors:
436 ampSide = sides[ampId][ampNeighbor]
437 interfaceWeight = (
438 1 if not self.config.applyWeights else self.interfaceLengthLookupBySide[ampSide]
439 )
440 edgeA = ampsEdges[ampId][ampSide]
441 edgeB = ampsEdges[ampNeighbor][(ampSide + 2) % 4]
442 if ampId < ampNeighbor:
443 (
444 interfaceId,
445 interfaceOffset,
446 interfaceOffsetOriginal,
447 ampEdgeGoodFrac,
448 minFracFail,
449 maxOffsetFail,
450 ) = self.getInterfaceOffset(ampsNames[ampId], ampsNames[ampNeighbor], edgeA, edgeB)
451 interfaceIds.append(interfaceId)
452 interfaceOffsetOriginals.append(interfaceOffsetOriginal)
453 ampEdgeGoodFracs.append(ampEdgeGoodFrac)
454 minFracFails.append(minFracFail)
455 maxOffsetFails.append(maxOffsetFail)
456 interfaceOffsetLookup[f"{ampId:02d}:{ampNeighbor:02d}"] = interfaceOffset
457 else:
458 interfaceOffset = -interfaceOffsetLookup[f"{ampNeighbor:02d}:{ampId:02d}"]
459 ampsOffsets[ampId] += interfaceWeight * interfaceOffset
460 if interfaceOffsetOriginals:
461 self.log.debug(
462 "Raw (uncapped) amp offset values for all interfaces: %s",
463 ", ".join(
464 [
465 f"{interfaceId}={interfaceOffset:0.2f}"
466 for interfaceId, interfaceOffset in zip(interfaceIds, interfaceOffsetOriginals)
467 ]
468 ),
469 )
470 quartile_summary = np.nanpercentile(interfaceOffsetOriginals, [0, 25, 50, 75, 100])
471 self.log.info(
472 "Raw amp offset quartile summary for all interfaces (min, Q1, Q2, Q3, max): "
473 "%.4f, %.4f, %.4f, %.4f, %.4f",
474 *quartile_summary,
475 )
476 log_fn = self.log.warning if self.config.doApplyAmpOffset else self.log.info
477 if any(minFracFails):
478 log_fn(
479 "The fraction of unmasked edge pixels for the following amp interfaces is below the "
480 "configured threshold (%s): %s",
481 self.config.ampEdgeMinFrac,
482 ", ".join(
483 [
484 f"{interfaceId} ({ampEdgeGoodFrac:0.2f})"
485 for interfaceId, ampEdgeGoodFrac, minFracFail in zip(
486 interfaceIds, ampEdgeGoodFracs, minFracFails
487 )
488 if minFracFail
489 ]
490 ),
491 )
492 if any(maxOffsetFails):
493 log_fn(
494 "Absolute amp offsets exceed the configured maximum (%s) and have been set to zero for the "
495 "following amp interfaces: %s",
496 self.config.ampEdgeMaxOffset,
497 ", ".join(
498 [
499 f"{interfaceId}={np.abs(interfaceOffset):0.2f}"
500 for interfaceId, interfaceOffset, maxOffsetFail in zip(
501 interfaceIds, interfaceOffsetOriginals, maxOffsetFails
502 )
503 if maxOffsetFail
504 ]
505 ),
506 )
507
508 # Pair each interface ID with its corresponding original offset.
509 interfaceOffsetDict = dict(zip(interfaceIds, interfaceOffsetOriginals))
510
511 return ampsOffsets, interfaceOffsetDict
512
513 def getAmpEdges(self, im, amps, ampSides):
514 """Calculate the amp edges for all amplifiers.
515
516 Parameters
517 ----------
518 im : `lsst.afw.image._image.ImageF`
519 Amplifier image to extract data from.
520 amps : `list` [`lsst.afw.cameraGeom.Amplifier`]
521 List of amplifier objects.
522 ampSides : `numpy.ndarray`
523 An N x N matrix containing amp side information, where N is the
524 number of amplifiers.
525
526 Returns
527 -------
528 ampEdges : `dict` [`int`, `dict` [`int`, `numpy.ndarray`]]
529 A dictionary containing amp edge(s) for each amplifier,
530 corresponding to one or more potential sides, where each edge is
531 associated with a side. The outer dictionary has integer keys
532 representing amplifier IDs, and the inner dictionary has integer
533 keys representing side IDs for each amplifier and values that are
534 1D arrays of floats representing the 1D medianified strips from the
535 amp image, referred to as "amp edge":
536 {ampID: {sideID: numpy.ndarray}, ...}
537 """
538 ampEdgeOuter = self.config.ampEdgeInset + self.config.ampEdgeWidth
539 ampEdges = {}
540 slice_map = {
541 0: (slice(-ampEdgeOuter, -self.config.ampEdgeInset), slice(None)),
542 1: (slice(None), slice(-ampEdgeOuter, -self.config.ampEdgeInset)),
543 2: (slice(self.config.ampEdgeInset, ampEdgeOuter), slice(None)),
544 3: (slice(None), slice(self.config.ampEdgeInset, ampEdgeOuter)),
545 }
546 for ampId, (amp, ampSides) in enumerate(zip(amps, ampSides)):
547 ampEdges[ampId] = {}
548 ampIm = im[amp.getBBox()].array
549 # Loop over identified sides.
550 for ampSide in ampSides:
551 if ampSide < 0:
552 continue
553 strip = ampIm[slice_map[ampSide]]
554 # Catch warnings to prevent all-NaN slice RuntimeWarning.
555 with warnings.catch_warnings():
556 warnings.filterwarnings("ignore", r"All-NaN (slice|axis) encountered")
557 ampEdges[ampId][ampSide] = np.nanmedian(strip, axis=ampSide % 2) # 1D medianified strip
558 return ampEdges
559
560 def getInterfaceOffset(self, ampNameA, ampNameB, edgeA, edgeB):
561 """Calculate the amp offset for a given interface between two
562 amplifiers.
563
564 Parameters
565 ----------
566 ampNameA : str
567 Name of the first amplifier.
568 ampNameB : str
569 Name of the second amplifier.
570 edgeA : numpy.ndarray
571 Amp edge for the first amplifier.
572 edgeB : numpy.ndarray
573 Amp edge for the second amplifier.
574
575 Returns
576 -------
577 interfaceOffset : float
578 The calculated amp offset value for the given interface between
579 amps A and B.
580 interfaceOffsetOriginal : float
581 The original calculated amp offset value for the given interface
582 between amps A and B.
583 ampEdgeGoodFrac : float
584 Fraction of viable pixel rows along the amp edge.
585 minFracFail : bool
586 True if the fraction of unmasked pixel rows is below the
587 ampEdgeMinFrac threshold.
588 maxOffsetFail : bool
589 True if the absolute offset value exceeds the ampEdgeMaxOffset
590 threshold.
591 """
592 interfaceId = f"{ampNameA}-{ampNameB}"
593 sctrl = StatisticsControl()
594 # NOTE: Taking the difference with the order below fixes the sign flip
595 # in the B matrix.
596 edgeDiff = edgeA - edgeB
597 if self.config.doWindowSmoothing:
598 # Compute rolling averages.
599 window = int(self.config.ampEdgeWindowFrac * len(edgeDiff))
600 edgeDiffSum = np.convolve(np.nan_to_num(edgeDiff), np.ones(window), "same")
601 edgeDiffNum = np.convolve(~np.isnan(edgeDiff), np.ones(window), "same")
602 edgeDiffAvg = edgeDiffSum / np.clip(edgeDiffNum, 1, None)
603 else:
604 # Directly use the difference.
605 edgeDiffAvg = edgeDiff.copy()
606 edgeDiffAvg[np.isnan(edgeDiff)] = np.nan
607 # Take clipped mean of rolling average data as amp offset value.
608 interfaceOffset = makeStatistics(edgeDiffAvg, MEANCLIP, sctrl).getValue()
609 interfaceOffsetOriginal = interfaceOffset
610 ampEdgeGoodFrac = 1 - (np.sum(np.isnan(edgeDiffAvg)) / len(edgeDiffAvg))
611
612 # Perform a couple of do-no-harm safety checks:
613 # a) The fraction of unmasked pixel rows is > ampEdgeMinFrac,
614 # b) The absolute offset ADU value is < ampEdgeMaxOffset.
615 minFracFail = ampEdgeGoodFrac < self.config.ampEdgeMinFrac
616 maxOffsetFail = np.abs(interfaceOffset) > self.config.ampEdgeMaxOffset
617 if minFracFail or maxOffsetFail:
618 interfaceOffset = 0
619 self.log.debug(
620 f"amp interface '{interfaceId}': "
621 f"viable edge difference frac = {ampEdgeGoodFrac:.2f}, "
622 f"amp offset = {interfaceOffsetOriginal:.3f}"
623 )
624 return (
625 interfaceId,
626 interfaceOffset,
627 interfaceOffsetOriginal,
628 ampEdgeGoodFrac,
629 minFracFail,
630 maxOffsetFail,
631 )
Pass parameters to a Statistics object.
Definition Statistics.h:83
getAmpOffsets(self, im, amps, associations, sides)
Definition ampOffset.py:398
getInterfaceOffset(self, ampNameA, ampNameB, edgeA, edgeB)
Definition ampOffset.py:560
__init__(self, *args, **kwargs)
Definition ampOffset.py:136
getAmpEdges(self, im, amps, ampSides)
Definition ampOffset.py:513
getNeighbors(self, ampIds, ampId)
Definition ampOffset.py:363