Coverage for python / lsst / ip / diffim / subtractImages.py: 20%

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 

22from astropy import units as u 

23from astropy.stats import gaussian_fwhm_to_sigma 

24import numpy as np 

25 

26import lsst.afw.detection as afwDetection 

27import lsst.afw.image 

28import lsst.afw.math 

29import lsst.geom 

30from lsst.ip.diffim.utils import (evaluateMeanPsfFwhm, getPsfFwhm, 

31 computeDifferenceImageMetrics, 

32 checkMask, setSourceFootprints) 

33from lsst.meas.algorithms import ScaleVarianceTask, ScienceSourceSelectorTask 

34import lsst.pex.config 

35import lsst.pipe.base 

36import lsst.pex.exceptions 

37from lsst.pipe.base import connectionTypes 

38from . import MakeKernelTask, DecorrelateALKernelTask 

39from lsst.utils.timer import timeMethod 

40 

41__all__ = ["AlardLuptonSubtractConfig", "AlardLuptonSubtractTask", 

42 "AlardLuptonPreconvolveSubtractConfig", "AlardLuptonPreconvolveSubtractTask", 

43 "SimplifiedSubtractConfig", "SimplifiedSubtractTask", 

44 "InsufficientKernelSourcesError"] 

45 

46_dimensions = ("instrument", "visit", "detector") 

47_defaultTemplates = {"coaddName": "deep", "fakesType": ""} 

48 

49 

50class InsufficientKernelSourcesError(lsst.pipe.base.AlgorithmError): 

51 """Raised when there are too few sources to calculate the PSF matching 

52 kernel. 

53 """ 

54 def __init__(self, *, nSources, nRequired): 

55 msg = (f"Only {nSources} sources were selected for PSF matching," 

56 f" but {nRequired} are required.") 

57 super().__init__(msg) 

58 self.nSources = nSources 

59 self.nRequired = nRequired 

60 

61 @property 

62 def metadata(self): 

63 return {"nSources": self.nSources, 

64 "nRequired": self.nRequired 

65 } 

66 

67 

68class SubtractInputConnections(lsst.pipe.base.PipelineTaskConnections, 

69 dimensions=_dimensions, 

70 defaultTemplates=_defaultTemplates): 

71 template = connectionTypes.Input( 

72 doc="Input warped template to subtract.", 

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

74 storageClass="ExposureF", 

75 name="{fakesType}{coaddName}Diff_templateExp" 

76 ) 

77 science = connectionTypes.Input( 

78 doc="Input science exposure to subtract from.", 

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

80 storageClass="ExposureF", 

81 name="{fakesType}calexp" 

82 ) 

83 sources = connectionTypes.Input( 

84 doc="Sources measured on the science exposure; " 

85 "used to select sources for making the matching kernel.", 

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

87 storageClass="SourceCatalog", 

88 name="{fakesType}src" 

89 ) 

90 visitSummary = connectionTypes.Input( 

91 doc=("Per-visit catalog with final calibration objects. " 

92 "These catalogs use the detector id for the catalog id, " 

93 "sorted on id for fast lookup."), 

94 dimensions=("instrument", "visit"), 

95 storageClass="ExposureCatalog", 

96 name="finalVisitSummary", 

97 ) 

98 

99 def __init__(self, *, config=None): 

100 super().__init__(config=config) 

101 if not config.doApplyExternalCalibrations: 

102 del self.visitSummary 

103 

104 

105class SubtractImageOutputConnections(lsst.pipe.base.PipelineTaskConnections, 

106 dimensions=_dimensions, 

107 defaultTemplates=_defaultTemplates): 

108 difference = connectionTypes.Output( 

109 doc="Result of subtracting convolved template from science image.", 

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

111 storageClass="ExposureF", 

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

113 ) 

114 matchedTemplate = connectionTypes.Output( 

115 doc="Warped and PSF-matched template used to create `subtractedExposure`.", 

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

117 storageClass="ExposureF", 

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

119 ) 

120 psfMatchingKernel = connectionTypes.Output( 

121 doc="Kernel used to PSF match the science and template images.", 

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

123 storageClass="MatchingKernel", 

124 name="{fakesType}{coaddName}Diff_psfMatchKernel", 

125 ) 

126 kernelSources = connectionTypes.Output( 

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

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

129 storageClass="SourceCatalog", 

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

131 ) 

132 

133 

134class SubtractScoreOutputConnections(lsst.pipe.base.PipelineTaskConnections, 

135 dimensions=_dimensions, 

136 defaultTemplates=_defaultTemplates): 

137 scoreExposure = connectionTypes.Output( 

138 doc="The maximum likelihood image, used for the detection of diaSources.", 

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

140 storageClass="ExposureF", 

141 name="{fakesType}{coaddName}Diff_scoreTempExp", 

142 ) 

