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

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__ = ("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 SourceDetectionConfig(pexConfig.Config): 

44 """Configuration parameters for the SourceDetectionTask 

45 """ 

46 minPixels = pexConfig.RangeField( 

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

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

49 ) 

50 isotropicGrow = pexConfig.Field( 

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

52 dtype=bool, default=True, 

53 ) 

54 combinedGrow = pexConfig.Field( 

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

56 dtype=bool, default=True, 

57 ) 

58 nSigmaToGrow = pexConfig.Field( 

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

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

61 ) 

62 returnOriginalFootprints = pexConfig.Field( 

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

64 dtype=bool, optional=False, default=False, 

65 ) 

66 thresholdValue = pexConfig.RangeField( 

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

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

69 ) 

70 includeThresholdMultiplier = pexConfig.RangeField( 

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

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

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

74 "including pixels >5 sigma.", 

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

76 ) 

77 thresholdType = pexConfig.ChoiceField( 

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

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

80 allowed={ 

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

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

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

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

85 }, 

86 ) 

87 thresholdPolarity = pexConfig.ChoiceField( 

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

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

90 allowed={ 

91 "positive": "detect only positive sources", 

92 "negative": "detect only negative sources", 

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

94 }, 

95 ) 

96 adjustBackground = pexConfig.Field( 

97 dtype=float, 

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

99 default=0.0, 

100 ) 

101 reEstimateBackground = pexConfig.Field( 

102 dtype=bool, 

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

104 default=True, optional=False, 

105 ) 

106 doApplyFlatBackgroundRatio = pexConfig.Field( 

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

108 "Only used if reEstimateBackground is True." 

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

110 dtype=bool, 

111 default=False, 

112 ) 

113 background = pexConfig.ConfigurableField( 

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

115 target=SubtractBackgroundTask, 

116 ) 

117 tempLocalBackground = pexConfig.ConfigurableField( 

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

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

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

121 target=SubtractBackgroundTask, 

122 ) 

123 doTempLocalBackground = pexConfig.Field( 

124 dtype=bool, 

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

126 default=True, 

127 ) 

128 tempWideBackground = pexConfig.ConfigurableField( 

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

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

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

132 target=SubtractBackgroundTask, 

133 ) 

134 doTempWideBackground = pexConfig.Field( 

135 dtype=bool, 

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

137 default=False, 

138 ) 

139 nPeaksMaxSimple = pexConfig.Field( 

140 dtype=int, 

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

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

143 default=1, 

144 ) 

145 nSigmaForKernel = pexConfig.Field( 

146 dtype=float, 

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

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

149 default=7.0, 

150 ) 

151 statsMask = pexConfig.ListField( 

152 dtype=str, 

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

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

155 ) 

156 excludeMaskPlanes = lsst.pex.config.ListField( 

157 dtype=str, 

158 default=[], 

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

160 ) 

161 

162 def setDefaults(self): 

163 self.tempLocalBackground.binSize = 64 

164 self.tempLocalBackground.algorithm = "AKIMA_SPLINE" 

165 self.tempLocalBackground.useApprox = False 

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

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

168 self.tempWideBackground.binSize = 512 

169 self.tempWideBackground.algorithm = "AKIMA_SPLINE" 

170 self.tempWideBackground.useApprox = False 

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

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

173 if maskPlane in self.tempWideBackground.ignoredPixelMask: 

174 self.tempWideBackground.ignoredPixelMask.remove(maskPlane) 

175 

176 

177class SourceDetectionTask(pipeBase.Task): 

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

179 

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

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

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

183 peaks, or other pathological outputs. 

184 

185 Parameters 

186 ---------- 

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

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

189 **kwds 

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

191 

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

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

194 negative threshold. 

195 

196 Notes 

197 ----- 

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

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

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

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

202 detected sources. 

203 """ 

204 ConfigClass = SourceDetectionConfig 

205 _DefaultName = "sourceDetection" 

206 

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

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

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

210 self.negativeFlagKey = schema.addField( 

211 "is_negative", type="Flag", 

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

213 ) 

214 else: 

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

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

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

218 self.negativeFlagKey = None 

219 if self.config.reEstimateBackground: 

220 self.makeSubtask("background") 

221 if self.config.doTempLocalBackground: 

222 self.makeSubtask("tempLocalBackground") 

223 if self.config.doTempWideBackground: 

224 self.makeSubtask("tempWideBackground") 

225 

226 @timeMethod 

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

228 background=None, backgroundToPhotometricRatio=None): 

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

230 

231 Parameters 

232 ---------- 

233 table : `lsst.afw.table.SourceTable` 

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

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

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

237 doSmooth : `bool`, optional 

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

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

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

241 sigma : `float`, optional 

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

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

244 clearMask : `bool`, optional 

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

246 expId : `int`, optional 

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

248 RNG seed by subclasses. 

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

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

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

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

253 Image to convert photometric-flattened image to 

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

255 exposure has been photometric-flattened. 

256 

257 Returns 

258 ------- 

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

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

261 

262 ``sources`` 

263 Detected sources on the exposure. 

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

265 ``positive`` 

266 Positive polarity footprints. 

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

268 ``negative`` 

269 Negative polarity footprints. 

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

271 ``numPos`` 

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

273 negative. (`int`) 

274 ``numNeg`` 

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

276 positive. (`int`) 

277 ``background`` 

278 Re-estimated background. `None` if 

279 ``reEstimateBackground==False``. 

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

281 ``factor`` 

282 Multiplication factor applied to the configured detection 

283 threshold. (`float`) 

284 

285 Raises 

286 ------ 

287 ValueError 

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

289 lsst.pipe.base.TaskError 

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

291 

292 Notes 

293 ----- 

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

295 `detectFootprints()` to just get the 

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

297 """ 

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

299 raise ValueError("Table has incorrect Schema") 

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

301 clearMask=clearMask, expId=expId, background=background, 

302 backgroundToPhotometricRatio=backgroundToPhotometricRatio) 

