Coverage for tests/test_extended_psf.py: 99%

346 statements  

« prev     ^ index     » next       coverage.py v7.15.2, created at 2026-07-23 09:30 +0000

1# This file is part of pipe_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 

22import logging 

23import unittest 

24from typing import cast 

25 

26import astropy.units as u 

27import lsst.images.tests 

28import lsst.utils.tests 

29import numpy as np 

30from astropy.table import Table 

31from lsst.afw.cameraGeom import FOCAL_PLANE, PIXELS 

32from lsst.afw.cameraGeom.testUtils import CameraWrapper 

33from lsst.afw.geom import makeCdMatrix, makeSkyWcs 

34from lsst.afw.image import ( 

35 ExposureF, 

36 ImageD, 

37 ImageF, 

38 MaskedImageF, 

39 VisitInfo, 

40 makePhotoCalibFromCalibZeroPoint, 

41) 

42from lsst.afw.math import ( 

43 REDUCE_INTERP_ORDER, 

44 ApproximateControl, 

45 BackgroundControl, 

46 BackgroundList, 

47 FixedKernel, 

48 Interpolate, 

49 makeBackground, 

50) 

51from lsst.geom import Box2I, Extent2I, Point2D, Point2I, SpherePoint, arcseconds, degrees 

52from lsst.images import Box, Image, Mask, MaskedImage, MaskSchema 

53from lsst.meas.algorithms import KernelPsf 

54from lsst.pipe.tasks.extended_psf import ( 

55 ExtendedPsfCandidate, 

56 ExtendedPsfCandidateInfo, 

57 ExtendedPsfCandidates, 

58 ExtendedPsfCutoutConfig, 

59 ExtendedPsfCutoutTask, 

60 ExtendedPsfFit, 

61 ExtendedPsfImage, 

62 ExtendedPsfImageInfo, 

63 ExtendedPsfMoffatFit, 

64 ExtendedPsfStackConfig, 

65 ExtendedPsfStackTask, 

66 ExtendedPsfSubtractConfig, 

67 ExtendedPsfSubtractTask, 

68) 

69 

70 

71class ExtendedPsfCandidatesTestCase(lsst.utils.tests.TestCase): 

72 """Test ExtendedPsfCandidate and ExtendedPsfCandidates data structures.""" 

73 

74 @classmethod 

75 def setUpClass(cls): 

76 rng = np.random.Generator(np.random.MT19937(seed=5)) 

77 cutout_size = (25, 25) 

78 

79 # Generate simulated stars 

80 extended_psf_candidates = [] 

81 cls.star_infos = [] 

82 mask_schema = MaskSchema([]) 

83 

84 for i in range(3): 

85 candidate_masked_image = MaskedImage( 

86 image=Image(rng.random(cutout_size).astype(np.float32)), 

87 mask=Mask(0, schema=mask_schema, shape=cutout_size), 

88 variance=Image(rng.random(cutout_size).astype(np.float32)), 

89 ) 

90 

91 star_info = ExtendedPsfCandidateInfo( 

92 visit=100 + i, 

93 detector=200 + i, 

94 ref_id=1000 + i, 

95 ref_mag=10.0 + i, 

96 position_x=float(rng.random()), 

97 position_y=float(rng.random()), 

98 focal_plane_radius=float(rng.random()) * u.mm, 

99 focal_plane_angle=float(rng.random()) * u.rad, 

100 ) 

101 cls.star_infos.append(star_info) 

102 

103 # Build a normalized 2-D Gaussian kernel image directly. 

104 yy, xx = np.mgrid[: cutout_size[0], : cutout_size[1]] 

105 cy = (cutout_size[0] - 1) / 2.0 

106 cx = (cutout_size[1] - 1) / 2.0 

107 sigma = 1.5 

108 kernel_array = np.exp(-((yy - cy) ** 2 + (xx - cx) ** 2) / (2.0 * sigma**2)) 

109 kernel_array /= np.sum(kernel_array) 

110 psf_kernel_image = Image(kernel_array.astype(np.float64)) 

111 

