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

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 dataclasses 

23import numpy as np 

24import requests 

25import os 

26 

27 

28import lsst.afw.detection as afwDetection 

29import lsst.afw.image as afwImage 

30import lsst.afw.math as afwMath 

31import lsst.afw.table as afwTable 

32import lsst.daf.base as dafBase 

33import lsst.geom 

34from lsst.images import Mask, get_legacy_difference_image_mask_planes 

35from lsst.images.fields import field_from_legacy_background 

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

37 populate_sattle_visit_cache, 

38 get_difference_image_provenance) 

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

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

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

42import lsst.meas.deblender 

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

44import lsst.meas.extensions.shapeHSM 

45import lsst.pex.config as pexConfig 

46from lsst.pex.exceptions import InvalidParameterError, LengthError 

47import lsst.pipe.base as pipeBase 

48import lsst.utils 

49from lsst.utils.timer import timeMethod 

50 

51from . import DipoleFitTask 

52 

53__all__ = ["DetectAndMeasureConfig", "DetectAndMeasureTask", 

54 "DetectAndMeasureScoreConfig", "DetectAndMeasureScoreTask"] 

55 

56 

57class BadSubtractionError(pipeBase.AlgorithmError): 

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

59 psf-matching kernel exceeds the configured maximum. 

60 """ 

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

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

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

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

65 f" threshold of {threshold}") 

66 super().__init__(msg) 

67 self.ratio = ratio 

68 self.threshold = threshold 

69 

70 @property 

71 def metadata(self): 

72 return {"ratio": self.ratio, 

73 "threshold": self.threshold 

74 } 

75 

76 

77class NoDiaSourcesError(pipeBase.AlgorithmError): 

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

79 """ 

80 def __init__(self): 

81 msg = ("No diaSources detected!") 

82 super().__init__(msg) 

83 

84 @property 

85 def metadata(self): 

86 return {} 

87 

88 

89class TooManyCosmicRays(pipeBase.AlgorithmError): 

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

91 

92 Parameters 

93 ---------- 

94 maxCosmicRays : `int` 

95 Maximum number of cosmic rays allowed. 

96 """ 

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

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

99 self.msg = msg 

100 self._metadata = kwargs 

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

102 self._metadata["maxCosmicRays"] = maxCosmicRays 

103 

104 def __str__(self): 

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

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

107 

108 @property 

109 def metadata(self): 

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

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

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

113 return self._metadata 

114 

115 

116class DetectAndMeasureConnections(pipeBase.PipelineTaskConnections, 

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

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

119 "warpTypeSuffix": "", 

120 "fakesType": ""}): 

121 science = pipeBase.connectionTypes.Input( 

122 doc="Input science exposure.", 

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

124 storageClass="ExposureF", 

125 name="{fakesType}calexp" 

126 ) 

127 matchedTemplate = pipeBase.connectionTypes.Input( 

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

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

130 storageClass="ExposureF", 

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

132 ) 

133 difference = pipeBase.connectionTypes.Input( 

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

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

136 storageClass="ExposureF", 

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

138 ) 

139 kernelSources = pipeBase.connectionTypes.Input( 

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

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

142 storageClass="SourceCatalog", 

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

144 ) 

145 outputSchema = pipeBase.connectionTypes.InitOutput( 

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

147 storageClass="SourceCatalog", 

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

149 ) 

150 diaSources = pipeBase.connectionTypes.Output( 

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

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

153 storageClass="SourceCatalog", 

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

155 ) 

156 subtractedMeasuredExposure = pipeBase.connectionTypes.Output( 

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

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

159 storageClass="ExposureF", 

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

161 ) 

162 differenceBackground = pipeBase.connectionTypes.Output( 

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

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

165 storageClass="Background", 

166 name="difference_background", 

167 ) 

168 maskedStreaks = pipeBase.connectionTypes.Output( 

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

170 storageClass="ArrowNumpyDict", 

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

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

173 ) 

174 glintTrailInfo = pipeBase.connectionTypes.Output( 

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

176 storageClass="ArrowNumpyDict", 

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

178 name="trailed_glints", 

179 ) 

180 

181 def __init__(self, *, config): 

182 super().__init__(config=config) 

183 if not (self.config.doCalculateResidualMetics): 183 ↛ 184line 183 didn't jump to line 184 because the condition on line 183 was never true

184 self.inputs.remove("kernelSources") 

185 if not (self.config.writeStreakInfo and self.config.doMaskStreaks): 185 ↛ 188line 185 didn't jump to line 188 because the condition on line 185 was always true

186 self.outputs.remove("maskedStreaks") 

187 # In future mode the background is recorded on the output images. 

188 if (self.config.image_type == "future" 

189 or not (self.config.doSubtractBackground and self.config.doWriteBackground)): 

190 self.outputs.remove("differenceBackground") 

191 if not (self.config.writeGlintInfo): 191 ↛ 192line 191 didn't jump to line 192 because the condition on line 191 was never true

192 self.outputs.remove("glintTrailInfo") 

193 if self.config.image_type == "future": 

194 self.difference = dataclasses.replace(self.difference, storageClass="DifferenceImage") 

195 self.subtractedMeasuredExposure = dataclasses.replace(self.subtractedMeasuredExposure, 

196 storageClass="DifferenceImage") 

197 

198 

199class DetectAndMeasureConfig(pipeBase.PipelineTaskConfig, 

200 pipelineConnections=DetectAndMeasureConnections): 

201 """Config for DetectAndMeasureTask 