303 sources = afwTable.SourceCatalog(table) 

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

305 if results.negative: 

306 results.negative.makeSources(sources) 

307 if self.negativeFlagKey: 

308 for record in sources: 

309 record.set(self.negativeFlagKey, True) 

310 if results.positive: 

311 results.positive.makeSources(sources) 

312 results.sources = sources 

313 return results 

314 

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

316 """Display detections if so configured 

317 

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

319 

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

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

322 

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

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

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

326 

327 Parameters 

328 ---------- 

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

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

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

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

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

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

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

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

337 Convolved image used for thresholding. 

338 """ 

339 try: 

340 import lsstDebug 

341 display = lsstDebug.Info(__name__).display 

342 except ImportError: 

343 try: 

344 display 

345 except NameError: 

346 display = False 

347 if not display: 

348 return 

349 

350 afwDisplay.setDefaultMaskTransparency(75) 

351 

352 disp0 = afwDisplay.Display(frame=0) 

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

354 

355 def plotPeaks(fps, ctype): 

356 if fps is None: 

357 return 

358 with disp0.Buffering(): 

359 for fp in fps.getFootprints(): 

360 for pp in fp.getPeaks(): 

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

362 plotPeaks(results.positive, "yellow") 

363 plotPeaks(results.negative, "red") 

364 

365 if convolvedImage and display > 1: 

366 disp1 = afwDisplay.Display(frame=1) 

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

368 

369 disp2 = afwDisplay.Display(frame=2) 

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

371 

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

373 """Apply a temporary local background subtraction 

374 

375 This temporary local background serves to suppress noise fluctuations 

376 in the wings of bright objects. 

377 

378 Peaks in the footprints will be updated. 

379 

380 Parameters 

381 ---------- 

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

383 Exposure for which to fit local background. 

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

385 Convolved image on which detection will be performed 

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

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

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

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

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

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

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

393 """ 

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

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

396 # it back in. 

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

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

399 self.tempLocalBackground.config.undersampleStyle) 

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

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

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

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

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

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

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

407 

408 def clearMask(self, mask): 

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

410 

411 Removes any previous detection mask in preparation for a new 

412 detection pass. 

413 

414 Parameters 

415 ---------- 

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

417 Mask to be cleared. 

418 """ 

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

420 

421 def calculateKernelSize(self, sigma): 

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

423 

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

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

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

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

428 

429 Parameters 

430 ---------- 

431 sigma : `float` 

432 Gaussian sigma of smoothing kernel. 

433 

434 Returns 

435 ------- 

436 size : `int` 

437 Size of the smoothing kernel. 

438 """ 

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

440 

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

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

443 

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

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

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

447 

448 Parameters 

449 ---------- 

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

451 Exposure from which to retrieve the PSF. 

452 sigma : `float`, optional 

453 Gaussian sigma to use if provided. 

454 

455 Returns 

456 ------- 

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

458 PSF to use for detection. 

459 

460 Raises 

461 ------ 

462 RuntimeError 

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

464 contain a ``Psf`` object. 

465 """ 

466 if sigma is None: 

467 psf = exposure.getPsf() 

468 if psf is None: 

469 raise RuntimeError("Unable to determine PSF to use for detection: no sigma provided") 

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

471 size = self.calculateKernelSize(sigma) 

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

473 return psf 

474 

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

476 """Convolve the image with the PSF. 

477 

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

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

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

481 Gaussian there's no difference. 

482 

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

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

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

486 because the kernel would extend off the image. 

487 

488 Parameters 

489 ---------- 

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

