Coverage for python / lsst / cp / pipe / cpDefects.py: 17%

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# 

22 

23__all__ = ['MeasureDefectsTaskConfig', 'MeasureDefectsTask', 

24 'MergeDefectsTaskConfig', 'MergeDefectsTask', 

25 'MeasureDefectsCombinedTaskConfig', 'MeasureDefectsCombinedTask', 

26 'MeasureDefectsCombinedWithFilterTaskConfig', 'MeasureDefectsCombinedWithFilterTask', 

27 'MergeDefectsCombinedTaskConfig', 'MergeDefectsCombinedTask', ] 

28 

29import numpy as np 

30 

31import lsst.pipe.base as pipeBase 

32import lsst.pipe.base.connectionTypes as cT 

33 

34from lsstDebug import getDebugFrame 

35import lsst.pex.config as pexConfig 

36 

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 

46 

47from .utils import bin_flat, FlatGradientFitter 

48from lsst.ip.isr import FlatGradient 

49 

50 

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 ) 

67 

68 outputDefects = cT.Output( 

69 name="singleExpDefects", 

70 doc="Output measured defects.", 

71 storageClass="Defects", 

72 dimensions=("instrument", "detector", "exposure"), 

73 ) 

74 

75 

76class MeasureDefectsTaskConfig(pipeBase.PipelineTaskConfig, 

77 pipelineConnections=MeasureDefectsConnections): 

78 """Configuration for measuring defects from a list of exposures 

79 """ 

80 

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 badOnAndOffPixelColumnThreshold = pexConfig.Field( 

157 dtype=int, 

158 doc=("If BPC is the set of all the bad pixels in a given column (not necessarily consecutive) " 

159 "and the size of BPC is at least 'badOnAndOffPixelColumnThreshold', all the pixels between the " 

160 "pixels that satisfy minY (BPC) and maxY (BPC) will be marked as bad, with 'Y' being the long " 

161 "axis of the amplifier (and 'X' the other axis, which for a column is a constant for all " 

162 "pixels in the set BPC). If there are more than 'goodPixelColumnGapThreshold' consecutive " 

163 "non-bad pixels in BPC, an exception to the above is made and those consecutive " 

164 "'goodPixelColumnGapThreshold' are not marked as bad."), 

165 default=50, 

166 ) 

167 goodPixelColumnGapThreshold = pexConfig.Field( 

168 dtype=int, 

169 doc=("Size, in pixels, of usable consecutive pixels in a column with on and off bad pixels (see " 

170 "'badOnAndOffPixelColumnThreshold')."), 

171 default=30, 

172 ) 

173 badPixelsToFillColumnThreshold = pexConfig.Field( 

174 dtype=float, 

175 doc=("If the number of bad pixels in an amplifier column is above this threshold " 

176 "then the full amplifier column will be marked bad. This operation is performed after " 

177 "any merging of blinking columns performed with badOnAndOffPixelColumnThreshold. If this" 

178 "value is less than 0 then no bad column filling will be performed."), 

179 default=-1, 

180 ) 

181 saturatedColumnMask = pexConfig.Field( 

182 dtype=str, 

183 default="SAT", 

184 doc="Saturated mask plane for dilation.", 

185 ) 

186 saturatedColumnDilationRadius = pexConfig.Field( 

187 dtype=int, 

188 doc=("Dilation radius (along rows) to use to expand saturated columns " 

189 "to mitigate glow."), 

190 default=0, 

191 ) 

192 saturatedPixelsToFillColumnThreshold = pexConfig.Field( 

193 dtype=int, 

194 doc=("If the number of saturated pixels in an amplifier column is above this threshold " 

195 "then the full amplifier column will be marked bad. If this value is less than 0" 

196 "then no saturated column filling will be performed."), 

197 default=-1, 

198 ) 

