Coverage for python / lsst / ip / diffim / detectAndMeasure.py: 17%

1# This file is part of ip_diffim. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

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 GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22import numpy as np 

23import requests 

24import os 

25 

26import lsst.afw.detection as afwDetection 

27import lsst.afw.image as afwImage 

28import lsst.afw.math as afwMath 

29import lsst.afw.table as afwTable 

30import lsst.daf.base as dafBase 

31import lsst.geom 

32from lsst.ip.diffim.utils import (evaluateMaskFraction, computeDifferenceImageMetrics, 

33 populate_sattle_visit_cache) 

34from lsst.meas.algorithms import SkyObjectsTask, SourceDetectionTask, SetPrimaryFlagsTask, MaskStreaksTask 

35from lsst.meas.algorithms import FindGlintTrailsTask, FindCosmicRaysConfig, findCosmicRays 

36from lsst.meas.base import ForcedMeasurementTask, ApplyApCorrTask, DetectorVisitIdGeneratorConfig 

37import lsst.meas.deblender 

38import lsst.meas.extensions.trailedSources # noqa: F401 

39import lsst.meas.extensions.shapeHSM 

40import lsst.pex.config as pexConfig 

41from lsst.pex.exceptions import InvalidParameterError, LengthError 

42import lsst.pipe.base as pipeBase 

43import lsst.utils 

44from lsst.utils.timer import timeMethod 

45 

46from . import DipoleFitTask 

47 

48__all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask", 

49 "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] 

50 

51 

52class BadSubtractionError(pipeBase.AlgorithmError): 

53 """Raised when the residuals in footprints of stars used to compute the 

54 psf-matching kernel exceeds the configured maximum. 

55 """ 

56 def __init__(self, *, ratio, threshold): 

57 msg = ("The ratio of residual power in source footprints on the" 

58 " difference image to the power in the footprints on the" 

59 f" science image was {ratio}, which exceeds the maximum" 

60 f" threshold of {threshold}") 

61 super().__init__(msg) 

62 self.ratio = ratio 

63 self.threshold = threshold 

64 

65 @property 

66 def metadata(self): 

67 return {"ratio": self.ratio, 

68 "threshold": self.threshold 

69 } 

70 

71 

72class NoDiaSourcesError(pipeBase.AlgorithmError): 

73 """Raised when there are no diaSources detected on an image difference. 

74 """ 

75 def __init__(self): 

76 msg = ("No diaSources detected!") 

77 super().__init__(msg) 

78 

79 @property 

80 def metadata(self): 

81 return {} 

82 

83 

84class TooManyCosmicRays(pipeBase.AlgorithmError): 

85 """Raised if the cosmic ray task fails with too many cosmics. 

86 

87 Parameters 

88 ---------- 

89 maxCosmicRays : `int` 

90 Maximum number of cosmic rays allowed. 

91 """ 

92 def __init__(self, maxCosmicRays, **kwargs): 

93 msg = f"Cosmic ray task found more than {maxCosmicRays} cosmics." 

94 self.msg = msg 

95 self._metadata = kwargs 

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

97 self._metadata["maxCosmicRays"] = maxCosmicRays 

98 

99 def __str__(self): 

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

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

102 

103 @property 

104 def metadata(self): 

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

106 if not (isinstance(value, int) or isinstance(value, float) or isinstance(value, str)): 

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

108 return self._metadata 

109 

110 

111class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections, 

112 dimensions=("instrument", "visit", "detector"), 

113 defaultTemplates={"coaddName": "deep", 

114 "warpTypeSuffix": "", 

115 "fakesType": ""}): 

116 science = pipeBase.connectionTypes.Input( 

117 doc="Input science exposure.", 

118 dimensions=("instrument", "visit", "detector"), 

119 storageClass="ExposureF", 

120 name="{fakesType}calexp" 

121 ) 

122 matchedTemplate = pipeBase.connectionTypes.Input( 

123 doc="Warped and PSF-matched template used to create the difference image.", 

124 dimensions=("instrument", "visit", "detector"), 

125 storageClass="ExposureF", 

126 name="{fakesType}{coaddName}Diff_matchedExp", 

127 ) 

128 difference = pipeBase.connectionTypes.Input( 

129 doc="Result of subtracting template from science.", 

130 dimensions=("instrument", "visit", "detector"), 

131 storageClass="ExposureF", 

132 name="{fakesType}{coaddName}Diff_differenceTempExp", 

133 ) 

134 kernelSources = pipeBase.connectionTypes.Input( 

135 doc="Final selection of sources used for psf matching.", 

136 dimensions=("instrument", "visit", "detector"), 

137 storageClass="SourceCatalog", 

138 name="{fakesType}{coaddName}Diff_psfMatchSources" 

139 ) 

