Coverage for python/lsst/meas/algorithms/dynamicDetection.py: 13%

2__all__ = [

3 "DynamicDetectionConfig",

4 "DynamicDetectionTask",

5 "InsufficientSourcesError",

8import numpy as np

10from lsst.pex.config import Field, ConfigurableField

12from .detection import SourceDetectionConfig, SourceDetectionTask

13from .skyObjects import SkyObjectsTask

15from lsst.afw.detection import FootprintSet

16from lsst.afw.geom import makeCdMatrix, makeSkyWcs, SpanSet

17from lsst.afw.table import SourceCatalog, SourceTable

18from lsst.meas.base import ForcedMeasurementTask

19from lsst.pipe.base import Struct

21import lsst.afw.image

22import lsst.afw.math

23import lsst.geom as geom

26class InsufficientSourcesError(Exception):

27 """Raised if an insufficient number of sky sources are found for

28 dynamic detection.

30 Parameters

31 ----------

32 msg : `str`

33 Error message.

34 nGoodPix : `int`

35 Number of good pixels (i.e. not NO_DATA or BAD).

36 nPix : `int`

37 Total number of pixels.

38 **kwargs : `dict`, optional

39 Additional keyword arguments to initialize the Exception base class.

40 """

41 def __init__(self, msg, nGoodPix, nPix, **kwargs):

42 self.msg = msg

43 self._metadata = kwargs

44 super().__init__(msg, **kwargs)

45 self._metadata["nGoodPix"] = int(nGoodPix)

46 self._metadata["nPix"] = int(nPix)

48 def __str__(self):

49 # Exception doesn't handle **kwargs, so we need a custom str.

50 return f"{self.msg}: {self.metadata}"

52 @property

53 def metadata(self):

54 for key, value in self._metadata.items():

55 if not isinstance(value, (int, float, str)):

56 raise TypeError(f"{key} is of type {type(value)}, but only (int, float, str) are allowed.")

57 return self._metadata

60class DynamicDetectionConfig(SourceDetectionConfig):

61 """Configuration for DynamicDetectionTask

62 """

63 prelimThresholdFactor = Field(dtype=float, default=0.5,

64 doc="Factor by which to multiply the main detection threshold "

65 "(thresholdValue) to use for first pass (to find sky objects).")

66 prelimNegMultiplier = Field(dtype=float, default=2.5,

67 doc="Multiplier for the negative (relative to positive) polarity "

68 "detections threshold to use for first pass (to find sky objects).")

69 skyObjects = ConfigurableField(target=SkyObjectsTask, doc="Generate sky objects.")

70 doBackgroundTweak = Field(dtype=bool, default=True,

71 doc="Tweak background level so median PSF flux of sky objects is zero?")

72 minFractionSources = Field(dtype=float, default=0.02,

73 doc="Minimum fraction of the requested number of sky sources for dynamic "

74 "detection to be considered a success. If the number of good sky sources "

75 "identified falls below this threshold, a NoWorkFound error is raised so "

76 "that this dataId is no longer considered in downstream processing.")

77 doBrightPrelimDetection = Field(dtype=bool, default=True,

78 doc="Do initial bright detection pass where footprints are grown "

79 "by brightGrowFactor?")

80 brightMultiplier = Field(dtype=float, default=2000.0,

81 doc="Multiplier to apply to the prelimThresholdFactor for the "

82 "\"bright\" detections stage (want this to be large to only "

83 "detect the brightest sources).")

84 brightNegFactor = Field(dtype=float, default=2.2,

85 doc="Factor by which to multiply the threshold for the negative polatiry "

86 "detections for the \"bright\" detections stage (this needs to be fairly "

87 "low given the nature of the negative polarity detections in the very "

88 "large positive polarity threshold).")

89 brightGrowFactor = Field(dtype=int, default=40,

90 doc="Factor by which to grow the footprints of sources detected in the "

91 "\"bright\" detections stage (want this to be large to mask wings of "

92 "bright sources).")

93 brightMaskFractionMax = Field(dtype=float, default=0.95,

94 doc="Maximum allowed fraction of masked pixes from the \"bright\" "

95 "detection stage (to mask regions unsuitable for sky sourcess). "

96 "If this fraction is exeeded, the detection threshold for this stage "

97 "will be increased by bisectFactor until the fraction of masked "

98 "pixels drops below this threshold.")

99 bisectFactor = Field(dtype=float, default=1.2,

100 doc="Factor by which to increase thresholds in brightMaskFractionMax loop.")

101 allowMaskErode = Field(dtype=bool, default=True,

102 doc="Crowded/large fill-factor fields make it difficult to find "

103 "suitable locations to lay down sky objects. To allow for best effort "

104 "sky source placement, if True, this allows for a slight erosion of "

105 "the detection masks.")

106 maxPeakToFootRatio = Field(dtype=float, default=150.0,

107 doc="Maximum ratio of peak per footprint in the detection mask. "

108 "This is to help prevent single contiguous footprints that nothing "

109 "can be done with (i.e. deblending will be skipped). If the current "

110 "detection plane does not satisfy this constraint, the detection "

111 "threshold is increased iteratively until it is. This behaviour is "

112 "intended to be an effective no-op for most \"typical\" scenes/standard "

113 "quality observations, but can avoid total meltdown in, e.g. very "

114 "crowded regions.")

115

116 def setDefaults(self):

117 SourceDetectionConfig.setDefaults(self)

118 self.skyObjects.nSources = 1000 # For good statistics

119 for maskStr in ["SAT"]:

120 if maskStr not in self.skyObjects.avoidMask:

121 self.skyObjects.avoidMask.append(maskStr)

122

123 def validate(self):

124 super().validate()

125

126 if self.doApplyFlatBackgroundRatio:

127 raise ValueError("DynamicDetectionTask does not support doApplyFlatBackgroundRatio.")

128

129

130class DynamicDetectionTask(SourceDetectionTask):

131 """Detection of sources on an image with a dynamic threshold

132

133 We first detect sources using a lower threshold than normal (see config

134 parameter ``prelimThresholdFactor``) in order to identify good sky regions

135 (configurable ``skyObjects``). Then we perform forced PSF photometry on

136 those sky regions. Using those PSF flux measurements and estimated errors,

137 we set the threshold so that the stdev of the measurements matches the

138 median estimated error.

139

140 Besides the usual initialisation of configurables, we also set up

141 the forced measurement which is deliberately not represented in

142 this Task's configuration parameters because we're using it as

143 part of the algorithm and we don't want to allow it to be modified.

144 """

145 ConfigClass = DynamicDetectionConfig

146 _DefaultName = "dynamicDetection"

147

148 def __init__(self, *args, **kwargs):

149

150 SourceDetectionTask.__init__(self, *args, **kwargs)

151 self.makeSubtask("skyObjects")

152

153 # Set up forced measurement.

154 config = ForcedMeasurementTask.ConfigClass()

155 config.plugins.names = ["base_TransformedCentroid", "base_PsfFlux"]

156 # We'll need the "centroid" and "psfFlux" slots

157 for slot in ("shape", "psfShape", "apFlux", "modelFlux", "gaussianFlux", "calibFlux"):

158 setattr(config.slots, slot, None)

159 config.copyColumns = {}

160 self.skySchema = SourceTable.makeMinimalSchema()

161 self.skyMeasurement = ForcedMeasurementTask(config=config, name="skyMeasurement", parentTask=self,

162 refSchema=self.skySchema)

163

164 def calculateThreshold(self, exposure, seed, sigma=None, minFractionSourcesFactor=1.0,

165 isBgTweak=False, nPixMaskErode=None, maxMaskErodeIter=10):

166 """Calculate new threshold

167

168 This is the main functional addition to the vanilla

169 `SourceDetectionTask`.

170

171 We identify sky objects and perform forced PSF photometry on

172 them. Using those PSF flux measurements and estimated errors,

173 we set the threshold so that the stdev of the measurements

174 matches the median estimated error.

175

176 Parameters

177 ----------

178 exposure : `lsst.afw.image.Exposure`

179 Exposure on which we're detecting sources.

180 seed : `int`

181 RNG seed to use for finding sky objects.

182 sigma : `float`, optional

183 Gaussian sigma of smoothing kernel; if not provided,

184 will be deduced from the exposure's PSF.

185 minFractionSourcesFactor : `float`

186 Change the fraction of required sky sources from that set in

187 ``self.config.minFractionSources`` by this factor. NOTE: this

188 is intended for use in the background tweak pass (the detection

189 threshold is much lower there, so many more pixels end up marked

190 as DETECTED or DETECTED_NEGATIVE, leaving less room for sky

191 object placement).

192 isBgTweak : `bool`

193 Set to ``True`` for the background tweak pass (for more helpful

194 log messages).

195 nPixMaskErode : `int`, optional

196 Number of pixels by which to erode the detection masks on each

197 iteration of best-effort sky object placement.

198 maxMaskErodeIter : `int`, optional

199 Maximum number of iterations for the detection mask erosion.

200

201 Returns

202 -------

203 result : `lsst.pipe.base.Struct`

204 Result struct with components:

205

206 ``multiplicative``

207 Multiplicative factor to be applied to the

208 configured detection threshold (`float`).

209 ``additive``

210 Additive factor to be applied to the background

211 level (`float`).

212

213 Raises

214 ------

215 NoWorkFound

216 Raised if the number of good sky sources found is less than the

217 minimum fraction

218 (``self.config.minFractionSources``*``minFractionSourcesFactor``)

219 of the number requested (``self.skyObjects.config.nSources``).

220 """

221 wcsIsNone = exposure.getWcs() is None

222 if wcsIsNone: # create a dummy WCS as needed by ForcedMeasurementTask

223 self.log.info("WCS for exposure is None. Setting a dummy WCS for dynamic detection.")

224 exposure.setWcs(makeSkyWcs(crpix=geom.Point2D(0, 0),

225 crval=geom.SpherePoint(0, 0, geom.degrees),

226 cdMatrix=makeCdMatrix(scale=1e-5*geom.degrees)))

227 minNumSources = int(self.config.minFractionSources*self.skyObjects.config.nSources)

228 # Reduce the number of sky sources required if requested, but ensure

229 # a minumum of 3.

230 if minFractionSourcesFactor != 1.0:

231 minNumSources = max(3, int(minNumSources*minFractionSourcesFactor))

232 fp = self.skyObjects.run(exposure.maskedImage.mask, seed)

233

234 if self.config.allowMaskErode:

235 detectedMaskPlanes = ["DETECTED", "DETECTED_NEGATIVE"]

236 mask = exposure.maskedImage.mask

237 for nIter in range(maxMaskErodeIter):

238 if nIter > 0:

239 fp = self.skyObjects.run(mask, seed)

240 if len(fp) < int(2*minNumSources): # Allow for measurement failures

241 self.log.info("Current number of sky sources is below 2*minimum required "

242 "(%d < %d, allowing for some subsequent measurement failures). "

243 "Allowing erosion of detected mask planes for sky placement "

244 "nIter: %d [of %d max]",

245 len(fp), 2*minNumSources, nIter, maxMaskErodeIter)

246 if nPixMaskErode is None:

247 if len(fp) == 0:

248 nPixMaskErode = 4

249 elif len(fp) < int(0.75*minNumSources):

250 nPixMaskErode = 2

251 else:

252 nPixMaskErode = 1

253 for maskName in detectedMaskPlanes:

254 # Compute the eroded detection mask plane using SpanSet

255 detectedMaskBit = mask.getPlaneBitMask(maskName)

256 detectedMaskSpanSet = SpanSet.fromMask(mask, detectedMaskBit)

257 detectedMaskSpanSet = detectedMaskSpanSet.eroded(nPixMaskErode)

258 # Clear the detected mask plane

259 detectedMask = mask.getMaskPlane(maskName)

260 mask.clearMaskPlane(detectedMask)

261 # Set the mask plane to the eroded one

262 detectedMaskSpanSet.setMask(mask, detectedMaskBit)

263 else:

264 break

265

266 skyFootprints = FootprintSet(exposure.getBBox())

267 skyFootprints.setFootprints(fp)

268 table = SourceTable.make(self.skyMeasurement.schema)

269 catalog = SourceCatalog(table)

270 catalog.reserve(len(skyFootprints.getFootprints()))

271 skyFootprints.makeSources(catalog)

272 key = catalog.getCentroidSlot().getMeasKey()

273 for source in catalog:

274 peaks = source.getFootprint().getPeaks()

275 assert len(peaks) == 1

276 source.set(key, peaks[0].getF())

277 # Coordinate covariance is not used, so don't bother calulating it.

278 source.updateCoord(exposure.getWcs(), include_covariance=False)

279

280 # Forced photometry on sky objects

281 self.skyMeasurement.run(catalog, exposure, catalog, exposure.getWcs())

282

283 # Calculate new threshold

284 fluxes = catalog["base_PsfFlux_instFlux"]

285 area = catalog["base_PsfFlux_area"]

286 good = (~catalog["base_PsfFlux_flag"] & np.isfinite(fluxes))

287

288 if good.sum() < minNumSources:

289 if not isBgTweak:

290 msg = (f"Insufficient good sky source flux measurements ({good.sum()} < "

291 f"{minNumSources}) for dynamic threshold calculation.")

292 else:

293 msg = (f"Insufficient good sky source flux measurements ({good.sum()} < "

294 f"{minNumSources}) for background tweak calculation.")

295

296 nPix = exposure.mask.array.size

297 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])

