Coverage for python / lsst / drp / tasks / tract_background.py: 0%

199 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-05-04 17:41 +0000

1# This file is part of drp_tasks. 

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 __future__ import annotations 

23 

24__all__ = ( 

25 "BinnedTractGeometry", 

26 "MeasureTractBackgroundConfig", 

27 "MeasureTractBackgroundConnections", 

28 "MeasureTractBackgroundTask", 

29 "SubtractTractBackgroundConfig", 

30 "SubtractTractBackgroundConnections", 

31 "SubtractTractBackgroundTask", 

32) 

33 

34import dataclasses 

35from collections.abc import Iterable, Mapping 

36from typing import ClassVar 

37 

38import numpy as np 

39 

40from lsst.afw.geom import SkyWcs, SpanSet, makeModifiedWcs, makeTransform 

41from lsst.afw.image import ExposureF, FilterLabel, MaskX 

42from lsst.afw.math import ( 

43 BackgroundList, 

44 binImage, 

45) 

46from lsst.geom import AffineTransform, Box2I, Extent2I, LinearTransform, Point2I 

47from lsst.meas.algorithms import SubtractBackgroundTask 

48from lsst.pex.config import ConfigurableField, Field, ListField 

49from lsst.pipe.base import ( 

50 InputQuantizedConnection, 

51 OutputQuantizedConnection, 

52 PipelineTask, 

53 PipelineTaskConfig, 

54 PipelineTaskConnections, 

55 QuantumContext, 

56 Struct, 

57) 

58from lsst.pipe.base import connectionTypes as cT 

59from lsst.skymap import BaseSkyMap, TractInfo 

60 

61 

62@dataclasses.dataclass 

63class BinnedTractGeometry: 

64 """Geometry for binned tract-level images.""" 

65 

66 tract_info: TractInfo 

67 """Definition of the tract.""" 

68 

69 original_bbox: Box2I 

70 """The original bounding box of the tract level image in true pixel 

71 coordinates. 

72 

73 Since this is constructed from a union of patch-level bounding boxes that 

74 are in general between the patch inner and outer boundin boxes, it is not 

75 the same as the tract's usual bounding box. 

76 """ 

77 

78 superpixel_size: Extent2I 

79 """The number of true pixels in each superpixel in each dimension.""" 

80 

81 binned_bbox: Box2I 

82 """The bounding box of the tract-level image in superpixel coordinates.""" 

83 

84 binned_wcs: SkyWcs 

85 """A WCS that maps the superpixel image to the sky.""" 

86 

87 def slices(self, child_bbox: Box2I) -> tuple[slice, slice]: 

88 """Return the array slices that select the given pixel-coordinate 

89 bounding box within the superpixel tract image. 

90 

91 Parameters 

92 ---------- 

93 child_bbox : `lsst.geom.Box2I` 

94 True-pixel bounding box to select. 

95 

96 Returns 

97 ------- 

98 y_slice : `slice` 

99 Superpixel image slice in the y dimension. 

100 x_slice : `slice` 

101 Superpixel image slice in the x dimension. 

102 """ 

103 assert child_bbox.y.begin % self.superpixel_size.y == 0, "superpixel size evenly divides bboxes" 

104 assert child_bbox.x.begin % self.superpixel_size.x == 0, "superpixel size evenly divides bboxes" 

105 assert child_bbox.y.end % self.superpixel_size.y == 0, "superpixel size evenly divides bboxes" 

106 assert child_bbox.x.end % self.superpixel_size.x == 0, "superpixel size evenly divides bboxes" 

107 return ( 

108 slice( 

109 (child_bbox.y.begin - self.original_bbox.y.begin) // self.superpixel_size.y, 

110 (child_bbox.y.end - self.original_bbox.y.begin) // self.superpixel_size.y, 

111 ), 

112 slice( 

113 (child_bbox.x.begin - self.original_bbox.x.begin) // self.superpixel_size.x, 

114 (child_bbox.x.end - self.original_bbox.x.begin) // self.superpixel_size.x, 

115 ), 

116 ) 

117 

118 