140 outputSchema = pipeBase.connectionTypes.InitOutput( 

141 doc="Schema (as an example catalog) for output DIASource catalog.", 

142 storageClass="SourceCatalog", 

143 name="{fakesType}{coaddName}Diff_diaSrc_schema", 

144 ) 

145 diaSources = pipeBase.connectionTypes.Output( 

146 doc="Detected diaSources on the difference image.", 

147 dimensions=("instrument", "visit", "detector"), 

148 storageClass="SourceCatalog", 

149 name="{fakesType}{coaddName}Diff_diaSrc", 

150 ) 

151 subtractedMeasuredExposure = pipeBase.connectionTypes.Output( 

152 doc="Difference image with detection mask plane filled in.", 

153 dimensions=("instrument", "visit", "detector"), 

154 storageClass="ExposureF", 

155 name="{fakesType}{coaddName}Diff_differenceExp", 

156 ) 

157 differenceBackground = pipeBase.connectionTypes.Output( 

158 doc="Background model that was subtracted from the difference image.", 

159 dimensions=("instrument", "visit", "detector"), 

160 storageClass="Background", 

161 name="difference_background", 

162 ) 

163 maskedStreaks = pipeBase.connectionTypes.Output( 

164 doc='Catalog of streak fit parameters for the difference image.', 

165 storageClass="ArrowNumpyDict", 

166 dimensions=("instrument", "visit", "detector"), 

167 name="{fakesType}{coaddName}Diff_streaks", 

168 ) 

169 glintTrailInfo = pipeBase.connectionTypes.Output( 

170 doc='Dict of fit parameters for glint trails in the catalog.', 

171 storageClass="ArrowNumpyDict", 

172 dimensions=("instrument", "visit", "detector"), 

173 name="trailed_glints", 

174 ) 

175 

176 def __init__(self, *, config): 

177 super().__init__(config=config) 

178 if not (self.config.doCalculateResidualMetics): 

179 self.inputs.remove("kernelSources") 

180 if not (self.config.writeStreakInfo and self.config.doMaskStreaks): 

181 self.outputs.remove("maskedStreaks") 

182 if not (self.config.doSubtractBackground and self.config.doWriteBackground): 

183 self.outputs.remove("differenceBackground") 

184 if not (self.config.writeGlintInfo): 

185 self.outputs.remove("glintTrailInfo") 

186 

187 

188class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, 

189 pipelineConnections=DetectAndMeasureConnections): 

190 """Config for DetectAndMeasureTask 

191 """ 

192 doMerge = pexConfig.Field( 

193 dtype=bool, 

194 default=True, 

195 doc="Merge positive and negative diaSources with grow radius " 

196 "set by growFootprint" 

197 ) 

198 doForcedMeasurement = pexConfig.Field( 

199 dtype=bool, 

200 default=True, 

201 doc="Force photometer diaSource locations on PVI?" 

202 ) 

203 doAddMetrics = pexConfig.Field( 

204 dtype=bool, 

205 default=False, 

206 doc="Add columns to the source table to hold analysis metrics?" 

207 ) 

208 doSubtractBackground = pexConfig.Field( 

209 dtype=bool, 

210 doc="Subtract a background model from the image before detection?", 

211 default=True, 

212 ) 

213 doWriteBackground = pexConfig.Field( 

214 dtype=bool, 

215 doc="Persist the fitted background model?", 

216 default=False, 

217 ) 

218 doCalculateResidualMetics = pexConfig.Field( 

219 dtype=bool, 

220 doc="Calculate metrics to assess image subtraction quality for the task" 

221 "metadata?", 

222 default=True, 

223 ) 

224 subtractInitialBackground = pexConfig.ConfigurableField( 

225 target=lsst.meas.algorithms.SubtractBackgroundTask, 

226 doc="Task to perform intial background subtraction, before first detection pass.", 

227 ) 

228 subtractFinalBackground = pexConfig.ConfigurableField( 

229 target=lsst.meas.algorithms.SubtractBackgroundTask, 

230 doc="Task to perform final background subtraction, after first detection pass.", 

231 ) 

232 doFindCosmicRays = pexConfig.Field( 

233 dtype=bool, 

234 doc="Detect and mask cosmic rays on the difference image?" 

235 "CRs can be interpolated over by setting cosmicray.keepCRs=False", 

236 default=True, 

237 ) 

238 cosmicray = pexConfig.ConfigField( 

239 dtype=FindCosmicRaysConfig, 

240 doc="Options for finding and masking cosmic rays", 

241 ) 

242 detection = pexConfig.ConfigurableField( 

243 target=SourceDetectionTask, 

244 doc="Final source detection for diaSource measurement", 

245 ) 

