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, growOverride=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 if growOverride is not None: 

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

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

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

636 factorNeg = factor if factorNeg is None else factorNeg 

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

638 fpSet = getattr(results, polarity) 

639 if fpSet is None: 

640 continue 

641 if self.config.nSigmaToGrow > 0: 

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

643 self.metadata["nGrow"] = nGrow 

644 if self.config.combinedGrow: 

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

646 else: 

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

648 afwGeom.Stencil.MANHATTAN) 

649 for fp in fpSet: 

650 fp.dilate(nGrow, stencil) 

651 fpSet.setMask(mask, maskName) 

652 if not self.config.returnOriginalFootprints: 

653 setattr(results, polarity, fpSet) 

654 

655 results.numPos = 0 

656 results.numPosPeaks = 0 

657 results.numNeg = 0 

658 results.numNegPeaks = 0 

659 positive = "" 

660 negative = "" 

661 

662 if results.positive is not None: 

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

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

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

666 if results.negative is not None: 

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

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

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

670 

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

672 positive, " and" if positive and negative else "", negative, 

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

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

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

676 

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

678 """Estimate the background after detection 

679 

680 Parameters 

681 ---------- 

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

683 Image on which to estimate the background. 

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

685 List of backgrounds; modified. 

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

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

688 background-flattened image. 

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

690 

691 Returns 

692 ------- 

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

694 Empirical background model. 

695 """ 

696 with backgroundFlatContext( 

697 maskedImage, 

698 self.config.doApplyFlatBackgroundRatio, 

699 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

700 ): 

701 bg = self.background.fitBackground(maskedImage) 

702 if self.config.adjustBackground: 

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

704 bg += self.config.adjustBackground 

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

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

707 self.background.config.undersampleStyle) 

708 

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

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

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

712 actrl.getOrderY(), actrl.getWeighting())) 

713 return bg 

714 

715 def clearUnwantedResults(self, mask, results): 

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

717 

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

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

720 

721 Parameters 

722 ---------- 

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

724 Mask image. 

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

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

727 modified. 

728 """ 

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

730 if self.config.reEstimateBackground: 

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

732 results.negative = None 

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

734 if self.config.reEstimateBackground: 

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

736 results.positive = None 

737 

738 @timeMethod 

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

740 background=None, backgroundToPhotometricRatio=None): 

741 """Detect footprints on an exposure. 

742 

743 Parameters 

744 ---------- 

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

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

747 set in-place. 

748 doSmooth : `bool`, optional 

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

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

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

752 plane. 

753 sigma : `float`, optional 

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

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

756 ``exposure``. 

757 clearMask : `bool`, optional 

758 Clear both DETECTED and DETECTED_NEGATIVE planes before running 

759 detection. 

760 expId : `dict`, optional 

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

762 RNG seed by subclasses. 

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

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

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

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

767 Image to convert photometric-flattened image to 

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

769 exposure has been photometric-flattened. 

770 

771 Returns 

772 ------- 

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

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

775 

776 ``positive`` 

777 Positive polarity footprints. 

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

779 ``negative`` 

780 Negative polarity footprints. 

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

782 ``numPos`` 

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

784 negative. (`int`) 

785 ``numNeg`` 

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

787 positive. (`int`) 

788 ``background`` 

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

790 if ``reEstimateBackground==False``. 

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

792 ``factor`` 

793 Multiplication factor applied to the configured detection 

794 threshold. (`float`) 

795 """ 

796 maskedImage = exposure.maskedImage 

797 

798 if clearMask: 

799 self.clearMask(maskedImage.getMask()) 

800 

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

802 with self.tempWideBackgroundContext(exposure): 

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

804 middle = convolveResults.middle 

805 sigma = convolveResults.sigma 

806 self.removeBadPixels(middle) 

807 

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

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

810 

811 if self.config.doTempLocalBackground: 

812 self.applyTempLocalBackground(exposure, middle, results) 

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

814 

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

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

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

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

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

820 negative=True) 

821 

822 if self.config.reEstimateBackground: 

823 self.reEstimateBackground( 

824 maskedImage, 

825 results.background, 

826 backgroundToPhotometricRatio=backgroundToPhotometricRatio, 

827 ) 

828 

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

830 

831 self.display(exposure, results, middle) 

832 

833 return results 

834 

835 def removeBadPixels(self, middle): 

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

837 

838 Parameters 

839 ---------- 

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

841 Score or maximum likelihood difference image. 

842 The image plane will be modified in place. 

843 """ 

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

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

846 middle.image.array[badPixels] = 0 

847 

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

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

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

851 

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

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

854 well-defined meaning in those cases. 

855 

856 Parameters 

857 ---------- 

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

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

860 local background-subtracted image. 

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

862 Footprints detected on the image. 

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

864 Threshold used to find footprints. 

865 negative : `bool`, optional 

866 Are we calculating for negative sources? 

867 """ 

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

869 return footprints 

870 polarity = -1 if negative else 1 

871 

872 # All incoming footprints have the same schema. 

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

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

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

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

877 

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

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

880 # significance field. 

881 newFootprints = afwDet.FootprintSet(footprints) 

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

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

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

885 new.getPeaks().clear() 

886 new.setPeakCatalog(newPeaks) 

887 

888 # Compute the significance values. 

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

890 for footprint in newFootprints.getFootprints(): 

891 footprint.updatePeakSignificance(exposure.variance, polarity) 

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

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

894 for footprint in newFootprints.getFootprints(): 

895 footprint.updatePeakSignificance(polarity*sigma) 

896 else: 

897 for footprint in newFootprints.getFootprints(): 

898 for peak in footprint.peaks: 

899 peak["significance"] = 0 

900 

901 return newFootprints 

902 

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

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

905 configuration and the statistics of the given image. 

906 

907 Parameters 

908 ---------- 

909 image : `afw.image.MaskedImage` 

910 Image to measure noise statistics from if needed. 

911 thresholdParity: `str`