Coverage for python/lsst/meas/algorithms/detection.py: 15%

1# 

2# LSST Data Management System 

3# 

4# Copyright 2008-2017 AURA/LSST. 

5# 

6# This product includes software developed by the 

7# LSST Project (http://www.lsst.org/). 

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 LSST License Statement and 

20# the GNU General Public License along with this program. If not, 

21# see <https://www.lsstcorp.org/LegalNotices/>. 

22# 

23 

24__all__ = ("PsfGenerationError", "SourceDetectionConfig", "SourceDetectionTask", "addExposures") 

25 

26from contextlib import contextmanager 

27 

28import numpy as np 

29 

30import lsst.geom 

31import lsst.afw.display as afwDisplay 

32import lsst.afw.detection as afwDet 

33import lsst.afw.geom as afwGeom 

34import lsst.afw.image as afwImage 

35import lsst.afw.math as afwMath 

36import lsst.afw.table as afwTable 

37import lsst.pex.config as pexConfig 

38import lsst.pipe.base as pipeBase 

39from lsst.utils.timer import timeMethod 

40from .subtractBackground import SubtractBackgroundTask, backgroundFlatContext 

41 

42 

43class PsfGenerationError(pipeBase.AlgorithmError): 

44 """Raised when we cannot generate a PSF for detection 

45 

46 Parameters 

47 ---------- 

48 msg : `str` 

49 Error message. 

50 **kwargs : `dict`, optional 

51 Additional keyword arguments to initialize the Exception base class. 

52 """ 

53 def __init__(self, msg, **kwargs): 

54 self.msg = msg 

55 self._metadata = kwargs 

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

57 

58 def __str__(self): 

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

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

61 

62 @property 

63 def metadata(self): 

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

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

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

67 return self._metadata 

68 

69 

70class SourceDetectionConfig(pexConfig.Config): 

71 """Configuration parameters for the SourceDetectionTask 

72 """ 

73 minPixels = pexConfig.RangeField( 

74 doc="detected sources with fewer than the specified number of pixels will be ignored", 

75 dtype=int, optional=False, default=1, min=0, 

76 ) 

77 isotropicGrow = pexConfig.Field( 

78 doc="Grow pixels as isotropically as possible? If False, use a Manhattan metric instead.", 

79 dtype=bool, default=True, 

80 ) 

81 combinedGrow = pexConfig.Field( 

82 doc="Grow all footprints at the same time? This allows disconnected footprints to merge.", 

83 dtype=bool, default=True, 

84 ) 

85 nSigmaToGrow = pexConfig.Field( 

86 doc="Grow detections by nSigmaToGrow * [PSF RMS width]; if 0 then do not grow", 

87 dtype=float, default=2.4, # 2.4 pixels/sigma is roughly one pixel/FWHM 

88 ) 

89 returnOriginalFootprints = pexConfig.Field( 

90 doc="Grow detections to set the image mask bits, but return the original (not-grown) footprints", 

91 dtype=bool, optional=False, default=False, 

92 ) 

93 thresholdValue = pexConfig.RangeField( 

94 doc="Threshold for detecting footprints; exact meaning and units depend on thresholdType.", 

95 dtype=float, optional=False, default=5.0, min=0.0, 

96 ) 

97 includeThresholdMultiplier = pexConfig.RangeField( 

98 doc="Multiplier on thresholdValue for whether a source is included in the output catalog." 

99 " For example, thresholdValue=5, includeThresholdMultiplier=10, thresholdType='pixel_stdev' " 

100 "results in a catalog of sources at >50 sigma with the detection mask and footprints " 

101 "including pixels >5 sigma.", 

102 dtype=float, default=1.0, min=0.0, 

103 ) 

104 thresholdType = pexConfig.ChoiceField( 

105 doc="Specifies the meaning of thresholdValue.", 

106 dtype=str, optional=False, default="pixel_stdev", 

107 allowed={ 

108 "variance": "threshold applied to image variance", 

109 "stdev": "threshold applied to image std deviation", 

110 "value": "threshold applied to image value", 

111 "pixel_stdev": "threshold applied to per-pixel std deviation", 

112 }, 

113 ) 

114 thresholdPolarity = pexConfig.ChoiceField( 

115 doc="Specifies whether to detect positive, or negative sources, or both.", 

116 dtype=str, optional=False, default="positive", 

117 allowed={ 

118 "positive": "detect only positive sources", 

119 "negative": "detect only negative sources", 

120 "both": "detect both positive and negative sources", 

121 }, 

122 ) 

123 adjustBackground = pexConfig.Field( 

124 dtype=float, 

125 doc="Fiddle factor to add to the background; debugging only", 

126 default=0.0, 

127 ) 