298 nGoodPix = np.sum(exposure.mask.array & badPixelMask == 0)

299 if nGoodPix/nPix > 0.2:

300 detectedPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["DETECTED", "DETECTED_NEGATIVE"])

301 nDetectedPix = np.sum(exposure.mask.array & detectedPixelMask != 0)

302 msg += (f" However, {nGoodPix}/{nPix} pixels are not marked NO_DATA or BAD, "

303 "so there should be sufficient area to locate suitable sky sources. "

304 f"Note that {nDetectedPix} of {nGoodPix} \"good\" pixels were marked "

305 "as DETECTED or DETECTED_NEGATIVE.")

306 raise InsufficientSourcesError(msg, nGoodPix, nPix)

307 raise InsufficientSourcesError(msg, nGoodPix, nPix)

308

309 if not isBgTweak:

310 self.log.info("Number of good sky sources used for dynamic detection: %d (of %d requested).",

311 good.sum(), self.skyObjects.config.nSources)

312 else:

313 self.log.info("Number of good sky sources used for dynamic detection background tweak:"

314 " %d (of %d requested).", good.sum(), self.skyObjects.config.nSources)

315

316 bgMedian = np.median((fluxes/area)[good])

317 lq, uq = np.percentile(fluxes[good], [25.0, 75.0])

