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` 

580 Multiplier for the configured threshold. 

581 factorNeg : `float` or `None` 

582 Multiplier for the configured threshold for negative detection polarity. 

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

584 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. 

599 (`float`) 

600 ``factorNeg`` 

601 Multiplier for the configured threshold for negative detection polarity. 

602 (`float`) 

603 """ 

604 if factorNeg is None: 

605 factorNeg = factor 

606 self.log.info("Setting factor for negative detections equal to that for positive " 

607 "detections: %f", factor) 

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` 

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

654 used here for logging purposes. 

655 factorNeg : `float` or `None` 

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 """ 

661 if growOverride is not None: 

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

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

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

665 factorNeg = factor if factorNeg is None else factorNeg 

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

667 fpSet = getattr(results, polarity) 

668 if fpSet is None: 

669 continue 

670 if self.config.nSigmaToGrow > 0: 

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

672 self.metadata["nGrow"] = nGrow 

673 if self.config.combinedGrow: 

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

675 else: 

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

677 afwGeom.Stencil.MANHATTAN) 

678 for fp in fpSet: 

679 fp.dilate(nGrow, stencil) 

680 fpSet.setMask(mask, maskName) 

681 if not self.config.returnOriginalFootprints: 

682 setattr(results, polarity, fpSet) 

683 

684 results.numPos = 0 

685 results.numPosPeaks = 0 

686 results.numNeg = 0 

687 results.numNegPeaks = 0 

688 positive = "" 

689 negative = "" 

690 

691 if results.positive is not None: 

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

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

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

695 if results.negative is not None: 

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

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

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

699 

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

701 positive, " and" if positive and negative else "", negative, 

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

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

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

705 

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

707 """Estimate the background after detection 

708 

709 Parameters 

710 ---------- 

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

712 Image on which to estimate the background. 

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

714 List of backgrounds; modified. 

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

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

717 background-flattened image. 

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

719 

720 Returns 

721 ------- 

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

723 Empirical background model. 

724 """ 

725 with backgroundFlatContext( 

726 maskedImage, 

727 self.config.doApplyFlatBackgroundRatio, 

728 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

729 ): 

730 bg = self.background.fitBackground(maskedImage) 

731 if self.config.adjustBackground: 

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

733 bg += self.config.adjustBackground 

734 self.log.info("Resubtracting the background after object detection") 

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

736 self.background.config.undersampleStyle) 

737 

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

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

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

741 actrl.getOrderY(), actrl.getWeighting())) 

742 return bg 

743 

744 def clearUnwantedResults(self, mask, results): 

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

746 

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

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

749 

750 Parameters 

751 ---------- 

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

753 Mask image. 

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

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

756 modified. 

757 """ 

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

759 if self.config.reEstimateBackground: 

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

761 results.negative = None 

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

763 if self.config.reEstimateBackground: 

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

765 results.positive = None 

766 

767 @timeMethod 

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

769 background=None, backgroundToPhotometricRatio=None): 

770 """Detect footprints on an exposure. 

771 

772 Parameters 

773 ---------- 

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

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

776 set in-place. 

777 doSmooth : `bool`, optional 

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

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

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

781 plane. 

782 sigma : `float`, optional 

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

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

785 ``exposure``. 

786 clearMask : `bool`, optional 

787 Clear both DETECTED and DETECTED_NEGATIVE planes before running 

788 detection. 

789 expId : `dict`, optional 

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

791 RNG seed by subclasses. 

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

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

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

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

796 Image to convert photometric-flattened image to 

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

798 exposure has been photometric-flattened. 

799 

800 Returns 

801 ------- 

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

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

804 

805 ``positive`` 

806 Positive polarity footprints. 

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

808 ``negative`` 

809 Negative polarity footprints. 

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

811 ``numPos`` 

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

813 negative. (`int`) 

814 ``numNeg`` 

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

816 positive. (`int`) 

817 ``background`` 

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

819 if ``reEstimateBackground==False``. 

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

821 ``factor`` 

822 Multiplication factor applied to the configured detection 

823 threshold. (`float`) 

824 """ 

825 maskedImage = exposure.maskedImage 

826 

827 if clearMask: 

828 self.clearMask(maskedImage.getMask()) 

829 

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

831 with self.tempWideBackgroundContext(exposure): 

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

833 middle = convolveResults.middle 

834 sigma = convolveResults.sigma 

835 self.removeBadPixels(middle) 

836 

837 results = self.applyThreshold(middle, maskedImage.getBBox()) 

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

839 

840 if self.config.doTempLocalBackground: 

841 self.applyTempLocalBackground(exposure, middle, results) 

842 self.finalizeFootprints(maskedImage.mask, results, sigma) 

843 

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

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

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

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

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

849 negative=True) 

850 

851 if self.config.reEstimateBackground: 

852 self.reEstimateBackground( 

853 maskedImage, 

854 results.background, 

855 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

856 ) 

857 

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

859 

860 self.display(exposure, results, middle) 

861 

862 return results 

863 

864 def removeBadPixels(self, middle): 

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

866 

867 Parameters 

868 ---------- 

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

870 Score or maximum likelihood difference image. 

871 The image plane will be modified in place. 

872 """ 

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

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

875 middle.image.array[badPixels] = 0 

876 

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

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

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

880 

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

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

883 well-defined meaning in those cases. 

884 

885 Parameters 

886 ---------- 

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

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

889 local background-subtracted image. 

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

891 Footprints detected on the image. 

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

893 Threshold used to find footprints. 

894 negative : `bool`, optional 

895 Are we calculating for negative sources? 

896 """ 

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

898 return footprints 

899 polarity = -1 if negative else 1 

900 

901 # All incoming footprints have the same schema. 

902 mapper = afwTable.SchemaMapper(footprints.getFootprints()[0].peaks.schema) 

903 mapper.addMinimalSchema(footprints.getFootprints()[0].peaks.schema) 

904 mapper.addOutputField("significance", type=float, 

905 doc="Ratio of peak value to configured standard deviation.") 

906 

907 # Copy the old peaks to the new ones with a significance field. 

908 # Do this independent of the threshold type, so we always have a 

909 # significance field. 

910 newFootprints = afwDet.FootprintSet(footprints) 

911 for old, new in zip(footprints.getFootprints(), newFootprints.getFootprints()): 

912 newPeaks = afwDet.PeakCatalog(mapper.getOutputSchema())