199 e2vMidlineBreakNRow = pexConfig.Field( 

200 dtype=int, 

201 doc="E2V midline break number of midline break rows at the bottom of the top amps " 

202 "or the top of the bottom amps to always mask or ignore. Number of rows will " 

203 "be twice this config value. Only used if detector is E2V type. Set to <=0 " 

204 "to treat E2V midline break rows like any other flat-field pixel.", 

205 default=1, 

206 ) 

207 e2vMidlineBreakOption = pexConfig.ChoiceField( 

208 dtype=str, 

209 doc="How should the E2V midline break be treated? Only used if e2vMidlineBreakNRow > 0.", 

210 default="NEVERMASK", 

211 allowed={ 

212 "NEVERMASK": "Never mask the E2V midline break, no matter the flat values.", 

213 "MASK": "Always mask the E2V midline break.", 

214 }, 

215 ) 

216 ampGradientBinFactor = pexConfig.Field( 

217 dtype=int, 

218 doc="Binning factor used for fitting per-amp focal plane gradient.", 

219 default=8, 

220 ) 

221 ampGradientBoundarySize = pexConfig.Field( 

222 dtype=int, 

223 doc="Amp boundary exclusion size for fitting per-amp focal plane gradient.", 

224 default=20, 

225 ) 

226 ampGradientNodes = pexConfig.Field( 

227 dtype=int, 

228 doc="Number of spline nodes for per-amp focal plane gradient.", 

229 default=4, 

230 ) 

231 

232 def validate(self): 

233 super().validate() 

234 if self.nSigmaBright < 0.0: 

235 raise ValueError("nSigmaBright must be above 0.0.") 

236 if self.nSigmaDark > 0.0: 

237 raise ValueError("nSigmaDark must be below 0.0.") 

238 

239 

240class MeasureDefectsTask(pipeBase.PipelineTask): 

241 """Measure the defects from one exposure. 

242 """ 

243 

244 ConfigClass = MeasureDefectsTaskConfig 

245 _DefaultName = 'cpDefectMeasure' 

246 

247 def run(self, inputExp, camera): 

248 """Measure one exposure for defects. 

249 

250 Parameters 

251 ---------- 

252 inputExp : `lsst.afw.image.Exposure` 

253 Exposure to examine. 

254 camera : `lsst.afw.cameraGeom.Camera` 

255 Camera to use for metadata. 

256 

257 Returns 

258 ------- 

259 results : `lsst.pipe.base.Struct` 

260 Results struct containing: 

261 

262 ``outputDefects`` 

263 The defects measured from this exposure 

264 (`lsst.ip.isr.Defects`). 

265 """ 

266 detector = inputExp.getDetector() 

267 try: 

268 filterName = inputExp.getFilter().physicalLabel 

269 except AttributeError: 

270 filterName = None 

271 

272 defects = self._findHotAndColdPixels(inputExp) 

273 

274 datasetType = inputExp.getMetadata().get('IMGTYPE', 'UNKNOWN') 

275 msg = "Found %s defects containing %s pixels in %s" 

276 self.log.info(msg, len(defects), self._nPixFromDefects(defects), datasetType) 

277 

278 defects.updateMetadataFromExposures([inputExp]) 

279 defects.updateMetadata(camera=camera, detector=detector, filterName=filterName, 

280 setCalibId=True, setDate=True, 

281 cpDefectGenImageType=datasetType) 

282 

283 return pipeBase.Struct( 

284 outputDefects=defects, 

285 ) 

286 

287 @staticmethod 

288 def _nPixFromDefects(defects): 

289 """Count pixels in a defect. 

290 

291 Parameters 

292 ---------- 

293 defects : `lsst.ip.isr.Defects` 

294 Defects to measure. 

295 

296 Returns 

297 ------- 

298 nPix : `int` 

299 Number of defect pixels. 

300 """ 

301 nPix = 0 

302 for defect in defects: 

303 nPix += defect.getBBox().getArea() 

304 return nPix 

305 

306 def getVampirePixels(self, ampImg): 