143 psfMatchingKernel = connectionTypes.Output( 

144 doc="Kernel used to PSF match the science and template images.", 

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

146 storageClass="MatchingKernel", 

147 name="{fakesType}{coaddName}Diff_psfScoreMatchKernel", 

148 ) 

149 kernelSources = connectionTypes.Output( 

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

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

152 storageClass="SourceCatalog", 

153 name="{fakesType}{coaddName}Diff_psfScoreMatchSources" 

154 ) 

155 

156 

157class AlardLuptonSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

158 pass 

159 

160 

161class SimplifiedSubtractConnections(SubtractInputConnections, SubtractImageOutputConnections): 

162 inputPsfMatchingKernel = connectionTypes.Input( 

163 doc="Kernel used to PSF match the science and template images.", 

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

165 storageClass="MatchingKernel", 

166 name="{fakesType}{coaddName}Diff_psfMatchKernel", 

167 ) 

168 

169 def __init__(self, *, config=None): 

170 super().__init__(config=config) 

171 del self.sources 

172 if config.useExistingKernel: 

173 del self.psfMatchingKernel 

174 del self.kernelSources 

175 else: 

176 del self.inputPsfMatchingKernel 

177 

178 

179class AlardLuptonSubtractBaseConfig(lsst.pex.config.Config): 

180 makeKernel = lsst.pex.config.ConfigurableField( 

181 target=MakeKernelTask, 

182 doc="Task to construct a matching kernel for convolution.", 

183 ) 

184 doDecorrelation = lsst.pex.config.Field( 

185 dtype=bool, 

186 default=True, 

187 doc="Perform diffim decorrelation to undo pixel correlation due to A&L " 

188 "kernel convolution? If True, also update the diffim PSF." 

189 ) 

190 decorrelate = lsst.pex.config.ConfigurableField( 

191 target=DecorrelateALKernelTask, 

192 doc="Task to decorrelate the image difference.", 

193 ) 

194 requiredTemplateFraction = lsst.pex.config.Field( 

195 dtype=float, 

196 default=0.1, 

197 doc="Raise NoWorkFound and do not attempt image subtraction if template covers less than this " 

198 " fraction of pixels. Setting to 0 will always attempt image subtraction." 

199 ) 

200 minTemplateFractionForExpectedSuccess = lsst.pex.config.Field( 

201 dtype=float, 

202 default=0.2, 

203 doc="Raise NoWorkFound if PSF-matching fails and template covers less than this fraction of pixels." 

204 " If the fraction of pixels covered by the template is less than this value (and greater than" 

205 " requiredTemplateFraction) this task is attempted but failure is anticipated and tolerated." 

206 ) 

207 doScaleVariance = lsst.pex.config.Field( 

208 dtype=bool, 

209 default=True, 

210 doc="Scale variance of the image difference?" 

211 ) 

212 scaleVariance = lsst.pex.config.ConfigurableField( 

213 target=ScaleVarianceTask, 

214 doc="Subtask to rescale the variance of the template to the statistically expected level." 

215 ) 

216 doSubtractBackground = lsst.pex.config.Field( 

217 doc="Subtract the background fit when solving the kernel? " 

218 "It is generally better to instead subtract the background in detectAndMeasure.", 

219 dtype=bool, 

220 default=False, 

221 ) 

222 doApplyExternalCalibrations = lsst.pex.config.Field( 

223 doc=( 

224 "Replace science Exposure's calibration objects with those" 

225 " in visitSummary. Ignored if `doApplyFinalizedPsf is True." 

226 ), 

227 dtype=bool, 

228 default=False, 

229 ) 

230 sourceSelector = lsst.pex.config.ConfigurableField( 

231 target=ScienceSourceSelectorTask, 

232 doc="Task to select sources to be used for PSF matching.", 

233 ) 

234 fallbackSourceSelector = lsst.pex.config.ConfigurableField( 

235 target=ScienceSourceSelectorTask, 

236 doc="Task to select sources to be used for PSF matching." 

237 "Used only if the kernel calculation fails and" 

238 "`allowKernelSourceDetection` is set. The fallback source detection" 

239 " will not include all of the same plugins as the original source " 

240 " detection, so not all of the same flags can be used.", 

241 ) 

242 detectionThreshold = lsst.pex.config.Field( 

243 dtype=float, 

244 default=10, 

245 doc="Minimum signal to noise ratio of detected sources " 

246 "to use for calculating the PSF matching kernel.", 

247 deprecated="No longer used. Will be removed after v30" 

248 ) 