246 streakDetection = pexConfig.ConfigurableField( 

247 target=SourceDetectionTask, 

248 doc="Separate source detection used only for streak masking", 

249 ) 

250 doDeblend = pexConfig.Field( 

251 dtype=bool, 

252 default=False, 

253 doc="Deblend DIASources after detection?" 

254 ) 

255 deblend = pexConfig.ConfigurableField( 

256 target=lsst.meas.deblender.SourceDeblendTask, 

257 doc="Task to split blended sources into their components." 

258 ) 

259 measurement = pexConfig.ConfigurableField( 

260 target=DipoleFitTask, 

261 doc="Task to measure sources on the difference image.", 

262 ) 

263 doApCorr = lsst.pex.config.Field( 

264 dtype=bool, 

265 default=True, 

266 doc="Run subtask to apply aperture corrections" 

267 ) 

268 applyApCorr = lsst.pex.config.ConfigurableField( 

269 target=ApplyApCorrTask, 

270 doc="Task to apply aperture corrections" 

271 ) 

272 forcedMeasurement = pexConfig.ConfigurableField( 

273 target=ForcedMeasurementTask, 

274 doc="Task to force photometer science image at diaSource locations.", 

275 ) 

276 growFootprint = pexConfig.Field( 

277 dtype=int, 

278 default=2, 

279 doc="Grow positive and negative footprints by this many pixels before merging" 

280 ) 

281 diaSourceMatchRadius = pexConfig.Field( 

282 dtype=float, 

283 default=0.5, 

284 doc="Match radius (in arcseconds) for DiaSource to Source association" 

285 ) 

286 doSkySources = pexConfig.Field( 

287 dtype=bool, 

288 default=False, 

289 doc="Generate sky sources?", 

290 ) 

291 skySources = pexConfig.ConfigurableField( 

292 target=SkyObjectsTask, 

293 doc="Generate sky sources", 

294 ) 

295 doMaskStreaks = pexConfig.Field( 

296 dtype=bool, 

297 default=True, 

298 doc="Turn on streak masking", 

299 ) 

300 maskStreaks = pexConfig.ConfigurableField( 

301 target=MaskStreaksTask, 

302 doc="Subtask for masking streaks. Only used if doMaskStreaks is True. " 

303 "Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.", 

304 ) 

305 streakBinFactor = pexConfig.Field( 

306 dtype=int, 

307 default=4, 

308 doc="Bin scale factor to use when rerunning detection for masking streaks. " 

309 "Only used if doMaskStreaks is True.", 

310 ) 

311 writeStreakInfo = pexConfig.Field( 

312 dtype=bool, 

313 default=False, 

314 doc="Record the parameters of any detected streaks. For LSST, this should be turned off except for " 

315 "development work." 

316 ) 

317 findGlints = pexConfig.ConfigurableField( 

318 target=FindGlintTrailsTask, 

319 doc="Subtask for finding glint trails, usually caused by satellites or debris." 

320 ) 

321 writeGlintInfo = pexConfig.Field( 

322 dtype=bool, 

323 default=True, 

324 doc="Record the parameters of any detected glint trails." 

325 ) 

326 setPrimaryFlags = pexConfig.ConfigurableField( 

327 target=SetPrimaryFlagsTask, 

328 doc="Task to add isPrimary and deblending-related flags to the catalog." 

329 ) 

330 badSourceFlags = lsst.pex.config.ListField( 

331 dtype=str, 

332 doc="Sources with any of these flags set are removed before writing the output catalog.", 

333 default=("base_PixelFlags_flag_offimage", 

334 "base_PixelFlags_flag_edge", 

335 "base_PixelFlags_flag_interpolatedCenterAll", 

336 "base_PixelFlags_flag_badCenter", 

337 "base_PixelFlags_flag_edgeCenter", 

338 "base_PixelFlags_flag_nodataCenter", 

339 "base_PixelFlags_flag_saturatedCenter", 

340 "base_PixelFlags_flag_saturated_templateCenter", 

341 "base_PixelFlags_flag_spikeCenter", 

342 ), 

343 ) 

344 clearMaskPlanes = lsst.pex.config.ListField( 

345 dtype=str, 

346 doc="Mask planes to clear before running detection.", 

347 default=("DETECTED", "DETECTED_NEGATIVE", "NOT_DEBLENDED", "STREAK"), 

348 ) 

349 doRejectBadMaskPlaneDetections = pexConfig.Field( 

350 dtype=bool, 

351 default=True, 

352 doc="Reject any peaks detected on ``badMaskPlanes`` before measurement." 

353 "These should all be rejected downstream by ``badSourceFlags``, " 

354 "but filtering earlier can save time." 

355 ) 