307 """Find vampire pixels (bright pixels in flats) and get footprint of 

308 extended area around them, 

309 

310 Parameters 

311 ---------- 

312 ampImg : `lsst.afw.image._maskedImage.MaskedImageF` 

313 The amplifier masked image to do the vampire pixels search on. 

314 

315 Returns 

316 ------- 

317 fs_grow : `lsst.afw.detection._detection.FootprintSet` 

318 The footprint set of areas around vampire pixels in the amplifier. 

319 """ 

320 

321 # Find bright pixels 

322 thresh = afwDetection.Threshold(self.config.thresholdVampirePixels) 

323 # Bright pixels footprint grown by a radius of radiusVampire pixels 

324 fs = afwDetection.FootprintSet(ampImg, thresh) 

325 fs_grow = afwDetection.FootprintSet(fs, rGrow=self.config.radiusVampirePixels, isotropic=True) 

326 

327 return fs_grow 

328 

329 def _setE2VMidline(self, amp, ampImg): 

330 """Set the E2V midline pixels on this amplifier. 

331 

332 If the task is configure to always mask it will set the midline rows to 

333 an illegal value; if it is to never mask it will set the midline rows 

334 to the median value. 

335 

336 Parameters 

337 ---------- 

338 amp : `lsst.afw.cameraGeom.Amplifier` 

339 Amplifier object. 

340 ampImg : `lsst.afw.image.ImageF` 

341 Image data for the amplifier. 

342 """ 

343 if self.config.e2vMidlineBreakNRow <= 0: 

344 return 

345 

346 if self.config.e2vMidlineBreakOption == "NEVERMASK": 

347 # Setting the midline values to the median value ensures 

348 # they will never trigger the defect finding threshold. 

349 value = np.nanmedian(ampImg.image.array) 

350 else: 

351 # Setting the midline values to a large negative value 

352 # ensures they will always trigger the defect finding 

353 # threshold. 

354 value = -1e30 

355 

356 if amp.getName().startswith("C0"): 

357 # This is a bottom amplifier, so change the top row(s). 

358 ampImg.image.array[-self.config.e2vMidlineBreakNRow:, :] = value 

359 else: 

360 # This is a top amplifier, so change the bottom row(s). 

361 ampImg.image.array[:self.config.e2vMidlineBreakNRow, :] = value 

362 

363 def _findHotAndColdPixels(self, exp): 

364 """Find hot and cold pixels in an image. 

365 

366 Using config-defined thresholds on a per-amp basis, mask 

367 pixels that are nSigma above threshold in dark frames (hot 

368 pixels), or nSigma away from the clipped mean in flats (hot & 

369 cold pixels). 

370 

371 Parameters 

372 ---------- 

373 exp : `lsst.afw.image.exposure.Exposure` 

374 The exposure in which to find defects. 

375 

376 Returns 

377 ------- 

378 defects : `lsst.ip.isr.Defects` 

379 The defects found in the image. 

380 """ 

381 self._setEdgeBits(exp) 

382 

383 # the detection polarity for afwDetection, True for positive, 

384 # False for negative, and therefore True for darks as they only have 

385 # bright pixels, and both for flats, as they have bright and dark pix 

386 footprintList = [] 

387 

388 hotPixelCount = {} 

389 coldPixelCount = {} 

390 

391 detector = exp.getDetector() 

392 

393 if self.config.fitAmpGradient: 

394 # 1. Bin flat. 

395 binnedExp = bin_flat( 

396 PhotonTransferCurveDataset(), 

397 exp, 

398 bin_factor=self.config.ampGradientBinFactor, 

399 amp_boundary=self.config.ampGradientBoundarySize, 

400 apply_gains=False, 

401 ) 

402 

403 # 2. Prepare transformations. 

404 transform = detector.getTransform( 

405 cameraGeom.PIXELS, 

406 cameraGeom.FOCAL_PLANE, 

407 ) 

408 

409 # 2a. Do transformation for binned coordinates. 

410 binnedXy = np.vstack((binnedExp["xd"], binnedExp["yd"])) 

411 binnedXf, binnedYf = np.vsplit( 

412 transform.getMapping().applyForward(binnedXy.astype(np.float64)), 

413 2, 

414 ) 