249 detectionThresholdMax = lsst.pex.config.Field( 

250 dtype=float, 

251 default=500, 

252 doc="Maximum signal to noise ratio of detected sources " 

253 "to use for calculating the PSF matching kernel.", 

254 deprecated="No longer used. Will be removed after v30" 

255 ) 

256 restrictKernelEdgeSources = lsst.pex.config.Field( 

257 dtype=bool, 

258 default=True, 

259 doc="Exclude sources close to the edge from the kernel calculation?" 

260 ) 

261 maxKernelSources = lsst.pex.config.Field( 

262 dtype=int, 

263 default=1000, 

264 doc="Maximum number of sources to use for calculating the PSF matching kernel." 

265 "Set to -1 to disable." 

266 ) 

267 minKernelSources = lsst.pex.config.Field( 

268 dtype=int, 

269 default=3, 

270 doc="Minimum number of sources needed for calculating the PSF matching kernel." 

271 ) 

272 excludeMaskPlanes = lsst.pex.config.ListField( 

273 dtype=str, 

274 default=("NO_DATA", "BAD", "SAT", "EDGE", "FAKE", "HIGH_VARIANCE"), 

275 doc="Template mask planes to exclude when selecting sources for PSF matching.", 

276 ) 

277 badMaskPlanes = lsst.pex.config.ListField( 

278 dtype=str, 

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

280 doc="Mask planes to interpolate over." 

281 ) 

282 preserveTemplateMask = lsst.pex.config.ListField( 

283 dtype=str, 

284 default=("NO_DATA", "BAD", "HIGH_VARIANCE"), 

285 doc="Mask planes from the template to propagate to the image difference." 

286 ) 

287 renameTemplateMask = lsst.pex.config.ListField( 

288 dtype=str, 

289 default=("SAT", "INJECTED", "INJECTED_CORE",), 

290 doc="Mask planes from the template to propagate to the image difference" 

291 "with '_TEMPLATE' appended to the name." 

292 ) 

293 allowKernelSourceDetection = lsst.pex.config.Field( 

294 dtype=bool, 

295 default=False, 

296 doc="Re-run source detection for kernel candidates if an error is" 

297 " encountered while calculating the matching kernel." 

298 ) 

299 

300 def setDefaults(self): 

301 self.makeKernel.kernel.name = "AL" 

302 # Always include background fitting in the kernel fit, 

303 # even if it is not subtracted 

304 self.makeKernel.kernel.active.fitForBackground = True 

305 self.makeKernel.kernel.active.spatialKernelOrder = 1 

306 self.makeKernel.kernel.active.spatialBgOrder = 2 

307 # Shared source selector settings 

308 doSkySources = False # Do not include sky sources 

309 doSignalToNoise = True # apply signal to noise filter 

310 doUnresolved = True # apply star-galaxy separation 

311 signalToNoiseMinimum = 10 

312 signalToNoiseMaximum = 500 

313 self.sourceSelector.doIsolated = True # apply isolated star selection 

314 self.sourceSelector.doRequirePrimary = True # apply primary flag selection 

315 self.sourceSelector.doUnresolved = doUnresolved 

316 self.sourceSelector.doSkySources = doSkySources 

317 self.sourceSelector.doSignalToNoise = doSignalToNoise 

318 self.sourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

319 self.sourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

320 # The following two configs should not be necessary to be turned on for 

321 # PSF-matching, and the fallback kernel source selection will fail if 

322 # they are set since it does not run deblending. 

323 self.fallbackSourceSelector.doIsolated = False # Do not apply isolated star selection 

324 self.fallbackSourceSelector.doRequirePrimary = False # Do not apply primary flag selection 

325 self.fallbackSourceSelector.doUnresolved = doUnresolved 

326 self.fallbackSourceSelector.doSkySources = doSkySources 

327 self.fallbackSourceSelector.doSignalToNoise = doSignalToNoise 

328 self.fallbackSourceSelector.signalToNoise.minimum = signalToNoiseMinimum 

329 self.fallbackSourceSelector.signalToNoise.maximum = signalToNoiseMaximum 

330 

331 

332class AlardLuptonSubtractConfig(AlardLuptonSubtractBaseConfig, lsst.pipe.base.PipelineTaskConfig, 

333 pipelineConnections=AlardLuptonSubtractConnections): 

334 mode = lsst.pex.config.ChoiceField( 

335 dtype=str, 

336 default="convolveTemplate", 

337 allowed={"auto": "Choose which image to convolve at runtime.", 

338 "convolveScience": "Only convolve the science image.", 

339 "convolveTemplate": "Only convolve the template image."}, 

340 doc="Choose which image to convolve at runtime, or require that a specific image is convolved." 

341 ) 

342 

343 

344class AlardLuptonSubtractTask(lsst.pipe.base.PipelineTask): 