112 extended_psf_candidates.append( 

113 ExtendedPsfCandidate( 

114 image=candidate_masked_image.image, 

115 mask=candidate_masked_image.mask, 

116 variance=candidate_masked_image.variance, 

117 psf_kernel_image=psf_kernel_image, 

118 star_info=star_info, 

119 ) 

120 ) 

121 

122 cls.global_metadata = {"CANDIDATES_TEST_KEY": "CANDIDATES_TEST_VALUE"} 

123 cls.extended_psf_candidates = ExtendedPsfCandidates(extended_psf_candidates, cls.global_metadata) 

124 

125 @classmethod 

126 def tearDownClass(cls): 

127 del cls.extended_psf_candidates 

128 del cls.star_infos 

129 del cls.global_metadata 

130 

131 def test_fits_roundtrip(self): 

132 """Test that ExtendedPsfCandidates can be serialized/deserialized.""" 

133 

134 with lsst.images.tests.RoundtripFits( 

135 self.extended_psf_candidates, storage_class="ExtendedPsfCandidates" 

136 ) as roundtrip: 

137 pass 

138 extended_psf_candidates = roundtrip.result 

139 

140 global_metadata = extended_psf_candidates.metadata 

141 self.assertEqual(self.global_metadata["CANDIDATES_TEST_KEY"], global_metadata["CANDIDATES_TEST_KEY"]) 

142 self.assertEqual(len(self.extended_psf_candidates), len(extended_psf_candidates)) 

143 for input_info, input_candidate, output_candidate in zip( 

144 self.star_infos, self.extended_psf_candidates, extended_psf_candidates 

145 ): 

146 self.assertEqual(input_candidate.metadata, {}) 

147 self.assertEqual(output_candidate.metadata, {}) 

148 

149 output_info = output_candidate.star_info 

150 self.assertEqual(input_info.visit, output_info.visit) 

151 self.assertEqual(input_info.detector, output_info.detector) 

152 self.assertEqual(input_info.ref_id, output_info.ref_id) 

153 self.assertAlmostEqual(input_info.ref_mag, output_info.ref_mag, places=10) 

154 self.assertAlmostEqual(input_info.position_x, output_info.position_x, places=10) 

155 self.assertAlmostEqual(input_info.position_y, output_info.position_y, places=10) 

156 

157 self.assertIsNotNone(input_info.focal_plane_radius) 

158 self.assertIsNotNone(output_info.focal_plane_radius) 

159 input_radius = cast(u.Quantity, input_info.focal_plane_radius) 

160 output_radius = cast(u.Quantity, output_info.focal_plane_radius) 

161 self.assertAlmostEqual(input_radius.to_value(u.mm), output_radius.to_value(u.mm), places=10) 

162 

163 self.assertIsNotNone(input_info.focal_plane_angle) 

164 self.assertIsNotNone(output_info.focal_plane_angle) 

165 input_angle = cast(u.Quantity, input_info.focal_plane_angle) 

166 output_angle = cast(u.Quantity, output_info.focal_plane_angle) 

167 self.assertAlmostEqual(input_angle.to_value(u.rad), output_angle.to_value(u.rad), places=10) 

168 

169 np.testing.assert_allclose( 

170 input_candidate.image.array, 

171 output_candidate.image.array, 

172 rtol=0.0, 

173 atol=1e-10, 

174 ) 

175 np.testing.assert_array_equal( 

176 input_candidate.mask.array, 

177 output_candidate.mask.array, 

178 ) 

179 np.testing.assert_allclose( 

180 input_candidate.variance.array, output_candidate.variance.array, rtol=0.0, atol=1e-10 

181 ) 

182 np.testing.assert_allclose( 

183 input_candidate.psf_kernel_image.array, 

184 output_candidate.psf_kernel_image.array, 

185 rtol=0.0, 

186 atol=1e-10, 

187 ) 

188 

189 def test_ref_id_map(self): 

190 """Test that ref_id_map maps reference catalog IDs to candidates.""" 

191 ref_id_map = self.extended_psf_candidates.ref_id_map 

192 self.assertEqual(len(ref_id_map), len(self.star_infos)) 