202 """ 

203 doMerge = pexConfig.Field( 

204 dtype=bool, 

205 default=True, 

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

207 "set by growFootprint" 

208 ) 

209 doForcedMeasurement = pexConfig.Field( 

210 dtype=bool, 

211 default=True, 

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

213 ) 

214 doAddMetrics = pexConfig.Field( 

215 dtype=bool, 

216 default=False, 

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

218 ) 

219 doSubtractBackground = pexConfig.Field( 

220 dtype=bool, 

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

222 default=True, 

223 ) 

224 doWriteBackground = pexConfig.Field( 

225 dtype=bool, 

226 doc="Persist the fitted background model? Ignored when" 

227 " image_type='future', which records it on the output images.", 

228 default=False, 

229 deprecated="Use image_type='future', which records the background on the output images." 

230 " Will be removed after v31.", 

231 ) 

232 doCalculateResidualMetics = pexConfig.Field( 

233 dtype=bool, 

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

235 "metadata?", 

236 default=True, 

237 ) 

238 subtractInitialBackground = pexConfig.ConfigurableField( 

239 target=lsst.meas.algorithms.SubtractBackgroundTask, 

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

241 ) 

242 subtractFinalBackground = pexConfig.ConfigurableField( 

243 target=lsst.meas.algorithms.SubtractBackgroundTask, 

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

245 ) 

246 doFindCosmicRays = pexConfig.Field( 

247 dtype=bool, 

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

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

250 default=True, 

251 ) 

252 cosmicray = pexConfig.ConfigField( 

253 dtype=FindCosmicRaysConfig, 

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

255 ) 

256 detection = pexConfig.ConfigurableField( 

257 target=SourceDetectionTask, 

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

259 ) 

260 streakDetection = pexConfig.ConfigurableField( 

261 target=SourceDetectionTask, 

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

263 ) 

264 doDeblend = pexConfig.Field( 

265 dtype=bool, 

266 default=False, 

267 doc="Deblend DIASources after detection?" 

268 ) 

269 deblend = pexConfig.ConfigurableField( 

270 target=lsst.meas.deblender.SourceDeblendTask, 

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

272 ) 

273 measurement = pexConfig.ConfigurableField( 

274 target=DipoleFitTask, 

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

276 ) 

277 doApCorr = lsst.pex.config.Field( 

278 dtype=bool, 

279 default=True, 

280 doc="Run subtask to apply aperture corrections" 

281 ) 

282 applyApCorr = lsst.pex.config.ConfigurableField( 

283 target=ApplyApCorrTask, 

284 doc="Task to apply aperture corrections" 

285 ) 

286 forcedMeasurement = pexConfig.ConfigurableField( 

287 target=ForcedMeasurementTask, 

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

289 ) 

290 growFootprint = pexConfig.Field( 

291 dtype=int, 

292 default=2, 

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

294 ) 

295 diaSourceMatchRadius = pexConfig.Field( 

296 dtype=float, 

297 default=0.5, 

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

299 ) 

300 doSkySources = pexConfig.Field( 

301 dtype=bool, 

302 default=False, 

303 doc="Generate sky sources?", 

304 ) 

305 skySources = pexConfig.ConfigurableField( 

306 target=SkyObjectsTask, 

307 doc="Generate sky sources", 

308 ) 

309 doMaskStreaks = pexConfig.Field( 

310 dtype=bool, 

311 default=True, 

312 doc="Turn on streak masking", 

313 ) 

314 maskStreaks = pexConfig.ConfigurableField( 

315 target=MaskStreaksTask, 

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

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

318 ) 

319 streakBinFactor = pexConfig.Field( 

320 dtype=int, 

321 default=4, 

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

323 "Only used if doMaskStreaks is True.", 

324 ) 

325 writeStreakInfo = pexConfig.Field( 

326 dtype=bool, 

327 default=False, 

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

329 "development work." 

330 ) 

331 findGlints = pexConfig.ConfigurableField( 

332 target=FindGlintTrailsTask, 

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

334 ) 

335 writeGlintInfo = pexConfig.Field( 

336 dtype=bool, 

337 default=True, 

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

339 ) 

340 setPrimaryFlags = pexConfig.ConfigurableField( 

341 target=SetPrimaryFlagsTask, 

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

343 ) 

344 badSourceFlags = lsst.pex.config.ListField( 

345 dtype=str, 

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

347 default=("base_PixelFlags_flag_offimage", 

348 "base_PixelFlags_flag_interpolatedCenterAll", 

349 "base_PixelFlags_flag_badCenter", 

350 "base_PixelFlags_flag_edgeCenter", 

351 "base_PixelFlags_flag_nodataCenter", 

352 "base_PixelFlags_flag_saturatedCenter", 

353 "base_PixelFlags_flag_sat_templateCenter", 

354 "base_PixelFlags_flag_spikeCenter", 

355 ), 

356 ) 

357 clearMaskPlanes = lsst.pex.config.ListField( 

358 dtype=str, 

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

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

361 ) 

362 doRejectBadMaskPlaneDetections = pexConfig.Field( 

363 dtype=bool, 

364 default=True, 

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

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

367 "but filtering earlier can save time." 

368 ) 

369 badMaskPlanes = lsst.pex.config.ListField( 

370 dtype=str, 

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

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

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

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

375 ) 

376 raiseOnBadSubtractionRatio = pexConfig.Field( 

377 dtype=bool, 

378 default=True, 

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

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

381 " image exceeds ``badSubtractionRatioThreshold``", 

382 ) 

383 badSubtractionRatioThreshold = pexConfig.Field( 

384 dtype=float, 

385 default=0.2, 

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

387 " the same footprints on the science image." 

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

389 ) 

390 badSubtractionVariationThreshold = pexConfig.Field( 

391 dtype=float, 

392 default=0.4, 

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

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

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

396 ) 

397 raiseOnNoDiaSources = pexConfig.Field( 

398 dtype=bool, 

399 default=True, 

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

401 ) 

402 run_sattle = pexConfig.Field( 

403 dtype=bool, 

404 default=False, 

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

406 "to the sattle service." 

407 ) 

408 sattle_historical = pexConfig.Field( 

409 dtype=bool, 

410 default=False, 

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

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

413 "closest in time to the exposure." 

414 ) 

415 idGenerator = DetectorVisitIdGeneratorConfig.make_field() 

416 image_type = pexConfig.ChoiceField[str]( 

417 "Image type of the input and measured difference images (and, for" 

418 " DetectAndMeasureScoreTask, the input and measured score images).", 

419 allowed={ 

420 "legacy": "Read and write lsst.afw.image.ExposureF.", 

421 "future": "Read and write lsst.images.DifferenceImage.", 

422 }, 

423 optional=False, 

424 default="legacy", 

425 ) 

426 

427 def setDefaults(self): 

428 # Background subtraction 

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

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

431 # faint diaSources 

432 self.subtractInitialBackground.binSize = 8 

433 self.subtractInitialBackground.useApprox = False 

434 self.subtractInitialBackground.statisticsProperty = "MEDIAN" 

435 self.subtractInitialBackground.doFilterSuperPixels = True 

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

437 "EDGE", 

438 "DETECTED", 

439 "DETECTED_NEGATIVE", 

440 "NO_DATA", 

441 ] 

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

443 # over-subtraction of bright objects. 

444 self.subtractFinalBackground.binSize = 40 

445 self.subtractFinalBackground.useApprox = False 

446 self.subtractFinalBackground.statisticsProperty = "MEDIAN" 

447 self.subtractFinalBackground.doFilterSuperPixels = True 

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

449 "EDGE", 

450 "DETECTED", 

451 "DETECTED_NEGATIVE", 

452 "NO_DATA", 

453 ] 

454 # Cosmic ray detection 

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

456 # DiaSource Detection 

457 self.detection.thresholdPolarity = "both" 

458 self.detection.thresholdValue = 5.0 

459 self.detection.reEstimateBackground = False 

460 self.detection.thresholdType = "pixel_stdev" 

461 self.detection.excludeMaskPlanes = [] 

462 

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

464 self.streakDetection.thresholdType = self.detection.thresholdType 

465 self.streakDetection.reEstimateBackground = False 

466 self.streakDetection.excludeMaskPlanes = self.detection.excludeMaskPlanes 

467 self.streakDetection.thresholdValue = self.detection.thresholdValue 

468 # Only detect positive streaks 

469 self.streakDetection.thresholdPolarity = "positive" 

470 # Do not grow detected mask for streaks 

471 self.streakDetection.nSigmaToGrow = 0 

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

473 # detected mask is set. 

474 self.maskStreaks.onlyMaskDetected = False 

475 # Restrict streak masking from growing too large 

476 self.maskStreaks.maxStreakWidth = 100 

477 # Restrict the number of iterations allowed for fitting streaks 

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

479 self.maskStreaks.maxFitIter = 10 

480 # Only mask to 2 sigma in width 

481 self.maskStreaks.nSigmaMask = 2 

482 # Threshold for including streaks after the Hough Transform. 

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

484 self.maskStreaks.absMinimumKernelHeight = 2 

485 

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

487 "base_LocalPhotoCalib", 

488 "base_LocalWcs", 

489 "ext_shapeHSM_HsmSourceMoments", 

490 "ext_shapeHSM_HsmPsfMoments", 

491 "base_ClassificationSizeExtendedness", 

492 ] 

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

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

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

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

497 self.forcedMeasurement.copyColumns = { 

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

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

500 self.forcedMeasurement.slots.shape = None 

501 

502 # Keep track of which footprints contain streaks 

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

504 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"] 

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

506 "STREAK", "INJECTED", "INJECTED_TEMPLATE", "HIGH_VARIANCE", "SAT_TEMPLATE", "SPIKE"] 

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

508 

509 def validate(self): 

510 super().validate() 

511 

512 if self.run_sattle: 

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

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

515 "Sattle requested but SATTLE_URI_BASE " 

516 "environment variable not set.") 

517 

518 

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

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

521 """ 