345 """Compute the image difference of a science and template image using 

346 the Alard & Lupton (1998) algorithm. 

347 """ 

348 ConfigClass = AlardLuptonSubtractConfig 

349 _DefaultName = "alardLuptonSubtract" 

350 usePreconvolution = False 

351 """Whether this task preconvolves the science image with its own PSF 

352 before kernel-matching. Subclasses that preconvolve override this to 

353 `True`.""" 

354 

355 def __init__(self, **kwargs): 

356 super().__init__(**kwargs) 

357 self.makeSubtask("decorrelate") 

358 self.makeSubtask("makeKernel") 

359 self.makeSubtask("sourceSelector") 

360 self.makeSubtask("fallbackSourceSelector") 

361 if self.config.doScaleVariance: 

362 self.makeSubtask("scaleVariance") 

363 

364 self.convolutionControl = lsst.afw.math.ConvolutionControl() 

365 # Normalization is an extra, unnecessary, calculation and will result 

366 # in mis-subtraction of the images if there are calibration errors. 

367 self.convolutionControl.setDoNormalize(False) 

368 self.convolutionControl.setDoCopyEdge(True) 

369 

370 def _applyExternalCalibrations(self, exposure, visitSummary): 

371 """Replace calibrations (psf, and ApCorrMap) on this exposure with 

372 external ones.". 

373 

374 Parameters 

375 ---------- 

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

377 Input exposure to adjust calibrations. 

378 visitSummary : `lsst.afw.table.ExposureCatalog` 

379 Exposure catalog with external calibrations to be applied. Catalog 

380 uses the detector id for the catalog id, sorted on id for fast 

381 lookup. 

382 

383 Returns 

384 ------- 

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

386 Exposure with adjusted calibrations. 

387 """ 

388 detectorId = exposure.info.getDetector().getId() 

389 

390 row = visitSummary.find(detectorId) 

391 if row is None: 

392 self.log.warning("Detector id %s not found in external calibrations catalog; " 

393 "Using original calibrations.", detectorId) 

394 else: 

395 psf = row.getPsf() 

396 apCorrMap = row.getApCorrMap() 

397 if psf is None: 

398 self.log.warning("Detector id %s has None for psf in " 

399 "external calibrations catalog; Using original psf and aperture correction.", 

400 detectorId) 

401 elif apCorrMap is None: 

402 self.log.warning("Detector id %s has None for apCorrMap in " 

403 "external calibrations catalog; Using original psf and aperture correction.", 

404 detectorId) 

405 else: 

406 exposure.setPsf(psf) 

407 exposure.info.setApCorrMap(apCorrMap) 

408 

409 return exposure 

410 

411 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

412 inputs = butlerQC.get(inputRefs) 

413 

414 try: 

415 results = self.run(**inputs) 

416 except lsst.pipe.base.AlgorithmError as e: 

417 error = lsst.pipe.base.AnnotatedPartialOutputsError.annotate(e, self, log=self.log) 

418 # No partial outputs for butler to put 

419 raise error from e 

420 

421 butlerQC.put(results, outputRefs) 

422 

423 @timeMethod 

424 def run(self, template, science, sources, visitSummary=None): 

425 """PSF match, subtract, and decorrelate two images. 

426 

427 Parameters 

428 ---------- 

429 template : `lsst.afw.image.ExposureF` 

430 Template exposure, warped to match the science exposure. 

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

432 Science exposure to subtract from the template. 

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

434 Identified sources on the science exposure. This catalog is used to 

435 select sources in order to perform the AL PSF matching on stamp 

436 images around them. 

437 visitSummary : `lsst.afw.table.ExposureCatalog`, optional 

438 Exposure catalog with external calibrations to be applied. Catalog 

439 uses the detector id for the catalog id, sorted on id for fast 

440 lookup. 

441 

442 Returns 

443 ------- 

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

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

446 Result of subtracting template and science. 

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

448 Warped and PSF-matched template exposure. 

449 ``backgroundModel`` : `lsst.afw.math.Function2D` 

450 Background model that was fit while solving for the 

451 PSF-matching kernel 

452 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

453 Kernel used to PSF-match the convolved image. 

454 ``kernelSources` : `lsst.afw.table.SourceCatalog` 

455 Sources from the input catalog that were used to construct the 

456 PSF-matching kernel. 

457 """ 

458 self._prepareInputs(template, science, visitSummary=visitSummary) 

459 

460 convolveTemplate = self.chooseConvolutionMethod(template, science) 

461 self.matchedPsfSize = self.sciencePsfSize if convolveTemplate else self.templatePsfSize 

462 

463 kernelResult = self.runMakeKernel(template, science, sources=sources, 

464 convolveTemplate=convolveTemplate, 

465 runSourceDetection=False) 