193 for i, star_info in enumerate(self.star_infos): 

194 self.assertIn(star_info.ref_id, ref_id_map) 

195 self.assertIs(ref_id_map[star_info.ref_id], self.extended_psf_candidates[i]) 

196 

197 def test_slice_preserves_subcollection_and_metadata(self): 

198 """Test that slicing candidates returns a sub-collection.""" 

199 sub = self.extended_psf_candidates[1:3] # length=2 

200 self.assertIsInstance(sub, ExtendedPsfCandidates) 

201 self.assertEqual(len(sub), 2) 

202 self.assertEqual(sub.metadata, self.global_metadata) 

203 self.assertIs(sub[0], self.extended_psf_candidates[1]) 

204 self.assertIs(sub[1], self.extended_psf_candidates[2]) 

205 

206 def test_spatial_box_slicing(self): 

207 """Test that candidates support spatial slicing with a Box.""" 

208 candidate = self.extended_psf_candidates[0] 

209 subbox = Box.from_shape((10, 10), start=(5, 5)) 

210 sub = candidate[subbox] 

211 self.assertIsInstance(sub, ExtendedPsfCandidate) 

212 self.assertEqual(sub.bbox, subbox) 

213 self.assertIs(sub.star_info, candidate.star_info) 

214 np.testing.assert_array_equal(sub.image.array, candidate.image[subbox].array) 

215 np.testing.assert_array_equal(sub.mask.array, candidate.mask[subbox].array) 

216 

217 def test_copy(self): 

218 """Test that copy() produces an independent deep copy.""" 

219 candidate = self.extended_psf_candidates[0] 

220 copied = candidate.copy() 

221 self.assertIsInstance(copied, ExtendedPsfCandidate) 

222 np.testing.assert_array_equal(copied.image.array, candidate.image.array) 

223 np.testing.assert_array_equal(copied.variance.array, candidate.variance.array) 

224 # Modifying the copy must not affect the original. 

225 copied.image.array[:] = 0.0 

226 self.assertFalse(np.all(candidate.image.array == 0.0)) 

227 

228 

229class ExtendedPsfImageTestCase(lsst.utils.tests.TestCase): 

230 """Test ExtendedPsfImage data model, properties, and operations.""" 

231 

232 def setUp(self): 

233 image = Image( 

234 np.arange(30, dtype=np.float32).reshape(5, 6), 

235 yx0=(10, -3), 

236 unit=u.nJy, 

237 ) 

238 variance = Image( 

239 np.full((5, 6), 2.5, dtype=np.float32), 

240 bbox=image.bbox, 

241 unit=u.nJy**2, 

242 ) 

243 info = ExtendedPsfImageInfo( 

244 n_stars=17, 

245 ) 

246 fit = ExtendedPsfMoffatFit( 

247 chi2=12.3, 

248 reduced_chi2=1.23, 

249 dof=10, 

250 amplitude=0.8, 

251 x_0=-0.2, 

252 y_0=0.4, 

253 gamma=2.3, 

254 alpha=4.5, 

255 ) 

256 self.extended_psf_image = ExtendedPsfImage( 

257 image=image, 

258 variance=variance, 

259 info=info, 

260 fit=fit, 

261 metadata={"EPSF_TEST_KEY": "EPSF_TEST VALUE"}, 

262 ) 

263 

264 def tearDown(self): 

265 del self.extended_psf_image 

266 

267 def test_fits_roundtrip(self): 

268 """Test that ExtendedPsfImage can be serialized/deserialized.""" 

269 with lsst.images.tests.RoundtripFits( 

270 self.extended_psf_image, 

271 storage_class="ExtendedPsfImage", 

272 ) as roundtrip: 

273 subbox = Box.from_shape((3, 3), start=(11, -1)) 

274 subimage = roundtrip.get(bbox=subbox) 

275 expected_subimage = self.extended_psf_image[subbox] 

276 

277 np.testing.assert_allclose( 

278 subimage.image.array, 

279 expected_subimage.image.array, 

280 rtol=0.0, 

281 atol=1e-8, 

282 ) 