356 badMaskPlanes = lsst.pex.config.ListField( 

357 dtype=str, 

358 doc="Detections whose footprint peak lies on a pixel with any of these" 

359 " mask planes set will be rejected before measurement." 

360 " Any missing mask planes will be silently ignored.", 

361 default=("NO_DATA", "BAD", "SAT", "EDGE"), 

362 ) 

363 raiseOnBadSubtractionRatio = pexConfig.Field( 

364 dtype=bool, 

365 default=True, 

366 doc="Raise an error if the ratio of power in detected footprints" 

367 " on the difference image to the power in footprints on the science" 

368 " image exceeds ``badSubtractionRatioThreshold``", 

369 ) 

370 badSubtractionRatioThreshold = pexConfig.Field( 

371 dtype=float, 

372 default=0.2, 

373 doc="Maximum ratio of power in footprints on the difference image to" 

374 " the same footprints on the science image." 

375 "Only used if ``raiseOnBadSubtractionRatio`` is set", 

376 ) 

377 badSubtractionVariationThreshold = pexConfig.Field( 

378 dtype=float, 

379 default=0.4, 

380 doc="Maximum standard deviation of the ratio of power in footprints on" 

381 " the difference image to the same footprints on the science image." 

382 "Only used if ``raiseOnBadSubtractionRatio`` is set", 

383 ) 

384 raiseOnNoDiaSources = pexConfig.Field( 

385 dtype=bool, 

386 default=True, 

387 doc="Raise an algorithm error if no diaSources are detected.", 

388 ) 

389 run_sattle = pexConfig.Field( 

390 dtype=bool, 

391 default=False, 

392 doc="If true, dia source bounding boxes will be sent for verification" 

393 "to the sattle service." 

394 ) 

395 sattle_historical = pexConfig.Field( 

396 dtype=bool, 

397 default=False, 

398 doc="If re-running a pipeline that requires sattle, this should be set " 

399 "to True. This will populate sattle's cache with the historic data " 

400 "closest in time to the exposure." 

401 ) 

402 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

403 

404 def setDefaults(self): 

405 # Background subtraction 

406 # Use a small binsize for the first pass to reduce detections on glints 

407 # and extended structures. Should not affect the detectability of 

408 # faint diaSources 

409 self.subtractInitialBackground.binSize = 8 

410 self.subtractInitialBackground.useApprox = False 

411 self.subtractInitialBackground.statisticsProperty = "MEDIAN" 

412 self.subtractInitialBackground.doFilterSuperPixels = True 

413 self.subtractInitialBackground.ignoredPixelMask = ["BAD", 

414 "EDGE", 

415 "DETECTED", 

416 "DETECTED_NEGATIVE", 

417 "NO_DATA", 

418 ] 

419 # Use a larger binsize for the final background subtraction, to reduce 

420 # over-subtraction of bright objects. 

421 self.subtractFinalBackground.binSize = 40 

422 self.subtractFinalBackground.useApprox = False 

423 self.subtractFinalBackground.statisticsProperty = "MEDIAN" 

424 self.subtractFinalBackground.doFilterSuperPixels = True 

425 self.subtractFinalBackground.ignoredPixelMask = ["BAD", 

426 "EDGE", 

427 "DETECTED", 

428 "DETECTED_NEGATIVE", 

429 "NO_DATA", 

430 ] 

431 # Cosmic ray detection 

432 self.cosmicray.keepCRs = True # do not interpolate over detected CRs 

433 # DiaSource Detection 

434 self.detection.thresholdPolarity = "both" 

435 self.detection.thresholdValue = 5.0 

436 self.detection.reEstimateBackground = False 

437 self.detection.thresholdType = "pixel_stdev" 

438 self.detection.excludeMaskPlanes = [] 

439 

440 # Copy configs for binned streak detection from the base detection task 

441 self.streakDetection.thresholdType = self.detection.thresholdType 

442 self.streakDetection.reEstimateBackground = False 

443 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes 

444 self.streakDetection.thresholdValue = self.detection.thresholdValue 

445 # Only detect positive streaks 

446 self.streakDetection.thresholdPolarity = "positive" 

447 # Do not grow detected mask for streaks 

448 self.streakDetection.nSigmaToGrow = 0 

449 # Set the streak mask along the entire fit line, not only where the 

450 # detected mask is set. 

451 self.maskStreaks.onlyMaskDetected = False 

452 # Restrict streak masking from growing too large 

453 self.maskStreaks.maxStreakWidth = 100 

454 # Restrict the number of iterations allowed for fitting streaks 

455 # When the fit is good it should solve quickly, and exit a bad fit quickly 