119class MeasureTractBackgroundConnections(PipelineTaskConnections, dimensions=["tract"]): 

120 sky_map = cT.Input( 

121 doc="Input definition of geometry/bbox and projection/wcs for warps.", 

122 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

123 storageClass="SkyMap", 

124 dimensions=("skymap",), 

125 ) 

126 hf_backgrounds = cT.Input( 

127 "deep_coadd_hf_background", 

128 doc=( 

129 "Original high-frequency background. Bins must exactly divide the inner patch region. " 

130 "This input can be produced by running DetectCoaddSourcesTask with forceExactBinning " 

131 "and writeEmptyBackgrounds both set to True." 

132 "" 

133 ), 

134 storageClass="Background", 

135 multiple=True, 

136 dimensions=["patch", "band"], 

137 ) 

138 lf_backgrounds = cT.Output( 

139 "deep_coadd_lf_background", 

140 doc=( 

141 "Remeasured low-frequency background, covering the same patches and area within those patches " 

142 "as the input high-frequency backgrounds." 

143 ), 

144 storageClass="Background", 

145 multiple=True, 

146 dimensions=["patch", "band"], 

147 ) 

148 hf_diagnostic_images = cT.Output( 

149 "deep_coadd_hf_background_image", 

150 doc="Tract-level image of the original high-frequency background superpixels.", 

151 storageClass="ExposureF", 

152 multiple=True, 

153 dimensions=["tract", "band"], 

154 ) 

155 lf_diagnostic_images = cT.Output( 

156 "deep_coadd_lf_background_image", 

157 doc="Tract-level image of the remeasured low-frequency background superpixels.", 

158 storageClass="ExposureF", 

159 multiple=True, 

160 dimensions=["tract", "band"], 

161 ) 

162 

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

164 assert isinstance(config, MeasureTractBackgroundConfig) 

165 if not config.write_diagnostic_images: 

166 del self.hf_diagnostic_images 

167 del self.lf_diagnostic_images 

168 

169 def adjustQuantum(self, inputs, outputs, label, data_id): 

170 # Make sure we have an input for every output, and vice-versa. 

171 hf_connection, hf_refs = inputs["hf_backgrounds"] 

172 lf_connection, lf_refs = outputs["lf_backgrounds"] 

173 hf_refs_by_data_id = {ref.dataId: ref for ref in hf_refs} 

174 lf_refs_by_data_id = {ref.dataId: ref for ref in lf_refs} 

175 common = sorted(hf_refs_by_data_id.keys() & lf_refs_by_data_id.keys()) 

176 inputs["hf_backgrounds"] = (hf_connection, [hf_refs_by_data_id[d] for d in common]) 

177 outputs["lf_backgrounds"] = (lf_connection, [lf_refs_by_data_id[d] for d in common]) 

178 return super().adjustQuantum(inputs, outputs, label, data_id) 

179 

180 

181class MeasureTractBackgroundConfig(PipelineTaskConfig, pipelineConnections=MeasureTractBackgroundConnections): 

182 background = ConfigurableField( 

183 "Configuration for the task used to subtract the background on the tract-level " 

184 "image of superpixels. The bin sizes in this configuration are in superpixel units, " 

185 "and the interpolation/approximation parameter control how the superpixel background is " 

186 "interpolated back onto the input high-frequency background bins. ", 

187 target=SubtractBackgroundTask, 

188 ) 

189 mask_grow_input_planes = ListField( 

190 "Mask planes to union across bands, grow by 'mask_grow_radius', and set as " 

191 "'mask_grow_output_plane'. Superpixel values that are NaN are masked as NO_DATA before this " 

192 "operation.", 

193 dtype=str, 

194 default=["NO_DATA"], 

195 ) 

196 mask_grow_radius = Field( 

197 "Radius to grow the background subtraction mask by (in tract-level superpixels, " 

198 "i.e. bins used in the high-frequency backgrounds).", 

199 dtype=int, 

200 default=8, 

201 ) 

202 mask_grow_output_plane = Field( 

203 "Name of the mask plane that will hold the grown union of 'mask_grow_input_planes'.", 

204 dtype=str, 

205 default="BG_IGNORE", 

206 ) 