491 Image to convolve. 

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

493 PSF to convolve with (actually with a Gaussian approximation 

494 to it). 

495 doSmooth : `bool` 

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

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

498 

499 Returns 

500 ------- 

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

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

503 

504 ``middle`` 

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

506 ``sigma`` 

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

508 """ 

509 self.metadata["doSmooth"] = doSmooth 

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

511 self.metadata["sigma"] = sigma 

512 

513 if not doSmooth: 

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

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

516 

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

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

519 kWidth = self.calculateKernelSize(sigma) 

520 self.metadata["smoothingKernelWidth"] = kWidth 

521 gaussFunc = afwMath.GaussianFunction1D(sigma) 

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

523 

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

525 

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

527 

528 # Only search psf-smoothed part of frame 

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

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

531 

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

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

534 

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

536 

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

538 r"""Apply thresholds to the convolved image 

539 

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

541 The threshold can be modified by the provided multiplication 

542 ``factor``. 

543 

544 Parameters 

545 ---------- 

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

547 Convolved image to threshold. 

548 bbox : `lsst.geom.Box2I` 

549 Bounding box of unconvolved image. 

550 factor : `float` 

551 Multiplier for the configured threshold. 

552 factorNeg : `float` or `None` 

553 Multiplier for the configured threshold for negative detection polarity. 

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

555 for positive detection polarity). 

556 

557 Returns 

558 ------- 

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

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

561 

562 ``positive`` 

563 Positive detection footprints, if configured. 

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

565 ``negative`` 

566 Negative detection footprints, if configured. 

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

568 ``factor`` 

569 Multiplier for the configured threshold. 

570 (`float`) 

571 ``factorNeg`` 

572 Multiplier for the configured threshold for negative detection polarity. 

573 (`float`) 

574 """ 

575 if factorNeg is None: 

576 factorNeg = factor 

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

578 "detections: %f", factor) 

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

580 positiveThreshold=None, negativeThreshold=None) 

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

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

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

584 results.positive = afwDet.FootprintSet( 

585 middle, 

586 results.positiveThreshold, 

587 "DETECTED", 

588 self.config.minPixels 

589 ) 

590 results.positive.setRegion(bbox) 

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

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

593 results.negative = afwDet.FootprintSet( 

594 middle, 

595 results.negativeThreshold, 

596 "DETECTED_NEGATIVE", 

597 self.config.minPixels 

598 ) 

599 results.negative.setRegion(bbox) 

600 

601 return results 

602 

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

604 """Finalize the detected footprints. 

605 

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

607 mask planes, and log the results. 

608 

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

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

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

612 ``results`` struct. 

613 

614 Parameters 

615 ---------- 

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

617 Mask image on which to flag detected pixels. 

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

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

620 ``negative`` entries; modified. 

621 sigma : `float` 

622 Gaussian sigma of PSF. 

623 factor : `float` 

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

625 used here for logging purposes. 

626 factorNeg : `float` or `None` 

627 Multiplier used for the negative detection polarity threshold. 

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

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

630 used here for logging purposes. 

631 """ 

632 factorNeg = factor if factorNeg is None else factorNeg 

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

634 fpSet = getattr(results, polarity) 

635 if fpSet is None: 

636 continue 

637 if self.config.nSigmaToGrow > 0: 

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

639 self.metadata["nGrow"] = nGrow 

640 if self.config.combinedGrow: 

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

642 else: 

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

644 afwGeom.Stencil.MANHATTAN) 

645 for fp in fpSet: 

646 fp.dilate(nGrow, stencil) 

647 fpSet.setMask(mask, maskName) 

648 if not self.config.returnOriginalFootprints: 

649 setattr(results, polarity, fpSet) 

650 

651 results.numPos = 0 

652 results.numPosPeaks = 0 

653 results.numNeg = 0 

654 results.numNegPeaks = 0 

655 positive = "" 

656 negative = "" 

657 

658 if results.positive is not None: 

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

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

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

662 if results.negative is not None: 

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

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

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

666 

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

668 positive, " and" if positive and negative else "", negative, 

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

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

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

672 

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

674 """Estimate the background after detection 

675 

676 Parameters 

677 ---------- 

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

679 Image on which to estimate the background. 

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

681 List of backgrounds; modified. 

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

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

684 background-flattened image. 

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

686 

687 Returns 

688 ------- 

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

690 Empirical background model. 

691 """ 

692 with backgroundFlatContext( 

693 maskedImage, 

694 self.config.doApplyFlatBackgroundRatio, 

695 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

696 ): 

697 bg = self.background.fitBackground(maskedImage) 

698 if self.config.adjustBackground: 

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

700 bg += self.config.adjustBackground 

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

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

703 self.background.config.undersampleStyle) 

704 

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

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

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

708 actrl.getOrderY(), actrl.getWeighting())) 

709 return bg 

710 

711 def clearUnwantedResults(self, mask, results): 

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

713 

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

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

716 

717 Parameters 

718 ---------- 

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

720 Mask image. 

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

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

723 modified. 

724 """ 

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

