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 

45from lsst.pex.exceptions import InvalidParameterError 

46 

47 

48class MeasureDefectsConnections(pipeBase.PipelineTaskConnections, 

49 dimensions=("instrument", "exposure", "detector")): 

50 inputExp = cT.Input( 

51 name="defectExps", 

52 doc="Input ISR-processed exposures to measure.", 

53 storageClass="Exposure", 

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

55 multiple=False 

56 ) 

57 camera = cT.PrerequisiteInput( 

58 name='camera', 

59 doc="Camera associated with this exposure.", 

60 storageClass="Camera", 

61 dimensions=("instrument", ), 

62 isCalibration=True, 

63 ) 

64 

65 outputDefects = cT.Output( 

66 name="singleExpDefects", 

67 doc="Output measured defects.", 

68 storageClass="Defects", 

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

70 ) 

71 

72 

73class MeasureDefectsTaskConfig(pipeBase.PipelineTaskConfig, 

74 pipelineConnections=MeasureDefectsConnections): 

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

76 """ 

77 

78 thresholdType = pexConfig.ChoiceField( 

79 dtype=str, 

80 doc=("Defects threshold type: ``STDEV`` or ``VALUE``. If ``VALUE``, cold pixels will be found " 

81 "in flats, and hot pixels in darks. If ``STDEV``, cold and hot pixels will be found " 

82 "in flats, and hot pixels in darks."), 

83 default='STDEV', 

84 allowed={'STDEV': "Use a multiple of the image standard deviation to determine detection threshold.", 

85 'VALUE': "Use pixel value to determine detection threshold."}, 

86 ) 

87 doVampirePixels = pexConfig.Field( 

88 dtype=bool, 

89 doc=("Search for vampire pixels (bright pixels surrounded by ring of low flux) in ComCam " 

90 "flatBootstrap and mask the area arount them."), 

91 default=False, 

92 ) 

93 thresholdVampirePixels = pexConfig.Field( 

94 dtype=float, 

95 doc=("Pixel value threshold to find bright pixels in ComCam flatBootstrap."), 

96 default=1.9, 

97 ) 

98 radiusVampirePixels = pexConfig.Field( 

99 dtype=int, 

100 doc=("Radius (in pixels) of the area to mask around ComCam flatBootstrap bright pixels."), 

101 default=8, 

102 ) 

103 darkCurrentThreshold = pexConfig.Field( 

104 dtype=float, 

105 doc=("If thresholdType=``VALUE``, dark current threshold (in e-/sec) to define " 

106 "hot/bright pixels in dark images. Unused if thresholdType==``STDEV``."), 

107 default=5, 

108 ) 

109 biasThreshold = pexConfig.Field( 

110 dtype=float, 

111 doc=("If thresholdType==``VALUE``, bias threshold (in ADU) to define " 

112 "hot/bright pixels in bias frame. Unused if thresholdType==``STDEV``."), 

113 default=1000.0, 

114 ) 

115 fracThresholdFlat = pexConfig.Field( 

116 dtype=float, 

117 doc=("If thresholdType=``VALUE``, fractional threshold to define cold/dark " 

118 "pixels in flat images (fraction of the mean value per amplifier)." 

119 "Unused if thresholdType==``STDEV``."), 

120 default=0.8, 

121 ) 

122 nSigmaBright = pexConfig.Field( 

123 dtype=float, 

124 doc=("If thresholdType=``STDEV``, number of sigma above mean for bright/hot " 

125 "pixel detection. The default value was found to be " 

126 "appropriate for some LSST sensors in DM-17490. " 

127 "Unused if thresholdType==``VALUE``"), 

128 default=4.8, 

129 ) 

130 nSigmaDark = pexConfig.Field( 

131 dtype=float, 

132 doc=("If thresholdType=``STDEV``, number of sigma below mean for dark/cold pixel " 

133 "detection. The default value was found to be " 

134 "appropriate for some LSST sensors in DM-17490. " 

135 "Unused if thresholdType==``VALUE``"), 

136 default=-5.0, 

137 ) 

138 nPixBorderUpDown = pexConfig.Field( 

139 dtype=int, 

140 doc="Number of pixels to exclude from top & bottom of image when looking for defects.", 

141 default=0, 

142 ) 

143 nPixBorderLeftRight = pexConfig.Field( 

144 dtype=int, 

145 doc="Number of pixels to exclude from left & right of image when looking for defects.", 

146 default=0, 

147 ) 

148 badOnAndOffPixelColumnThreshold = pexConfig.Field( 

149 dtype=int, 

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

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

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

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

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

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

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

157 default=50, 

158 ) 

159 goodPixelColumnGapThreshold = pexConfig.Field( 

160 dtype=int, 

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

162 "'badOnAndOffPixelColumnThreshold')."), 

163 default=30, 

164 ) 

165 badPixelsToFillColumnThreshold = pexConfig.Field( 

166 dtype=float, 

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

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

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

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

171 default=-1, 

172 ) 

173 saturatedColumnMask = pexConfig.Field( 

174 dtype=str, 

175 default="SAT", 

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

177 ) 

178 saturatedColumnDilationRadius = pexConfig.Field( 

179 dtype=int, 

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

181 "to mitigate glow."), 

182 default=0, 

183 ) 

184 saturatedPixelsToFillColumnThreshold = pexConfig.Field( 

185 dtype=int, 

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

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

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

189 default=-1, 

190 ) 

191 e2vMidlineBreakNRow = pexConfig.Field( 

192 dtype=int, 

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

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

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

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

197 default=1, 

198 ) 

199 e2vMidlineBreakOption = pexConfig.ChoiceField( 

200 dtype=str, 

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

202 default="NEVERMASK", 

203 allowed={ 

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

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

206 }, 

207 ) 

208 

209 def validate(self): 

210 super().validate() 

211 if self.nSigmaBright < 0.0: 

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

213 if self.nSigmaDark > 0.0: 

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

215 

216 

217class MeasureDefectsTask(pipeBase.PipelineTask): 

218 """Measure the defects from one exposure. 