207 write_diagnostic_images = Field( 

208 "Whether to write diagnostic image outputs.", 

209 dtype=bool, 

210 default=True, 

211 ) 

212 

213 def setDefaults(self): 

214 super().setDefaults() 

215 self.background.binSize = 8 

216 self.background.useApprox = False 

217 self.background.ignoredPixelMask.append("BG_IGNORE") 

218 

219 def validate(self): 

220 super().validate() 

221 if self.mask_grow_output_plane not in self.background.ignoredPixelMask: 

222 raise ValueError( 

223 f"mask_grow_output_plane={self.mask_grow_output_plane} " 

224 "is not used in 'background.ignoredMaskPixels'." 

225 ) 

226 

227 

228class MeasureTractBackgroundTask(PipelineTask): 

229 """A task that fits a (relatively) low-frequency background to tract-level 

230 superpixels, using high-frequency patch-level backgrounds as inputs. 

231 

232 Notes 

233 ----- 

234 This task starts from an initial round of high-frequency background 

235 estimation done in a previous task (e.g. DetectCoaddSourcesTask), combining 

236 the per-bin background estimates into a tract-level image in which each 

237 superpixel corresponds to a single bin from an input background object. 

238 The input background's configuration for how to interpolate or approximate 

239 those bin values is ignored at this stage; we use backgrounds as inputs 

240 rather than simple binned images because each background bin value is 

241 already a robust mean of non-detected pixels, which does a passable job 

242 at avoiding contamination from sources. 

243 

244 For tracts with very large objects, it is expected that the tract-level 

245 image will have many missing superpixels (i.e. with NaN values). These 

246 are used to form a mask, which is then unioned across all bands and grown 

247 by a configurable amount. That mask is then input to a round of background 

248 subtraction on the superpixel image to form the "low frequency" output 

249 background models. This background is propagated to patch-level output 

250 `lsst.afw.math.BackgroundList` objects by replacing the bin values in the 

251 input objects with the interpolated low-frequency model. The original 

252 input interpolation/approximation configuration is then used in the 

253 outputs. 

254 

255 This task runs over multiple bands at once only to enable that cross-band 

256 union of masks, which is nevertheless very important (at least in the 

257 context of making RGB images from the background-subtracted coadds, but 

258 avoiding unphysical colors is good for science productions, too). 

259 """ 

260 

261 _DefaultName: ClassVar[str] = "measureTractBackground" 

262 ConfigClass: ClassVar[type[MeasureTractBackgroundConfig]] = MeasureTractBackgroundConfig 

263 config: MeasureTractBackgroundConfig 

264 

265 def __init__(self, *, config=None, log=None, initInputs=None, **kwargs): 

266 super().__init__(config=config, log=log, initInputs=initInputs, **kwargs) 

267 self.makeSubtask("background") 

268 

269 def runQuantum( 

270 self, 

271 butlerQC: QuantumContext, 

272 inputRefs: InputQuantizedConnection, 

273 outputRefs: OutputQuantizedConnection, 

274 ): 

275 sky_map = butlerQC.get(inputRefs.sky_map) 

276 tract_info = sky_map[butlerQC.quantum.dataId["tract"]] 

277 hf_backgrounds = { 

278 (ref.dataId["band"], ref.dataId["patch"]): butlerQC.get(ref) for ref in inputRefs.hf_backgrounds 

279 } 

280 results = self.run(hf_backgrounds=hf_backgrounds, tract_info=tract_info) 

281 for ref in outputRefs.lf_backgrounds: 

282 butlerQC.put(results.lf_backgrounds[ref.dataId["band"], ref.dataId["patch"]], ref) 

283 if self.config.write_diagnostic_images: 

284 for ref in outputRefs.hf_diagnostic_images: 

285 butlerQC.put(results.hf_diagnostic_images[ref.dataId["band"]], ref) 

286 for ref in outputRefs.lf_diagnostic_images: 

287 butlerQC.put(results.lf_diagnostic_images[ref.dataId["band"]], ref) 

288 