128 reEstimateBackground = pexConfig.Field( 

129 dtype=bool, 

130 doc="Estimate the background again after final source detection?", 

131 default=True, optional=False, 

132 ) 

133 doApplyFlatBackgroundRatio = pexConfig.Field( 

134 doc="Convert from a photometrically flat image to one suitable for background subtraction? " 

135 "Only used if reEstimateBackground is True." 

136 "If True, then a backgroundToPhotometricRatio must be supplied to the task run method.", 

137 dtype=bool, 

138 default=False, 

139 ) 

140 background = pexConfig.ConfigurableField( 

141 doc="Background re-estimation; ignored if reEstimateBackground false", 

142 target=SubtractBackgroundTask, 

143 ) 

144 tempLocalBackground = pexConfig.ConfigurableField( 

145 doc=("A local (small-scale), temporary background estimation step run between " 

146 "detecting above-threshold regions and detecting the peaks within " 

147 "them; used to avoid detecting spuerious peaks in the wings."), 

148 target=SubtractBackgroundTask, 

149 ) 

150 doTempLocalBackground = pexConfig.Field( 

151 dtype=bool, 

152 doc="Enable temporary local background subtraction? (see tempLocalBackground)", 

153 default=True, 

154 ) 

155 tempWideBackground = pexConfig.ConfigurableField( 

156 doc=("A wide (large-scale) background estimation and removal before footprint and peak detection. " 

157 "It is added back into the image after detection. The purpose is to suppress very large " 

158 "footprints (e.g., from large artifacts) that the deblender may choke on."), 

159 target=SubtractBackgroundTask, 

160 ) 

161 doTempWideBackground = pexConfig.Field( 

162 dtype=bool, 

163 doc="Do temporary wide (large-scale) background subtraction before footprint detection?", 

164 default=False, 

165 ) 

166 nPeaksMaxSimple = pexConfig.Field( 

167 dtype=int, 

168 doc=("The maximum number of peaks in a Footprint before trying to " 

169 "replace its peaks using the temporary local background"), 

170 default=1, 

171 ) 

172 nSigmaForKernel = pexConfig.Field( 

173 dtype=float, 

174 doc=("Multiple of PSF RMS size to use for convolution kernel bounding box size; " 

175 "note that this is not a half-size. The size will be rounded up to the nearest odd integer"), 

176 default=7.0, 

177 ) 

178 statsMask = pexConfig.ListField( 

179 dtype=str, 

180 doc="Mask planes to ignore when calculating statistics of image (for thresholdType=stdev)", 

181 default=['BAD', 'SAT', 'EDGE', 'NO_DATA'], 

182 ) 

183 excludeMaskPlanes = lsst.pex.config.ListField( 

184 dtype=str, 

185 default=[], 

186 doc="Mask planes to exclude when detecting sources." 

187 ) 

188 

189 def setDefaults(self): 

190 self.tempLocalBackground.binSize = 64 

191 self.tempLocalBackground.algorithm = "AKIMA_SPLINE" 

192 self.tempLocalBackground.useApprox = False 

193 # Background subtraction to remove a large-scale background (e.g., scattered light); restored later. 

194 # Want to keep it from exceeding the deblender size limit of 1 Mpix, so half that is reasonable. 

195 self.tempWideBackground.binSize = 512 

196 self.tempWideBackground.algorithm = "AKIMA_SPLINE" 

197 self.tempWideBackground.useApprox = False 

198 # Ensure we can remove even bright scattered light that is DETECTED 

199 for maskPlane in ("DETECTED", "DETECTED_NEGATIVE"): 

200 if maskPlane in self.tempWideBackground.ignoredPixelMask: 

201 self.tempWideBackground.ignoredPixelMask.remove(maskPlane) 

202 

203 

204class SourceDetectionTask(pipeBase.Task): 

205 """Detect peaks and footprints of sources in an image. 

206 

207 This task expects the image to have been background subtracted first. 

208 Running detection on images with a non-zero-centered background may result 

209 in a single source detected on the entire image containing thousands of 

210 peaks, or other pathological outputs. 

211 

212 Parameters 

213 ---------- 

214 schema : `lsst.afw.table.Schema` 

215 Schema object used to create the output `lsst.afw.table.SourceCatalog` 

216 **kwds 

217 Keyword arguments passed to `lsst.pipe.base.Task.__init__` 

218 

219 If schema is not None and configured for 'both' detections, 

220 a 'flags.negative' field will be added to label detections made with a 

221 negative threshold. 

222 

223 Notes 

224 ----- 

225 This task convolves the image with a Gaussian approximation to the PSF, 

226 matched to the sigma of the input exposure, because this is separable and 

227 fast. The PSF would have to be very non-Gaussian or non-circular for this 

228 approximation to have a significant impact on the signal-to-noise of the 

229 detected sources. 

230 """ 