283 np.testing.assert_allclose( 

284 subimage.variance.array, 

285 expected_subimage.variance.array, 

286 rtol=0.0, 

287 atol=1e-8, 

288 ) 

289 self.assertEqual(subimage.info, self.extended_psf_image.info) 

290 self.assertEqual(subimage.fit, self.extended_psf_image.fit) 

291 

292 roundtripped = roundtrip.result 

293 np.testing.assert_allclose( 

294 roundtripped.image.array, 

295 self.extended_psf_image.image.array, 

296 rtol=0.0, 

297 atol=1e-8, 

298 ) 

299 np.testing.assert_allclose( 

300 roundtripped.variance.array, 

301 self.extended_psf_image.variance.array, 

302 rtol=0.0, 

303 atol=1e-8, 

304 ) 

305 self.assertEqual(roundtripped.info, self.extended_psf_image.info) 

306 self.assertEqual(roundtripped.fit, self.extended_psf_image.fit) 

307 self.assertEqual(roundtripped.metadata["EPSF_TEST_KEY"], "EPSF_TEST VALUE") 

308 

309 def test_unit_sky_projection_bbox_properties(self): 

310 """Test ExtendedPsfImage properties: unit, sky_projection, and bbox.""" 

311 self.assertEqual(self.extended_psf_image.unit, u.nJy) 

312 self.assertIsNone(self.extended_psf_image.sky_projection) 

313 self.assertEqual(self.extended_psf_image.bbox, self.extended_psf_image.image.bbox) 

314 

315 def test_copy_independence(self): 

316 """Test that copy() produces an independent deep copy.""" 

317 copied = self.extended_psf_image.copy() 

318 np.testing.assert_array_equal(copied.image.array, self.extended_psf_image.image.array) 

319 np.testing.assert_array_equal(copied.variance.array, self.extended_psf_image.variance.array) 

320 self.assertEqual(copied.info, self.extended_psf_image.info) 

321 self.assertEqual(copied.fit, self.extended_psf_image.fit) 

322 # Modifying the copy must not affect the original. 

323 copied.image.array[:] = 0.0 

324 self.assertFalse(np.all(self.extended_psf_image.image.array == 0.0)) 

325 

326 def test_setitem_subregion_assignment(self): 

327 """Test __setitem__ writes image and variance into a subregion.""" 

328 target = ExtendedPsfImage(Image(np.zeros((5, 6), dtype=np.float32), yx0=(10, -3), unit=u.nJy)) 

329 subbox = Box.from_shape((3, 3), start=(11, -1)) 

330 target[subbox] = self.extended_psf_image[subbox] 

331 np.testing.assert_array_equal( 

332 target[subbox].image.array, 

333 self.extended_psf_image[subbox].image.array, 

334 ) 

335 np.testing.assert_array_equal( 

336 target[subbox].variance.array, 

337 self.extended_psf_image[subbox].variance.array, 

338 ) 

339 

340 def test_constructor_invalid_inputs(self): 

341 """Test that ValueError is raised for invalid constructor inputs.""" 

342 # Mismatched bboxes 

343 image = Image(np.ones((5, 6), dtype=np.float32)) 

344 variance = Image(np.ones((4, 6), dtype=np.float32)) 

345 with self.assertRaises(ValueError): 

346 ExtendedPsfImage(image=image, variance=variance) 

347 

348 # Variance has wrong units (nJy instead of nJy**2). 

349 image = Image(np.ones((5, 6), dtype=np.float32), unit=u.nJy) 

350 variance_wrong_units = Image(np.ones((5, 6), dtype=np.float32), unit=u.nJy) 

351 with self.assertRaises(ValueError): 

352 ExtendedPsfImage(image=image, variance=variance_wrong_units) 

353 

354 # Image has no units but variance does. 

355 image_no_units = Image(np.ones((5, 6), dtype=np.float32)) 

356 variance_with_units = Image(np.ones((5, 6), dtype=np.float32), unit=u.nJy**2) 

357 with self.assertRaises(ValueError): 

358 ExtendedPsfImage(image=image_no_units, variance=variance_with_units) 