318 stdevMeas = 0.741*(uq - lq)

319 medianError = np.median(catalog["base_PsfFlux_instFluxErr"][good])

320 if wcsIsNone:

321 exposure.setWcs(None)

322 return Struct(multiplicative=medianError/stdevMeas, additive=bgMedian)

323

324 def detectFootprints(self, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None,

325 background=None, backgroundToPhotometricRatio=None):

326 """Detect footprints with a dynamic threshold

327

328 This varies from the vanilla ``detectFootprints`` method because we

329 do detection three times: first with a high threshold to detect

330 "bright" (both positive and negative, the latter to identify very

331 over-subtracted regions) sources for which we grow the DETECTED and

332 DETECTED_NEGATIVE masks significantly to account for wings. Second,

333 with a low threshold to mask all non-empty regions of the image. These

334 two masks are combined and used to identify regions of sky

335 uncontaminated by objects. A final round of detection is then done

336 with the new calculated threshold.

337

338 Parameters

339 ----------

340 exposure : `lsst.afw.image.Exposure`

341 Exposure to process; DETECTED{,_NEGATIVE} mask plane will be

342 set in-place.

343 doSmooth : `bool`, optional

344 If True, smooth the image before detection using a Gaussian

345 of width ``sigma``.

346 sigma : `float`, optional

347 Gaussian Sigma of PSF (pixels); used for smoothing and to grow

348 detections; if `None` then measure the sigma of the PSF of the

349 ``exposure``.

350 clearMask : `bool`, optional

351 Clear both DETECTED and DETECTED_NEGATIVE planes before running

352 detection.

353 expId : `int`, optional

354 Exposure identifier, used as a seed for the random number

355 generator. If absent, the seed will be the sum of the image.

356 background : `lsst.afw.math.BackgroundList`, optional

357 Background that was already subtracted from the exposure; will be

358 modified in-place if ``reEstimateBackground=True``.

359 backgroundToPhotometricRatio : `lsst.afw.image.Image`, optional

360 Unused; if set will Raise.

361

362 Returns

363 -------

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

365 The results `~lsst.pipe.base.Struct` contains:

366

367 ``positive``

368 Positive polarity footprints.

369 (`lsst.afw.detection.FootprintSet` or `None`)

370 ``negative``

371 Negative polarity footprints.

372 (`lsst.afw.detection.FootprintSet` or `None`)

373 ``numPos``

374 Number of footprints in positive or 0 if detection polarity was

375 negative. (`int`)

376 ``numNeg``

377 Number of footprints in negative or 0 if detection polarity was

378 positive. (`int`)

379 ``background``

380 Re-estimated background. `None` or the input ``background``

381 if ``reEstimateBackground==False``.

382 (`lsst.afw.math.BackgroundList`)

383 ``factor``

384 Multiplication factor applied to the configured detection

385 threshold. (`float`)

386 ``prelim``

387 Results from preliminary detection pass.

388 (`lsst.pipe.base.Struct`)

389 """