231 ConfigClass = SourceDetectionConfig 

232 _DefaultName = "sourceDetection" 

233 

234 def __init__(self, schema=None, **kwds): 

235 pipeBase.Task.__init__(self, **kwds) 

236 if schema is not None and self.config.thresholdPolarity == "both": 

237 self.negativeFlagKey = schema.addField( 

238 "is_negative", type="Flag", 

239 doc="Set if source peak was detected as negative." 

240 ) 

241 else: 

242 if self.config.thresholdPolarity == "both": 

243 self.log.warning("Detection polarity set to 'both', but no flag will be " 

244 "set to distinguish between positive and negative detections") 

245 self.negativeFlagKey = None 

246 if self.config.reEstimateBackground: 

247 self.makeSubtask("background") 

248 if self.config.doTempLocalBackground: 

249 self.makeSubtask("tempLocalBackground") 

250 if self.config.doTempWideBackground: 

251 self.makeSubtask("tempWideBackground") 

252 

253 @timeMethod 

254 def run(self, table, exposure, doSmooth=True, sigma=None, clearMask=True, expId=None, 

255 background=None, backgroundToPhotometricRatio=None): 

256 r"""Detect sources and return catalog(s) of detections. 

257 

258 Parameters 

259 ---------- 

260 table : `lsst.afw.table.SourceTable` 

261 Table object that will be used to create the SourceCatalog. 

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

263 Exposure to process; DETECTED mask plane will be set in-place. 

264 doSmooth : `bool`, optional 

265 If True, smooth the image before detection using a Gaussian of width 

266 ``sigma``, or the measured PSF width. Set to False when running on 

267 e.g. a pre-convolved image, or a mask plane. 

268 sigma : `float`, optional 

269 Sigma of PSF (pixels); used for smoothing and to grow detections; 

270 if None then measure the sigma of the PSF of the exposure 

271 clearMask : `bool`, optional 

272 Clear DETECTED{,_NEGATIVE} planes before running detection. 

273 expId : `int`, optional 

274 Exposure identifier; unused by this implementation, but used for 

275 RNG seed by subclasses. 

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

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

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

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

280 Image to convert photometric-flattened image to 

281 background-flattened image if ``reEstimateBackground=True`` and 

282 exposure has been photometric-flattened. 

283 

284 Returns 

285 ------- 

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

287 The `~lsst.pipe.base.Struct` contains: 

288 

289 ``sources`` 

290 Detected sources on the exposure. 

291 (`lsst.afw.table.SourceCatalog`) 

292 ``positive`` 

293 Positive polarity footprints. 

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

295 ``negative`` 

296 Negative polarity footprints. 

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

298 ``numPos`` 

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

300 negative. (`int`) 

301 ``numNeg`` 

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

303 positive. (`int`) 

304 ``background`` 

305 Re-estimated background. `None` if 

306 ``reEstimateBackground==False``. 

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

308 ``factor`` 

309 Multiplication factor applied to the configured detection 

310 threshold. (`float`) 

311 

312 Raises 

313 ------ 

314 ValueError 

315 Raised if flags.negative is needed, but isn't in table's schema. 

316 lsst.pipe.base.TaskError 

317 Raised if sigma=None, doSmooth=True and the exposure has no PSF. 

318 

319 Notes 

320 ----- 

321 If you want to avoid dealing with Sources and Tables, you can use 

322 `detectFootprints()` to just get the 

323 `~lsst.afw.detection.FootprintSet`\s. 

324 """ 

325 if self.negativeFlagKey is not None and self.negativeFlagKey not in table.getSchema(): 

326 raise ValueError("Table has incorrect Schema") 

327 results = self.detectFootprints(exposure=exposure, doSmooth=doSmooth, sigma=sigma, 

328 clearMask=clearMask, expId=expId, background=background, 

329 backgroundToPhotometricRatio=backgroundToPhotometricRatio) 

330 sources = afwTable.SourceCatalog(table) 

331 sources.reserve(results.numPos + results.numNeg) 

332 if results.negative: 

333 results.negative.makeSources(sources) 

334 if self.negativeFlagKey: 

335 for record in sources: 

336 record.set(self.negativeFlagKey, True) 

337 if results.positive: 

338 results.positive.makeSources(sources) 

339 results.sources = sources 

340 return results 

341 

342 def display(self, exposure, results, convolvedImage=None): 