289 def run( 

290 self, *, hf_backgrounds: Mapping[tuple[str, int], BackgroundList], tract_info: TractInfo 

291 ) -> Struct: 

292 """Estimate a tract-level background model. 

293 

294 Parameters 

295 ---------- 

296 hf_backgrounds : `~collections.abc.Mapping` 

297 Mapping of input `lsst.afw.math.BackgroundList` instances, keyed 

298 by tuples of ``(band, patch)``. Each background must have a single 

299 layer with bins that exactly divide the patch inner region (but may 

300 have bins of the same size outside the patch inner region). 

301 tract_info : `lsst.skymap.TractInfo` 

302 Geometry information for the tract. 

303 

304 Returns 

305 ------- 

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

307 Results struct with the following attributes: 

308 

309 - ``lf_backgrounds``: mapping of output 

310 `lsst.afw.math.BackgroundList` instances, with the same keys as 

311 ``hf_backgrounds``. 

312 - ``hf_diagnostic_images``: mapping of tract-level superpixel 

313 images with the input backgrounds, keyed by band. 

314 - ``lf_diagnostic_images``: mapping of tract-level superpixel 

315 images with the output backgrounds, keyed by band. 

316 - ``geometry``: `BinnedTractGeometry`: struct with the geometry 

317 of the bins. 

318 """ 

319 geometry = self.compute_geometry(hf_backgrounds=hf_backgrounds.values(), tract_info=tract_info) 

320 hf_bg_exposures = self.make_hf_background_exposures(geometry=geometry, hf_backgrounds=hf_backgrounds) 

321 common_mask = MaskX(geometry.binned_bbox) 

322 for hf_bg_exposure in hf_bg_exposures.values(): 

323 common_mask |= hf_bg_exposure.mask 

324 grow_spans = ( 

325 SpanSet.fromMask(common_mask, common_mask.getPlaneBitMask(self.config.mask_grow_input_planes)) 

326 .dilated(self.config.mask_grow_radius) 

327 .clippedTo(geometry.binned_bbox) 

328 ) 

329 grow_spans.setMask(common_mask, common_mask.getPlaneBitMask(self.config.mask_grow_output_plane)) 

330 lf_bg_exposures = {} 

331 for band, hf_bg_exposure in hf_bg_exposures.items(): 

332 hf_bg_exposure.mask = common_mask 

333 lf_bg_list = self.background.run(hf_bg_exposure).background 

334 lf_bg_exposure = hf_bg_exposure.clone() 

335 lf_bg_exposure.image = lf_bg_list.getImage() 

336 lf_bg_exposures[band] = lf_bg_exposure 

337 lf_backgrounds = self.make_lf_background_lists( 

338 lf_bg_exposures=lf_bg_exposures, geometry=geometry, hf_backgrounds=hf_backgrounds 

339 ) 

340 return Struct( 

341 lf_backgrounds=lf_backgrounds, 

342 hf_diagnostic_images=hf_bg_exposures, 

343 lf_diagnostic_images=lf_bg_exposures, 

344 geometry=geometry, 

345 ) 

346 

347 def compute_geometry( 

348 self, *, hf_backgrounds: Iterable[BackgroundList], tract_info: TractInfo 

349 ) -> BinnedTractGeometry: 

350 """Work out the geometry of a binned tract-level superpixel image. 

351 

352 Parameters 

353 ---------- 

354 hf_backgrounds : `~collections.abc.Iterable` 

355 Iterable of all input `lsst.afw.math.BackgroundList` instances that 

356 cover the tract. 

357 tract_info : `lsst.skymap.TractInfo` 

358 Geometry information for the tract. 

359 

360 Returns 

361 ------- 

362 geometry: `BinnedTractGeometry` 

363 A struct describing the bounding boxes and binning of the 

364 tract-level superpixel image. 

365 """ 

366 original_bbox = Box2I() 

367 superpixel_size: Extent2I | None = None 

368 for hf_bg_list in hf_backgrounds: 

369 patch_bbox: Box2I | None = None 

370 for hf_bg, *_ in hf_bg_list: 

371 if patch_bbox is None: 

372 patch_bbox = hf_bg.getImageBBox() 

373 else: 