466 

467 if self.config.doSubtractBackground: 

468 backgroundModel = kernelResult.backgroundModel 

469 else: 

470 backgroundModel = None 

471 if convolveTemplate: 

472 subtractResults = self.runConvolveTemplate(template, science, kernelResult.psfMatchingKernel, 

473 backgroundModel=backgroundModel) 

474 else: 

475 subtractResults = self.runConvolveScience(template, science, kernelResult.psfMatchingKernel, 

476 backgroundModel=backgroundModel) 

477 subtractResults.kernelSources = kernelResult.kernelSources 

478 

479 metrics = computeDifferenceImageMetrics(science, subtractResults.difference, sources) 

480 

481 self.metadata["differenceFootprintRatioMean"] = metrics.differenceFootprintRatioMean 

482 self.metadata["differenceFootprintRatioStdev"] = metrics.differenceFootprintRatioStdev 

483 self.metadata["differenceFootprintSkyRatioMean"] = metrics.differenceFootprintSkyRatioMean 

484 self.metadata["differenceFootprintSkyRatioStdev"] = metrics.differenceFootprintSkyRatioStdev 

485 self.log.info("Mean, stdev of ratio of difference to science " 

486 "pixels in star footprints: %5.4f, %5.4f", 

487 self.metadata["differenceFootprintRatioMean"], 

488 self.metadata["differenceFootprintRatioStdev"]) 

489 

490 return subtractResults 

491 

492 def chooseConvolutionMethod(self, template, science): 

493 """Determine whether the template should be convolved with the PSF 

494 matching kernel. 

495 

496 Parameters 

497 ---------- 

498 template : `lsst.afw.image.ExposureF` 

499 Template exposure, warped to match the science exposure. 

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

501 Science exposure to subtract from the template. 

502 

503 Returns 

504 ------- 

505 convolveTemplate : `bool` 

506 Convolve the template to match the two images? 

507 

508 Raises 

509 ------ 

510 RuntimeError 

511 If an unsupported convolution mode is supplied. 

512 """ 

513 if self.usePreconvolution: 

514 raise RuntimeError("Choosing a convolution method is incompatible with preconvolution!") 

515 if self.config.mode == "auto": 

516 convolveTemplate = _shapeTest(template, 

517 science, 

518 fwhmExposureBuffer=self.config.makeKernel.fwhmExposureBuffer, 

519 fwhmExposureGrid=self.config.makeKernel.fwhmExposureGrid) 

520 if convolveTemplate: 

521 if self.sciencePsfSize < self.templatePsfSize: 

522 self.log.info("Average template PSF size is greater, " 

523 "but science PSF greater in one dimension: convolving template image.") 

524 else: 

525 self.log.info("Science PSF size is greater: convolving template image.") 

526 else: 

527 self.log.info("Template PSF size is greater: convolving science image.") 

528 elif self.config.mode == "convolveTemplate": 

529 self.log.info("`convolveTemplate` is set: convolving template image.") 

530 convolveTemplate = True 

531 elif self.config.mode == "convolveScience": 

532 self.log.info("`convolveScience` is set: convolving science image.") 

533 convolveTemplate = False 

534 else: 

535 raise RuntimeError(f"Cannot handle AlardLuptonSubtract mode: {self.config.mode}") 

536 return convolveTemplate 

537 

538 def runMakeKernel(self, template, science, sources=None, convolveTemplate=True, runSourceDetection=False): 

539 """Construct the PSF-matching kernel. Not used for preconvolution. 

540 

541 Parameters 

542 ---------- 

543 template : `lsst.afw.image.ExposureF` 

544 Template exposure, warped to match the science exposure. 

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

546 Science exposure to subtract from the template. 

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

548 Identified sources on the science exposure. This catalog is used to 

549 select sources in order to perform the AL PSF matching on stamp 

550 images around them. 

551 Not used if ``runSourceDetection`` is set. 

552 convolveTemplate : `bool`, optional 

553 Construct the matching kernel to convolve the template? 

554 runSourceDetection : `bool`, optional 

555 Run a minimal version of source detection to determine kernel 

556 candidates? If False, a source list to select kernel candidates 

557 from must be supplied. 

558 

559 Returns 

560 ------- 

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

562 ``backgroundModel`` : `lsst.afw.math.Function2D` 

563 Background model that was fit while solving for the 

564 PSF-matching kernel 

565 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

566 Kernel used to PSF-match the convolved image. 

567 ``kernelSources` : `lsst.afw.table.SourceCatalog` 

568 Sources from the input catalog that were used to construct the 

569 PSF-matching kernel. 

570 """ 

571 if self.usePreconvolution: 