343 """Display detections if so configured 

344 

345 Displays the ``exposure`` in frame 0, overlays the detection peaks. 

346 

347 Requires that ``lsstDebug`` has been set up correctly, so that 

348 ``lsstDebug.Info("lsst.meas.algorithms.detection")`` evaluates `True`. 

349 

350 If the ``convolvedImage`` is non-`None` and 

351 ``lsstDebug.Info("lsst.meas.algorithms.detection") > 1``, the 

352 ``convolvedImage`` will be displayed in frame 1. 

353 

354 Parameters 

355 ---------- 

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

357 Exposure to display, on which will be plotted the detections. 

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

359 Results of the 'detectFootprints' method, containing positive and 

360 negative footprints (which contain the peak positions that we will 

361 plot). This is a `Struct` with ``positive`` and ``negative`` 

362 elements that are of type `lsst.afw.detection.FootprintSet`. 

363 convolvedImage : `lsst.afw.image.Image`, optional 

364 Convolved image used for thresholding. 

365 """ 

366 try: 

367 import lsstDebug 

368 display = lsstDebug.Info(__name__).display 

369 except ImportError: 

370 try: 

371 display 

372 except NameError: 

373 display = False 

374 if not display: 

375 return 

376 

377 afwDisplay.setDefaultMaskTransparency(75) 

378 

379 disp0 = afwDisplay.Display(frame=0) 

380 disp0.mtv(exposure, title="detection") 

381 

382 def plotPeaks(fps, ctype): 

383 if fps is None: 

384 return 

385 with disp0.Buffering(): 

386 for fp in fps.getFootprints(): 

387 for pp in fp.getPeaks(): 

388 disp0.dot("+", pp.getFx(), pp.getFy(), ctype=ctype) 

389 plotPeaks(results.positive, "yellow") 

390 plotPeaks(results.negative, "red") 

391 

392 if convolvedImage and display > 1: 

393 disp1 = afwDisplay.Display(frame=1) 

394 disp1.mtv(convolvedImage, title="PSF smoothed") 

395 

396 disp2 = afwDisplay.Display(frame=2) 

397 disp2.mtv(afwImage.ImageF(np.sqrt(exposure.variance.array)), title="stddev") 

398 

399 def applyTempLocalBackground(self, exposure, middle, results): 

400 """Apply a temporary local background subtraction 

401 

402 This temporary local background serves to suppress noise fluctuations 

403 in the wings of bright objects. 

404 

405 Peaks in the footprints will be updated. 

406 

407 Parameters 

408 ---------- 

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

410 Exposure for which to fit local background. 

411 middle : `lsst.afw.image.MaskedImage` 

412 Convolved image on which detection will be performed 

413 (typically smaller than ``exposure`` because the 

414 half-kernel has been removed around the edges). 

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

416 Results of the 'detectFootprints' method, containing positive and 

417 negative footprints (which contain the peak positions that we will 

418 plot). This is a `Struct` with ``positive`` and ``negative`` 

419 elements that are of type `lsst.afw.detection.FootprintSet`. 

420 """ 

421 # Subtract the local background from the smoothed image. Since we 

422 # never use the smoothed again we don't need to worry about adding 

423 # it back in. 

424 bg = self.tempLocalBackground.fitBackground(exposure.getMaskedImage()) 

425 bgImage = bg.getImageF(self.tempLocalBackground.config.algorithm, 

426 self.tempLocalBackground.config.undersampleStyle) 

427 middle -= bgImage.Factory(bgImage, middle.getBBox()) 

428 if self.config.thresholdPolarity != "negative": 

429 results.positiveThreshold = self.makeThreshold(middle, "positive") 

430 self.updatePeaks(results.positive, middle, results.positiveThreshold) 

431 if self.config.thresholdPolarity != "positive": 

432 results.negativeThreshold = self.makeThreshold(middle, "negative") 

433 self.updatePeaks(results.negative, middle, results.negativeThreshold) 

434 

435 def clearMask(self, mask): 

436 """Clear the DETECTED and DETECTED_NEGATIVE mask planes. 

437 

438 Removes any previous detection mask in preparation for a new 

439 detection pass. 

440 

441 Parameters 

442 ---------- 

443 mask : `lsst.afw.image.Mask` 

444 Mask to be cleared. 

445 """ 

446 mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE")) 

447 

448 def calculateKernelSize(self, sigma): 

449 """Calculate the size of the smoothing kernel. 

450 

451 Uses the ``nSigmaForKernel`` configuration parameter. Note 

452 that that is the full width of the kernel bounding box 

453 (so a value of 7 means 3.5 sigma on either side of center). 

454 The value will be rounded up to the nearest odd integer. 

455 

456 Parameters 

457 ---------- 

458 sigma : `float` 

459 Gaussian sigma of smoothing kernel. 

460 

461 Returns 

462 ------- 

463 size : `int` 

464 Size of the smoothing kernel. 

465 """ 