359 

360 

361class ExtendedPsfFitTestCase(lsst.utils.tests.TestCase): 

362 """Test ExtendedPsfFit.""" 

363 

364 def test_construction(self): 

365 """Test construction with and without optional fields.""" 

366 # With all fields 

367 fit = ExtendedPsfFit(chi2=10.5, dof=10, reduced_chi2=1.05) 

368 self.assertEqual(fit.chi2, 10.5) 

369 self.assertEqual(fit.dof, 10) 

370 self.assertEqual(fit.reduced_chi2, 1.05) 

371 

372 # With only the required field 

373 fit_minimal = ExtendedPsfFit(chi2=10.5) 

374 self.assertEqual(fit_minimal.chi2, 10.5) 

375 self.assertIsNone(fit_minimal.dof) 

376 self.assertIsNone(fit_minimal.reduced_chi2) 

377 

378 

379class ExtendedPsfSubtractionTestCase(lsst.utils.tests.TestCase): 

380 """Integration tests for all extended PSF subtraction pipeline tasks.""" 

381 

382 @classmethod 

383 def setUpClass(cls): 

384 rng = np.random.default_rng(42) 

385 

386 # Background coefficients 

387 sigma = 60.0 # noise 

388 pedestal = 3210.1 

389 coef_x = 1e-2 

390 coef_y = 2e-2 

391 coef_x2 = 1e-5 

392 coef_xy = 2e-5 

393 coef_y2 = 3e-5 

394 

395 # Make an input exposure 

396 wcs = makeSkyWcs( 

397 crpix=Point2D(0, 0), 

398 crval=SpherePoint(0, 0, degrees), 

399 cdMatrix=makeCdMatrix(scale=0.2 * arcseconds, flipX=True), 

400 ) 

401 cls.exposure = ExposureF(1001, 1001, wcs) 

402 cls.exposure.setPhotoCalib(makePhotoCalibFromCalibZeroPoint(10 ** (0.4 * 30), 1.0)) 

403 ny, nx = cls.exposure.image.array.shape 