572 raise RuntimeError("Incorrect matching kernel calculation configured. " 

573 "`runMakeKernel` can't be called if `usePreconvolution` is set.") 

574 if convolveTemplate: 

575 reference = template 

576 target = science 

577 referenceFwhmPix = self.templatePsfSize 

578 targetFwhmPix = self.sciencePsfSize 

579 else: 

580 reference = science 

581 target = template 

582 referenceFwhmPix = self.sciencePsfSize 

583 targetFwhmPix = self.templatePsfSize 

584 try: 

585 if runSourceDetection: 

586 kernelSources = self.runKernelSourceDetection(template, science) 

587 else: 

588 kernelSources = self._sourceSelector(template, science, sources) 

589 kernelResult = self.makeKernel.run(reference, target, kernelSources, 

590 preconvolved=False, 

591 templateFwhmPix=referenceFwhmPix, 

592 scienceFwhmPix=targetFwhmPix) 

593 except (RuntimeError, lsst.pex.exceptions.Exception) as e: 

594 self.log.warning("Failed to match template. Checking coverage") 

595 # Raise NoWorkFound if template fraction is insufficient 

596 checkTemplateIsSufficient(template[science.getBBox()], science, self.log, 

597 self.config.minTemplateFractionForExpectedSuccess, 

598 exceptionMessage="Template coverage lower than expected to succeed." 

599 f" Failure is tolerable: {e}") 

600 # checkTemplateIsSufficient did not raise NoWorkFound, so raise original exception 

601 raise e 

602 

603 return lsst.pipe.base.Struct(backgroundModel=kernelResult.backgroundModel, 

604 psfMatchingKernel=kernelResult.psfMatchingKernel, 

605 kernelSources=kernelSources) 

606 

607 def runKernelSourceDetection(self, template, science): 

608 """Run detection on the science image and use the template mask plane 

609 to reject candidate sources. 

610 

611 Parameters 

612 ---------- 

613 template : `lsst.afw.image.ExposureF` 

614 Template exposure, warped to match the science exposure. 

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

616 Science exposure to subtract from the template. 

617 

618 Returns 

619 ------- 

620 kernelSources : `lsst.afw.table.SourceCatalog` 

621 Sources from the input catalog to use to construct the 

622 PSF-matching kernel. 

623 """ 

624 kernelSize = self.makeKernel.makeKernelBasisList( 

625 self.templatePsfSize, self.matchedPsfSize)[0].getWidth() 

626 sources = self.makeKernel.makeCandidateList(template, science, kernelSize, 

627 candidateList=None, 

628 sigma=gaussian_fwhm_to_sigma*self.sciencePsfSize) 

629 return self._sourceSelector(template, science, sources, fallback=True) 

630 

631 def runConvolveTemplate(self, template, science, psfMatchingKernel, backgroundModel=None): 

632 """Convolve the template image with a PSF-matching kernel and subtract 

633 from the science image. 

634 

635 Parameters 

636 ---------- 

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

638 Template exposure, warped to match the science exposure. 

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

640 Science exposure to subtract from the template. 

641 psfMatchingKernel : `lsst.afw.math.Kernel` 

642 Kernel to be used to PSF-match the science image to the template. 

643 backgroundModel : `lsst.afw.math.Function2D`, optional 

644 Background model that was fit while solving for the PSF-matching 

645 kernel. 

646 

647 Returns 

648 ------- 

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

650 

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

652 Result of subtracting template and science. 

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

654 Warped and PSF-matched template exposure. 

655 ``backgroundModel`` : `lsst.afw.math.Function2D` 

656 Background model that was fit while solving for the PSF-matching kernel 

657 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

658 Kernel used to PSF-match the template to the science image. 

659 """ 

660 self.metadata["convolvedExposure"] = "Template" 

661 

662 matchedTemplate = self._convolveExposure(template, psfMatchingKernel, 

663 self.convolutionControl, 

664 bbox=science.getBBox(), 

665 psf=science.psf, 

666 photoCalib=science.photoCalib) 

667 

668 difference = _subtractImages(science, matchedTemplate, backgroundModel=backgroundModel) 

669 correctedExposure = self.finalize(template, science, difference, 

670 psfMatchingKernel, 

671 templateMatched=True) 

672 

673 return lsst.pipe.base.Struct(difference=correctedExposure, 

674 matchedTemplate=matchedTemplate, 

675 matchedScience=science, 

676 backgroundModel=backgroundModel, 

677 psfMatchingKernel=psfMatchingKernel) 

678 

679 def runConvolveScience(self, template, science, psfMatchingKernel, backgroundModel=None): 