466 return (int(sigma * self.config.nSigmaForKernel + 0.5)//2)*2 + 1 # make sure it is odd 

467 

468 def getPsf(self, exposure, sigma=None): 

469 """Create a single Gaussian PSF for an exposure. 

470 

471 If ``sigma`` is provided, we make a `~lsst.afw.detection.GaussianPsf` 

472 with that, otherwise use the sigma from the psf of the ``exposure`` to 

473 make the `~lsst.afw.detection.GaussianPsf`. 

474 

475 Parameters 

476 ---------- 

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

478 Exposure from which to retrieve the PSF. 

479 sigma : `float`, optional 

480 Gaussian sigma to use if provided. 

481 

482 Returns 

483 ------- 

484 psf : `lsst.afw.detection.GaussianPsf` 

485 PSF to use for detection. 

486 

487 Raises 

488 ------ 

489 RuntimeError 

490 Raised if ``sigma`` is not provided and ``exposure`` does not 

491 contain a ``Psf`` object. 

492 """ 

493 if sigma is None: 

494 psf = exposure.getPsf() 

495 if psf is None: 

496 raise PsfGenerationError("Unable to determine PSF to use for detection: no sigma provided") 

497 sigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() 

498 if not np.isfinite(sigma) or sigma <= 0.0: 

499 raise PsfGenerationError("Invalid sigma=%s for PSF used in detection" % (sigma,)) 

500 size = self.calculateKernelSize(sigma) 

501 psf = afwDet.GaussianPsf(size, size, sigma) 

502 return psf 

503 

504 def convolveImage(self, maskedImage, psf, doSmooth=True): 

505 """Convolve the image with the PSF. 

506 

507 We convolve the image with a Gaussian approximation to the PSF, 

508 because this is separable and therefore fast. It's technically a 

509 correlation rather than a convolution, but since we use a symmetric 

510 Gaussian there's no difference. 

511 

512 The convolution can be disabled with ``doSmooth=False``. If we do 

513 convolve, we mask the edges as ``EDGE`` and return the convolved image 

514 with the edges removed. This is because we can't convolve the edges 

515 because the kernel would extend off the image. 

516 

517 Parameters 

518 ---------- 

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

520 Image to convolve. 

521 psf : `lsst.afw.detection.Psf` 

522 PSF to convolve with (actually with a Gaussian approximation 

523 to it). 

524 doSmooth : `bool` 

525 Actually do the convolution? Set to False when running on 

526 e.g. a pre-convolved image, or a mask plane. 

527 

528 Returns 

529 ------- 

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

531 The `~lsst.pipe.base.Struct` contains: 

532 

533 ``middle`` 

534 Convolved image, without the edges. (`lsst.afw.image.MaskedImage`) 

535 ``sigma`` 

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

537 """ 

538 self.metadata["doSmooth"] = doSmooth 

539 sigma = psf.computeShape(psf.getAveragePosition()).getDeterminantRadius() 

540 self.metadata["sigma"] = sigma 

541 

542 if not doSmooth: 

543 middle = maskedImage.Factory(maskedImage, deep=True) 

544 return pipeBase.Struct(middle=middle, sigma=sigma) 

545 

546 # Smooth using a Gaussian (which is separable, hence fast) of width sigma 

547 # Make a SingleGaussian (separable) kernel with the 'sigma' 

548 kWidth = self.calculateKernelSize(sigma) 

549 self.metadata["smoothingKernelWidth"] = kWidth 

550 gaussFunc = afwMath.GaussianFunction1D(sigma) 

551 gaussKernel = afwMath.SeparableKernel(kWidth, kWidth, gaussFunc, gaussFunc) 

552 

553 convolvedImage = maskedImage.Factory(maskedImage.getBBox()) 

554 

555 afwMath.convolve(convolvedImage, maskedImage, gaussKernel, afwMath.ConvolutionControl()) 

556 

557 # Only search psf-smoothed part of frame 

558 goodBBox = gaussKernel.shrinkBBox(convolvedImage.getBBox()) 

559 middle = convolvedImage.Factory(convolvedImage, goodBBox, afwImage.PARENT, False) 

560 

561 # Mark the parts of the image outside goodBBox as EDGE 

562 self.setEdgeBits(maskedImage, goodBBox, maskedImage.getMask().getPlaneBitMask("EDGE")) 

563 

564 return pipeBase.Struct(middle=middle, sigma=sigma) 

565 

566 def applyThreshold(self, middle, bbox, factor=1.0, factorNeg=None): 

567 r"""Apply thresholds to the convolved image 

568 

569 Identifies `~lsst.afw.detection.Footprint`\s, both positive and negative. 

570 The threshold can be modified by the provided multiplication 

571 ``factor``. 

572 

573 Parameters 

574 ---------- 

575 middle : `lsst.afw.image.MaskedImage` 

576 Convolved image to threshold. 

577 bbox : `lsst.geom.Box2I` 

578 Bounding box of unconvolved image. 

579 factor : `float`, optional 

580 Multiplier for the configured threshold. 

581 factorNeg : `float` or `None`, optional 

582 Multiplier for the configured threshold for negative detection 

583 polarity. If `None`, will be set equal to ``factor`` (i.e. equal 

584 to the factor used for positive detection polarity). 

585 

586 Returns 

587 ------- 

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

589 The `~lsst.pipe.base.Struct` contains: 

590 

591 ``positive`` 

592 Positive detection footprints, if configured. 

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

594 ``negative`` 

595 Negative detection footprints, if configured. 

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

597 ``factor`` 

598 Multiplier for the configured threshold for positive detection 

599 polarity. (`float`) 

600 ``factorNeg`` 

601 Multiplier for the configured threshold for negative detection 

602 polarity. (`float`) 

603 """ 

604 if factorNeg is None: 

605 factorNeg = factor 

606 self.log.info("Threshold scaling factor for positive detections is: %.3f. For negative " 

607 "detections it is: %.3f", factor, factorNeg) 

608 results = pipeBase.Struct(positive=None, negative=None, factor=factor, factorNeg=factorNeg, 

609 positiveThreshold=None, negativeThreshold=None) 

610 # Detect the Footprints (peaks may be replaced if doTempLocalBackground) 

611 if self.config.reEstimateBackground or self.config.thresholdPolarity != "negative": 

612 results.positiveThreshold = self.makeThreshold(middle, "positive", factor=factor) 

613 results.positive = afwDet.FootprintSet( 

614 middle, 

615 results.positiveThreshold, 

616 "DETECTED", 

617 self.config.minPixels 

618 ) 

619 results.positive.setRegion(bbox) 

620 if self.config.reEstimateBackground or self.config.thresholdPolarity != "positive": 

621 results.negativeThreshold = self.makeThreshold(middle, "negative", factor=factorNeg) 

622 results.negative = afwDet.FootprintSet( 

623 middle, 

624 results.negativeThreshold, 

625 "DETECTED_NEGATIVE", 

626 self.config.minPixels 

627 ) 

628 results.negative.setRegion(bbox) 

629 

630 return results 

631 

632 def finalizeFootprints(self, mask, results, sigma, factor=1.0, factorNeg=None, growOverride=None): 

633 """Finalize the detected footprints. 

634 

635 Grow the footprints, set the ``DETECTED`` and ``DETECTED_NEGATIVE`` 

636 mask planes, and log the results. 

637 

638 ``numPos`` (number of positive footprints), ``numPosPeaks`` (number 

639 of positive peaks), ``numNeg`` (number of negative footprints), 

640 ``numNegPeaks`` (number of negative peaks) entries are added to the 

641 ``results`` struct. 

642 

643 Parameters 

644 ---------- 

645 mask : `lsst.afw.image.Mask` 

646 Mask image on which to flag detected pixels. 

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

648 Struct of detection results, including ``positive`` and 

649 ``negative`` entries; modified. 

650 sigma : `float` 

651 Gaussian sigma of PSF. 

652 factor : `float`, optional 

653 Multiplier for the configured threshold. Note that this is only 

654 used here for logging purposes. 

655 factorNeg : `float` or `None`, optional 

656 Multiplier used for the negative detection polarity threshold. 

657 If `None`, a factor equal to ``factor`` (i.e. equal to the one used 

658 for positive detection polarity) is assumed. Note that this is only 

659 used here for logging purposes. 

660 growOverride : `float` or `None`, optional 

661 Override to use for ``nSigmaToGrow``, regardless of the value set 

662 in ``config.nSigmaToGrow``. 

663 """ 

664 if growOverride is not None: 

665 self.log.warning("config.nSigmaToGrow is set to %.2f, but the caller has set " 

666 "growOverride to %.2f, so the footprints will be grown by " 

667 "%.2f sigma.", self.config.nSigmaToGrow, growOverride, growOverride) 

668 nSigmaToGrow = growOverride 

669 else: 

670 nSigmaToGrow = self.config.nSigmaToGrow 

671 factorNeg = factor if factorNeg is None else factorNeg 

672 for polarity, maskName in (("positive", "DETECTED"), ("negative", "DETECTED_NEGATIVE")): 

673 fpSet = getattr(results, polarity) 

674 if fpSet is None: 

675 continue 

676 if nSigmaToGrow > 0: 

677 nGrow = int((self.config.nSigmaToGrow * sigma) + 0.5) 

678 self.metadata["nGrow"] = nGrow 

679 if self.config.combinedGrow: 

680 fpSet = afwDet.FootprintSet(fpSet, nGrow, self.config.isotropicGrow) 

681 else: 

682 stencil = (afwGeom.Stencil.CIRCLE if self.config.isotropicGrow else 

683 afwGeom.Stencil.MANHATTAN) 

684 for fp in fpSet: 

685 fp.dilate(nGrow, stencil) 

686 fpSet.setMask(mask, maskName) 

687 if not self.config.returnOriginalFootprints: 

688 setattr(results, polarity, fpSet) 

689 

690 results.numPos = 0 

691 results.numPosPeaks = 0 

692 results.numNeg = 0 

693 results.numNegPeaks = 0 

694 positive = "" 

695 negative = "" 

696 

697 if results.positive is not None: 

698 results.numPos = len(results.positive.getFootprints()) 

699 results.numPosPeaks = sum(len(fp.getPeaks()) for fp in results.positive.getFootprints()) 

700 positive = " %d positive peaks in %d footprints" % (results.numPosPeaks, results.numPos) 

701 if results.negative is not None: 

702 results.numNeg = len(results.negative.getFootprints()) 

703 results.numNegPeaks = sum(len(fp.getPeaks()) for fp in results.negative.getFootprints()) 

704 negative = " %d negative peaks in %d footprints" % (results.numNegPeaks, results.numNeg) 

705 

706 self.log.info("Detected%s%s%s to %g +ve and %g -ve %s", 

707 positive, " and" if positive and negative else "", negative, 

708 self.config.thresholdValue*self.config.includeThresholdMultiplier*factor, 

709 self.config.thresholdValue*self.config.includeThresholdMultiplier*factorNeg, 

710 "DN" if self.config.thresholdType == "value" else "sigma") 

711 

712 def reEstimateBackground(self, maskedImage, backgrounds, backgroundToPhotometricRatio=None): 

713 """Estimate the background after detection 

714 

715 Parameters 

716 ---------- 

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

718 Image on which to estimate the background. 

719 backgrounds : `lsst.afw.math.BackgroundList` 

720 List of backgrounds; modified. 

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

722 Image to multiply a photometrically-flattened image by to obtain a 

723 background-flattened image. 

724 Only used if ``config.doApplyFlatBackgroundRatio`` is ``True``. 

725 

726 Returns 

727 ------- 

728 bg : `lsst.afw.math.backgroundMI` 

729 Empirical background model. 

730 """ 

731 with backgroundFlatContext( 

732 maskedImage, 

733 self.config.doApplyFlatBackgroundRatio, 

734 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

735 ): 

736 bg = self.background.fitBackground(maskedImage) 

737 if self.config.adjustBackground: 

738 self.log.warning("Fiddling the background by %g", self.config.adjustBackground) 

739 bg += self.config.adjustBackground 

740 self.log.info("Resubtracting the background after object detection (median background " 

741 "value = %.2f)", np.median(bg.getImageF().array)) 

742 maskedImage -= bg.getImageF(self.background.config.algorithm, 

743 self.background.config.undersampleStyle) 

744 

745 actrl = bg.getBackgroundControl().getApproximateControl() 

746 backgrounds.append((bg, getattr(afwMath.Interpolate, self.background.config.algorithm), 

747 bg.getAsUsedUndersampleStyle(), actrl.getStyle(), actrl.getOrderX(), 

748 actrl.getOrderY(), actrl.getWeighting())) 

749 return bg 

750 

751 def clearUnwantedResults(self, mask, results): 

752 """Clear unwanted results from the Struct of results 

753 

754 If we specifically want only positive or only negative detections, 

755 drop the ones we don't want, and its associated mask plane. 

756 

757 Parameters 

758 ---------- 

759 mask : `lsst.afw.image.Mask` 

760 Mask image. 

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

762 Detection results, with ``positive`` and ``negative`` elements; 

763 modified. 

764 """ 

765 if self.config.thresholdPolarity == "positive": 

766 if self.config.reEstimateBackground: 

767 mask &= ~mask.getPlaneBitMask("DETECTED_NEGATIVE") 

768 results.negative = None 

769 elif self.config.thresholdPolarity == "negative": 

770 if self.config.reEstimateBackground: 

771 mask &= ~mask.getPlaneBitMask("DETECTED") 

772 results.positive = None 

773 

774 @timeMethod 

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

776 background=None, backgroundToPhotometricRatio=None, factor=1.0, factorNeg=None): 