415 binnedXf = binnedXf.ravel() 

416 binnedYf = binnedYf.ravel() 

417 binnedRadius = np.sqrt(binnedXf**2. + binnedYf**2.) 

418 

419 # 2b. Do transformation for full detector coordinates. 

420 xx = np.arange(exp.image.array.shape[1], dtype=np.int64) 

421 yy = np.arange(exp.image.array.shape[0], dtype=np.int64) 

422 x, y = np.meshgrid(xx, yy) 

423 x = x.ravel() 

424 y = y.ravel() 

425 

426 xy = np.vstack((x, y)) 

427 xf, yf = np.vsplit(transform.getMapping().applyForward(xy.astype(np.float64)), 2) 

428 xf = xf.ravel() 

429 yf = yf.ravel() 

430 

431 # Measure the focal plane gradient on the binned image and remove 

432 # on the full image, amp by amp. 

433 for i, amp in enumerate(detector): 

434 # Fit the focal plane radial gradient in the binned image. 

435 binnedInAmp = (binnedExp["amp_index"] == i) 

436 

437 nodes = np.linspace( 

438 np.min(binnedRadius[binnedInAmp]), 

439 np.max(binnedRadius[binnedInAmp]), 

440 self.config.ampGradientNodes, 

441 ) 

442 

443 norm = np.nanpercentile(binnedExp["value"][binnedInAmp], 95.0) 

444 

445 fitter = FlatGradientFitter( 

446 nodes, 

447 binnedXf[binnedInAmp], 

448 binnedYf[binnedInAmp], 

449 binnedExp["value"][binnedInAmp]/norm, 

450 np.array([]), 

451 constrain_zero=False, 

452 ) 

453 p0 = fitter.compute_p0() 

454 pars = fitter.fit(p0) 

455 

456 gradient = FlatGradient() 

457 gradient.setParameters( 

458 radialSplineNodes=nodes, 

459 radialSplineValues=pars[fitter.indices["spline"]], 

460 ) 

461 

462 # Apply the focal plane gradient to the full image. 

463 pixelsInAmp = amp.getBBox().contains(x, y) 

464 

465 model = gradient.computeRadialSplineModelXY(xf[pixelsInAmp], yf[pixelsInAmp]) 

466 

467 exp.image.array[y[pixelsInAmp], x[pixelsInAmp]] /= model 

468 

469 maskedIm = exp.maskedImage 

470 

471 for amp in exp.getDetector(): 

472 ampName = amp.getName() 

473 

474 hotPixelCount[ampName] = 0 

475 coldPixelCount[ampName] = 0 

476 

477 ampImg = maskedIm[amp.getBBox()].clone() 

478 

479 # crop ampImage depending on where the amp lies in the image 

480 if self.config.nPixBorderLeftRight: 

481 if ampImg.getX0() == 0: 

482 ampImg = ampImg[self.config.nPixBorderLeftRight:, :, afwImage.LOCAL] 

483 else: 

484 ampImg = ampImg[:-self.config.nPixBorderLeftRight, :, afwImage.LOCAL] 

485 if self.config.nPixBorderUpDown: 

486 if ampImg.getY0() == 0: 

487 ampImg = ampImg[:, self.config.nPixBorderUpDown:, afwImage.LOCAL] 

488 else: 

489 ampImg = ampImg[:, :-self.config.nPixBorderUpDown, afwImage.LOCAL] 

490 

491 if self._getNumGoodPixels(ampImg) == 0: # amp contains no usable pixels 

492 continue 

493 

494 if self.config.doVampirePixels: 

495 # This is only applied in LSSTComCam flatBootstrap pipeline 

496 footprintSet_VampirePixel = self.getVampirePixels(ampImg) 

497 footprintSet_VampirePixel.setMask(maskedIm.mask, ("BAD")) 

498 

499 # Remove a background estimate 

500 meanClip = afwMath.makeStatistics(ampImg, afwMath.MEANCLIP, ).getValue() 