390 if backgroundToPhotometricRatio is not None:

391 raise RuntimeError("DynamicDetectionTask does not support backgroundToPhotometricRatio.")

392 maskedImage = exposure.maskedImage

393

394 if clearMask:

395 self.clearMask(maskedImage.mask)

396 else:

397 oldDetected = maskedImage.mask.array & maskedImage.mask.getPlaneBitMask(["DETECTED",

398 "DETECTED_NEGATIVE"])

399 nPix = maskedImage.mask.array.size

400 badPixelMask = lsst.afw.image.Mask.getPlaneBitMask(["NO_DATA", "BAD"])

401 nGoodPix = np.sum(maskedImage.mask.array & badPixelMask == 0)

402 self.log.info("Number of good data pixels (i.e. not NO_DATA or BAD): {} ({:.1f}% of total)".

403 format(nGoodPix, 100*nGoodPix/nPix))

404

405 with self.tempWideBackgroundContext(exposure):

406 # Could potentially smooth with a wider kernel than the PSF in

407 # order to better pick up the wings of stars and galaxies, but for

408 # now sticking with the PSF as that's more simple.

409 psf = self.getPsf(exposure, sigma=sigma)

410 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)

411

412 if self.config.doBrightPrelimDetection:

413 brightDetectedMask = self._computeBrightDetectionMask(maskedImage, convolveResults)