522 ConfigClass = DetectAndMeasureConfig 

523 _DefaultName = "detectAndMeasure" 

524 

525 def __init__(self, **kwargs): 

526 super().__init__(**kwargs) 

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

528 

529 self.algMetadata = dafBase.PropertyList() 

530 if self.config.doSubtractBackground: 

531 self.makeSubtask("subtractInitialBackground") 

532 self.makeSubtask("subtractFinalBackground") 

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

534 if self.config.doDeblend: 

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

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

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

538 algMetadata=self.algMetadata) 

539 if self.config.doApCorr: 539 ↛ 541line 539 didn't jump to line 541 because the condition on line 539 was always true

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

541 if self.config.doForcedMeasurement: 

542 self.schema.addField( 

543 "ip_diffim_forced_PsfFlux_instFlux", "D", 

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

545 units="count") 

546 self.schema.addField( 

547 "ip_diffim_forced_PsfFlux_instFluxErr", "D", 

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

549 units="count") 

550 self.schema.addField( 

551 "ip_diffim_forced_PsfFlux_area", "F", 

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

553 units="pixel") 

554 self.schema.addField( 

555 "ip_diffim_forced_PsfFlux_flag", "Flag", 

556 "Forced PSF flux general failure flag.") 

557 self.schema.addField( 

558 "ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag", 

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

560 self.schema.addField( 

561 "ip_diffim_forced_PsfFlux_flag_edge", "Flag", 

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

563 self.schema.addField( 

564 "ip_diffim_forced_template_PsfFlux_instFlux", "D", 

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

566 units="count") 

567 self.schema.addField( 

568 "ip_diffim_forced_template_PsfFlux_instFluxErr", "D", 

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

570 units="count") 

571 self.schema.addField( 

572 "ip_diffim_forced_template_PsfFlux_area", "F", 

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

574 units="pixel") 

575 self.schema.addField( 

576 "ip_diffim_forced_template_PsfFlux_flag", "Flag", 

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

578 self.schema.addField( 

579 "ip_diffim_forced_template_PsfFlux_flag_noGoodPixels", "Flag", 

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

581 self.schema.addField( 

582 "ip_diffim_forced_template_PsfFlux_flag_edge", "Flag", 

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

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

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

586 

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

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

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

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

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

592 if self.config.doMaskStreaks: 592 ↛ 595line 592 didn't jump to line 595 because the condition on line 592 was always true

593 self.makeSubtask("maskStreaks") 

594 self.makeSubtask("streakDetection") 

595 self.makeSubtask("findGlints") 

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

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

598 self.schema.addField("reliabilityVersion", type=str, size=7, doc="Version of the reliability model") 

599 

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

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

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

603 

604 # Check that the schema and config are consistent 

605 for flag in self.config.badSourceFlags: 

606 if flag not in self.schema: 

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

608 

609 # initialize InitOutputs 

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

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

612 

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

614 inputRefs: pipeBase.InputQuantizedConnection, 

615 outputRefs: pipeBase.OutputQuantizedConnection): 

616 inputs = butlerQC.get(inputRefs) 

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

618 idFactory = idGenerator.make_table_id_factory() 

619 

620 # Input image that each output is measured from, keyed by output name. 

621 inputImages = {} 

622 if self.config.image_type == "future": 

623 dataId = butlerQC.quantum.dataId 

624 for output, connection in (("subtractedMeasuredExposure", "difference"), 

625 ("scoreMeasuredExposure", "scoreExposure")): 

626 if connection not in inputs: 

627 continue 

628 image = inputs[connection] 

629 # The output is written back into this image, so it must 

630 # already carry the kernel and the template records. 

631 get_difference_image_provenance(image, "kernel", connection, dataId) 

632 get_difference_image_provenance(image, "templates", connection, dataId) 

633 inputImages[output] = image 

634 inputs[connection] = image.to_legacy() 

635 

636 # Specify the fields that `annotate` and `convert_outputs_to_future` 

637 # need below, to ensure they exist, even as None. 

638 measurementResults = pipeBase.Struct( 

639 subtractedMeasuredExposure=None, 

640 scoreMeasuredExposure=None, 

641 diaSources=None, 

642 maskedStreaks=None, 

643 differenceBackground=None, 

644 ) 

645 try: 

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

647 except pipeBase.AlgorithmError as e: 

648 # Annotate before converting 

649 error = pipeBase.AnnotatedPartialOutputsError.annotate( 

650 e, 

651 self, 

652 measurementResults.subtractedMeasuredExposure, 

653 measurementResults.scoreMeasuredExposure, 

654 measurementResults.diaSources, 

655 measurementResults.maskedStreaks, 

656 log=self.log 

657 ) 

658 if self.config.image_type == "future": 

659 self.convert_outputs_to_future(measurementResults, inputImages) 

660 butlerQC.put(measurementResults, outputRefs) 

661 raise error from e 

662 if self.config.image_type == "future": 

663 self.convert_outputs_to_future(measurementResults, inputImages) 

664 butlerQC.put(measurementResults, outputRefs) 

665 

666 def convert_outputs_to_future(self, results, inputImages): 

667 """Write the output images in a result struct back into the 

668 `lsst.images` images they were measured from. 

669 

670 The image, variance, and mask planes of each input image are replaced 

671 by those of the corresponding legacy output, and the background 

672 subtracted by this task is attached. The kernel, template records, 

673 and observation metadata of the input are kept. Outputs that are 

674 absent or `None` are skipped. 

675 

676 Parameters 

677 ---------- 

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

679 Output struct to read and modify in place. 

680 inputImages : `dict` [`str`, `lsst.images.DifferenceImage`] 

681 Input image each output was measured from, keyed by the name of 

682 the output. Each is modified in place. 

683 

684 Notes 

685 ----- 

686 This uses the default legacy mask plane map of 

687 `lsst.images.get_legacy_difference_image_mask_planes` and adds any 

688 planes from `~lsst.images.get_legacy_optional_mask_planes` that are set 

689 for any pixels in the image. 

690 """ 

691 background = getattr(results, "differenceBackground", None) 

692 for name, image in inputImages.items(): 

693 exposure = getattr(results, name, None) 

694 if exposure is None: 694 ↛ 695line 694 didn't jump to line 695 because the condition on line 694 was never true

695 continue 

696 # These are no-ops when the legacy output is a view of the input, 

697 # as returned by `~lsst.images.DifferenceImage.to_legacy`. 

698 image.image.array[...] = exposure.image.array 

699 image.variance.array[...] = exposure.variance.array 

700 image.mask = Mask.from_legacy(exposure.mask, get_legacy_difference_image_mask_planes(), 

701 sky_projection=image.sky_projection) 

702 if background is not None and len(background) > 0: 

703 # The same background is subtracted from the difference image 

704 # and the score image, so it applies to both. 

705 image.backgrounds.add( 

706 "subtracted", 

707 field_from_legacy_background(background, bounds=image.bbox, unit=image.unit), 

708 description="Background subtracted from this image before detection.", 

709 is_subtracted=True, 

710 ) 

711 setattr(results, name, image) 

712 if background is not None and len(background) == 0: 

713 # An empty BackgroundList cannot be converted to a field, and means 

714 # that no background was subtracted. 

715 self.log.debug("No background model to attach to the output images.") 

716 

717 @timeMethod 

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

719 idFactory=None, measurementResults=None): 

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

721 

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

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

724 Close positive and negative detections will optionally be merged into 

725 dipole diaSources. 

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

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

728 detections. 

729 

730 Parameters 

731 ---------- 

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

733 Science exposure that the template was subtracted from. 

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

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

736 difference image. 

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

738 Result of subtracting template from the science image. 

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

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

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

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

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

744 deblending and merging positive/negative peaks. 

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

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

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

748 this is also returned. 

749 

750 Returns 

751 ------- 

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

753 

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

755 Subtracted exposure with detection mask applied. 

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

757 The catalog of detected sources. 

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

759 Background that was subtracted from the difference image. 

760 """ 

761 if measurementResults is None: 

762 measurementResults = pipeBase.Struct() 

763 if idFactory is None: 

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

765 

766 # Check image properties and clear detection mask planes. 

767 self._prepareInputs(difference) 

768 if self.config.doSubtractBackground: 

769 background = self._fitAndSubtractBackground(differenceExposure=difference) 

770 else: 

771 background = afwMath.BackgroundList() 

772 

773 if self.config.doFindCosmicRays: 773 ↛ 779line 773 didn't jump to line 779 because the condition on line 773 was always true

774 self.findAndMaskCosmicRays(difference) 

775 

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

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

778 # know about past ids). 

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

780 results = self.detection.run( 

781 table=table, 

782 exposure=difference, 

783 doSmooth=True, 

784 background=background, 

785 clearMask=True, 

786 ) 

787 measurementResults.differenceBackground = background 

788 

789 if self.config.doDeblend: 

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

791 results.positive, 

792 results.negative) 

793 

794 else: 

795 positives = afwTable.SourceCatalog(self.schema) 

796 results.positive.makeSources(positives) 

797 negatives = afwTable.SourceCatalog(self.schema) 

798 results.negative.makeSources(negatives) 

799 sources = results.sources 

800 

801 self.processResults(science, matchedTemplate, difference, 

802 sources, idFactory, kernelSources, 

803 positives=positives, 

804 negatives=negatives, 

805 measurementResults=measurementResults) 

806 return measurementResults 

807 

808 def _prepareInputs(self, difference): 

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

810 

811 Parameters 

812 ---------- 

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

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

815 The mask plane will be modified in place. 

816 

817 Raises 

818 ------ 

819 lsst.pipe.base.UpstreamFailureNoWorkFound 

820 If the PSF is not usable for measurement. 

821 """ 

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

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

824 if np.isnan(sigma): 

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

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

827 mask = difference.mask 

828 for mp in self.config.clearMaskPlanes: 

829 if mp not in mask.getMaskPlaneDict(): 

830 mask.addMaskPlane(mp) 

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

832 

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

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

835 

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

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

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

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

840 scheme is what suppresses spurious detections from glints and 

841 extended structures: the small-binsize initial fit aggressively 

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

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

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

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

846 detection runs on a properly background-subtracted score. 

847 

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

849 leaves on the clone is used downstream. 

850 

851 Parameters 

852 ---------- 

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

854 Difference image. The final background is subtracted from its 

855 image in place but its mask is not modified.