501 ampImg -= meanClip 

502 

503 # Determine thresholds 

504 stDev = afwMath.makeStatistics(ampImg, afwMath.STDEVCLIP, ).getValue() 

505 expTime = exp.getInfo().getVisitInfo().getExposureTime() 

506 datasetType = exp.getMetadata().get('IMGTYPE', 'UNKNOWN') 

507 if np.isnan(expTime): 

508 self.log.warning("expTime=%s for AMP %s in %s. Setting expTime to 1 second", 

509 expTime, ampName, datasetType) 

510 expTime = 1. 

511 thresholdType = self.config.thresholdType 

512 if thresholdType == 'VALUE': 

513 # LCA-128 and eoTest: bright/hot pixels in dark images are 

514 # defined as any pixel with more than 5 e-/s of dark current. 

515 # We scale by the exposure time. 

516 if datasetType.lower() == 'dark': 

517 # hot pixel threshold 

518 valueThreshold = self.config.darkCurrentThreshold*expTime/amp.getGain() 

519 elif datasetType.lower() == 'bias': 

520 # hot pixel threshold, no exposure time. 

521 valueThreshold = self.config.biasThreshold 

522 else: 

523 # LCA-128 and eoTest: dark/cold pixels in flat images as 

524 # defined as any pixel with photoresponse <80% of 

525 # the mean (at 500nm). 

526 

527 # We subtracted the mean above, so the threshold will be 

528 # negative cold pixel threshold. 

529 valueThreshold = (self.config.fracThresholdFlat-1)*meanClip 

530 # Find equivalent sigma values. 

531 if stDev == 0.0: 

532 self.log.warning("stDev=%s for AMP %s in %s. Setting nSigma to inf.", 

533 stDev, ampName, datasetType) 

534 nSigmaList = [np.inf] 

535 else: 

536 nSigmaList = [valueThreshold/stDev] 

537 else: 

538 hotPixelThreshold = self.config.nSigmaBright 

539 coldPixelThreshold = self.config.nSigmaDark 

540 if datasetType.lower() == 'dark': 

541 nSigmaList = [hotPixelThreshold] 

542 valueThreshold = stDev*hotPixelThreshold 

543 elif datasetType.lower() == 'bias': 

544 self.log.warning( 

545 "Bias frame detected, but thresholdType == STDEV; not looking for defects.", 

546 ) 

547 return Defects.fromFootprintList([]) 

548 else: 

549 nSigmaList = [hotPixelThreshold, coldPixelThreshold] 

550 valueThreshold = [x*stDev for x in nSigmaList] 

551 

552 self.log.info("Image type: %s. Amp: %s. Threshold Type: %s. Sigma values and Pixel" 

553 "Values (hot and cold pixels thresholds): %s, %s", 

554 datasetType, ampName, thresholdType, nSigmaList, valueThreshold) 

555 

556 if datasetType.lower() == "flat" and exp.getDetector().getPhysicalType() == "E2V": 

557 self._setE2VMidline(amp, ampImg) 

558 

559 mergedSet = None 

560 for sigma in nSigmaList: 

561 nSig = np.abs(sigma) 

562 self.debugHistogram('ampFlux', ampImg, nSig, exp) 

563 polarity = {-1: False, 1: True}[np.sign(sigma)] 

564 

565 threshold = afwDetection.createThreshold(nSig, 'stdev', polarity=polarity) 

566 

567 try: 

568 footprintSet = afwDetection.FootprintSet(ampImg, threshold) 

569 except InvalidParameterError: 

570 # This occurs if the image sigma value is 0.0. 

571 # Let's mask the whole area. 

572 minValue = np.nanmin(ampImg.image.array) - 1.0 

573 threshold = afwDetection.createThreshold(minValue, 'value', polarity=True) 

574 footprintSet = afwDetection.FootprintSet(ampImg, threshold) 

575 

576 footprintSet.setMask(maskedIm.mask, ("DETECTED" if polarity else "DETECTED_NEGATIVE")) 