456 self.maskStreaks.maxFitIter = 10 

457 # Only mask to 2 sigma in width 

458 self.maskStreaks.nSigmaMask = 2 

459 # Threshold for including streaks after the Hough Transform. 

460 # A lower value will detect more features that are less linear. 

461 self.maskStreaks.absMinimumKernelHeight = 2 

462 

463 self.measurement.plugins.names |= ["ext_trailedSources_Naive", 

464 "base_LocalPhotoCalib", 

465 "base_LocalWcs", 

466 "ext_shapeHSM_HsmSourceMoments", 

467 "ext_shapeHSM_HsmPsfMoments", 

468 "base_ClassificationSizeExtendedness", 

469 ] 

470 self.measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments" 

471 self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments" 

472 self.measurement.plugins["base_SdssCentroid"].maxDistToPeak = 5.0 

473 self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"] 

474 self.forcedMeasurement.copyColumns = { 

475 "id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"} 

476 self.forcedMeasurement.slots.centroid = "base_TransformedCentroid" 

477 self.forcedMeasurement.slots.shape = None 

478 

479 # Keep track of which footprints contain streaks 

480 self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [ 

481 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SATURATED_TEMPLATE", "SPIKE"] 

482 self.measurement.plugins["base_PixelFlags"].masksFpCenter = [ 

483 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SATURATED_TEMPLATE", "SPIKE"] 

484 self.skySources.avoidMask = ["DETECTED", "DETECTED_NEGATIVE", "BAD", "NO_DATA", "EDGE"] 

485 

486 def validate(self): 

487 super().validate() 

488 

489 if self.run_sattle: 

490 if not os.getenv("SATTLE_URI_BASE"): 

491 raise pexConfig.FieldValidationError(DetectAndMeasureConfig.run_sattle, self, 

492 "Sattle requested but SATTLE_URI_BASE " 

493 "environment variable not set.") 

494 

495 

496class DetectAndMeasureTask(lsst.pipe.base.PipelineTask): 

497 """Detect and measure sources on a difference image. 

498 """ 

499 ConfigClass = DetectAndMeasureConfig 

500 _DefaultName = "detectAndMeasure" 

501 

502 def __init__(self, **kwargs): 

503 super().__init__(**kwargs) 

504 self.schema = afwTable.SourceTable.makeMinimalSchema() 

505 

506 self.algMetadata = dafBase.PropertyList() 

507 if self.config.doSubtractBackground: 

508 self.makeSubtask("subtractInitialBackground") 

509 self.makeSubtask("subtractFinalBackground") 

510 self.makeSubtask("detection", schema=self.schema) 

511 if self.config.doDeblend: 

512 self.makeSubtask("deblend", schema=self.schema) 

513 self.makeSubtask("setPrimaryFlags", schema=self.schema, isSingleFrame=True) 

514 self.makeSubtask("measurement", schema=self.schema, 

515 algMetadata=self.algMetadata) 

516 if self.config.doApCorr: 

517 self.makeSubtask("applyApCorr", schema=self.measurement.schema) 

518 if self.config.doForcedMeasurement: 

519 self.schema.addField( 

520 "ip_diffim_forced_PsfFlux_instFlux", "D", 

521 "Forced PSF flux measured on the direct image.", 

522 units="count") 

523 self.schema.addField( 

524 "ip_diffim_forced_PsfFlux_instFluxErr", "D", 

525 "Forced PSF flux error measured on the direct image.", 

526 units="count") 

527 self.schema.addField( 

528 "ip_diffim_forced_PsfFlux_area", "F", 

529 "Forced PSF flux effective area of PSF.", 

530 units="pixel") 

531 self.schema.addField( 

532 "ip_diffim_forced_PsfFlux_flag", "Flag", 

533 "Forced PSF flux general failure flag.") 

534 self.schema.addField( 

535 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag", 

536 "Forced PSF flux not enough non-rejected pixels in data to attempt the fit.") 

537 self.schema.addField( 

538 "ip_diffim_forced_PsfFlux_flag_edge", "Flag", 

539 "Forced PSF flux object was too close to the edge of the image to use the full PSF model.") 

540 self.schema.addField( 

541 "ip_diffim_forced_template_PsfFlux_instFlux", "D", 

542 "Forced PSF flux measured on the template image.", 

543 units="count") 

544 self.schema.addField( 

545 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D", 

546 "Forced PSF flux error measured on the template image.", 

547 units="count") 

548 self.schema.addField( 

549 "ip_diffim_forced_template_PsfFlux_area", "F", 

550 "Forced template PSF flux effective area of PSF.", 

551 units="pixel") 

552 self.schema.addField( 

553 "ip_diffim_forced_template_PsfFlux_flag", "Flag", 

554 "Forced template PSF flux general failure flag.") 

555 self.schema.addField( 

556 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag", 

557 "Forced template PSF flux not enough non-rejected pixels in data to attempt the fit.") 

558 self.schema.addField( 

559 "ip_diffim_forced_template_PsfFlux_flag_edge", "Flag", 

560 """Forced template PSF flux object was too close to the edge of the image """ 

561 """to use the full PSF model.""") 

562 self.makeSubtask("forcedMeasurement", refSchema=self.schema) 

563 

564 self.schema.addField("refMatchId", "L", "unique id of reference catalog match") 

565 self.schema.addField("srcMatchId", "L", "unique id of source match") 

566 # Create the sky source task for use by metrics, 

567 # even if sky sources are not added to the diaSource catalog 

568 self.makeSubtask("skySources", schema=self.schema) 

569 if self.config.doMaskStreaks: 

570 self.makeSubtask("maskStreaks") 

571 self.makeSubtask("streakDetection") 

572 self.makeSubtask("findGlints") 

573 self.schema.addField("glint_trail", "Flag", "DiaSource is part of a glint trail.") 

574 self.schema.addField("reliability", type="F", doc="Reliability score of the DiaSource") 

575 

576 # To get the "merge_*" fields in the schema; have to re-initialize 

577 # this later, once we have a peak schema post-detection. 

578 lsst.afw.detection.FootprintMergeList(self.schema, ["positive", "negative"]) 

579 

580 # Check that the schema and config are consistent 

581 for flag in self.config.badSourceFlags: 

582 if flag not in self.schema: 

583 raise pipeBase.InvalidQuantumError("Field %s not in schema" % flag) 

584 

585 # initialize InitOutputs 

586 self.outputSchema = afwTable.SourceCatalog(self.schema) 

587 self.outputSchema.getTable().setMetadata(self.algMetadata) 

588 

589 def runQuantum(self, butlerQC: pipeBase.QuantumContext, 

590 inputRefs: pipeBase.InputQuantizedConnection, 

591 outputRefs: pipeBase.OutputQuantizedConnection): 

592 inputs = butlerQC.get(inputRefs) 

593 idGenerator = self.config.idGenerator.apply(butlerQC.quantum.dataId) 

594 idFactory = idGenerator.make_table_id_factory() 

595 # Specify the fields that `annotate` needs below, to ensure they 

596 # exist, even as None. 

597 measurementResults = pipeBase.Struct( 

598 subtractedMeasuredExposure=None, 

599 diaSources=None, 

600 maskedStreaks=None, 

601 differenceBackground=None, 

602 ) 

603 try: 

604 self.run(**inputs, idFactory=idFactory, measurementResults=measurementResults) 

605 except pipeBase.AlgorithmError as e: 

606 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

607 e, 

608 self, 

609 measurementResults.subtractedMeasuredExposure, 

610 measurementResults.diaSources, 

611 measurementResults.maskedStreaks, 

612 log=self.log 

613 ) 

614 butlerQC.put(measurementResults, outputRefs) 

615 raise error from e 

616 butlerQC.put(measurementResults, outputRefs) 

617 

618 @timeMethod 

619 def run(self, science, matchedTemplate, difference, kernelSources=None, 

620 idFactory=None, measurementResults=None): 

621 """Detect and measure sources on a difference image. 

622 

623 The difference image will be convolved with a gaussian approximation of 

624 the PSF to form a maximum likelihood image for detection. 

625 Close positive and negative detections will optionally be merged into 

626 dipole diaSources. 

627 Sky sources, or forced detections in background regions, will optionally 

628 be added, and the configured measurement algorithm will be run on all 

629 detections. 

630 

631 Parameters 

632 ---------- 

633 science : `lsst.afw.image.ExposureF` 

634 Science exposure that the template was subtracted from. 

635 matchedTemplate : `lsst.afw.image.ExposureF` 

636 Warped and PSF-matched template that was used produce the 

637 difference image. 

638 difference : `lsst.afw.image.ExposureF` 

639 Result of subtracting template from the science image. 

640 kernelSources : `lsst.afw.table.SourceCatalog`, optional 

641 Final selection of sources that was used for psf matching. 

642 idFactory : `lsst.afw.table.IdFactory`, optional 

643 Generator object used to assign ids to detected sources in the 

644 difference image. Ids from this generator are not set until after 

645 deblending and merging positive/negative peaks. 

646 measurementResults : `lsst.pipe.base.Struct`, optional 

647 Result struct that is modified to allow saving of partial outputs 

648 for some failure conditions. If the task completes successfully, 

649 this is also returned. 

650 

651 Returns 

652 ------- 

653 measurementResults : `lsst.pipe.base.Struct` 

654 

655 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF` 

656 Subtracted exposure with detection mask applied. 

657 ``diaSources`` : `lsst.afw.table.SourceCatalog` 

658 The catalog of detected sources. 

659 ``differenceBackground`` : `lsst.afw.math.BackgroundList` 

660 Background that was subtracted from the difference image. 

661 """ 

662 if measurementResults is None: 

663 measurementResults = pipeBase.Struct() 

664 if idFactory is None: 

665 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory() 

666 

667 # Check image properties and clear detection mask planes. 

668 self._prepareInputs(difference) 

669 if self.config.doSubtractBackground: 

670 background = self._fitAndSubtractBackground(differenceExposure=difference) 

671 else: 

672 background = afwMath.BackgroundList() 

673 

674 if self.config.doFindCosmicRays: 

675 self.findAndMaskCosmicRays(difference) 

676 

677 # Don't use the idFactory until after deblend+merge, so that we aren't 

678 # generating ids that just get thrown away (footprint merge doesn't 

679 # know about past ids). 

680 table = afwTable.SourceTable.make(self.schema) 

681 results = self.detection.run( 

682 table=table, 

683 exposure=difference, 

684 doSmooth=True, 

685 background=background, 

686 clearMask=True, 

687 ) 

688 measurementResults.differenceBackground = background 

689 

690 if self.config.doDeblend: 

691 sources, positives, negatives = self._deblend(difference, 

692 results.positive, 

693 results.negative) 

694 

695 else: 

696 positives = afwTable.SourceCatalog(self.schema) 

697 results.positive.makeSources(positives) 

698 negatives = afwTable.SourceCatalog(self.schema) 

699 results.negative.makeSources(negatives) 

700 sources = results.sources 

701 

702 self.processResults(science, matchedTemplate, difference, 

703 sources, idFactory, kernelSources, 

704 positives=positives, 

705 negatives=negatives, 

706 measurementResults=measurementResults) 

707 return measurementResults 

708 

709 def _prepareInputs(self, difference): 

710 """Ensure that we start with an empty detection and deblended mask. 

711 

712 Parameters 

713 ---------- 

714 difference : `lsst.afw.image.ExposureF` 

715 The difference image that will be used for detecting diaSources. 

716 The mask plane will be modified in place. 

717 

718 Raises 

719 ------ 

720 lsst.pipe.base.UpstreamFailureNoWorkFound 

721 If the PSF is not usable for measurement. 

722 """ 

723 # Check that we have a valid PSF now before we do more work 

724 sigma = difference.psf.computeShape(difference.psf.getAveragePosition()).getDeterminantRadius() 

725 if np.isnan(sigma): 

726 raise pipeBase.UpstreamFailureNoWorkFound("Invalid PSF detected! PSF width evaluates to NaN.") 

727 # Ensure that we start with an empty detection and deblended mask. 

728 mask = difference.mask 

729 for mp in self.config.clearMaskPlanes: 

730 if mp not in mask.getMaskPlaneDict(): 

731 mask.addMaskPlane(mp) 

732 mask &= ~mask.getPlaneBitMask(self.config.clearMaskPlanes) 

733 

734 def _fitAndSubtractBackground(self, differenceExposure, scoreExposure=None): 

735 """Fit the final background to ``differenceExposure``. 

736 

737 A rough background is first fit on the difference and an internal 

738 first-pass detection is run on the (background-subtracted) clone; the 

739 resulting detection mask is then transferred to ``differenceExposure`` 

740 so the final background fit can mask out real sources. This two-pass 

741 scheme is what suppresses spurious detections from glints and 

742 extended structures: the small-binsize initial fit aggressively 

743 over-subtracts those features, and the resulting peak mask lets the 

744 large-binsize final fit produce a clean background. When 

745 ``scoreExposure`` is supplied (score-image variant), the final 

746 background is also subtracted from its image so the caller's final 

747 detection runs on a properly background-subtracted score. 

748 

749 The first-pass detection's `Struct` is discarded; only the mask it 

750 leaves on the clone is used downstream. 

751 

752 Parameters 

753 ---------- 

754 differenceExposure : `lsst.afw.image.ExposureF` 

755 Difference image. The final background is subtracted from its 

756 image in place but its mask is not modified. 

757 scoreExposure : `lsst.afw.image.ExposureF`, optional 

758 Maximum-likelihood score image. When supplied, the final 

759 background is also subtracted from its image. 

760 

761 Returns 

762 ------- 

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

764 The fit final background. 

765 """ 

766 if scoreExposure is None: 

767 detectionExposure = differenceExposure.clone() 

768 background = self.subtractInitialBackground.run(detectionExposure).background 

769 doSmooth = True 

770 else: 

771 # Use a clone of differenceExposure because the background is 

772 # subtracted in place. 

773 background = self.subtractInitialBackground.run(differenceExposure.clone()).background 

774 detectionExposure = scoreExposure.clone() 

775 detectionExposure.image -= background.getImage() 

776 doSmooth = False 

777 table = afwTable.SourceTable.make(self.schema) 

778 self.detection.run( 

779 table=table, 

780 exposure=detectionExposure, 

781 doSmooth=doSmooth, 

782 background=background, 

783 clearMask=True, 

784 ) 

785 # Use the temporary detection mask for the final background subtraction. 

786 # The detection mask planes will be cleared before the final detection 

787 # step, so it is OK if they get set for differenceExposure. 

788 differenceExposure.setMask(detectionExposure.mask) 

789 background = self.subtractFinalBackground.run(differenceExposure).background 

790 if scoreExposure is not None: 

791 # The preconvolution kernel is normalized to 1, so the same 

792 # background level applies to the difference and score images. 

793 scoreExposure.image -= background.getImage() 

794 return background 

795 

796 def processResults(self, science, matchedTemplate, difference, sources, idFactory, 

797 kernelSources=None, positives=None, negatives=None, measurementResults=None): 

798 """Measure and process the results of source detection. 

799 

800 Parameters 

801 ---------- 

802 science : `lsst.afw.image.ExposureF` 

803 Science exposure that the template was subtracted from. 

804 matchedTemplate : `lsst.afw.image.ExposureF` 

805 Warped and PSF-matched template that was used produce the 

806 difference image. 

807 difference : `lsst.afw.image.ExposureF` 

808 Result of subtracting template from the science image. 

809 sources : `lsst.afw.table.SourceCatalog` 

810 Detected sources on the difference exposure. 

811 idFactory : `lsst.afw.table.IdFactory` 

812 Generator object used to assign ids to detected sources in the 

813 difference image. 

814 kernelSources : `lsst.afw.table.SourceCatalog`, optional 

815 Final selection of sources that was used for psf matching. 

816 positives : `lsst.afw.table.SourceCatalog`, optional 

817 Positive polarity footprints. 

818 negatives : `lsst.afw.table.SourceCatalog`, optional 

819 Negative polarity footprints. 

820 measurementResults : `lsst.pipe.base.Struct`, optional 

821 Result struct that is modified to allow saving of partial outputs 

822 for some failure conditions. If the task completes successfully, 

823 this is also returned. 

824 

825 Returns 

826 ------- 

827 measurementResults : `lsst.pipe.base.Struct` 

828 

829 ``subtractedMeasuredExposure`` : `lsst.afw.image.ExposureF` 

830 Subtracted exposure with detection mask applied. 

831 ``diaSources`` : `lsst.afw.table.SourceCatalog` 

832 The catalog of detected sources. 

833 """ 

834 if measurementResults is None: 

835 measurementResults = pipeBase.Struct() 

836 self.metadata["nUnmergedDiaSources"] = len(sources) 

837 if self.config.doMerge: 

838 # preserve peak schema, if there are any footprints 

839 if len(positives) > 0: 

840 peakSchema = positives[0].getFootprint().peaks.schema 

841 elif len(negatives) > 0: 

842 peakSchema = negatives[0].getFootprint().peaks.schema 

843 else: 

844 peakSchema = afwDetection.PeakTable.makeMinimalSchema() 

845 mergeList = afwDetection.FootprintMergeList(self.schema, 

846 ["positive", "negative"], peakSchema) 

847 initialDiaSources = afwTable.SourceCatalog(self.schema) 

848 # Start with positive, as FootprintMergeList will self-merge the 

849 # subsequent added catalogs, and we want to try to preserve 

850 # deblended positive sources. 

851 mergeList.addCatalog(initialDiaSources.table, positives, "positive", minNewPeakDist=0) 

852 mergeList.addCatalog(initialDiaSources.table, negatives, "negative", minNewPeakDist=0) 

853 mergeList.getFinalSources(initialDiaSources) 

854 # Flag as negative those sources that *only* came from the negative 

855 # footprint set. 

856 initialDiaSources["is_negative"] = initialDiaSources["merge_footprint_negative"] & \ 

857 ~initialDiaSources["merge_footprint_positive"] 

858 self.log.info("Merging detections into %d sources", len(initialDiaSources)) 

859 else: 

860 initialDiaSources = sources 

861 

862 # Assign source ids at the end: deblend/merge mean that we don't keep 

863 # track of parents and children, we only care about the final ids. 

864 for source in initialDiaSources: 

865 source.setId(idFactory())