680 """Convolve the science image with a PSF-matching kernel and subtract 

681 the template image. 

682 

683 Parameters 

684 ---------- 

685 template : `lsst.afw.image.ExposureF` 

686 Template exposure, warped to match the science exposure. 

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

688 Science exposure to subtract from the template. 

689 psfMatchingKernel : `lsst.afw.math.Kernel` 

690 Kernel to be used to PSF-match the science image to the template. 

691 backgroundModel : `lsst.afw.math.Function2D`, optional 

692 Background model that was fit while solving for the PSF-matching 

693 kernel. 

694 

695 Returns 

696 ------- 

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

698 

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

700 Result of subtracting template and science. 

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

702 Warped template exposure. Note that in this case, the template 

703 is not PSF-matched to the science image. 

704 ``backgroundModel`` : `lsst.afw.math.Function2D` 

705 Background model that was fit while solving for the PSF-matching kernel 

706 ``psfMatchingKernel`` : `lsst.afw.math.Kernel` 

707 Kernel used to PSF-match the science image to the template. 

708 """ 

709 self.metadata["convolvedExposure"] = "Science" 

710 bbox = science.getBBox() 

711 

712 kernelImage = lsst.afw.image.ImageD(psfMatchingKernel.getDimensions()) 

713 xcen, ycen = bbox.getCenter() 

714 norm = psfMatchingKernel.computeImage(kernelImage, doNormalize=False, x=xcen, y=ycen) 

715 

716 matchedScience = self._convolveExposure(science, psfMatchingKernel, 

717 self.convolutionControl, 

718 psf=template.psf) 

719 

720 # Place back on native photometric scale 

721 matchedScience.maskedImage /= norm 

722 matchedTemplate = template.clone()[bbox] 

723 matchedTemplate.setPhotoCalib(science.photoCalib) 

724 

725 if backgroundModel is not None: 

726 # We must invert the background model if the matching kernel is solved for the science image. 

727 invertedBackground = backgroundModel.clone() 

728 invertedBackground.setParameters([-p for p in backgroundModel.getParameters()]) 

729 backgroundModel = invertedBackground 

730 

731 difference = _subtractImages(matchedScience, matchedTemplate, backgroundModel=backgroundModel) 

732 

733 correctedExposure = self.finalize(template, science, difference, 

734 psfMatchingKernel, 

735 templateMatched=False) 

736 

737 return lsst.pipe.base.Struct(difference=correctedExposure, 

738 matchedTemplate=matchedTemplate, 

739 matchedScience=matchedScience, 

740 backgroundModel=backgroundModel, 

741 psfMatchingKernel=psfMatchingKernel) 

742 

743 def finalize(self, template, science, difference, kernel, 

744 templateMatched=True, 

745 preConvMode=False, 

746 preConvKernel=None, 

747 spatiallyVarying=False): 

748 """Decorrelate the difference image to undo the noise correlations 

749 caused by convolution. 

750 

751 Parameters 

752 ---------- 

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

754 Template exposure, warped to match the science exposure. 

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

756 Science exposure to subtract from the template. 

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

758 Result of subtracting template and science. 

759 kernel : `lsst.afw.math.Kernel` 

760 An (optionally spatially-varying) PSF matching kernel 

761 templateMatched : `bool`, optional 

762 Was the template PSF-matched to the science image? 

763 preConvMode : `bool`, optional 

764 Was the science image preconvolved with its own PSF 

765 before PSF matching the template? 

766 preConvKernel : `lsst.afw.detection.Psf`, optional 

767 If not `None`, then the science image was pre-convolved with 

768 (the reflection of) this kernel. Must be normalized to sum to 1. 

769 spatiallyVarying : `bool`, optional 

770 Compute the decorrelation kernel spatially varying across the image? 

771 

772 Returns 

773 ------- 

774 correctedExposure : `lsst.afw.image.ExposureF` 

775 The decorrelated image difference. 

776 """ 

777 if self.config.doDecorrelation: 

778 self.log.info("Decorrelating image difference.") 

779 # We have cleared the template mask plane, so copy the mask plane of 

780 # the image difference so that we can calculate correct statistics 

781 # during decorrelation 

782 correctedExposure = self.decorrelate.run(science, template[science.getBBox()], difference, kernel, 

783 templateMatched=templateMatched, 

784 preConvMode=preConvMode, 

785 preConvKernel=preConvKernel, 

786 spatiallyVarying=spatiallyVarying).correctedExposure 

787 else: 

788 self.log.info("NOT decorrelating image difference.") 

789 correctedExposure = difference 

790 return correctedExposure 

791 

792 def _calculateMagLim(self, exposure, nsigma=5.0, fallbackPsfSize=None): 