577 

578 if mergedSet is None: 

579 mergedSet = footprintSet 

580 else: 

581 mergedSet.merge(footprintSet) 

582 

583 if polarity: 

584 # hot pixels 

585 for fp in footprintSet.getFootprints(): 

586 hotPixelCount[ampName] += fp.getArea() 

587 else: 

588 # cold pixels 

589 for fp in footprintSet.getFootprints(): 

590 coldPixelCount[ampName] += fp.getArea() 

591 

592 if self.config.doVampirePixels: 

593 # Count the number of pixels masked 

594 vampirePixelCount = 0 

595 for fp in footprintSet_VampirePixel.getFootprints(): 

596 vampirePixelCount += fp.getArea() 

597 self.log.info("%s Vampire pixels are masked", vampirePixelCount) 

598 # Add vampire pixels to footprint set 

599 mergedSet.merge(footprintSet_VampirePixel) 

600 

601 footprintList += mergedSet.getFootprints() 

602 

603 self.debugView('defectMap', ampImg, 

604 Defects.fromFootprintList(mergedSet.getFootprints()), exp.getDetector()) 

605 

606 defects = Defects.fromFootprintList(footprintList) 

607 defects = self.dilateSaturatedColumns(exp, defects) 

608 defects, _ = self.maskBlocksIfIntermitentBadPixelsInColumn(defects) 

609 defects, count = self.maskBadColumns(exp, defects) 

610 # We want this to reflect the number of completely bad columns. 

611 defects.updateCounters(columns=count, hot=hotPixelCount, cold=coldPixelCount) 

612 

613 return defects 

614 

615 @staticmethod 

616 def _getNumGoodPixels(maskedIm, badMaskString="NO_DATA"): 

617 """Return the number of non-bad pixels in the image.""" 

618 nPixels = maskedIm.mask.array.size 

619 nBad = countMaskedPixels(maskedIm, badMaskString) 

620 return nPixels - nBad 

621 

622 def _setEdgeBits(self, exposureOrMaskedImage, maskplaneToSet='EDGE'): 

623 """Set edge bits on an exposure or maskedImage. 

624 

625 Raises 

626 ------ 

627 TypeError 

628 Raised if parameter ``exposureOrMaskedImage`` is an invalid type. 

629 """ 

630 if isinstance(exposureOrMaskedImage, afwImage.Exposure): 

631 mi = exposureOrMaskedImage.maskedImage 

632 elif isinstance(exposureOrMaskedImage, afwImage.MaskedImage): 

633 mi = exposureOrMaskedImage 

634 else: 

635 t = type(exposureOrMaskedImage) 

636 raise TypeError(f"Function supports exposure or maskedImage but not {t}") 

637 

638 MASKBIT = mi.mask.getPlaneBitMask(maskplaneToSet) 

639 if self.config.nPixBorderLeftRight: 

640 mi.mask[: self.config.nPixBorderLeftRight, :, afwImage.LOCAL] |= MASKBIT 

641 mi.mask[-self.config.nPixBorderLeftRight:, :, afwImage.LOCAL] |= MASKBIT 

642 if self.config.nPixBorderUpDown: 

643 mi.mask[:, : self.config.nPixBorderUpDown, afwImage.LOCAL] |= MASKBIT 

644 mi.mask[:, -self.config.nPixBorderUpDown:, afwImage.LOCAL] |= MASKBIT 

645 

646 def maskBlocksIfIntermitentBadPixelsInColumn(self, defects): 

647 """Mask blocks in a column if there are on-and-off bad pixels 

648 

649 If there's a column with on and off bad pixels, mask all the 

650 pixels in between, except if there is a large enough gap of 

651 consecutive good pixels between two bad pixels in the column. 

652 

653 Parameters 

654 ---------- 

655 defects : `lsst.ip.isr.Defects` 

656 The defects found in the image so far 

657 

658 Returns 

659 ------- 

660 defects : `lsst.ip.isr.Defects` 

661 If the number of bad pixels in a column is not larger or 

662 equal than self.config.badPixelColumnThreshold, the input 

663 list is returned. Otherwise, the defects list returned 

664 will include boxes that mask blocks of on-and-of pixels. 

665 badColumnCount : `int` 

666 Number of bad columns partially masked. 

667 """ 