219 """ 

220 

221 ConfigClass = MeasureDefectsTaskConfig 

222 _DefaultName = 'cpDefectMeasure' 

223 

224 def run(self, inputExp, camera): 

225 """Measure one exposure for defects. 

226 

227 Parameters 

228 ---------- 

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

230 Exposure to examine. 

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

232 Camera to use for metadata. 

233 

234 Returns 

235 ------- 

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

237 Results struct containing: 

238 

239 ``outputDefects`` 

240 The defects measured from this exposure 

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

242 """ 

243 detector = inputExp.getDetector() 

244 try: 

245 filterName = inputExp.getFilter().physicalLabel 

246 except AttributeError: 

247 filterName = None 

248 

249 defects = self._findHotAndColdPixels(inputExp) 

250 

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

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

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

254 

255 defects.updateMetadataFromExposures([inputExp]) 

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

257 setCalibId=True, setDate=True, 

258 cpDefectGenImageType=datasetType) 

259 

260 return pipeBase.Struct( 

261 outputDefects=defects, 

262 ) 

263 

264 @staticmethod 

265 def _nPixFromDefects(defects): 

266 """Count pixels in a defect. 

267 

268 Parameters 

269 ---------- 

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

271 Defects to measure. 

272 

273 Returns 

274 ------- 

275 nPix : `int` 

276 Number of defect pixels. 

277 """ 

278 nPix = 0 

279 for defect in defects: 

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

281 return nPix 

282 

283 def getVampirePixels(self, ampImg): 

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

285 extended area around them, 

286 

287 Parameters 

288 ---------- 

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

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

291 

292 Returns 

293 ------- 

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

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

296 """ 

297 

298 # Find bright pixels 

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

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

301 fs = afwDetection.FootprintSet(ampImg, thresh) 

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

303 

304 return fs_grow 

305 

306 def _setE2VMidline(self, amp, ampImg): 

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

308 

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

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

311 to the median value. 

312 

313 Parameters 

314 ---------- 

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

316 Amplifier object. 

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

318 Image data for the amplifier. 

319 """ 

320 if self.config.e2vMidlineBreakNRow <= 0: 

321 return 

322 

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

324 # Setting the midline values to the median value ensures 

325 # they will never trigger the defect finding threshold. 

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

327 else: 

328 # Setting the midline values to a large negative value 

329 # ensures they will always trigger the defect finding 

330 # threshold. 

331 value = -1e30 

332 

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

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

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

336 else: 

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

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

339 

340 def _findHotAndColdPixels(self, exp): 

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

342 

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

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

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

346 cold pixels). 

347 

348 Parameters 

349 ---------- 

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

351 The exposure in which to find defects. 

352 

353 Returns 

354 ------- 

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

356 The defects found in the image. 

357 """ 

358 self._setEdgeBits(exp) 

359 maskedIm = exp.maskedImage 

360 

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

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

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

364 footprintList = [] 

365 

366 hotPixelCount = {} 

367 coldPixelCount = {} 