793 """Calculate an exposure's limiting magnitude. 

794 

795 This method uses the photometric zeropoint together with the 

796 PSF size from the average position of the exposure. 

797 

798 Parameters 

799 ---------- 

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

801 The target exposure to calculate the limiting magnitude for. 

802 nsigma : `float`, optional 

803 The detection threshold in sigma. 

804 fallbackPsfSize : `float`, optional 

805 PSF FWHM to use in the event the exposure PSF cannot be retrieved. 

806 

807 Returns 

808 ------- 

809 maglim : `astropy.units.Quantity` 

810 The limiting magnitude of the exposure, or np.nan. 

811 """ 

812 if exposure.photoCalib is None: 

813 return np.nan 

814 try: 

815 psf = exposure.getPsf() 

816 psf_shape = psf.computeShape(psf.getAveragePosition()) 

817 except (lsst.pex.exceptions.InvalidParameterError, 

818 afwDetection.InvalidPsfError, 

819 lsst.pex.exceptions.RangeError): 

820 if fallbackPsfSize is None: 

821 self.log.info("Unable to evaluate PSF, setting maglim to nan") 

822 return np.nan 

823 self.log.info("Unable to evaluate PSF, using fallback FWHM %f", fallbackPsfSize) 

824 psf_area = np.pi*(fallbackPsfSize/2)**2 

825 else: 

826 # Get a more accurate area than `psf_shape.getArea()` via moments 

827 psf_area = np.pi*np.sqrt(psf_shape.getIxx()*psf_shape.getIyy()) 

828 

829 zeropoint = exposure.photoCalib.instFluxToMagnitude(1) 

830 return zeropoint - 2.5*np.log10(nsigma*np.sqrt(psf_area)) 

831 

832 @staticmethod 

833 def _validateExposures(template, science): 

834 """Check that the WCS of the two Exposures match, the template bbox 

835 contains the science bbox, and that the bands match. 

836 

837 Parameters 

838 ---------- 

839 template : `lsst.afw.image.ExposureF` 

840 Template exposure, warped to match the science exposure. 

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

842 Science exposure to subtract from the template. 

843 

844 Raises 

845 ------ 

846 AssertionError 

847 Raised if the WCS of the template is not equal to the science WCS, 

848 if the science image is not fully contained in the template 

849 bounding box, or if the bands do not match. 

850 """ 

851 assert template.wcs == science.wcs, \ 

852 "Template and science exposure WCS are not identical." 

853 templateBBox = template.getBBox() 

854 scienceBBox = science.getBBox() 

855 assert science.filter.bandLabel == template.filter.bandLabel, \ 

856 "Science and template exposures have different bands: %s, %s" % \ 

857 (science.filter, template.filter) 

858 

859 assert templateBBox.contains(scienceBBox), \ 

860 "Template bbox does not contain all of the science image." 

861 

862 def _convolveExposure(self, exposure, kernel, convolutionControl, 

863 bbox=None, 

864 psf=None, 

865 photoCalib=None, 

866 interpolateBadMaskPlanes=False, 

867 ): 

868 """Convolve an exposure with the given kernel. 

869 

870 Parameters 

871 ---------- 

872 exposure : `lsst.afw.Exposure` 

873 exposure to convolve. 

874 kernel : `lsst.afw.math.LinearCombinationKernel` 

875 PSF matching kernel computed in the ``makeKernel`` subtask. 

876 convolutionControl : `lsst.afw.math.ConvolutionControl` 

877 Configuration for convolve algorithm. 

878 bbox : `lsst.geom.Box2I`, optional 

879 Bounding box to trim the convolved exposure to. 

880 psf : `lsst.afw.detection.Psf`, optional 

881 Point spread function (PSF) to set for the convolved exposure. 

882 photoCalib : `lsst.afw.image.PhotoCalib`, optional 

883 Photometric calibration of the convolved exposure. 

884 interpolateBadMaskPlanes : `bool`, optional 

885 If set, interpolate over mask planes specified in 

886 ``config.badMaskPlanes`` before convolving the image. 

887 

888 Returns 

889 ------- 

890 convolvedExp : `lsst.afw.Exposure` 

891 The convolved image. 

892 """ 

893 convolvedExposure = exposure.clone() 

894 if psf is not None: 

895 convolvedExposure.setPsf(psf) 

896 if photoCalib is not None: 

897 convolvedExposure.setPhotoCalib(photoCalib) 

898 if interpolateBadMaskPlanes and self.config.badMaskPlanes is not None: 

899 nInterp = _interpolateImage(convolvedExposure.maskedImage, 

900 self.config.badMaskPlanes) 

901 self.metadata["nInterpolated"] = nInterp 

902 convolvedImage = lsst.afw.image.MaskedImageF(convolvedExposure.getBBox()) 

903 lsst.afw.math.convolve(convolvedImage, convolvedExposure.maskedImage, kernel, convolutionControl)