668 badColumnCount = 0 

669 # Get the (x, y) values of each bad pixel in amp. 

670 coordinates = [] 

671 for defect in defects: 

672 bbox = defect.getBBox() 

673 x0, y0 = bbox.getMinX(), bbox.getMinY() 

674 deltaX0, deltaY0 = bbox.getDimensions() 

675 for j in np.arange(y0, y0+deltaY0): 

676 for i in np.arange(x0, x0 + deltaX0): 

677 coordinates.append((i, j)) 

678 

679 x, y = [], [] 

680 for coordinatePair in coordinates: 

681 x.append(coordinatePair[0]) 

682 y.append(coordinatePair[1]) 

683 

684 x = np.array(x) 

685 y = np.array(y) 

686 # Find the defects with same "x" (vertical) coordinate (column). 

687 unique, counts = np.unique(x, return_counts=True) 

688 multipleX = [] 

689 for (a, b) in zip(unique, counts): 

690 if b >= self.config.badOnAndOffPixelColumnThreshold: 

691 multipleX.append(a) 

692 if len(multipleX) != 0: 

693 defects = self._markBlocksInBadColumn(x, y, multipleX, defects) 

694 badColumnCount += 1 

695 

696 return defects, badColumnCount 

697 

698 def dilateSaturatedColumns(self, exp, defects): 

699 """Dilate saturated columns by a configurable amount. 

700 

701 Parameters 

702 ---------- 

703 exp : `lsst.afw.image.exposure.Exposure` 

704 The exposure in which to find defects. 

705 defects : `lsst.ip.isr.Defects` 

706 The defects found in the image so far 

707 

708 Returns 

709 ------- 

710 defects : `lsst.ip.isr.Defects` 

711 The expanded defects. 

712 """ 

713 if self.config.saturatedColumnDilationRadius <= 0: 

714 # This is a no-op. 

715 return defects 

716 

717 mask = afwImage.Mask.getPlaneBitMask(self.config.saturatedColumnMask) 

718 

719 satY, satX = np.where((exp.mask.array & mask) > 0) 

720 

721 if len(satX) == 0: 

722 # No saturated pixels, nothing to do. 

723 return defects 

724 

725 radius = self.config.saturatedColumnDilationRadius 

726 

727 with defects.bulk_update(): 

728 for index in range(len(satX)): 

729 minX = np.clip(satX[index] - radius, 0, None) 

730 maxX = np.clip(satX[index] + radius, None, exp.image.array.shape[1] - 1) 

731 s = Box2I(minimum=Point2I(minX, satY[index]), 

732 maximum=Point2I(maxX, satY[index])) 

733 defects.append(s) 

734 

735 return defects 

736 

737 def maskBadColumns(self, exp, defects): 

738 """Mask full amplifier columns if they are sufficiently bad. 

739 

740 Parameters 

741 ---------- 

742 defects : `lsst.ip.isr.Defects` 

743 The defects found in the image so far 

744 

745 Returns 

746 ------- 

747 exp : `lsst.afw.image.exposure.Exposure` 

748 The exposure in which to find defects. 

749 defects : `lsst.ip.isr.Defects` 

750 If the number of bad pixels in a column is not larger or 

751 equal than self.config.badPixelColumnThreshold, the input 

752 list is returned. Otherwise, the defects list returned 

753 will include boxes that mask blocks of on-and-of pixels. 

754 badColumnCount : `int` 

755 Number of bad columns masked. 

756 """ 

757 # Render the defects into an image. 

758 defectImage = afwImage.ImageI(exp.getBBox()) 

759 

760 for defect in defects: 

761 defectImage[defect.getBBox()] = 1 

762 