368 

369 for amp in exp.getDetector(): 

370 ampName = amp.getName() 

371 

372 hotPixelCount[ampName] = 0 

373 coldPixelCount[ampName] = 0 

374 

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

376 

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

378 if self.config.nPixBorderLeftRight: 

379 if ampImg.getX0() == 0: 

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

381 else: 

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

383 if self.config.nPixBorderUpDown: 

384 if ampImg.getY0() == 0: 

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

386 else: 

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

388 

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

390 continue 

391 

392 if self.config.doVampirePixels: 

393 # This is only applied in LSSTComCam flatBootstrap pipeline 

394 footprintSet_VampirePixel = self.getVampirePixels(ampImg) 

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

396 

397 # Remove a background estimate 

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

399 ampImg -= meanClip 

400 

401 # Determine thresholds 

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

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

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

405 if np.isnan(expTime): 

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

407 expTime, ampName, datasetType) 

408 expTime = 1. 

409 thresholdType = self.config.thresholdType 

410 if thresholdType == 'VALUE': 

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

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

413 # We scale by the exposure time. 

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

415 # hot pixel threshold 

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

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

418 # hot pixel threshold, no exposure time. 

419 valueThreshold = self.config.biasThreshold 

420 else: 

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

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

423 # the mean (at 500nm). 

424 

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

426 # negative cold pixel threshold. 

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

428 # Find equivalent sigma values. 

429 if stDev == 0.0: 

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

431 stDev, ampName, datasetType) 

432 nSigmaList = [np.inf] 

433 else: 

434 nSigmaList = [valueThreshold/stDev] 

435 else: 

436 hotPixelThreshold = self.config.nSigmaBright 

437 coldPixelThreshold = self.config.nSigmaDark 

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

439 nSigmaList = [hotPixelThreshold] 

440 valueThreshold = stDev*hotPixelThreshold 

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

442 self.log.warning( 

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

444 ) 

445 return Defects.fromFootprintList([]) 

446 else: 

447 nSigmaList = [hotPixelThreshold, coldPixelThreshold] 

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

449 

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

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

452 datasetType, ampName, thresholdType, nSigmaList, valueThreshold) 

453 

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

455 self._setE2VMidline(amp, ampImg) 

456 

457 mergedSet = None 

458 for sigma in nSigmaList: 

459 nSig = np.abs(sigma) 

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

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

462 

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

464 

465 try: 

466 footprintSet = afwDetection.FootprintSet(ampImg, threshold) 

467 except InvalidParameterError: 

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

469 # Let's mask the whole area. 

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

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

472 footprintSet = afwDetection.FootprintSet(ampImg, threshold) 

473 

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

475 

476 if mergedSet is None: 

477 mergedSet = footprintSet 

478 else: 

479 mergedSet.merge(footprintSet) 

480 

481 if polarity: 

482 # hot pixels 

483 for fp in footprintSet.getFootprints(): 

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

485 else: 

486 # cold pixels 

487 for fp in footprintSet.getFootprints(): 

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

489 

490 if self.config.doVampirePixels: 

491 # Count the number of pixels masked 

492 vampirePixelCount = 0 

493 for fp in footprintSet_VampirePixel.getFootprints(): 

494 vampirePixelCount += fp.getArea() 

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

496 # Add vampire pixels to footprint set 

497 mergedSet.merge(footprintSet_VampirePixel) 

498 

499 footprintList += mergedSet.getFootprints() 

500 