374 assert patch_bbox == hf_bg.getImageBBox() 

375 stats_image_bbox = hf_bg.getStatsImage().getBBox() 

376 if superpixel_size is None: 

377 superpixel_size = Extent2I( 

378 patch_bbox.width // stats_image_bbox.width, 

379 patch_bbox.height // stats_image_bbox.height, 

380 ) 

381 assert superpixel_size.x * stats_image_bbox.width == patch_bbox.width 

382 assert superpixel_size.y * stats_image_bbox.height == patch_bbox.height 

383 original_bbox.include(patch_bbox) 

384 binned_bbox = Box2I( 

385 Point2I(original_bbox.x.begin // superpixel_size.x, original_bbox.y.begin // superpixel_size.y), 

386 Extent2I(original_bbox.width // superpixel_size.x, original_bbox.height // superpixel_size.y), 

387 ) 

388 return BinnedTractGeometry( 

389 tract_info=tract_info, 

390 original_bbox=original_bbox, 

391 superpixel_size=superpixel_size, 

392 binned_bbox=binned_bbox, 

393 binned_wcs=_make_binned_wcs(tract_info.getWcs(), superpixel_size.x, superpixel_size.y), 

394 ) 

395 

396 def make_hf_background_exposures( 

397 self, *, geometry: BinnedTractGeometry, hf_backgrounds: Mapping[tuple[str, int], BackgroundList] 

398 ) -> dict[str, ExposureF]: 

399 """Build tract-level superpixel images from the bin values in input 

400 backgrounds. 

401 

402 Parameters 

403 ---------- 

404 hf_backgrounds : `~collections.abc.Mapping` 

405 Mapping of input `lsst.afw.math.BackgroundList` instances, keyed 

406 by tuples of ``(band, patch)``. Each background must have a single 

407 layer with bins that exactly divide the patch inner region (but may 

408 have bins of the same size outside the patch inner region). 

409 geometry: `BinnedTractGeometry` 

410 A struct describing the bounding boxes and binning of the 

411 tract-level superpixel image. 

412 

413 Returns 

414 ------- 

415 hf_images : `dict` 

416 Mapping of tract-level superpixel images 

417 (`lsst.afw.image.ExposureF` objects) with the input background bin 

418 values, keyed by band. 

419 """ 

420 result = {} 

421 no_data_bitmask = MaskX.getPlaneBitMask("NO_DATA") 

422 for (band, _), hf_bg_list in hf_backgrounds.items(): 

423 if (exposure := result.get(band)) is None: 

424 exposure = ExposureF(geometry.binned_bbox) 

425 exposure.setWcs(geometry.binned_wcs) 

426 exposure.setFilter(FilterLabel.fromBand(band)) 

427 exposure.image.array[...] = np.nan 

428 exposure.mask.array[...] = no_data_bitmask 

429 exposure.variance.array[...] = 0.0 

430 exposure.mask.addMaskPlane(self.config.mask_grow_output_plane) 

431 result[band] = exposure 

432 if len(hf_bg_list) != 1: 

433 raise RuntimeError("Input BackgroundList must have exactly one layer.") 

434 hf_bg = hf_bg_list[0][0] 

435 tract_slices = geometry.slices(hf_bg.getImageBBox()) 

436 stats_image = hf_bg.getStatsImage() 

437 exposure.image.array[tract_slices] = stats_image.image.array 

438 exposure.mask.array[tract_slices] = stats_image.mask.array 

439 exposure.mask.array[tract_slices] |= no_data_bitmask * np.isnan(stats_image.image.array) 

440 exposure.variance.array[tract_slices] = stats_image.variance.array 

441 return result 

442 

443 def make_lf_background_lists( 

444 self, 

445 *, 

446 lf_bg_exposures: Mapping[str, ExposureF], 

447 geometry: BinnedTractGeometry, 

448 hf_backgrounds: Mapping[tuple[str, int], BackgroundList], 

449 ) -> dict[tuple[str, int], BackgroundList]: 

450 """Make output `lsst.afw.math.BackgroundList` objects for each patch 

451 and band. 

452 

453 Parameters 

454 ---------- 

455 lf_bg_exposures : `~collections.abc.Mapping` 

456 Mapping from band to binned tract-level superpixel images 

457 (`lsst.afw.image.ExposureF` instance) that hold the low-frequency 

458 background model fit to the tract-level superpixel image for that 

459 band. 

460 geometry: `BinnedTractGeometry` 

461 A struct describing the bounding boxes and binning of the 

462 tract-level superpixel image. 

463 hf_backgrounds : `~collections.abc.Mapping` 

464 Mapping of input `lsst.afw.math.BackgroundList` instances, keyed 

465 by tuples of ``(band, patch)``. THESE ARE MODIFIED IN PLACE AND 

466 SHOULD BE CONSIDERED CONSUMED. 

467 

468 Returns 

469 ------- 

470 lf_backgrounds : `dict` 

471 A dictionary of output low-frequency `lsst.afw.math.BackgroundList` 

472 instances, keyed by tuples of ``(band, patch)``. These use the 

473 same interpolation/approximation parameters as the input 

474 high-frequency background objects, but have their bin values 

475 replaced by evaluating the low-frequency background model for the 

476 tract. 

477 """ 

478 result = {} 

479 for (band, patch_id), hf_bg_list in hf_backgrounds.items(): 

480 lf_bg_list = BackgroundList() 

481 lf_bg_list.append(hf_bg_list[0]) 

482 lf_bg, *_ = lf_bg_list[0] 

483 stats_image = lf_bg.getStatsImage() 

484 slices = geometry.slices(lf_bg.getImageBBox()) 

485 stats_image.image.array[...] = lf_bg_exposures[band].image.array[slices] 

486 stats_image.mask.array[...] = 0 

487 stats_image.variance.array[...] = 1.0 

488 result[band, patch_id] = lf_bg_list 

489 return result 

490 

491 

492class SubtractTractBackgroundConnections(PipelineTaskConnections, dimensions=["patch", "band"]): 

493 input_coadd = cT.Input( 

494 "deep_coadd_predetection", 

495 doc="Original coadd with only large-scale visit-level background subtraction.", 

496 storageClass="ExposureF", 

497 dimensions=["patch", "band"], 

498 ) 

499 lf_background = cT.Input( 

500 "deep_coadd_lf_background", 

501 doc="Remeasured low-frequency background, covering inner patch regions only.", 

502 storageClass="Background", 

503 dimensions=["patch", "band"], 

504 ) 

505 output_coadd = cT.Output( 

506 "deep_coadd_lf_subtracted", 

507 doc="Output coadd with tract-level backgrounds subtracted.", 

508 storageClass="ExposureF", 

509 dimensions=["patch", "band"], 

510 ) 

511 input_coadd_binned = cT.Output( 

512 "deep_coadd_predetection_binned", 

513 doc="Binned version of 'input_coadd`, with bin size set to the superpixel size.", 

514 storageClass="ExposureF", 

515 dimensions=["patch", "band"], 

516 ) 

517 output_coadd_binned = cT.Output( 

518 "deep_coadd_lf_subtracted_binned", 

519 doc="Binned version of 'output_coadd`, with bin size set to the superpixel size.", 

520 storageClass="ExposureF", 

521 dimensions=["patch", "band"], 

522 ) 

523 

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

525 assert isinstance(config, SubtractTractBackgroundConfig) 

526 if not config.write_binned_images: 

527 del self.input_coadd_binned 

528 del self.output_coadd_binned 

529 

530 

531class SubtractTractBackgroundConfig( 

532 PipelineTaskConfig, pipelineConnections=SubtractTractBackgroundConnections 

533): 

534 write_binned_images = Field( 

535 "Whether to write out a binned versions of the input and output coadds.", 

536 dtype=bool, 

537 default=True, 

538 ) 

539 

540 

541class SubtractTractBackgroundTask(PipelineTask): 

542 """A task that subtracts a previously-fit background model from a coadd. 

543 

544 Notes 

545 ----- 

546 This task is designed to apply the background models fit by 

547 `MeasureTractBackgroundTask`. It also produces diagnostic outputs that are 

548 useful for developing and configuring that task. 

549 """ 

550 

551 _DefaultName: ClassVar[str] = "subtractTractBackground" 

552 ConfigClass: ClassVar[type[SubtractTractBackgroundConfig]] = SubtractTractBackgroundConfig 

553 config: SubtractTractBackgroundConfig 

554 

555 def run(self, *, input_coadd: ExposureF, lf_background: BackgroundList) -> Struct: 

556 """Subtract a background model from a coadd and optionally write 

557 binned versions of the original and background-subtracted images. 

558 

559 Parameters 

560 ---------- 

561 input_coadd : `lsst.afw.image.ExposureF` 

562 Input coadd to subtract a background from. 

563 lf_background : `lsst.afw.math.BackgroundList` 

564 Background model to subtract. May only cover a subset of the given 

565 coadd. 

566 

567 Returns 

568 ------- 

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

570 Output struct containing: 

571 

572 - ``output_coadd`` (`lsst.afw.image.ExposureF`): 

573 The output background-subtracted coadd. If ``lf_background`` 

574 only covers a subset of ``input_coadd``, this will be cropped 

575 to just that subset. 

576 - ``input_coadd_binned`` (``lsst.afw.image.ExposureF``): 

577 A binned version of the input coadd. Only produced if 

578 `~SubtractTractBackgroundConfig.write_binned_images` is `True`. 

579 - ``output_coadd_binned`` (``lsst.afw.image.ExposureF``): 

580 A binned version of the output coadd. Only produced if 

581 `~SubtractTractBackgroundConfig.write_binned_images` is `True`. 

582 """ 

583 for bg, *_ in lf_background: 

584 bg_bbox = bg.getImageBBox() 

585 bg_binned_bbox = bg.getStatsImage().getBBox() 

586 bin_size_x = bg_bbox.width // bg_binned_bbox.width 

587 bin_size_y = bg_bbox.height // bg_binned_bbox.height 

588 break 

589 else: 

590 raise AssertionError("No background in background list.") 

591 input_coadd = input_coadd[bg_bbox] 

592 result = Struct() 

593 if self.config.write_binned_images: 

594 result.input_coadd_binned = ExposureF(binImage(input_coadd.maskedImage, bin_size_x, bin_size_y)) 

595 result.input_coadd_binned.setXY0( 

596 Point2I(input_coadd.getX0() // bin_size_x, input_coadd.getY0() // bin_size_y) 

597 ) 

598 result.input_coadd_binned.setWcs(_make_binned_wcs(input_coadd.getWcs(), bin_size_x, bin_size_y)) 

599 result.input_coadd_binned.setFilter(input_coadd.getFilter()) 

600 bg_image = lf_background.getImage() 

601 result.output_coadd = input_coadd.clone() 

602 result.output_coadd.image -= bg_image 

603 if self.config.write_binned_images: 

604 result.output_coadd_binned = result.input_coadd_binned.clone() 

605 result.output_coadd_binned.image = binImage(result.output_coadd.image, bin_size_x, bin_size_y) 

606 return result 

607 

608 

609def _make_binned_wcs(original: SkyWcs, bin_size_x: int, bin_size_y: int) -> SkyWcs: 

610 """Make a WCS appropriate for a binned version of an image. 

611 

612 Parameters 

613 ---------- 

614 original : `lsst.afw.geom.SkyWcs` 

615 Original WCS for the unbinned image. 

616 bin_size_x : `int` 

617 Size of each bin in the X dimension. 

618 bin_size_y : `int` 

619 Size of each bin in the Y dimension. 

620 

621 Returns 

622 ------- 

623 binned_wcs : `lsst.afw.geom.SkyWcs` 

624 WCS for the binned image. 

625 

626 Notes 

627 ----- 

628 This function assumes that the original image's XY0 has been divided by 

629 the bin size in each dimension (or is ``(0, 0)`` so this is irrelevant). 

630 """ 

631 binned_to_original = makeTransform(AffineTransform(LinearTransform.makeScaling(bin_size_x, bin_size_y))) 

632 return makeModifiedWcs(binned_to_original, original, False)