414

415 middle = convolveResults.middle

416 sigma = convolveResults.sigma

417 prelim = self.applyThreshold(

418 middle, maskedImage.getBBox(), factor=self.config.prelimThresholdFactor,

419 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor

420 )

421 self.finalizeFootprints(

422 maskedImage.mask, prelim, sigma, factor=self.config.prelimThresholdFactor,

423 factorNeg=self.config.prelimNegMultiplier*self.config.prelimThresholdFactor

424 )

425 if self.config.doBrightPrelimDetection:

426 # Combine prelim and bright detection masks for multiplier

427 # determination.

428 maskedImage.mask.array |= brightDetectedMask

429

430 # Calculate the proper threshold

431 # seed needs to fit in a C++ 'int' so pybind doesn't choke on it

432 seed = (expId if expId is not None else int(maskedImage.image.array.sum())) % (2**31 - 1)

433 threshResults = self.calculateThreshold(exposure, seed, sigma=sigma)

434 minMultiplicative = 0.5

435 if threshResults.multiplicative < minMultiplicative:

436 self.log.warning("threshResults.multiplicative = %.2f is less than minimum value (%.2f). "

437 "Setting to %.2f.", threshResults.multiplicative, minMultiplicative,

438 minMultiplicative)

439 factor = max(minMultiplicative, threshResults.multiplicative)

440 self.log.info("Modifying configured detection threshold by factor %.2f to %.2f",

441 factor, factor*self.config.thresholdValue)

442

443 growOverride = None

444 inFinalize = True

445 while inFinalize:

446 inFinalize = False

447 # Blow away preliminary (low threshold) detection mask

448 self.clearMask(maskedImage.mask)

449 if not clearMask:

450 maskedImage.mask.array |= oldDetected

451

452 # Rinse and repeat thresholding with new calculated threshold

453 results = self.applyThreshold(middle, maskedImage.getBBox(), factor)

454 results.prelim = prelim

455 results.background = background if background is not None else lsst.afw.math.BackgroundList()

456 if self.config.doTempLocalBackground:

457 self.applyTempLocalBackground(exposure, middle, results)

458 self.finalizeFootprints(maskedImage.mask, results, sigma, factor=factor,

459 growOverride=growOverride)

460 self.log.warning("nPeaks/nFootprint = %.2f (max is %.1f)",

461 results.numPosPeaks/results.numPos,

462 self.config.maxPeakToFootRatio)