777 """Detect footprints on an exposure. 

778 

779 Parameters 

780 ---------- 

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

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

783 set in-place. 

784 doSmooth : `bool`, optional 

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

786 of width ``sigma``, or the measured PSF width of ``exposure``. 

787 Set to False when running on e.g. a pre-convolved image, or a mask 

788 plane. 

789 sigma : `float`, optional 

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

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

792 ``exposure``. 

793 clearMask : `bool`, optional 

794 Clear both DETECTED and DETECTED_NEGATIVE planes before running 

795 detection. 

796 expId : `dict`, optional 

797 Exposure identifier; unused by this implementation, but used for 

798 RNG seed by subclasses. 

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

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

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

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

803 Image to convert photometric-flattened image to 

804 background-flattened image if ``reEstimateBackground=True`` and 

805 exposure has been photometric-flattened. 

806 factor : `float`, optional 

807 Multiplier for the configured threshold for positive detection 

808 polarity. 

809 factorNeg : `float` or `None`, optional 

810 Multiplier for the configured threshold for negative detection 

811 polarity. If `None`, will be set equal to ``factor`` (i.e. equal 

812 to the factor used for positive detection polarity). 

813 

814 Returns 

815 ------- 

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

817 A `~lsst.pipe.base.Struct` containing: 

818 

819 ``positive`` 

820 Positive polarity footprints. 

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

822 ``negative`` 

823 Negative polarity footprints. 

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

825 ``numPos`` 

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

827 negative. (`int`) 

828 ``numNeg`` 

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

830 positive. (`int`) 

831 ``background`` 

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

833 if ``reEstimateBackground==False``. 

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

835 ``factor`` 

836 Multiplication factor applied to the configured threshold 

837 for positive detection polarity. (`float`) 

838 ``factorNeg`` 

839 Multiplication factor applied to the configured threshold 

840 for negative detection polarity. (`float`) 

841 """ 

