Coverage for python/lsst/cp/pipe/cpDefects.py: 71%
548 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-16 10:36 +0000
1# This file is part of cp_pipe.
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#
23__all__ = ['MeasureDefectsTaskConfig', 'MeasureDefectsTask',
24 'MergeDefectsTaskConfig', 'MergeDefectsTask',
25 'MeasureDefectsCombinedTaskConfig', 'MeasureDefectsCombinedTask',
26 'MeasureDefectsCombinedWithFilterTaskConfig', 'MeasureDefectsCombinedWithFilterTask',
27 'MergeDefectsCombinedTaskConfig', 'MergeDefectsCombinedTask', ]
29import numpy as np
31import lsst.pipe.base as pipeBase
32import lsst.pipe.base.connectionTypes as cT
34from lsstDebug import getDebugFrame
35import lsst.pex.config as pexConfig
37import lsst.afw.image as afwImage
38import lsst.afw.math as afwMath
39import lsst.afw.detection as afwDetection
40import lsst.afw.display as afwDisplay
41from lsst.afw import cameraGeom
42from lsst.geom import Box2I, Point2I, Extent2I
43from lsst.meas.algorithms import SourceDetectionTask
44from lsst.ip.isr import Defects, countMaskedPixels, PhotonTransferCurveDataset
45from lsst.pex.exceptions import InvalidParameterError
47from .utils import bin_flat, FlatGradientFitter
48from lsst.ip.isr import FlatGradient
51class MeasureDefectsConnections(pipeBase.PipelineTaskConnections,
52 dimensions=("instrument", "exposure", "detector")):
53 inputExp = cT.Input(
54 name="defectExps",
55 doc="Input ISR-processed exposures to measure.",
56 storageClass="Exposure",
57 dimensions=("instrument", "detector", "exposure"),
58 multiple=False
59 )
60 camera = cT.PrerequisiteInput(
61 name='camera',
62 doc="Camera associated with this exposure.",
63 storageClass="Camera",
64 dimensions=("instrument", ),
65 isCalibration=True,
66 )
68 outputDefects = cT.Output(
69 name="singleExpDefects",
70 doc="Output measured defects.",
71 storageClass="Defects",
72 dimensions=("instrument", "detector", "exposure"),
73 )
76class MeasureDefectsTaskConfig(pipeBase.PipelineTaskConfig,
77 pipelineConnections=MeasureDefectsConnections):
78 """Configuration for measuring defects from a list of exposures
79 """
81 thresholdType = pexConfig.ChoiceField(
82 dtype=str,
83 doc=("Defects threshold type: ``STDEV`` or ``VALUE``. If ``VALUE``, cold pixels will be found "
84 "in flats, and hot pixels in darks. If ``STDEV``, cold and hot pixels will be found "
85 "in flats, and hot pixels in darks."),
86 default='STDEV',
87 allowed={'STDEV': "Use a multiple of the image standard deviation to determine detection threshold.",
88 'VALUE': "Use pixel value to determine detection threshold."},
89 )
90 fitAmpGradient = pexConfig.Field(
91 dtype=bool,
92 doc="Fit out the focal plane radial gradient per amplifier.",
93 default=False,
94 )
95 doVampirePixels = pexConfig.Field(
96 dtype=bool,
97 doc=("Search for vampire pixels (bright pixels surrounded by ring of low flux) in ComCam "
98 "flatBootstrap and mask the area arount them."),
99 default=False,
100 )
101 thresholdVampirePixels = pexConfig.Field(
102 dtype=float,
103 doc=("Pixel value threshold to find bright pixels in ComCam flatBootstrap."),
104 default=1.9,
105 )
106 radiusVampirePixels = pexConfig.Field(
107 dtype=int,
108 doc=("Radius (in pixels) of the area to mask around ComCam flatBootstrap bright pixels."),
109 default=8,
110 )
111 darkCurrentThreshold = pexConfig.Field(
112 dtype=float,
113 doc=("If thresholdType=``VALUE``, dark current threshold (in e-/sec) to define "
114 "hot/bright pixels in dark images. Unused if thresholdType==``STDEV``."),
115 default=5,
116 )
117 biasThreshold = pexConfig.Field(
118 dtype=float,
119 doc=("If thresholdType==``VALUE``, bias threshold (in ADU) to define "
120 "hot/bright pixels in bias frame. Unused if thresholdType==``STDEV``."),
121 default=1000.0,
122 )
123 fracThresholdFlat = pexConfig.Field(
124 dtype=float,
125 doc=("If thresholdType=``VALUE``, fractional threshold to define cold/dark "
126 "pixels in flat images (fraction of the mean value per amplifier)."
127 "Unused if thresholdType==``STDEV``."),
128 default=0.8,
129 )
130 nSigmaBright = pexConfig.Field(
131 dtype=float,
132 doc=("If thresholdType=``STDEV``, number of sigma above mean for bright/hot "
133 "pixel detection. The default value was found to be "
134 "appropriate for some LSST sensors in DM-17490. "
135 "Unused if thresholdType==``VALUE``"),
136 default=4.8,
137 )
138 nSigmaDark = pexConfig.Field(
139 dtype=float,
140 doc=("If thresholdType=``STDEV``, number of sigma below mean for dark/cold pixel "
141 "detection. The default value was found to be "
142 "appropriate for some LSST sensors in DM-17490. "
143 "Unused if thresholdType==``VALUE``"),
144 default=-5.0,
145 )
146 nPixBorderUpDown = pexConfig.Field(
147 dtype=int,
148 doc="Number of pixels to exclude from top & bottom of image when looking for defects.",
149 default=0,
150 )
151 nPixBorderLeftRight = pexConfig.Field(
152 dtype=int,
153 doc="Number of pixels to exclude from left & right of image when looking for defects.",
154 default=0,
155 )
156 nPixBorderUpDownITL = pexConfig.Field(
157 dtype=int,
158 doc="Number of pixels to exclude from up & down of image when looking for defects in ITL.",
159 default=0,
160 )
161 nPixBorderLeftRightITL = pexConfig.Field(
162 dtype=int,
163 doc="Number of pixels to exclude from left & right of image when looking for defects in ITL.",
164 default=0,
165 )
166 nPixBorderUpDownE2V = pexConfig.Field(
167 dtype=int,
168 doc="Number of pixels to exclude from up & down of image when looking for defects in E2V.",
169 default=0,
170 )
171 nPixBorderLeftRightE2V = pexConfig.Field(
172 dtype=int,
173 doc="Number of pixels to exclude from left & right of image when looking for defects in E2V.",
174 default=0,
175 )
176 badOnAndOffPixelColumnThreshold = pexConfig.Field(
177 dtype=int,
178 doc=("If BPC is the set of all the bad pixels in a given column (not necessarily consecutive) "
179 "and the size of BPC is at least 'badOnAndOffPixelColumnThreshold', all the pixels between the "
180 "pixels that satisfy minY (BPC) and maxY (BPC) will be marked as bad, with 'Y' being the long "
181 "axis of the amplifier (and 'X' the other axis, which for a column is a constant for all "
182 "pixels in the set BPC). If there are more than 'goodPixelColumnGapThreshold' consecutive "
183 "non-bad pixels in BPC, an exception to the above is made and those consecutive "
184 "'goodPixelColumnGapThreshold' are not marked as bad."),
185 default=50,
186 )
187 goodPixelColumnGapThreshold = pexConfig.Field(
188 dtype=int,
189 doc=("Size, in pixels, of usable consecutive pixels in a column with on and off bad pixels (see "
190 "'badOnAndOffPixelColumnThreshold')."),
191 default=30,
192 )
193 badPixelsToFillColumnThreshold = pexConfig.Field(
194 dtype=float,
195 doc=("If the number of bad pixels in an amplifier column is above this threshold "
196 "then the full amplifier column will be marked bad. This operation is performed after "
197 "any merging of blinking columns performed with badOnAndOffPixelColumnThreshold. If this"
198 "value is less than 0 then no bad column filling will be performed."),
199 default=-1,
200 )
201 saturatedColumnMask = pexConfig.Field(
202 dtype=str,
203 default="SAT",
204 doc="Saturated mask plane for dilation.",
205 )
206 saturatedColumnDilationRadius = pexConfig.Field(
207 dtype=int,
208 doc=("Dilation radius (along rows) to use to expand saturated columns "
209 "to mitigate glow."),
210 default=0,
211 )
212 saturatedPixelsToFillColumnThreshold = pexConfig.Field(
213 dtype=int,
214 doc=("If the number of saturated pixels in an amplifier column is above this threshold "
215 "then the full amplifier column will be marked bad. If this value is less than 0"
216 "then no saturated column filling will be performed."),
217 default=-1,
218 )
219 e2vMidlineBreakNRow = pexConfig.Field(
220 dtype=int,
221 doc="E2V midline break number of midline break rows at the bottom of the top amps "
222 "or the top of the bottom amps to always mask or ignore. Number of rows will "
223 "be twice this config value. Only used if detector is E2V type. Set to <=0 "
224 "to treat E2V midline break rows like any other flat-field pixel.",
225 default=1,
226 )
227 e2vMidlineBreakOption = pexConfig.ChoiceField(
228 dtype=str,
229 doc="How should the E2V midline break be treated? Only used if e2vMidlineBreakNRow > 0.",
230 default="NEVERMASK",
231 allowed={
232 "NEVERMASK": "Never mask the E2V midline break, no matter the flat values.",
233 "MASK": "Always mask the E2V midline break.",
234 },
235 )
236 ampGradientBinFactor = pexConfig.Field(
237 dtype=int,
238 doc="Binning factor used for fitting per-amp focal plane gradient.",
239 default=8,
240 )
241 ampGradientBoundarySize = pexConfig.Field(
242 dtype=int,
243 doc="Amp boundary exclusion size for fitting per-amp focal plane gradient.",
244 default=20,
245 )
246 ampGradientNodes = pexConfig.Field(
247 dtype=int,
248 doc="Number of spline nodes for per-amp focal plane gradient.",
249 default=4,
250 )
252 def validate(self):
253 super().validate()
254 if self.nSigmaBright < 0.0: 254 ↛ 255line 254 didn't jump to line 255 because the condition on line 254 was never true
255 raise ValueError("nSigmaBright must be above 0.0.")
256 if self.nSigmaDark > 0.0: 256 ↛ 257line 256 didn't jump to line 257 because the condition on line 256 was never true
257 raise ValueError("nSigmaDark must be below 0.0.")
260class MeasureDefectsTask(pipeBase.PipelineTask):
261 """Measure the defects from one exposure.
262 """
264 ConfigClass = MeasureDefectsTaskConfig
265 _DefaultName = 'cpDefectMeasure'
267 def run(self, inputExp, camera):
268 """Measure one exposure for defects.
270 Parameters
271 ----------
272 inputExp : `lsst.afw.image.Exposure`
273 Exposure to examine.
274 camera : `lsst.afw.cameraGeom.Camera`
275 Camera to use for metadata.
277 Returns
278 -------
279 results : `lsst.pipe.base.Struct`
280 Results struct containing:
282 ``outputDefects``
283 The defects measured from this exposure
284 (`lsst.ip.isr.Defects`).
285 """
286 detector = inputExp.getDetector()
287 try:
288 filterName = inputExp.getFilter().physicalLabel
289 except AttributeError:
290 filterName = None
292 defects = self._findHotAndColdPixels(inputExp)
294 datasetType = inputExp.getMetadata().get('IMGTYPE', 'UNKNOWN')
295 msg = "Found %s defects containing %s pixels in %s"
296 self.log.info(msg, len(defects), self._nPixFromDefects(defects), datasetType)
298 defects.updateMetadataFromExposures([inputExp])
299 defects.updateMetadata(camera=camera, detector=detector, filterName=filterName,
300 setCalibId=True, setDate=True,
301 cpDefectGenImageType=datasetType)
303 return pipeBase.Struct(
304 outputDefects=defects,
305 )
307 @staticmethod
308 def _nPixFromDefects(defects):
309 """Count pixels in a defect.
311 Parameters
312 ----------
313 defects : `lsst.ip.isr.Defects`
314 Defects to measure.
316 Returns
317 -------
318 nPix : `int`
319 Number of defect pixels.
320 """
321 nPix = 0
322 for defect in defects:
323 nPix += defect.getBBox().getArea()
324 return nPix
326 def getVampirePixels(self, ampImg):
327 """Find vampire pixels (bright pixels in flats) and get footprint of
328 extended area around them,
330 Parameters
331 ----------
332 ampImg : `lsst.afw.image._maskedImage.MaskedImageF`
333 The amplifier masked image to do the vampire pixels search on.
335 Returns
336 -------
337 fs_grow : `lsst.afw.detection._detection.FootprintSet`
338 The footprint set of areas around vampire pixels in the amplifier.
339 """
341 # Find bright pixels
342 thresh = afwDetection.Threshold(self.config.thresholdVampirePixels)
343 # Bright pixels footprint grown by a radius of radiusVampire pixels
344 fs = afwDetection.FootprintSet(ampImg, thresh)
345 fs_grow = afwDetection.FootprintSet(fs, rGrow=self.config.radiusVampirePixels, isotropic=True)
347 return fs_grow
349 def _setE2VMidline(self, amp, ampImg):
350 """Set the E2V midline pixels on this amplifier.
352 If the task is configure to always mask it will set the midline rows to
353 an illegal value; if it is to never mask it will set the midline rows
354 to the median value.
356 Parameters
357 ----------
358 amp : `lsst.afw.cameraGeom.Amplifier`
359 Amplifier object.
360 ampImg : `lsst.afw.image.ImageF`
361 Image data for the amplifier.
362 """
363 if self.config.e2vMidlineBreakNRow <= 0:
364 return
366 if self.config.e2vMidlineBreakOption == "NEVERMASK":
367 # Setting the midline values to the median value ensures
368 # they will never trigger the defect finding threshold.
369 value = np.nanmedian(ampImg.image.array)
370 else:
371 # Setting the midline values to a large negative value
372 # ensures they will always trigger the defect finding
373 # threshold.
374 value = -1e30
376 if amp.getName().startswith("C0"):
377 # This is a bottom amplifier, so change the top row(s).
378 ampImg.image.array[-self.config.e2vMidlineBreakNRow:, :] = value
379 else:
380 # This is a top amplifier, so change the bottom row(s).
381 ampImg.image.array[:self.config.e2vMidlineBreakNRow, :] = value
383 def _findHotAndColdPixels(self, exp):
384 """Find hot and cold pixels in an image.
386 Using config-defined thresholds on a per-amp basis, mask
387 pixels that are nSigma above threshold in dark frames (hot
388 pixels), or nSigma away from the clipped mean in flats (hot &
389 cold pixels).
391 Parameters
392 ----------
393 exp : `lsst.afw.image.exposure.Exposure`
394 The exposure in which to find defects.
396 Returns
397 -------
398 defects : `lsst.ip.isr.Defects`
399 The defects found in the image.
400 """
402 # the detection polarity for afwDetection, True for positive,
403 # False for negative, and therefore True for darks as they only have
404 # bright pixels, and both for flats, as they have bright and dark pix
405 footprintList = []
407 hotPixelCount = {}
408 coldPixelCount = {}
410 detector = exp.getDetector()
411 detectorType = detector.getPhysicalType()
413 self._setEdgeBits(exp, detectorType=detectorType)
415 if self.config.fitAmpGradient:
416 # 1. Bin flat.
417 binnedExp = bin_flat(
418 PhotonTransferCurveDataset(),
419 exp,
420 bin_factor=self.config.ampGradientBinFactor,
421 amp_boundary=self.config.ampGradientBoundarySize,
422 apply_gains=False,
423 )
425 # 2. Prepare transformations.
426 transform = detector.getTransform(
427 cameraGeom.PIXELS,
428 cameraGeom.FOCAL_PLANE,
429 )
431 # 2a. Do transformation for binned coordinates.
432 binnedXy = np.vstack((binnedExp["xd"], binnedExp["yd"]))
433 binnedXf, binnedYf = np.vsplit(
434 transform.getMapping().applyForward(binnedXy.astype(np.float64)),
435 2,
436 )
437 binnedXf = binnedXf.ravel()
438 binnedYf = binnedYf.ravel()
439 binnedRadius = np.sqrt(binnedXf**2. + binnedYf**2.)
441 # 2b. Do transformation for full detector coordinates.
442 xx = np.arange(exp.image.array.shape[1], dtype=np.int64)
443 yy = np.arange(exp.image.array.shape[0], dtype=np.int64)
444 x, y = np.meshgrid(xx, yy)
445 x = x.ravel()
446 y = y.ravel()
448 xy = np.vstack((x, y))
449 xf, yf = np.vsplit(transform.getMapping().applyForward(xy.astype(np.float64)), 2)
450 xf = xf.ravel()
451 yf = yf.ravel()
453 # Measure the focal plane gradient on the binned image and remove
454 # on the full image, amp by amp.
455 for i, amp in enumerate(detector):
456 # Fit the focal plane radial gradient in the binned image.
457 binnedInAmp = (binnedExp["amp_index"] == i)
459 nodes = np.linspace(
460 np.min(binnedRadius[binnedInAmp]),
461 np.max(binnedRadius[binnedInAmp]),
462 self.config.ampGradientNodes,
463 )
465 norm = np.nanpercentile(binnedExp["value"][binnedInAmp], 95.0)
467 fitter = FlatGradientFitter(
468 nodes,
469 binnedXf[binnedInAmp],
470 binnedYf[binnedInAmp],
471 binnedExp["value"][binnedInAmp]/norm,
472 np.array([]),
473 constrain_zero=False,
474 )
475 p0 = fitter.compute_p0()
476 pars = fitter.fit(p0)
478 gradient = FlatGradient()
479 gradient.setParameters(
480 radialSplineNodes=nodes,
481 radialSplineValues=pars[fitter.indices["spline"]],
482 )
484 # Apply the focal plane gradient to the full image.
485 pixelsInAmp = amp.getBBox().contains(x, y)
487 model = gradient.computeRadialSplineModelXY(xf[pixelsInAmp], yf[pixelsInAmp])
489 exp.image.array[y[pixelsInAmp], x[pixelsInAmp]] /= model
491 maskedIm = exp.maskedImage
493 for amp in exp.getDetector():
494 ampName = amp.getName()
496 hotPixelCount[ampName] = 0
497 coldPixelCount[ampName] = 0
499 ampImg = maskedIm[amp.getBBox()].clone()
501 # crop ampImage depending on where the amp lies in the image
502 # and depending on detector type to mask picture frame effect
503 if detectorType == 'E2V':
504 nPixBorderLeftRight = self.config.nPixBorderLeftRightE2V
505 elif 'ITL' in detectorType:
506 nPixBorderLeftRight = self.config.nPixBorderLeftRightITL
507 else:
508 nPixBorderLeftRight = self.config.nPixBorderLeftRight
510 if detectorType == 'E2V':
511 nPixBorderUpDown = self.config.nPixBorderUpDownE2V
512 elif 'ITL' in detectorType:
513 nPixBorderUpDown = self.config.nPixBorderUpDownITL
514 else:
515 nPixBorderUpDown = self.config.nPixBorderUpDown
517 if nPixBorderLeftRight:
518 if ampImg.getBBox().getMinX() == 0:
519 ampImg = ampImg[nPixBorderLeftRight:, :, afwImage.LOCAL]
520 elif ampImg.getBBox().getMaxX() == exp.getBBox().getMaxX(): 520 ↛ 523line 520 didn't jump to line 523 because the condition on line 520 was always true
521 ampImg = ampImg[:-nPixBorderLeftRight, :, afwImage.LOCAL]
523 if nPixBorderUpDown:
524 if ampImg.getBBox().getMinY() == 0:
525 ampImg = ampImg[:, nPixBorderUpDown:, afwImage.LOCAL]
526 elif ampImg.getBBox().getMaxY() == exp.getBBox().getMaxY():
527 ampImg = ampImg[:, :-nPixBorderUpDown, afwImage.LOCAL]
529 if self._getNumGoodPixels(ampImg) == 0: # amp contains no usable pixels 529 ↛ 530line 529 didn't jump to line 530 because the condition on line 529 was never true
530 continue
532 if self.config.doVampirePixels:
533 # This is only applied in LSSTComCam flatBootstrap pipeline
534 footprintSet_VampirePixel = self.getVampirePixels(ampImg)
535 footprintSet_VampirePixel.setMask(maskedIm.mask, ("BAD"))
537 # Remove a background estimate
538 meanClip = afwMath.makeStatistics(ampImg, afwMath.MEANCLIP, ).getValue()
539 ampImg -= meanClip
541 # Determine thresholds
542 stDev = afwMath.makeStatistics(ampImg, afwMath.STDEVCLIP, ).getValue()
543 expTime = exp.getInfo().getVisitInfo().getExposureTime()
544 datasetType = exp.getMetadata().get('IMGTYPE', 'UNKNOWN')
545 if np.isnan(expTime): 545 ↛ 546line 545 didn't jump to line 546 because the condition on line 545 was never true
546 self.log.warning("expTime=%s for AMP %s in %s. Setting expTime to 1 second",
547 expTime, ampName, datasetType)
548 expTime = 1.
549 thresholdType = self.config.thresholdType
550 if thresholdType == 'VALUE':
551 # LCA-128 and eoTest: bright/hot pixels in dark images are
552 # defined as any pixel with more than 5 e-/s of dark current.
553 # We scale by the exposure time.
554 if datasetType.lower() == 'dark':
555 # hot pixel threshold
556 valueThreshold = self.config.darkCurrentThreshold*expTime/amp.getGain()
557 elif datasetType.lower() == 'bias': 557 ↛ 559line 557 didn't jump to line 559 because the condition on line 557 was never true
558 # hot pixel threshold, no exposure time.
559 valueThreshold = self.config.biasThreshold
560 else:
561 # LCA-128 and eoTest: dark/cold pixels in flat images as
562 # defined as any pixel with photoresponse <80% of
563 # the mean (at 500nm).
565 # We subtracted the mean above, so the threshold will be
566 # negative cold pixel threshold.
567 valueThreshold = (self.config.fracThresholdFlat-1)*meanClip
568 # Find equivalent sigma values.
569 if stDev == 0.0:
570 self.log.warning("stDev=%s for AMP %s in %s. Setting nSigma to inf.",
571 stDev, ampName, datasetType)
572 nSigmaList = [np.inf]
573 else:
574 nSigmaList = [valueThreshold/stDev]
575 else:
576 hotPixelThreshold = self.config.nSigmaBright
577 coldPixelThreshold = self.config.nSigmaDark
578 if datasetType.lower() == 'dark': 578 ↛ 579line 578 didn't jump to line 579 because the condition on line 578 was never true
579 nSigmaList = [hotPixelThreshold]
580 valueThreshold = stDev*hotPixelThreshold
581 elif datasetType.lower() == 'bias': 581 ↛ 582line 581 didn't jump to line 582 because the condition on line 581 was never true
582 self.log.warning(
583 "Bias frame detected, but thresholdType == STDEV; not looking for defects.",
584 )
585 return Defects.fromFootprintList([])
586 else:
587 nSigmaList = [hotPixelThreshold, coldPixelThreshold]
588 valueThreshold = [x*stDev for x in nSigmaList]
590 self.log.info("Image type: %s. Amp: %s. Threshold Type: %s. Sigma values and Pixel"
591 "Values (hot and cold pixels thresholds): %s, %s",
592 datasetType, ampName, thresholdType, nSigmaList, valueThreshold)
594 if datasetType.lower() == "flat" and exp.getDetector().getPhysicalType() == "E2V":
595 self._setE2VMidline(amp, ampImg)
597 mergedSet = None
598 for sigma in nSigmaList:
599 nSig = np.abs(sigma)
600 self.debugHistogram('ampFlux', ampImg, nSig, exp)
601 polarity = {-1: False, 1: True}[np.sign(sigma)]
603 threshold = afwDetection.createThreshold(nSig, 'stdev', polarity=polarity)
605 try:
606 footprintSet = afwDetection.FootprintSet(ampImg, threshold)
607 except InvalidParameterError:
608 # This occurs if the image sigma value is 0.0.
609 # Let's mask the whole area.
610 minValue = np.nanmin(ampImg.image.array) - 1.0
611 threshold = afwDetection.createThreshold(minValue, 'value', polarity=True)
612 footprintSet = afwDetection.FootprintSet(ampImg, threshold)
614 footprintSet.setMask(maskedIm.mask, ("DETECTED" if polarity else "DETECTED_NEGATIVE"))
616 if mergedSet is None:
617 mergedSet = footprintSet
618 else:
619 mergedSet.merge(footprintSet)
621 if polarity:
622 # hot pixels
623 for fp in footprintSet.getFootprints():
624 hotPixelCount[ampName] += fp.getArea()
625 else:
626 # cold pixels
627 for fp in footprintSet.getFootprints():
628 coldPixelCount[ampName] += fp.getArea()
630 if self.config.doVampirePixels:
631 # Count the number of pixels masked
632 vampirePixelCount = 0
633 for fp in footprintSet_VampirePixel.getFootprints():
634 vampirePixelCount += fp.getArea()
635 self.log.info("%s Vampire pixels are masked", vampirePixelCount)
636 # Add vampire pixels to footprint set
637 mergedSet.merge(footprintSet_VampirePixel)
639 footprintList += mergedSet.getFootprints()
641 self.debugView('defectMap', ampImg,
642 Defects.fromFootprintList(mergedSet.getFootprints()), exp.getDetector())
644 defects = Defects.fromFootprintList(footprintList)
645 defects = self.dilateSaturatedColumns(exp, defects)
646 defects, _ = self.maskBlocksIfIntermitentBadPixelsInColumn(defects)
647 defects, count = self.maskBadColumns(exp, defects)
648 # We want this to reflect the number of completely bad columns.
649 defects.updateCounters(columns=count, hot=hotPixelCount, cold=coldPixelCount)
651 return defects
653 @staticmethod
654 def _getNumGoodPixels(maskedIm, badMaskString="NO_DATA"):
655 """Return the number of non-bad pixels in the image."""
656 nPixels = maskedIm.mask.array.size
657 nBad = countMaskedPixels(maskedIm, badMaskString)
658 return nPixels - nBad
660 def _setEdgeBits(self, exposureOrMaskedImage, maskplaneToSet='EDGE', detectorType='CCD'):
661 """Set edge bits on an exposure or maskedImage.
663 Parameters
664 ----------
665 exposureOrMaskedImage : `lsst.afw.image.exposure.Exposure`
666 or `lsst.afw.image.MaskedImage`
667 The exposure or masked image in which to find defects.
668 maskplaneToSet : `str`
669 Name of mask plane the edges are set to.
670 detectorType : `str`
671 Name of the detector type.
673 Raises
674 ------
675 TypeError
676 Raised if parameter ``exposureOrMaskedImage`` is an invalid type.
677 """
678 if isinstance(exposureOrMaskedImage, afwImage.Exposure):
679 mi = exposureOrMaskedImage.maskedImage
680 elif isinstance(exposureOrMaskedImage, afwImage.MaskedImage): 680 ↛ 683line 680 didn't jump to line 683 because the condition on line 680 was always true
681 mi = exposureOrMaskedImage
682 else:
683 t = type(exposureOrMaskedImage)
684 raise TypeError(f"Function supports exposure or maskedImage but not {t}")
686 MASKBIT = mi.mask.getPlaneBitMask(maskplaneToSet)
688 if detectorType == 'E2V':
689 nPixBorderLeftRight = self.config.nPixBorderLeftRightE2V
690 elif 'ITL' in detectorType:
691 nPixBorderLeftRight = self.config.nPixBorderLeftRightITL
692 else:
693 nPixBorderLeftRight = self.config.nPixBorderLeftRight
695 if detectorType == 'E2V':
696 nPixBorderUpDown = self.config.nPixBorderUpDownE2V
697 elif 'ITL' in detectorType:
698 nPixBorderUpDown = self.config.nPixBorderUpDownITL
699 else:
700 nPixBorderUpDown = self.config.nPixBorderUpDown
702 if nPixBorderLeftRight:
703 mi.mask[: nPixBorderLeftRight, :, afwImage.LOCAL] |= MASKBIT
704 mi.mask[-nPixBorderLeftRight:, :, afwImage.LOCAL] |= MASKBIT
705 if nPixBorderUpDown:
706 mi.mask[:, : nPixBorderUpDown, afwImage.LOCAL] |= MASKBIT
707 mi.mask[:, -nPixBorderUpDown:, afwImage.LOCAL] |= MASKBIT
709 def maskBlocksIfIntermitentBadPixelsInColumn(self, defects):
710 """Mask blocks in a column if there are on-and-off bad pixels
712 If there's a column with on and off bad pixels, mask all the
713 pixels in between, except if there is a large enough gap of
714 consecutive good pixels between two bad pixels in the column.
716 Parameters
717 ----------
718 defects : `lsst.ip.isr.Defects`
719 The defects found in the image so far
721 Returns
722 -------
723 defects : `lsst.ip.isr.Defects`
724 If the number of bad pixels in a column is not larger or
725 equal than self.config.badPixelColumnThreshold, the input
726 list is returned. Otherwise, the defects list returned
727 will include boxes that mask blocks of on-and-of pixels.
728 badColumnCount : `int`
729 Number of bad columns partially masked.
730 """
731 badColumnCount = 0
732 # Get the (x, y) values of each bad pixel in amp.
733 coordinates = []
734 for defect in defects:
735 bbox = defect.getBBox()
736 x0, y0 = bbox.getMinX(), bbox.getMinY()
737 deltaX0, deltaY0 = bbox.getDimensions()
738 for j in np.arange(y0, y0+deltaY0):
739 for i in np.arange(x0, x0 + deltaX0):
740 coordinates.append((i, j))
742 x, y = [], []
743 for coordinatePair in coordinates:
744 x.append(coordinatePair[0])
745 y.append(coordinatePair[1])
747 x = np.array(x)
748 y = np.array(y)
749 # Find the defects with same "x" (vertical) coordinate (column).
750 unique, counts = np.unique(x, return_counts=True)
751 multipleX = []
752 for (a, b) in zip(unique, counts):
753 if b >= self.config.badOnAndOffPixelColumnThreshold:
754 multipleX.append(a)
755 if len(multipleX) != 0:
756 defects = self._markBlocksInBadColumn(x, y, multipleX, defects)
757 badColumnCount += 1
759 return defects, badColumnCount
761 def dilateSaturatedColumns(self, exp, defects):
762 """Dilate saturated columns by a configurable amount.
764 Parameters
765 ----------
766 exp : `lsst.afw.image.exposure.Exposure`
767 The exposure in which to find defects.
768 defects : `lsst.ip.isr.Defects`
769 The defects found in the image so far
771 Returns
772 -------
773 defects : `lsst.ip.isr.Defects`
774 The expanded defects.
775 """
776 if self.config.saturatedColumnDilationRadius <= 0:
777 # This is a no-op.
778 return defects
780 mask = afwImage.Mask.getPlaneBitMask(self.config.saturatedColumnMask)
782 satY, satX = np.where((exp.mask.array & mask) > 0)
784 if len(satX) == 0:
785 # No saturated pixels, nothing to do.
786 return defects
788 radius = self.config.saturatedColumnDilationRadius
790 with defects.bulk_update():
791 for index in range(len(satX)):
792 minX = np.clip(satX[index] - radius, 0, None)
793 maxX = np.clip(satX[index] + radius, None, exp.image.array.shape[1] - 1)
794 s = Box2I(minimum=Point2I(minX, satY[index]),
795 maximum=Point2I(maxX, satY[index]))
796 defects.append(s)
798 return defects
800 def maskBadColumns(self, exp, defects):
801 """Mask full amplifier columns if they are sufficiently bad.
803 Parameters
804 ----------
805 defects : `lsst.ip.isr.Defects`
806 The defects found in the image so far
808 Returns
809 -------
810 exp : `lsst.afw.image.exposure.Exposure`
811 The exposure in which to find defects.
812 defects : `lsst.ip.isr.Defects`
813 If the number of bad pixels in a column is not larger or
814 equal than self.config.badPixelColumnThreshold, the input
815 list is returned. Otherwise, the defects list returned
816 will include boxes that mask blocks of on-and-of pixels.
817 badColumnCount : `int`
818 Number of bad columns masked.
819 """
820 # Render the defects into an image.
821 defectImage = afwImage.ImageI(exp.getBBox())
823 for defect in defects:
824 defectImage[defect.getBBox()] = 1
826 badColumnCount = 0
828 if self.config.badPixelsToFillColumnThreshold > 0:
829 with defects.bulk_update():
830 for amp in exp.getDetector():
831 subImage = defectImage[amp.getBBox()].array
832 nInCol = np.sum(subImage, axis=0)
834 badColIndices, = (nInCol >= self.config.badPixelsToFillColumnThreshold).nonzero()
835 badColumns = badColIndices + amp.getBBox().getMinX()
837 for badColumn in badColumns:
838 s = Box2I(minimum=Point2I(badColumn, amp.getBBox().getMinY()),
839 maximum=Point2I(badColumn, amp.getBBox().getMaxY()))
840 defects.append(s)
842 badColumnCount += len(badColIndices)
844 if self.config.saturatedPixelsToFillColumnThreshold > 0:
845 mask = afwImage.Mask.getPlaneBitMask(self.config.saturatedColumnMask)
847 with defects.bulk_update():
848 for amp in exp.getDetector():
849 subMask = exp.mask[amp.getBBox()].array
850 # Turn all the SAT bits into 1s
851 subMask &= mask
852 subMask[subMask > 0] = 1
854 nInCol = np.sum(subMask, axis=0)
856 badColIndices, = (nInCol >= self.config.saturatedPixelsToFillColumnThreshold).nonzero()
857 badColumns = badColIndices + amp.getBBox().getMinX()
859 for badColumn in badColumns:
860 s = Box2I(minimum=Point2I(badColumn, amp.getBBox().getMinY()),
861 maximum=Point2I(badColumn, amp.getBBox().getMaxY()))
862 defects.append(s)
864 badColumnCount += len(badColIndices)
866 return defects, badColumnCount
868 def _markBlocksInBadColumn(self, x, y, multipleX, defects):
869 """Mask blocks in a column if number of on-and-off bad pixels is above
870 threshold.
872 This function is called if the number of on-and-off bad pixels
873 in a column is larger or equal than
874 self.config.badOnAndOffPixelColumnThreshold.
876 Parameters
877 ---------
878 x : `list`
879 Lower left x coordinate of defect box. x coordinate is
880 along the short axis if amp.
881 y : `list`
882 Lower left y coordinate of defect box. x coordinate is
883 along the long axis if amp.
884 multipleX : list
885 List of x coordinates in amp. with multiple bad pixels
886 (i.e., columns with defects).
887 defects : `lsst.ip.isr.Defects`
888 The defcts found in the image so far
890 Returns
891 -------
892 defects : `lsst.ip.isr.Defects`
893 The defects list returned that will include boxes that
894 mask blocks of on-and-of pixels.
895 """
896 with defects.bulk_update():
897 goodPixelColumnGapThreshold = self.config.goodPixelColumnGapThreshold
898 for x0 in multipleX:
899 index = np.where(x == x0)
900 multipleY = y[index] # multipleY and multipleX are in 1-1 correspondence.
901 multipleY.sort() # Ensure that the y values are sorted to look for gaps.
902 minY, maxY = np.min(multipleY), np.max(multipleY)
903 # Next few lines: don't mask pixels in column if gap
904 # of good pixels between two consecutive bad pixels is
905 # larger or equal than 'goodPixelColumnGapThreshold'.
906 diffIndex = np.where(np.diff(multipleY) >= goodPixelColumnGapThreshold)[0]
907 if len(diffIndex) != 0:
908 limits = [minY] # put the minimum first
909 for gapIndex in diffIndex:
910 limits.append(multipleY[gapIndex])
911 limits.append(multipleY[gapIndex+1])
912 limits.append(maxY) # maximum last
913 for i in np.arange(0, len(limits)-1, 2):
914 s = Box2I(minimum=Point2I(x0, limits[i]), maximum=Point2I(x0, limits[i+1]))
915 defects.append(s)
916 else: # No gap is large enough
917 s = Box2I(minimum=Point2I(x0, minY), maximum=Point2I(x0, maxY))
918 defects.append(s)
919 return defects
921 def debugView(self, stepname, ampImage, defects, detector): # pragma: no cover
922 """Plot the defects found by the task.
924 Parameters
925 ----------
926 stepname : `str`
927 Debug frame to request.
928 ampImage : `lsst.afw.image.MaskedImage`
929 Amplifier image to display.
930 defects : `lsst.ip.isr.Defects`
931 The defects to plot.
932 detector : `lsst.afw.cameraGeom.Detector`
933 Detector holding camera geometry.
934 """
935 frame = getDebugFrame(self._display, stepname)
936 if frame:
937 disp = afwDisplay.Display(frame=frame)
938 disp.scale('asinh', 'zscale')
939 disp.setMaskTransparency(80)
940 disp.setMaskPlaneColor("BAD", afwDisplay.RED)
942 maskedIm = ampImage.clone()
943 defects.maskPixels(maskedIm, "BAD")
945 mpDict = maskedIm.mask.getMaskPlaneDict()
946 for plane in mpDict.keys():
947 if plane in ['BAD']:
948 continue
949 disp.setMaskPlaneColor(plane, afwDisplay.IGNORE)
951 disp.setImageColormap('gray')
952 disp.mtv(maskedIm)
953 cameraGeom.utils.overlayCcdBoxes(detector, isTrimmed=True, display=disp)
954 prompt = "Press Enter to continue [c]... "
955 while True:
956 ans = input(prompt).lower()
957 if ans in ('', 'c', ):
958 break
960 def debugHistogram(self, stepname, ampImage, nSigmaUsed, exp):
961 """Make a histogram of the distribution of pixel values for
962 each amp.
964 The main image data histogram is plotted in blue. Edge
965 pixels, if masked, are in red. Note that masked edge pixels
966 do not contribute to the underflow and overflow numbers.
968 Note that this currently only supports the 16-amp LSST
969 detectors.
971 Parameters
972 ----------
973 stepname : `str`
974 Debug frame to request.
975 ampImage : `lsst.afw.image.MaskedImage`
976 Amplifier image to display.
977 nSigmaUsed : `float`
978 The number of sigma used for detection
979 exp : `lsst.afw.image.exposure.Exposure`
980 The exposure in which the defects were found.
981 """
982 frame = getDebugFrame(self._display, stepname)
983 if frame: 983 ↛ 984line 983 didn't jump to line 984 because the condition on line 983 was never true
984 import matplotlib.pyplot as plt
986 detector = exp.getDetector()
987 nX = np.floor(np.sqrt(len(detector)))
988 nY = len(detector) // nX
989 fig, ax = plt.subplots(nrows=int(nY), ncols=int(nX), sharex='col', sharey='row', figsize=(13, 10))
991 expTime = exp.getInfo().getVisitInfo().getExposureTime()
993 for (amp, a) in zip(reversed(detector), ax.flatten()):
994 mi = exp.maskedImage[amp.getBBox()]
996 # normalize by expTime as we plot in ADU/s and don't
997 # always work with master calibs
998 mi.image.array /= expTime
999 stats = afwMath.makeStatistics(mi, afwMath.MEANCLIP | afwMath.STDEVCLIP)
1000 mean, sigma = stats.getValue(afwMath.MEANCLIP), stats.getValue(afwMath.STDEVCLIP)
1001 # Get array of pixels
1002 EDGEBIT = exp.maskedImage.mask.getPlaneBitMask("EDGE")
1003 imgData = mi.image.array[(mi.mask.array & EDGEBIT) == 0].flatten()
1004 edgeData = mi.image.array[(mi.mask.array & EDGEBIT) != 0].flatten()
1006 thrUpper = mean + nSigmaUsed*sigma
1007 thrLower = mean - nSigmaUsed*sigma
1009 nRight = len(imgData[imgData > thrUpper])
1010 nLeft = len(imgData[imgData < thrLower])
1012 nsig = nSigmaUsed + 1.2 # add something small so the edge of the plot is out from level used
1013 leftEdge = mean - nsig * nSigmaUsed*sigma
1014 rightEdge = mean + nsig * nSigmaUsed*sigma
1015 nbins = np.linspace(leftEdge, rightEdge, 1000)
1016 ey, bin_borders, patches = a.hist(edgeData, histtype='step', bins=nbins,
1017 lw=1, edgecolor='red')
1018 y, bin_borders, patches = a.hist(imgData, histtype='step', bins=nbins,
1019 lw=3, edgecolor='blue')
1021 # Report number of entries in over- and under-flow
1022 # bins, i.e. off the edges of the histogram
1023 nOverflow = len(imgData[imgData > rightEdge])
1024 nUnderflow = len(imgData[imgData < leftEdge])
1026 # Put v-lines and textboxes in
1027 a.axvline(thrUpper, c='k')
1028 a.axvline(thrLower, c='k')
1029 msg = f"{amp.getName()}\nmean:{mean: .2f}\n$\\sigma$:{sigma: .2f}"
1030 a.text(0.65, 0.6, msg, transform=a.transAxes, fontsize=11)
1031 msg = f"nLeft:{nLeft}\nnRight:{nRight}\nnOverflow:{nOverflow}\nnUnderflow:{nUnderflow}"
1032 a.text(0.03, 0.6, msg, transform=a.transAxes, fontsize=11.5)
1034 # set axis limits and scales
1035 a.set_ylim([1., 1.7*np.max(y)])
1036 lPlot, rPlot = a.get_xlim()
1037 a.set_xlim(np.array([lPlot, rPlot]))
1038 a.set_yscale('log')
1039 a.set_xlabel("ADU/s")
1040 fig.show()
1041 prompt = "Press Enter or c to continue [chp]..."
1042 while True:
1043 ans = input(prompt).lower()
1044 if ans in ("", " ", "c",):
1045 break
1046 elif ans in ("p", ):
1047 import pdb
1048 pdb.set_trace()
1049 elif ans in ("h", ):
1050 print("[h]elp [c]ontinue [p]db")
1051 plt.close()
1054class MeasureDefectsCombinedConnections(pipeBase.PipelineTaskConnections,
1055 dimensions=("instrument", "detector")):
1056 inputExp = cT.Input(
1057 name="dark",
1058 doc="Input ISR-processed combined exposure to measure.",
1059 storageClass="ExposureF",
1060 dimensions=("instrument", "detector"),
1061 multiple=False,
1062 isCalibration=True,
1063 )
1064 camera = cT.PrerequisiteInput(
1065 name='camera',
1066 doc="Camera associated with this exposure.",
1067 storageClass="Camera",
1068 dimensions=("instrument", ),
1069 isCalibration=True,
1070 )
1072 outputDefects = cT.Output(
1073 name="cpDefectsFromDark",
1074 doc="Output measured defects.",
1075 storageClass="Defects",
1076 dimensions=("instrument", "detector"),
1077 )
1080class MeasureDefectsCombinedTaskConfig(MeasureDefectsTaskConfig,
1081 pipelineConnections=MeasureDefectsCombinedConnections):
1082 """Configuration for measuring defects from combined exposures.
1083 """
1084 pass
1087class MeasureDefectsCombinedTask(MeasureDefectsTask):
1088 """Task to measure defects in combined images."""
1090 ConfigClass = MeasureDefectsCombinedTaskConfig
1091 _DefaultName = "cpDefectMeasureCombined"
1094class MeasureDefectsCombinedWithFilterConnections(pipeBase.PipelineTaskConnections,
1095 dimensions=("instrument", "detector", "physical_filter")):
1096 """Task to measure defects in combined flats under a certain filter."""
1097 inputExp = cT.Input(
1098 name="flat",
1099 doc="Input ISR-processed combined exposure to measure.",
1100 storageClass="ExposureF",
1101 dimensions=("instrument", "detector", "physical_filter"),
1102 multiple=False,
1103 isCalibration=True,
1104 )
1105 camera = cT.PrerequisiteInput(
1106 name='camera',
1107 doc="Camera associated with this exposure.",
1108 storageClass="Camera",
1109 dimensions=("instrument", ),
1110 isCalibration=True,
1111 )
1113 outputDefects = cT.Output(
1114 name="cpDefectsFromFlat",
1115 doc="Output measured defects.",
1116 storageClass="Defects",
1117 dimensions=("instrument", "detector", "physical_filter"),
1118 )
1121class MeasureDefectsCombinedWithFilterTaskConfig(
1122 MeasureDefectsTaskConfig,
1123 pipelineConnections=MeasureDefectsCombinedWithFilterConnections):
1124 """Configuration for measuring defects from combined exposures.
1125 """
1126 pass
1129class MeasureDefectsCombinedWithFilterTask(MeasureDefectsTask):
1130 """Task to measure defects in combined images."""
1132 ConfigClass = MeasureDefectsCombinedWithFilterTaskConfig
1133 _DefaultName = "cpDefectMeasureWithFilterCombined"
1136class MergeDefectsConnections(pipeBase.PipelineTaskConnections,
1137 dimensions=("instrument", "detector")):
1138 inputDefects = cT.Input(
1139 name="singleExpDefects",
1140 doc="Measured defect lists.",
1141 storageClass="Defects",
1142 dimensions=("instrument", "detector", "exposure",),
1143 multiple=True,
1144 )
1145 camera = cT.PrerequisiteInput(
1146 name='camera',
1147 doc="Camera associated with these defects.",
1148 storageClass="Camera",
1149 dimensions=("instrument", ),
1150 isCalibration=True,
1151 )
1153 mergedDefects = cT.Output(
1154 name="defects",
1155 doc="Final merged defects.",
1156 storageClass="Defects",
1157 dimensions=("instrument", "detector"),
1158 multiple=False,
1159 isCalibration=True,
1160 )
1163class MergeDefectsTaskConfig(pipeBase.PipelineTaskConfig,
1164 pipelineConnections=MergeDefectsConnections):
1165 """Configuration for merging single exposure defects.
1166 """
1168 assertSameRun = pexConfig.Field(
1169 dtype=bool,
1170 doc=("Ensure that all visits are from the same run? Raises if this is not the case, or "
1171 "if the run key isn't found."),
1172 default=False, # false because most obs_packages don't have runs. obs_lsst/ts8 overrides this.
1173 )
1174 ignoreFilters = pexConfig.Field(
1175 dtype=bool,
1176 doc=("Set the filters used in the CALIB_ID to NONE regardless of the filters on the input"
1177 " images. Allows mixing of filters in the input flats. Set to False if you think"
1178 " your defects might be chromatic and want to have registry support for varying"
1179 " defects with respect to filter."),
1180 default=True,
1181 )
1182 nullFilterName = pexConfig.Field(
1183 dtype=str,
1184 doc=("The name of the null filter if ignoreFilters is True. Usually something like NONE or EMPTY"),
1185 default="NONE",
1186 )
1187 combinationMode = pexConfig.ChoiceField(
1188 doc="Which types of defects to identify",
1189 dtype=str,
1190 default="FRACTION",
1191 allowed={
1192 "AND": "Logical AND the pixels found in each visit to form set ",
1193 "OR": "Logical OR the pixels found in each visit to form set ",
1194 "FRACTION": "Use pixels found in more than config.combinationFraction of visits ",
1195 }
1196 )
1197 combinationFraction = pexConfig.RangeField(
1198 dtype=float,
1199 doc=("The fraction (0..1) of visits in which a pixel was found to be defective across"
1200 " the visit list in order to be marked as a defect. Note, upper bound is exclusive, so use"
1201 " mode AND to require pixel to appear in all images."),
1202 default=0.7,
1203 min=0,
1204 max=1,
1205 )
1206 nPixBorderUpDown = pexConfig.Field(
1207 dtype=int,
1208 doc=("Width (in pixels) of CCD top and bottom edges set as defects"
1209 "if edgesAsDefects is True."),
1210 default=5,
1211 )
1212 nPixBorderLeftRight = pexConfig.Field(
1213 dtype=int,
1214 doc=("Width (in pixels) of CCD left and right edges set as defects"
1215 "if edgesAsDefects is True."),
1216 default=5,
1217 )
1218 nPixBorderUpDownITL = pexConfig.Field(
1219 dtype=int,
1220 doc=("Width (in pixels) of ITL CCD top and bottom edges set as defects"
1221 "if edgesAsDefects is True."),
1222 default=5,
1223 )
1224 nPixBorderLeftRightITL = pexConfig.Field(
1225 dtype=int,
1226 doc=("Width (in pixels) of ITL CCD left and right edges set as defects"
1227 "if edgesAsDefects is True."),
1228 default=5,
1229 )
1230 nPixBorderUpDownE2V = pexConfig.Field(
1231 dtype=int,
1232 doc=("Width (in pixels) of E2V CCD top and bottom edges set as defects"
1233 "if edgesAsDefects is True."),
1234 default=5,
1235 )
1236 nPixBorderLeftRightE2V = pexConfig.Field(
1237 dtype=int,
1238 doc=("Width (in pixels) of E2V CCD left and right edges set as defects"
1239 "if edgesAsDefects is True."),
1240 default=5,
1241 )
1242 edgesAsDefects = pexConfig.Field(
1243 dtype=bool,
1244 doc="Mark all edge pixels, as defined by nPixBorder[UpDown, LeftRight], as defects.",
1245 default=False,
1246 )
1249class MergeDefectsTask(pipeBase.PipelineTask):
1250 """Merge the defects from multiple exposures.
1251 """
1253 ConfigClass = MergeDefectsTaskConfig
1254 _DefaultName = 'cpDefectMerge'
1256 def run(self, inputDefects, camera):
1257 """Merge a list of single defects to find the common defect regions.
1259 Parameters
1260 ----------
1261 inputDefects : `list` [`lsst.ip.isr.Defects`]
1262 Partial defects from a single exposure.
1263 camera : `lsst.afw.cameraGeom.Camera`
1264 Camera to use for metadata.
1266 Returns
1267 -------
1268 results : `lsst.pipe.base.Struct`
1269 Results struct containing:
1271 ``mergedDefects``
1272 The defects merged from the input lists
1273 (`lsst.ip.isr.Defects`).
1274 """
1275 detectorId = inputDefects[0].getMetadata().get('DETECTOR', None)
1276 if detectorId is None:
1277 raise RuntimeError("Cannot identify detector id.")
1278 detector = camera[detectorId]
1280 imageTypes = set()
1281 for inDefect in inputDefects:
1282 imageType = inDefect.getMetadata().get('cpDefectGenImageType', 'UNKNOWN')
1283 imageTypes.add(imageType)
1285 # Determine common defect pixels separately for each input image type.
1286 splitDefects = list()
1287 for imageType in imageTypes:
1288 sumImage = afwImage.MaskedImageF(detector.getBBox())
1289 count = 0
1290 for inDefect in inputDefects:
1291 if imageType == inDefect.getMetadata().get('cpDefectGenImageType', 'UNKNOWN'):
1292 count += 1
1293 for defect in inDefect:
1294 sumImage.image[defect.getBBox()] += 1.0
1295 sumImage /= count
1296 nDetected = len(np.where(sumImage.getImage().getArray() > 0)[0])
1297 self.log.info("Pre-merge %s pixels with non-zero detections for %s", nDetected, imageType)
1299 if self.config.combinationMode == 'AND':
1300 threshold = 1.0
1301 elif self.config.combinationMode == 'OR':
1302 threshold = 0.0
1303 elif self.config.combinationMode == 'FRACTION':
1304 threshold = self.config.combinationFraction
1305 else:
1306 raise RuntimeError(f"Got unsupported combinationMode {self.config.combinationMode}")
1307 indices = np.where(sumImage.getImage().getArray() > threshold)
1308 BADBIT = sumImage.getMask().getPlaneBitMask('BAD')
1309 sumImage.getMask().getArray()[indices] |= BADBIT
1310 self.log.info("Post-merge %s pixels marked as defects for %s", len(indices[0]), imageType)
1311 partialDefect = Defects.fromMask(sumImage, 'BAD')
1312 splitDefects.append(partialDefect)
1314 # Do final combination of separate image types
1315 finalImage = afwImage.MaskedImageF(detector.getBBox())
1316 for inDefect in splitDefects:
1317 for defect in inDefect:
1318 finalImage.image[defect.getBBox()] += 1
1319 finalImage /= len(splitDefects)
1320 nDetected = len(np.where(finalImage.getImage().getArray() > 0)[0])
1321 self.log.info("Pre-final merge %s pixels with non-zero detections", nDetected)
1323 # This combination is the OR of all image types
1324 threshold = 0.0
1325 indices = np.where(finalImage.getImage().getArray() > threshold)
1326 BADBIT = finalImage.getMask().getPlaneBitMask('BAD')
1327 finalImage.getMask().getArray()[indices] |= BADBIT
1328 self.log.info("Post-final merge %s pixels marked as defects", len(indices[0]))
1330 if self.config.edgesAsDefects:
1331 self.log.info("Masking edge pixels as defects.")
1332 detectorType = detector.getPhysicalType()
1333 if detectorType == 'E2V':
1334 nPixBorderLeftRight = self.config.nPixBorderLeftRightE2V
1335 elif 'ITL' in detectorType:
1336 nPixBorderLeftRight = self.config.nPixBorderLeftRightITL
1337 else:
1338 nPixBorderLeftRight = self.config.nPixBorderLeftRight
1340 if detectorType == 'E2V':
1341 nPixBorderUpDown = self.config.nPixBorderUpDownE2V
1342 elif 'ITL' in detectorType:
1343 nPixBorderUpDown = self.config.nPixBorderUpDownITL
1344 else:
1345 nPixBorderUpDown = self.config.nPixBorderUpDown
1347 # This code follows the pattern from isrTask.maskEdges().
1348 if nPixBorderLeftRight > 0:
1349 box = detector.getBBox()
1350 subImage = finalImage[box]
1351 box.grow(Extent2I(-nPixBorderLeftRight, 0))
1352 SourceDetectionTask.setEdgeBits(subImage, box, BADBIT)
1353 if nPixBorderUpDown > 0:
1354 box = detector.getBBox()
1355 subImage = finalImage[box]
1356 box.grow(Extent2I(0, -nPixBorderUpDown))
1357 SourceDetectionTask.setEdgeBits(subImage, box, BADBIT)
1359 merged = Defects.fromMask(finalImage, 'BAD')
1360 merged.updateMetadataFromExposures(inputDefects)
1361 merged.updateMetadata(camera=camera, detector=detector, filterName=None,
1362 setCalibId=True, setDate=True)
1364 return pipeBase.Struct(
1365 mergedDefects=merged,
1366 )
1368# Subclass the MergeDefects task to reduce the input dimensions
1369# from ("instrument", "detector", "exposure") to
1370# ("instrument", "detector").
1373class MergeDefectsCombinedConnections(pipeBase.PipelineTaskConnections,
1374 dimensions=("instrument", "detector")):
1375 inputDarkDefects = cT.Input(
1376 name="cpDefectsFromDark",
1377 doc="Measured defect lists.",
1378 storageClass="Defects",
1379 dimensions=("instrument", "detector",),
1380 multiple=True,
1381 )
1382 inputBiasDefects = cT.Input(
1383 name="cpDefectsFromBias",
1384 doc="Additional measured defect lists.",
1385 storageClass="Defects",
1386 dimensions=("instrument", "detector",),
1387 multiple=True,
1388 )
1389 inputFlatDefects = cT.Input(
1390 name="cpDefectsFromFlat",
1391 doc="Additional measured defect lists.",
1392 storageClass="Defects",
1393 dimensions=("instrument", "detector", "physical_filter"),
1394 multiple=True,
1395 )
1396 inputManualDefects = cT.Input(
1397 name="cpManualDefects",
1398 doc="Additional manual defects.",
1399 storageClass="Defects",
1400 dimensions=("instrument", "detector"),
1401 multiple=True,
1402 isCalibration=True,
1403 )
1404 camera = cT.PrerequisiteInput(
1405 name='camera',
1406 doc="Camera associated with these defects.",
1407 storageClass="Camera",
1408 dimensions=("instrument", ),
1409 isCalibration=True,
1410 )
1412 mergedDefects = cT.Output(
1413 name="defects",
1414 doc="Final merged defects.",
1415 storageClass="Defects",
1416 dimensions=("instrument", "detector"),
1417 multiple=False,
1418 isCalibration=True,
1419 )
1421 def __init__(self, *, config=None):
1422 super().__init__(config=config)
1424 if config.doManualDefects is not True:
1425 del self.inputManualDefects
1428class MergeDefectsCombinedTaskConfig(MergeDefectsTaskConfig,
1429 pipelineConnections=MergeDefectsCombinedConnections):
1430 """Configuration for merging defects from combined exposure.
1431 """
1432 doManualDefects = pexConfig.Field(
1433 dtype=bool,
1434 doc="Apply manual defects?",
1435 default=False,
1436 )
1438 def validate(self):
1439 super().validate()
1440 if self.combinationMode != 'OR': 1440 ↛ 1441line 1440 didn't jump to line 1441 because the condition on line 1440 was never true
1441 raise ValueError("combinationMode must be 'OR'")
1444class MergeDefectsCombinedTask(MergeDefectsTask):
1445 """Task to measure defects in combined images."""
1447 ConfigClass = MergeDefectsCombinedTaskConfig
1448 _DefaultName = "cpMergeDefectsCombined"
1450 @staticmethod
1451 def chooseBest(inputs):
1452 """Select the input with the most exposures used."""
1453 best = 0
1454 if len(inputs) > 1:
1455 nInput = 0
1456 for num, exp in enumerate(inputs):
1457 # This technically overcounts by a factor of 3.
1458 N = len([k for k, v in exp.getMetadata().toDict().items() if "CPP_INPUT_" in k])
1459 if N > nInput:
1460 best = num
1461 nInput = N
1462 return inputs[best]
1464 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1465 inputs = butlerQC.get(inputRefs)
1466 # Turn inputFlatDefects and inputDarkDefects into a list which
1467 # is what MergeDefectsTask expects. If there are multiple,
1468 # use the one with the most inputs.
1469 tempList = [self.chooseBest(inputs['inputFlatDefects']),
1470 self.chooseBest(inputs['inputDarkDefects']),
1471 self.chooseBest(inputs['inputBiasDefects'])]
1473 if "inputManualDefects" in inputs.keys():
1474 tempList.extend(inputs["inputManualDefects"])
1476 # Rename inputDefects
1477 inputsCombined = {'inputDefects': tempList, 'camera': inputs['camera']}
1479 outputs = super().run(**inputsCombined)
1480 butlerQC.put(outputs, outputRefs)