404 grid_y, grid_x = np.mgrid[(-ny + 1) // 2 : ny // 2 + 1, (-nx + 1) // 2 : nx // 2 + 1] 

405 cls.exposure.image.array[:] += rng.normal(scale=sigma, size=cls.exposure.image.array.shape) 

406 cls.background_array = np.zeros_like(cls.exposure.image.array) 

407 cls.background_array += pedestal 

408 cls.background_array += coef_x * grid_x 

409 cls.background_array += coef_y * grid_y 

410 cls.background_array += coef_x2 * grid_x * grid_x 

411 cls.background_array += coef_xy * grid_x * grid_y 

412 cls.background_array += coef_y2 * grid_y * grid_y 

413 cls.exposure.image.array += cls.background_array 

414 cls.exposure.info.setVisitInfo(VisitInfo(id=12345)) 

415 camera = CameraWrapper().camera 

416 detector = camera[10] 

417 cls.exposure.setDetector(detector) 

418 for mask_plane in [ 

419 "BAD", 

420 "CR", 

421 "CROSSTALK", 

422 "EDGE", 

423 "NO_DATA", 

424 "SAT", 

425 "SUSPECT", 

426 "UNMASKEDNAN", 

427 "NEIGHBOR", 

428 ]: 

429 _ = cls.exposure.mask.addMaskPlane(mask_plane) 

430 cls.exposure.variance.array.fill(1.0) 

431 

432 # Make a background list 

433 background_image = ImageF(cls.exposure.getBBox()) 

434 background_image.array[:] = cls.background_array 

435 background_control = BackgroundControl(Interpolate.NATURAL_SPLINE) 

436 background_control.setNxSample(background_image.getWidth()) 

437 background_control.setNySample(background_image.getHeight()) 

438 background_model = makeBackground(background_image, background_control) 

439 cls.background = BackgroundList() 

440 cls.background.append( 

441 ( 

442 background_model, 

443 Interpolate.NATURAL_SPLINE, 

444 REDUCE_INTERP_ORDER, 

445 ApproximateControl.UNKNOWN, 

446 0, 

447 0, 

448 False, 

449 ) 

450 ) 

451 

452 # Make a table of extended PSF candidate stars 

453 corners = cls.exposure.wcs.pixelToSky([Point2D(x) for x in cls.exposure.getBBox().getCorners()]) 

454 corner_ras = [corner.getRa().asDegrees() for corner in corners] 

455 corner_decs = [corner.getDec().asDegrees() for corner in corners] 

456 num_stars = 10 

457 ras = rng.uniform(np.min(corner_ras), np.max(corner_ras), num_stars) 

458 decs = rng.uniform(np.min(corner_decs), np.max(corner_decs), num_stars) 

459 sky_coords = [SpherePoint(ra, dec, degrees) for ra, dec in zip(ras, decs)] 

460 pixel_coords = cls.exposure.wcs.skyToPixel(sky_coords) 

461 pixel_x = [coord.getX() for coord in pixel_coords] 

462 pixel_y = [coord.getY() for coord in pixel_coords] 

463 mags = rng.uniform(10, 20, num_stars) 

464 fluxes = [cls.exposure.photoCalib.magnitudeToInstFlux(mag) for mag in mags] 

465 mm_coords = detector.transform(pixel_coords, PIXELS, FOCAL_PLANE) 

466 mm_coords_x = np.array([mm_coord.x for mm_coord in mm_coords]) 

467 mm_coords_y = np.array([mm_coord.y for mm_coord in mm_coords]) 

468 radius_mm = np.sqrt(mm_coords_x**2 + mm_coords_y**2) 

469 theta_radians = np.arctan2(mm_coords_y, mm_coords_x) 

470 cls.extended_psf_candidate_table = Table( 

471 { 

472 "id": np.arange(num_stars), 

473 "coord_ra": ras, 

474 "coord_dec": decs, 

475 "phot_g_mean_flux": fluxes, 

476 "mag": mags, 

477 "pixel_x": pixel_x, 

478 "pixel_y": pixel_y, 

479 "radius_mm": radius_mm, 

480 "angle_radians": theta_radians, 

481 } 

482 ) 

483 

484 # Make a synthetic star 

485 cutout_radius = 25 

486 grid_y, grid_x = np.mgrid[-cutout_radius : cutout_radius + 1, -cutout_radius : cutout_radius + 1] 

487 dist_from_center = np.sqrt(grid_x**2 + grid_y**2) 

488 sigma = 1.5 

489 psf_array = np.exp(-(dist_from_center**2) / (2 * sigma**2)) 

490 psf_array /= np.sum(psf_array) 

491 fixed_kernel = FixedKernel(ImageD(psf_array)) 

492 kernel_psf = KernelPsf(fixed_kernel) 

493 cls.exposure.setPsf(kernel_psf) 

494 psf = kernel_psf.computeKernelImage(kernel_psf.getAveragePosition()) 

495 

496 # Add synthetic stars to the exposure 

497 footprints = ImageF(cls.exposure.getBBox()) 

498 cls.injected_stars_array = np.zeros_like(cls.exposure.image.array) 

499 for candidate_id, candidate in enumerate(cls.extended_psf_candidate_table): 

500 bbox_star = Box2I(Point2I(candidate["pixel_x"], candidate["pixel_y"]), Extent2I(1, 1)).dilatedBy( 

501 cutout_radius 

502 ) 

503 bbox_star_clipped = bbox_star.clippedTo(cls.exposure.getBBox()) 

504 candidate_im = MaskedImageF(bbox_star) 

505 candidate_im.image.array = candidate["phot_g_mean_flux"] * psf.getArray() 

506 candidate_im = candidate_im[bbox_star_clipped] 

507 detection_threshold = cls.exposure.getPhotoCalib().magnitudeToInstFlux(25) 

508 detected = candidate_im.image.array > detection_threshold 

509 footprints[bbox_star_clipped].array = detected * (candidate_id + 1) 

510 _ = candidate_im.mask.addMaskPlane("DETECTED") 

511 candidate_im.mask.array[detected] |= candidate_im.mask.getPlaneBitMask("DETECTED") 

512 candidate_im.variance.array.fill(1.0) 

513 y_slice = slice(bbox_star_clipped.getMinY(), bbox_star_clipped.getMaxY() + 1) 

514 x_slice = slice(bbox_star_clipped.getMinX(), bbox_star_clipped.getMaxX() + 1) 

515 cls.injected_stars_array[y_slice, x_slice] += candidate_im.image.array 

516 cls.exposure.maskedImage[bbox_star_clipped] += candidate_im 

517 cls.footprints = footprints.array 

518 

519 # Run the cutout task 

520 extendedPsfCutoutConfig = ExtendedPsfCutoutConfig() 

521 extendedPsfCutoutTask = ExtendedPsfCutoutTask(config=extendedPsfCutoutConfig) 

522 cls.extended_psf_candidates = extendedPsfCutoutTask._get_extended_psf_candidates( 

523 input_exposure=cls.exposure, 

524 input_background=None, 

525 footprints=cls.footprints, 

526 extended_psf_candidate_table=cls.extended_psf_candidate_table, 

527 ) 

528 

529 # Run the stack task 

530 extendedPsfStackConfig = ExtendedPsfStackConfig() 

531 extendedPsfStackTask = ExtendedPsfStackTask(config=extendedPsfStackConfig) 

532 stack_result = extendedPsfStackTask.run(extended_psf_candidates=cls.extended_psf_candidates) 

533 cls.extended_psf = stack_result.extended_psf if stack_result is not None else None 

534 

535 @classmethod 

536 def tearDownClass(cls): 

537 del cls.exposure 

538 del cls.background 

539 del cls.background_array 

540 del cls.injected_stars_array 

541 del cls.extended_psf_candidate_table 

542 del cls.footprints 

543 del cls.extended_psf_candidates 

544 del cls.extended_psf 

545 

546 def test_cutout_task_candidate_extraction(self): 

547 """Test ExtendedPsfCutoutTask candidate extraction.""" 

548 assert self.extended_psf_candidates is not None 

549 self.assertAlmostEqual( 

550 float(self.extended_psf_candidates.metadata["FOCAL_PLANE_RADIUS_MM_MIN"]), 5.22408977, 7 

551 ) 

552 self.assertAlmostEqual( 

553 float(self.extended_psf_candidates.metadata["FOCAL_PLANE_RADIUS_MM_MAX"]), 14.6045832, 7 

554 ) 

555 self.assertEqual(len(self.extended_psf_candidates), len(self.extended_psf_candidate_table)) 

556 self.assertEqual(self.extended_psf_candidates[0].star_info.visit, 12345) 

557 self.assertEqual(self.extended_psf_candidates[0].star_info.detector, 10) 

558 

559 for candidate, candidate_entry in zip( 

560 self.extended_psf_candidates, self.extended_psf_candidate_table 

561 ): 

562 self.assertEqual(candidate.star_info.ref_id, candidate_entry["id"]) 

563 self.assertEqual(candidate.star_info.ref_mag, candidate_entry["mag"]) 

564 self.assertEqual(candidate.star_info.position_x, candidate_entry["pixel_x"]) 

565 self.assertEqual(candidate.star_info.position_y, candidate_entry["pixel_y"]) 

566 self.assertIsInstance(candidate.psf_kernel_image, Image) 

567 self.assertEqual(candidate.psf_kernel_image.array.ndim, 2) 

568 self.assertGreater(candidate.psf_kernel_image.array.shape[0], 0) 

569 self.assertGreater(candidate.psf_kernel_image.array.shape[1], 0) 

570 self.assertTrue(np.isfinite(candidate.psf_kernel_image.array).all()) 

571 self.assertGreater(candidate.psf_kernel_image.array.sum(), 0.0) 

572 focal_plane_radius = cast(u.Quantity, candidate.star_info.focal_plane_radius) 

573 focal_plane_angle = cast(u.Quantity, candidate.star_info.focal_plane_angle) 

574 self.assertEqual(focal_plane_radius.to_value(u.mm), candidate_entry["radius_mm"]) 

575 self.assertEqual(focal_plane_angle.to_value(u.rad), candidate_entry["angle_radians"]) 

576 

577 def test_stack_task_moffat_fitting(self): 

578 """Test Moffat fitting.""" 

579 assert self.extended_psf is not None 

580 self.assertAlmostEqual(np.sum(self.extended_psf.image.array), 0.8233417, places=2) 

581 self.assertAlmostEqual(np.sum(self.extended_psf.variance.array), 0.0075618913, places=4) 

582 fit = self.extended_psf.fit 

583 self.assertAlmostEqual(fit.chi2, 107652.97393353, delta=5) 

584 self.assertEqual(fit.dof, 62996) 

585 self.assertAlmostEqual(fit.reduced_chi2, 1.7088858647141, places=2) 

586 self.assertAlmostEqual(fit.amplitude, 0.078900464260488, places=2) 

587 self.assertAlmostEqual(fit.x_0, -0.68834523633912, places=2) 

588 self.assertAlmostEqual(fit.y_0, -0.069005412739451, places=2) 

589 self.assertAlmostEqual(fit.gamma, 8.0966823485900, places=2) 

590 self.assertAlmostEqual(fit.alpha, 16.048683662812, places=2) 

591 

592 def test_stack_task_no_candidates(self): 

593 """Test that None returned when no candidates pass the radius check.""" 

594 config = ExtendedPsfStackConfig() 

595 config.min_focal_plane_radius = 1e6 # Excludes all candidates. 

596 task = ExtendedPsfStackTask(config=config) 

597 with self.assertLogs(level=logging.INFO): 

598 result = task.run(extended_psf_candidates=self.extended_psf_candidates) 

599 self.assertIsNone(result) 

600 

601 def test_subtract_task(self): 

602 """Test subtraction task on all synthetic stars. 

603 

604 This test validates subtraction against known injected-star flux 

605 while exercising default background restore and re-estimation. 

606 """ 

607 

608 assert self.extended_psf is not None 

609 

610 config = ExtendedPsfSubtractConfig() 

611 config.bad_mask_planes = [] 

612 task = ExtendedPsfSubtractTask(config=config) 

613 

614 preliminary_visit_image = self.exposure.clone() 

615 preliminary_visit_image.image.array -= self.background_array 

616 restored_input_image = self.exposure.image.array.copy() 

617 

618 star_table = self.extended_psf_candidate_table 

619 

620 with unittest.mock.patch.object(task, "_get_subtraction_star_table", return_value=star_table): 

621 result = task.run( 

622 preliminary_visit_image=preliminary_visit_image, 

623 preliminary_visit_image_background=self.background, 

624 extended_psf=self.extended_psf, 

625 ref_obj_loader=object(), 

626 ) 

627 

628 output_exposure = result.preliminary_visit_image_extended_psf_subtracted 

629 output_background = result.preliminary_visit_image_extended_psf_subtracted_background 

630 self.assertIsNotNone(output_background) 

631 

632 restored_output_image = output_exposure.image.array + output_background.getImage().array 

633 removed_image = restored_input_image - restored_output_image 

634 expected_flux = float(np.sum(self.injected_stars_array)) 

635 removed_flux = float(np.sum(removed_image)) 

636 self.assertFloatsAlmostEqual(removed_flux, expected_flux, rtol=0.25) 

637 

638 metadata = output_exposure.getMetadata() 

639 n_stars = len(star_table) 

640 self.assertEqual(metadata.getScalar("EPSFSUB_ATTEMPTED"), n_stars) 

641 self.assertEqual(metadata.getScalar("EPSFSUB_SUBTRACTED"), n_stars) 

642 

643 

644class MemoryTestCase(lsst.utils.tests.MemoryTestCase): 

645 pass 

646 

647 

648def setup_module(module): 

649 lsst.utils.tests.init() 

650 

651 

652if __name__ == "__main__": 652 ↛ 653line 652 didn't jump to line 653 because the condition on line 652 was never true

653 lsst.utils.tests.init() 

654 unittest.main()