842 maskedImage = exposure.maskedImage 

843 

844 if clearMask: 

845 self.clearMask(maskedImage.getMask()) 

846 

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

848 with self.tempWideBackgroundContext(exposure): 

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

850 middle = convolveResults.middle 

851 sigma = convolveResults.sigma 

852 self.removeBadPixels(middle) 

853 

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

855 results.background = background if background is not None else afwMath.BackgroundList() 

856 

857 if self.config.doTempLocalBackground: 

858 self.applyTempLocalBackground(exposure, middle, results) 

859 self.finalizeFootprints(maskedImage.mask, results, sigma, factor=factor, factorNeg=factorNeg) 

860 

861 # Compute the significance of peaks after the peaks have been 

862 # finalized and after local background correction/updatePeaks, so 

863 # that the significance represents the "final" detection S/N. 

864 results.positive = self.setPeakSignificance(middle, results.positive, results.positiveThreshold) 

865 results.negative = self.setPeakSignificance(middle, results.negative, results.negativeThreshold, 

866 negative=True) 

867 

868 if self.config.reEstimateBackground: 

869 self.reEstimateBackground( 

870 maskedImage, 

871 results.background, 

872 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

873 ) 

874 

875 self.clearUnwantedResults(maskedImage.getMask(), results) 

876 