463 if results.numPosPeaks/results.numPos > self.config.maxPeakToFootRatio:

464 if results.numPosPeaks/results.numPos > 3*self.config.maxPeakToFootRatio:

465 factor *= 1.4

466 else:

467 factor *= 1.2

468 if factor > 2.0:

469 if growOverride is None:

470 growOverride = 0.75*self.config.nSigmaToGrow

471 else:

472 growOverride *= 0.75

473 self.log.warning("Decreasing nSigmaToGrow to %.2f", growOverride)

474 if factor >= 5:

475 self.log.warning("New theshold value would be > 5 times the initially requested "

476 "one (%.2f > %.2f). Leaving inFinalize iteration without "

477 "getting the number of peaks per footprint below %.1f",

478 factor*self.config.thresholdValue, self.config.thresholdValue,

479 self.config.maxPeakToFootRatio)

480 inFinalize = False

481 else:

482 inFinalize = True

483 self.log.warning("numPosPeaks/numPos (%d) > maxPeakPerFootprint (%.1f). "

484 "Increasing threshold factor to %.2f and re-running,",

485 results.numPosPeaks/results.numPos, self.config.maxPeakToFootRatio,

486 factor)

487

488 self.clearUnwantedResults(maskedImage.mask, results)

489

490 if self.config.reEstimateBackground:

491 self.reEstimateBackground(maskedImage, results.background)

492

493 self.display(exposure, results, middle)

494

495 # Re-do the background tweak after any temporary backgrounds have

496 # been restored.

497 #

498 # But we want to keep any large-scale background (e.g., scattered

499 # light from bright stars) from being selected for sky objects in

500 # the calculation, so do another detection pass without either the

501 # local or wide temporary background subtraction; the DETECTED

502 # pixels will mark the area to ignore.

503

504 # The following if/else is to workaround the fact that it is

505 # currently not possible to persist an empty BackgroundList, so

506 # we instead set the value of the backround tweak to 0.0 if

507 # doBackgroundTweak is False and call the tweakBackground function

508 # regardless to get at least one background into the list (do we

509 # need a TODO here?).

510 if self.config.doBackgroundTweak:

511 originalMask = maskedImage.mask.array.copy()

512 try:

513 self.clearMask(exposure.mask)

514 convolveResults = self.convolveImage(maskedImage, psf, doSmooth=doSmooth)

515 tweakDetResults = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(), factor)

516 self.finalizeFootprints(maskedImage.mask, tweakDetResults, sigma, factor=factor)

517 bgLevel = self.calculateThreshold(exposure, seed, sigma=sigma, minFractionSourcesFactor=0.5,

518 isBgTweak=True).additive

519 finally:

520 maskedImage.mask.array[:] = originalMask

521 else:

522 bgLevel = 0.0

523 self.tweakBackground(exposure, bgLevel, results.background)

524

525 return results

526

527 def tweakBackground(self, exposure, bgLevel, bgList=None):

528 """Modify the background by a constant value

529

530 Parameters

531 ----------

532 exposure : `lsst.afw.image.Exposure`

533 Exposure for which to tweak background.

534 bgLevel : `float`

535 Background level to remove

536 bgList : `lsst.afw.math.BackgroundList`, optional

537 List of backgrounds to append to.

538

539 Returns

540 -------

541 bg : `lsst.afw.math.BackgroundMI`

542 Constant background model.

543 """

544 if bgLevel != 0.0:

545 self.log.info("Tweaking background by %f to match sky photometry", bgLevel)

546 exposure.image -= bgLevel

547 bgStats = lsst.afw.image.MaskedImageF(1, 1)

548 bgStats.set(bgLevel, 0, bgLevel)

549 bg = lsst.afw.math.BackgroundMI(exposure.getBBox(), bgStats)

550 bgData = (bg, lsst.afw.math.Interpolate.LINEAR, lsst.afw.math.REDUCE_INTERP_ORDER,

551 lsst.afw.math.ApproximateControl.UNKNOWN, 0, 0, False)

552 if bgList is not None:

553 bgList.append(bgData)

554 return bg

555

556 def _computeBrightDetectionMask(self, maskedImage, convolveResults):

