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

335 "base_PixelFlags_flag_badCenter", 

336 "base_PixelFlags_flag_edgeCenter", 

337 "base_PixelFlags_flag_nodataCenter", 

338 "base_PixelFlags_flag_saturatedCenter", 

339 "base_PixelFlags_flag_saturated_templateCenter", 

340 "base_PixelFlags_flag_spikeCenter", 

341 ), 

342 ) 

343 clearMaskPlanes = lsst.pex.config.ListField( 

344 dtype=str, 

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

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

347 ) 

348 doRejectBadMaskPlaneDetections = pexConfig.Field( 

349 dtype=bool, 

350 default=True, 

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

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

353 "but filtering earlier can save time." 

354 ) 

355 badMaskPlanes = lsst.pex.config.ListField( 

356 dtype=str, 

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

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

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

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

361 ) 

362 raiseOnBadSubtractionRatio = pexConfig.Field( 

363 dtype=bool, 

364 default=True, 

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

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

367 " image exceeds ``badSubtractionRatioThreshold``", 

368 ) 

369 badSubtractionRatioThreshold = pexConfig.Field( 

370 dtype=float, 

371 default=0.2, 

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

373 " the same footprints on the science image." 

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

375 ) 

376 badSubtractionVariationThreshold = pexConfig.Field( 

377 dtype=float, 

378 default=0.4, 

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

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

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

382 ) 

383 raiseOnNoDiaSources = pexConfig.Field( 

384 dtype=bool, 

385 default=True, 

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

387 ) 

388 run_sattle = pexConfig.Field( 

389 dtype=bool, 

390 default=False, 

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

392 "to the sattle service." 

393 ) 

394 sattle_historical = pexConfig.Field( 

395 dtype=bool, 

396 default=False, 

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

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

399 "closest in time to the exposure." 

400 ) 

401 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

402 

403 def setDefaults(self): 

404 # Background subtraction 

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

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

407 # faint diaSources 

408 self.subtractInitialBackground.binSize = 8 

409 self.subtractInitialBackground.useApprox = False 

410 self.subtractInitialBackground.statisticsProperty = "MEDIAN" 

411 self.subtractInitialBackground.doFilterSuperPixels = True 

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

413 "EDGE", 

414 "DETECTED", 

415 "DETECTED_NEGATIVE", 

416 "NO_DATA", 

417 ] 

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

419 # over-subtraction of bright objects. 

420 self.subtractFinalBackground.binSize = 40 

421 self.subtractFinalBackground.useApprox = False 

422 self.subtractFinalBackground.statisticsProperty = "MEDIAN" 

423 self.subtractFinalBackground.doFilterSuperPixels = True 

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

425 "EDGE", 

426 "DETECTED", 

427 "DETECTED_NEGATIVE", 

428 "NO_DATA", 

429 ] 

430 # Cosmic ray detection 

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

432 # DiaSource Detection 

433 self.detection.thresholdPolarity = "both" 

434 self.detection.thresholdValue = 5.0 

435 self.detection.reEstimateBackground = False 

436 self.detection.thresholdType = "pixel_stdev" 

437 self.detection.excludeMaskPlanes = [] 

438 

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

440 self.streakDetection.thresholdType = self.detection.thresholdType 

441 self.streakDetection.reEstimateBackground = False 

442 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes 

443 self.streakDetection.thresholdValue = self.detection.thresholdValue 

444 # Only detect positive streaks 

445 self.streakDetection.thresholdPolarity = "positive" 

446 # Do not grow detected mask for streaks 

447 self.streakDetection.nSigmaToGrow = 0 

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

449 # detected mask is set. 

450 self.maskStreaks.onlyMaskDetected = False 

451 # Restrict streak masking from growing too large 

452 self.maskStreaks.maxStreakWidth = 100 

453 # Restrict the number of iterations allowed for fitting streaks 

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

455 self.maskStreaks.maxFitIter = 10 

456 # Only mask to 2 sigma in width 

457 self.maskStreaks.nSigmaMask = 2 

458 # Threshold for including streaks after the Hough Transform. 

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

460 self.maskStreaks.absMinimumKernelHeight = 2 

461 

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

463 "base_LocalPhotoCalib", 

464 "base_LocalWcs", 

465 "ext_shapeHSM_HsmSourceMoments", 

466 "ext_shapeHSM_HsmPsfMoments", 

467 "base_ClassificationSizeExtendedness", 

468 ] 

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

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

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

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

473 self.forcedMeasurement.copyColumns = { 

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

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

476 self.forcedMeasurement.slots.shape = None 

477 

478 # Keep track of which footprints contain streaks 

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

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

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

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

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

484 

485 def validate(self): 

486 super().validate() 

487 

488 if self.run_sattle: 

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

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

491 "Sattle requested but SATTLE_URI_BASE " 

492 "environment variable not set.") 