726 if self.config.reEstimateBackground: 

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

728 results.negative = None 

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

730 if self.config.reEstimateBackground: 

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

732 results.positive = None 

733 

734 @timeMethod 

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

736 background=None, backgroundToPhotometricRatio=None): 

737 """Detect footprints on an exposure. 

738 

739 Parameters 

740 ---------- 

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

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

743 set in-place. 

744 doSmooth : `bool`, optional 

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

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

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

748 plane. 

749 sigma : `float`, optional 

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

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

752 ``exposure``. 

753 clearMask : `bool`, optional 

754 Clear both DETECTED and DETECTED_NEGATIVE planes before running 

755 detection. 

756 expId : `dict`, optional 

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

758 RNG seed by subclasses. 

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

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

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

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

763 Image to convert photometric-flattened image to 

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

765 exposure has been photometric-flattened. 

766 

767 Returns 

768 ------- 

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

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

771 

772 ``positive`` 

773 Positive polarity footprints. 

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

775 ``negative`` 

776 Negative polarity footprints. 

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

778 ``numPos`` 

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

780 negative. (`int`) 

781 ``numNeg`` 

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

783 positive. (`int`) 

784 ``background`` 

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

786 if ``reEstimateBackground==False``. 

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

788 ``factor`` 

789 Multiplication factor applied to the configured detection 

790 threshold. (`float`) 

791 """ 

792 maskedImage = exposure.maskedImage 

793 

794 if clearMask: 

795 self.clearMask(maskedImage.getMask()) 

796 

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

798 with self.tempWideBackgroundContext(exposure): 

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

800 middle = convolveResults.middle 

801 sigma = convolveResults.sigma 

802 self.removeBadPixels(middle) 

803 

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

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

806 

807 if self.config.doTempLocalBackground: 

808 self.applyTempLocalBackground(exposure, middle, results) 

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

810 

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

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

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

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

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

816 negative=True) 

817 

818 if self.config.reEstimateBackground: 

819 self.reEstimateBackground( 

820 maskedImage, 

821 results.background, 

822 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

823 ) 

824 

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

826 

827 self.display(exposure, results, middle) 

828 

829 return results 

830 

831 def removeBadPixels(self, middle): 

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

833 

834 Parameters 

835 ---------- 

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

837 Score or maximum likelihood difference image. 

838 The image plane will be modified in place. 

839 """ 

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

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

842 middle.image.array[badPixels] = 0 

843 

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

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

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

847 

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

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

850 well-defined meaning in those cases. 

851 

852 Parameters 

853 ---------- 

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

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

856 local background-subtracted image. 

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

858 Footprints detected on the image. 

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

860 Threshold used to find footprints. 

861 negative : `bool`, optional 

862 Are we calculating for negative sources? 

863 """ 

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

865 return footprints 

866 polarity = -1 if negative else 1 

867 

868 # All incoming footprints have the same schema. 

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

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

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

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

873 

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

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

876 # significance field. 

877 newFootprints = afwDet.FootprintSet(footprints) 

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

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

880 newPeaks.extend(old.peaks, mapper=mapper) 

881 new.getPeaks().clear() 

882 new.setPeakCatalog(newPeaks) 

883 

884 # Compute the significance values. 

885 if self.config.thresholdType == "pixel_stdev": 

886 for footprint in newFootprints.getFootprints(): 

887 footprint.updatePeakSignificance(exposure.variance, polarity) 

888 elif self.config.thresholdType == "stdev": 

889 sigma = threshold.getValue() / self.config.thresholdValue 

890 for footprint in newFootprints.getFootprints(): 

891 footprint.updatePeakSignificance(polarity*sigma) 

892 else: 

893 for footprint in newFootprints.getFootprints(): 

894 for peak in footprint.peaks: 

895 peak["significance"] = 0 

896 

897 return newFootprints 

898 

899 def makeThreshold(self, image, thresholdParity, factor=1.0): 

900 """Make an afw.detection.Threshold object corresponding to the task's 

901 configuration and the statistics of the given image. 

902 

903 Parameters 

904 ---------- 

905 image : `afw.image.MaskedImage` 

906 Image to measure noise statistics from if needed. 

907 thresholdParity: `str` 

908 One of "positive" or "negative", to set the kind of fluctuations 

909 the Threshold will detect. 

910 factor : `float` 

911 Factor by which to multiply the configured detection threshold. 

912 This is useful for tweaking the detection threshold slightly. 

913 

914 Returns 

915 ------- 

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