763 badColumnCount = 0 

764 

765 if self.config.badPixelsToFillColumnThreshold > 0: 

766 with defects.bulk_update(): 

767 for amp in exp.getDetector(): 

768 subImage = defectImage[amp.getBBox()].array 

769 nInCol = np.sum(subImage, axis=0) 

770 

771 badColIndices, = (nInCol >= self.config.badPixelsToFillColumnThreshold).nonzero() 

772 badColumns = badColIndices + amp.getBBox().getMinX() 

773 

774 for badColumn in badColumns: 

775 s = Box2I(minimum=Point2I(badColumn, amp.getBBox().getMinY()), 

776 maximum=Point2I(badColumn, amp.getBBox().getMaxY())) 

777 defects.append(s) 

778 

779 badColumnCount += len(badColIndices) 

780 

781 if self.config.saturatedPixelsToFillColumnThreshold > 0: 

782 mask = afwImage.Mask.getPlaneBitMask(self.config.saturatedColumnMask) 

783 

784 with defects.bulk_update(): 

785 for amp in exp.getDetector(): 

786 subMask = exp.mask[amp.getBBox()].array 

787 # Turn all the SAT bits into 1s 

788 subMask &= mask 

789 subMask[subMask > 0] = 1 

790 

791 nInCol = np.sum(subMask, axis=0) 

792 

793 badColIndices, = (nInCol >= self.config.saturatedPixelsToFillColumnThreshold).nonzero() 

794 badColumns = badColIndices + amp.getBBox().getMinX() 

795 

796 for badColumn in badColumns: 

797 s = Box2I(minimum=Point2I(badColumn, amp.getBBox().getMinY()), 

798 maximum=Point2I(badColumn, amp.getBBox().getMaxY())) 

799 defects.append(s) 

800 

801 badColumnCount += len(badColIndices) 

802 

803 return defects, badColumnCount 

804 

805 def _markBlocksInBadColumn(self, x, y, multipleX, defects): 

806 """Mask blocks in a column if number of on-and-off bad pixels is above 

807 threshold. 

808 

809 This function is called if the number of on-and-off bad pixels 

810 in a column is larger or equal than 

811 self.config.badOnAndOffPixelColumnThreshold. 

812 

813 Parameters 

814 --------- 

815 x : `list` 

816 Lower left x coordinate of defect box. x coordinate is 

817 along the short axis if amp. 

818 y : `list` 

819 Lower left y coordinate of defect box. x coordinate is 

820 along the long axis if amp. 

821 multipleX : list 

822 List of x coordinates in amp. with multiple bad pixels 

823 (i.e., columns with defects). 

824 defects : `lsst.ip.isr.Defects` 

825 The defcts found in the image so far 

826 

827 Returns 

828 ------- 

829 defects : `lsst.ip.isr.Defects` 

830 The defects list returned that will include boxes that 

831 mask blocks of on-and-of pixels. 

832 """ 

833 with defects.bulk_update(): 

834 goodPixelColumnGapThreshold = self.config.goodPixelColumnGapThreshold 

835 for x0 in multipleX: 

836 index = np.where(x == x0) 

837 multipleY = y[index] # multipleY and multipleX are in 1-1 correspondence. 

838 multipleY.sort() # Ensure that the y values are sorted to look for gaps. 

839 minY, maxY = np.min(multipleY), np.max(multipleY) 

840 # Next few lines: don't mask pixels in column if gap 

841 # of good pixels between two consecutive bad pixels is 

842 # larger or equal than 'goodPixelColumnGapThreshold'. 

843 diffIndex = np.where(np.diff(multipleY) >= goodPixelColumnGapThreshold)[0] 

844 if len(diffIndex) != 0: 

845 limits = [minY] # put the minimum first 

846 for gapIndex in diffIndex: 

847 limits.append(multipleY[gapIndex]) 

848 limits.append(multipleY[gapIndex+1]) 

849 limits.append(maxY) # maximum last 

850 for i in np.arange(0, len(limits)-1, 2):