557 """Perform an initial bright source detection pass.

558

559 Perform an initial bright object detection pass using a high detection

560 threshold. The footprints in this pass are grown significantly more

561 than is typical to account for wings around bright sources. The

562 negative polarity detections in this pass help in masking severely

563 over-subtracted regions.

564

565 A maximum fraction of masked pixel from this pass is ensured via

566 the config ``brightMaskFractionMax``. If the masked pixel fraction is

567 above this value, the detection thresholds here are increased by

568 ``bisectFactor`` in a while loop until the detected masked fraction

569 falls below this value.

570

571 Parameters

572 ----------

573 maskedImage : `lsst.afw.image.MaskedImage`

574 Masked image on which to run the detection.

575 convolveResults : `lsst.pipe.base.Struct`

576 The results of the self.convolveImage function with attributes:

577

578 ``middle``

579 Convolved image, without the edges

580 (`lsst.afw.image.MaskedImage`).

581 ``sigma``

582 Gaussian sigma used for the convolution (`float`).

583

584 Returns

585 -------

586 brightDetectedMask : `numpy.ndarray`

587 Boolean array representing the union of the bright detection pass

588 DETECTED and DETECTED_NEGATIVE masks.

589 """

590 # Initialize some parameters.

591 brightPosFactor = (

592 self.config.prelimThresholdFactor*self.config.brightMultiplier/self.config.bisectFactor

593 )

594 brightNegFactor = self.config.brightNegFactor/self.config.bisectFactor

595 nPix = 1

596 nPixDet = 1

597 nPixDetNeg = 1

598 brightMaskFractionMax = self.config.brightMaskFractionMax

599

600 # Loop until masked fraction is smaller than

601 # brightMaskFractionMax, increasing the thresholds by

602 # config.bisectFactor on each iteration (rarely necessary

603 # for current defaults).

604 while nPixDetNeg/nPix > brightMaskFractionMax or nPixDet/nPix > brightMaskFractionMax:

605 self.clearMask(maskedImage.mask)

606 brightPosFactor *= self.config.bisectFactor

607 brightNegFactor *= self.config.bisectFactor

608 prelimBright = self.applyThreshold(convolveResults.middle, maskedImage.getBBox(),

609 factor=brightPosFactor, factorNeg=brightNegFactor)

610 self.finalizeFootprints(

611 maskedImage.mask, prelimBright, convolveResults.sigma*self.config.brightGrowFactor,

612 factor=brightPosFactor, factorNeg=brightNegFactor

613 )

614 # Check that not too many pixels got masked.

615 nPix = maskedImage.mask.array.size

616 nPixDet = countMaskedPixels(maskedImage, "DETECTED")

617 self.log.info("Number (%) of bright DETECTED pix: {} ({:.1f}%)".

618 format(nPixDet, 100*nPixDet/nPix))

619 nPixDetNeg = countMaskedPixels(maskedImage, "DETECTED_NEGATIVE")

620 self.log.info("Number (%) of bright DETECTED_NEGATIVE pix: {} ({:.1f}%)".

621 format(nPixDetNeg, 100*nPixDetNeg/nPix))

622 if nPixDetNeg/nPix > brightMaskFractionMax or nPixDet/nPix > brightMaskFractionMax:

623 self.log.warn("Too high a fraction (%.1f > %.1f) of pixels were masked with current "

624 "\"bright\" detection round thresholds. Increasing by a factor of %.2f "

625 "and trying again.", max(nPixDetNeg, nPixDet)/nPix,

626 brightMaskFractionMax, self.config.bisectFactor)

627

628 # Save the mask planes from the "bright" detection round, then

629 # clear them before moving on to the "prelim" detection phase.

630 brightDetectedMask = (maskedImage.mask.array

631 & maskedImage.mask.getPlaneBitMask(["DETECTED", "DETECTED_NEGATIVE"]))

632 self.clearMask(maskedImage.mask)

633 return brightDetectedMask

634

635

636def countMaskedPixels(maskedIm, maskPlane):

637 """Count the number of pixels in a given mask plane.

638

639 Parameters

640 ----------

641 maskedIm : `lsst.afw.image.MaskedImage`

642 Masked image to examine.

643 maskPlane : `str`

644 Name of the mask plane to examine.

645

646 Returns

647 -------

648 nPixMasked : `int`

649 Number of pixels with ``maskPlane`` bit set.

650 """

651 maskBit = maskedIm.mask.getPlaneBitMask(maskPlane)

652 nPixMasked = np.sum(np.bitwise_and(maskedIm.mask.array, maskBit))/maskBit

653 return nPixMasked