493 

494 

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

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

497 """ 

498 ConfigClass = DetectAndMeasureConfig 

499 _DefaultName = "detectAndMeasure" 

500 

501 def __init__(self, **kwargs): 

502 super().__init__(**kwargs) 

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

504 

505 self.algMetadata = dafBase.PropertyList() 

506 if self.config.doSubtractBackground: 

507 self.makeSubtask("subtractInitialBackground") 

508 self.makeSubtask("subtractFinalBackground") 

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

510 if self.config.doDeblend: 

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

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

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

514 algMetadata=self.algMetadata) 

515 if self.config.doApCorr: 

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

517 if self.config.doForcedMeasurement: 

518 self.schema.addField( 

519 "ip_diffim_forced_PsfFlux_instFlux", "D", 

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

521 units="count") 

522 self.schema.addField( 

523 "ip_diffim_forced_PsfFlux_instFluxErr", "D", 

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

525 units="count") 

526 self.schema.addField( 

527 "ip_diffim_forced_PsfFlux_area", "F", 

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

529 units="pixel") 

530 self.schema.addField( 

531 "ip_diffim_forced_PsfFlux_flag", "Flag", 

532 "Forced PSF flux general failure flag.") 

533 self.schema.addField( 

534 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag", 

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

536 self.schema.addField( 

537 "ip_diffim_forced_PsfFlux_flag_edge", "Flag", 

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

539 self.schema.addField( 

540 "ip_diffim_forced_template_PsfFlux_instFlux", "D", 

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

542 units="count") 

543 self.schema.addField( 

544 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D", 

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

546 units="count") 

547 self.schema.addField( 

548 "ip_diffim_forced_template_PsfFlux_area", "F", 

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

550 units="pixel") 

551 self.schema.addField( 

552 "ip_diffim_forced_template_PsfFlux_flag", "Flag", 

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

554 self.schema.addField( 

555 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag", 

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

557 self.schema.addField( 

558 "ip_diffim_forced_template_PsfFlux_flag_edge", "Flag", 

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

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

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

562 

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

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

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

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

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

568 if self.config.doMaskStreaks: 

569 self.makeSubtask("maskStreaks") 

570 self.makeSubtask("streakDetection") 

571 self.makeSubtask("findGlints") 

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

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

574 

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

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

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

578 

579 # Check that the schema and config are consistent 

580 for flag in self.config.badSourceFlags: 

581 if flag not in self.schema: 

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

583 

584 # initialize InitOutputs 

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

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

587 

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

589 inputRefs: pipeBase.InputQuantizedConnection, 

590 outputRefs: pipeBase.OutputQuantizedConnection): 

591 inputs = butlerQC.get(inputRefs) 

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

593 idFactory = idGenerator.make_table_id_factory() 

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

595 # exist, even as None. 

596 measurementResults = pipeBase.Struct( 

597 subtractedMeasuredExposure=None, 

598 diaSources=None, 

599 maskedStreaks=None, 

600 differenceBackground=None, 

601 ) 

602 try: 

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

604 except pipeBase.AlgorithmError as e: 

605 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

606 e, 

607 self, 

608 measurementResults.subtractedMeasuredExposure, 

609 measurementResults.diaSources, 

610 measurementResults.maskedStreaks, 

611 log=self.log 

612 ) 

613 butlerQC.put(measurementResults, outputRefs) 

614 raise error from e 

615 butlerQC.put(measurementResults, outputRefs) 

616 

617 @timeMethod 

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

619 idFactory=None, measurementResults=None): 

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

621 

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

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

624 Close positive and negative detections will optionally be merged into 

625 dipole diaSources. 

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

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

628 detections. 

629 

630 Parameters 

631 ---------- 

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

633 Science exposure that the template was subtracted from. 

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

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

636 difference image. 

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

638 Result of subtracting template from the science image. 

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

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

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

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

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

644 deblending and merging positive/negative peaks. 

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

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

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

648 this is also returned. 

649 

650 Returns 

651 ------- 

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

653 

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

655 Subtracted exposure with detection mask applied. 

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

657 The catalog of detected sources. 

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

659 Background that was subtracted from the difference image. 

660 """ 

661 if measurementResults is None: 

662 measurementResults = pipeBase.Struct() 

663 if idFactory is None: 

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

665 

666 # Check image properties and clear detection mask planes. 

667 self._prepareInputs(difference) 

668 if self.config.doSubtractBackground: 

669 background = self._fitAndSubtractBackground(differenceExposure=difference) 

670 else: 