501 self.debugView('defectMap', ampImg, 

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

503 

504 defects = Defects.fromFootprintList(footprintList) 

505 defects = self.dilateSaturatedColumns(exp, defects) 

506 defects, _ = self.maskBlocksIfIntermitentBadPixelsInColumn(defects) 

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

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

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

510 

511 return defects 

512 

513 @staticmethod 

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

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

516 nPixels = maskedIm.mask.array.size 

517 nBad = countMaskedPixels(maskedIm, badMaskString) 

518 return nPixels - nBad 

519 

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

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

522 

523 Raises 

524 ------ 

525 TypeError 

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

527 """ 

528 if isinstance(exposureOrMaskedImage, afwImage.Exposure): 

529 mi = exposureOrMaskedImage.maskedImage 

530 elif isinstance(exposureOrMaskedImage, afwImage.MaskedImage): 

531 mi = exposureOrMaskedImage 

532 else: 

533 t = type(exposureOrMaskedImage) 

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

535 

536 MASKBIT = mi.mask.getPlaneBitMask(maskplaneToSet) 

537 if self.config.nPixBorderLeftRight: 

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

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

540 if self.config.nPixBorderUpDown: 

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

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

543 

544 def maskBlocksIfIntermitentBadPixelsInColumn(self, defects): 

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

546 

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

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

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

550 

551 Parameters 

552 ---------- 

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

554 The defects found in the image so far 

555 

556 Returns 

557 ------- 

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

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

560 equal than self.config.badPixelColumnThreshold, the input 

561 list is returned. Otherwise, the defects list returned 

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

563 badColumnCount : `int` 

564 Number of bad columns partially masked. 

565 """ 

566 badColumnCount = 0 

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

568 coordinates = [] 

569 for defect in defects: 

570 bbox = defect.getBBox() 

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

572 deltaX0, deltaY0 = bbox.getDimensions() 

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

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

575 coordinates.append((i, j)) 

576 

577 x, y = [], [] 

578 for coordinatePair in coordinates: 

579 x.append(coordinatePair[0]) 

580 y.append(coordinatePair[1]) 

581 

582 x = np.array(x) 

583 y = np.array(y) 

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

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

586 multipleX = [] 

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

588 if b >= self.config.badOnAndOffPixelColumnThreshold: 

589 multipleX.append(a) 

590 if len(multipleX) != 0: 

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

592 badColumnCount += 1 

593 

594 return defects, badColumnCount 

595 

596 def dilateSaturatedColumns(self, exp, defects): 

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

598 

599 Parameters 

600 ---------- 

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

602 The exposure in which to find defects. 

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

604 The defects found in the image so far 

605 

606 Returns 

607 ------- 

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

609 The expanded defects. 

610 """ 

611 if self.config.saturatedColumnDilationRadius <= 0: 

612 # This is a no-op. 

613 return defects 

614 

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

616 

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

618 

619 if len(satX) == 0: 

620 # No saturated pixels, nothing to do. 

621 return defects 

622 

623 radius = self.config.saturatedColumnDilationRadius 

624 

625 with defects.bulk_update(): 

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

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

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

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

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

631 defects.append(s) 

632 

633 return defects 

634 

635 def maskBadColumns(self, exp, defects): 

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

637 

638 Parameters 

639 ---------- 

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

641 The defects found in the image so far 

642 

643 Returns 

644 ------- 

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

646 The exposure in which to find defects. 

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

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

649 equal than self.config.badPixelColumnThreshold, the input 

650 list is returned. Otherwise, the defects list returned 

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

652 badColumnCount : `int` 

653 Number of bad columns masked. 

654 """ 

655 # Render the defects into an image. 

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

657 

658 for defect in defects: 

659 defectImage[defect.getBBox()] = 1 

660 

661 badColumnCount = 0 

662 

663 if self.config.badPixelsToFillColumnThreshold > 0: 

664 with defects.bulk_update(): 

665 for amp in exp.getDetector(): 

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

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

668 

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

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

671 

672 for badColumn in badColumns: 

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

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

675 defects.append(s) 

676 

677 badColumnCount += len(badColIndices) 

678 

679 if self.config.saturatedPixelsToFillColumnThreshold > 0: 

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

681 

682 with defects.bulk_update(): 

683 for amp in exp.getDetector(): 

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

685 # Turn all the SAT bits into 1s 

686 subMask &= mask 

687 subMask[subMask > 0] = 1 

688 

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

690 

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

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

693 

694 for badColumn in badColumns: 

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

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

697 defects.append(s) 

698 

699 badColumnCount += len(badColIndices) 

700 

701 return defects, badColumnCount 

702 

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

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

705 threshold. 

706 

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

708 in a column is larger or equal than 

709 self.config.badOnAndOffPixelColumnThreshold. 

710 

711 Parameters 

712 --------- 

713 x : `list` 

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

715 along the short axis if amp. 

716 y : `list` 

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

718 along the long axis if amp. 

719 multipleX : list 

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

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

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

723 The defcts found in the image so far 

724 

725 Returns 

726 ------- 

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

728 The defects list returned that will include boxes that 

729 mask blocks of on-and-of pixels. 

730 """ 

731 with defects.bulk_update(): 

732 goodPixelColumnGapThreshold = self.config.goodPixelColumnGapThreshold 

733 for x0 in multipleX: 

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

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

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

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

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

739 # of good pixels between two consecutive bad pixels is 

740 # larger or equal than 'goodPixelColumnGapThreshold'. 

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

742 if len(diffIndex) != 0: 

743 limits = [minY] # put the minimum first 

744 for gapIndex in diffIndex: 

745 limits.append(multipleY[gapIndex]) 

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

747 limits.append(maxY) # maximum last 

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

749 s = Box2I(minimum=Point2I(x0, limits[i]), maximum=Point2I(x0, limits[i+1])) 

750 defects.append(s) 

751 else: # No gap is large enough 

752 s = Box2I(minimum=Point2I(x0, minY), maximum=Point2I(x0, maxY)) 

753 defects.append(s) 

754 return defects 

755 

756 def debugView(self, stepname, ampImage, defects, detector): # pragma: no cover 

757 """Plot the defects found by the task. 

758 

759 Parameters 

760 ---------- 

761 stepname : `str` 

762 Debug frame to request. 

763 ampImage : `lsst.afw.image.MaskedImage` 

764 Amplifier image to display. 

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

766 The defects to plot. 

767 detector : `lsst.afw.cameraGeom.Detector` 

768 Detector holding camera geometry. 

769 """ 

770 frame = getDebugFrame(self._display, stepname) 

771 if frame: 

772 disp = afwDisplay.Display(frame=frame) 

773 disp.scale('asinh', 'zscale') 

774 disp.setMaskTransparency(80) 

775 disp.setMaskPlaneColor("BAD", afwDisplay.RED) 

776 

777 maskedIm = ampImage.clone() 

778 defects.maskPixels(maskedIm, "BAD") 

779 

780 mpDict = maskedIm.mask.getMaskPlaneDict() 

781 for plane in mpDict.keys(): 

782 if plane in ['BAD']: 

783 continue 

784 disp.setMaskPlaneColor(plane, afwDisplay.IGNORE) 

785 

786 disp.setImageColormap('gray') 

787 disp.mtv(maskedIm) 

788 cameraGeom.utils.overlayCcdBoxes(detector, isTrimmed=True, display=disp) 

789 prompt = "Press Enter to continue [c]... " 

790 while True: 

791 ans = input(prompt).lower() 

792 if ans in ('', 'c', ): 

793 break 

794 

795 def debugHistogram(self, stepname, ampImage, nSigmaUsed, exp): 

796 """Make a histogram of the distribution of pixel values for 

797 each amp. 

798 

799 The main image data histogram is plotted in blue. Edge 

800 pixels, if masked, are in red. Note that masked edge pixels 

801 do not contribute to the underflow and overflow numbers. 

802 

803 Note that this currently only supports the 16-amp LSST 

804 detectors. 

805 

806 Parameters 

807 ---------- 

808 stepname : `str` 

809 Debug frame to request. 

810 ampImage : `lsst.afw.image.MaskedImage` 

811 Amplifier image to display. 

812 nSigmaUsed : `float` 

813 The number of sigma used for detection 

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

815 The exposure in which the defects were found. 

816 """ 

817 frame = getDebugFrame(self._display, stepname) 

818 if frame: 

819 import matplotlib.pyplot as plt 

820 

821 detector = exp.getDetector() 

822 nX = np.floor(np.sqrt(len(detector))) 

823 nY = len(detector) // nX 

824 fig, ax = plt.subplots(nrows=int(nY), ncols=int(nX), sharex='col', sharey='row', figsize=(13, 10)) 

825 

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

827 

828 for (amp, a) in zip(reversed(detector), ax.flatten()): 

829 mi = exp.maskedImage[amp.getBBox()] 

830 

831 # normalize by expTime as we plot in ADU/s and don't 

832 # always work with master calibs 

833 mi.image.array /= expTime 

834 stats = afwMath.makeStatistics(mi, afwMath.MEANCLIP | afwMath.STDEVCLIP) 

835 mean, sigma = stats.getValue(afwMath.MEANCLIP), stats.getValue(afwMath.STDEVCLIP) 

836 # Get array of pixels 

837 EDGEBIT = exp.maskedImage.mask.getPlaneBitMask("EDGE") 

838 imgData = mi.image.array[(mi.mask.array & EDGEBIT) == 0].flatten() 

839 edgeData = mi.image.array[(mi.mask.array & EDGEBIT) != 0].flatten() 

840 

841 thrUpper = mean + nSigmaUsed*sigma 

842 thrLower = mean - nSigmaUsed*sigma 

843 

844 nRight = len(imgData[imgData > thrUpper]) 

845 nLeft = len(imgData[imgData < thrLower]) 

846 

847 nsig = nSigmaUsed + 1.2 # add something small so the edge of the plot is out from level used 

848 leftEdge = mean - nsig * nSigmaUsed*sigma