877 self.display(exposure, results, middle) 

878 

879 return results 

880 

881 def removeBadPixels(self, middle): 

882 """Set the significance of flagged pixels to zero. 

883 

884 Parameters 

885 ---------- 

886 middle : `lsst.afw.image.Exposure` 

887 Score or maximum likelihood difference image. 

888 The image plane will be modified in place. 

889 """ 

890 badPixelMask = middle.mask.getPlaneBitMask(self.config.excludeMaskPlanes) 

891 badPixels = middle.mask.array & badPixelMask > 0 

892 middle.image.array[badPixels] = 0 

893 

894 def setPeakSignificance(self, exposure, footprints, threshold, negative=False): 

895 """Set the significance of each detected peak to the pixel value divided 

896 by the appropriate standard-deviation for ``config.thresholdType``. 

897 

898 Only sets significance for "stdev" and "pixel_stdev" thresholdTypes; 

899 we leave it undefined for "value" and "variance" as it does not have a 

900 well-defined meaning in those cases. 

901 

902 Parameters 

903 ---------- 

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

905 Exposure that footprints were detected on, likely the convolved, 

906 local background-subtracted image. 

907 footprints : `lsst.afw.detection.FootprintSet` 

908 Footprints detected on the image. 

909 threshold : `lsst.afw.detection.Threshold` 

910 Threshold used to find footprints. 

911 negative : `bool`, optional 

912 Are we calculating for negative sources? 

913 """ 

914 if footprints is None or footprints.getFootprints() == []: 

915 return footprints 

916 polarity = -1 if negative else 1