671 background = afwMath.BackgroundList() 

672 

673 if self.config.doFindCosmicRays: 

674 self.findAndMaskCosmicRays(difference) 

675 

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

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

678 # know about past ids). 

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

680 results = self.detection.run( 

681 table=table, 

682 exposure=difference, 

683 doSmooth=True, 

684 background=background, 

685 clearMask=True, 

686 ) 

687 measurementResults.differenceBackground = background 

688 

689 if self.config.doDeblend: 

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

691 results.positive, 

692 results.negative) 

693 

694 else: 

695 positives = afwTable.SourceCatalog(self.schema) 

696 results.positive.makeSources(positives) 

697 negatives = afwTable.SourceCatalog(self.schema) 

698 results.negative.makeSources(negatives) 

699 sources = results.sources 

700 

701 self.processResults(science, matchedTemplate, difference, 

702 sources, idFactory, kernelSources, 

703 positives=positives, 

704 negatives=negatives, 

705 measurementResults=measurementResults) 

706 return measurementResults 

707 

708 def _prepareInputs(self, difference): 

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

710 

711 Parameters 

712 ---------- 

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

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

715 The mask plane will be modified in place. 

716 

717 Raises 

718 ------ 

719 lsst.pipe.base.UpstreamFailureNoWorkFound 

720 If the PSF is not usable for measurement. 

721 """ 

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

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

724 if np.isnan(sigma): 

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

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

727 mask = difference.mask 

728 for mp in self.config.clearMaskPlanes: 

729 if mp not in mask.getMaskPlaneDict(): 

730 mask.addMaskPlane(mp) 

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

732 

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

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

735 

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

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

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

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

740 scheme is what suppresses spurious detections from glints and 

741 extended structures: the small-binsize initial fit aggressively 

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

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

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

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

746 detection runs on a properly background-subtracted score. 

747 

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

749 leaves on the clone is used downstream. 

750 

751 Parameters 

752 ---------- 

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

754 Difference image. The final background is subtracted from its 

755 image in place but its mask is not modified. 

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

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

758 background is also subtracted from its image. 

759 

760 Returns 

761 ------- 

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

763 The fit final background. 

764 """ 

765 if scoreExposure is None: 

766 detectionExposure = differenceExposure.clone() 

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

768 doSmooth = True 

769 else: 

770 # Use a clone of differenceExposure because the background is 

771 # subtracted in place. 

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

773 detectionExposure = scoreExposure.clone() 

774 detectionExposure.image -= background.getImage() 

775 doSmooth = False 

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

777 self.detection.run( 

778 table=table, 

779 exposure=detectionExposure, 

780 doSmooth=doSmooth, 

781 background=background, 

782 clearMask=True, 

783 ) 

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

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

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

787 detectedBit = differenceExposure.mask.getPlaneBitMask(["DETECTED"]) 

788 detectedPix = detectionExposure.mask.array & detectedBit > 0 

789 detectedNegativeBit = differenceExposure.mask.getPlaneBitMask(["DETECTED_NEGATIVE"]) 

790 detectedNegativePix = detectionExposure.mask.array & detectedNegativeBit > 0 

791 differenceExposure.mask.array[detectedPix] |= detectedBit 

792 differenceExposure.mask.array[detectedNegativePix] |= detectedNegativeBit 

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

794 if scoreExposure is not None: 

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

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

797 scoreExposure.image -= background.getImage() 

798 return background 

799 

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

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

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

803 

804 Parameters 

805 ---------- 

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

807 Science exposure that the template was subtracted from. 

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

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

810 difference image. 

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

812 Result of subtracting template from the science image. 

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

814 Detected sources on the difference exposure. 

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

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

817 difference image. 

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

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

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

821 Positive polarity footprints. 

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

823 Negative polarity footprints. 

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

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

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

827 this is also returned. 

828 

829 Returns 

830 ------- 

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

832 

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

834 Subtracted exposure with detection mask applied. 

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

836 The catalog of detected sources. 

837 """ 

838 if measurementResults is None: 

839 measurementResults = pipeBase.Struct() 

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

841 if self.config.doMerge: 

842 # preserve peak schema, if there are any footprints 

843 if len(positives) > 0: 

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

845 elif len(negatives) > 0: 

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

847 else: 

848 peakSchema = afwDetection.PeakTable.makeMinimalSchema() 

849 mergeList = afwDetection.FootprintMergeList(self.schema, 

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

851 initialDiaSources = afwTable.SourceCatalog(self.schema) 

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

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

854 # deblended positive sources. 

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

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

857 mergeList.getFinalSources(initialDiaSources) 

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

859 # footprint set. 

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

861 ~initialDiaSources["merge_footprint_positive"]