Coverage for tests / test_extended_psf_subtraction.py: 18%
119 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-25 08:39 +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/>.
22import unittest
23from typing import cast
25import astropy.units as u
26import lsst.utils.tests
27import numpy as np
28from astropy.table import Table
29from lsst.afw.cameraGeom import FOCAL_PLANE, PIXELS
30from lsst.afw.cameraGeom.testUtils import CameraWrapper
31from lsst.afw.geom import makeCdMatrix, makeSkyWcs
32from lsst.afw.image import (
33 ExposureF,
34 ImageD,
35 ImageF,
36 MaskedImageF,
37 VisitInfo,
38 makePhotoCalibFromCalibZeroPoint,
39)
40from lsst.afw.math import FixedKernel
41from lsst.geom import Box2I, Extent2I, Point2D, Point2I, SpherePoint, arcseconds, degrees
42from lsst.images import Image
43from lsst.meas.algorithms import KernelPsf
44from lsst.pipe.tasks.extended_psf import ExtendedPsfCutoutConfig, ExtendedPsfCutoutTask
47class ExtendedPsfSubtractionTestCase(lsst.utils.tests.TestCase):
48 def setUp(self):
49 rng = np.random.default_rng(42)
51 # Background coefficients
52 sigma = 60.0 # noise
53 pedestal = 3210.1
54 coef_x = 1e-2
55 coef_y = 2e-2
56 coef_x2 = 1e-5
57 coef_xy = 2e-5
58 coef_y2 = 3e-5
60 # Make an input exposure
61 wcs = makeSkyWcs(
62 crpix=Point2D(0, 0),
63 crval=SpherePoint(0, 0, degrees),
64 cdMatrix=makeCdMatrix(scale=0.2 * arcseconds, flipX=True),
65 )
66 self.exposure = ExposureF(1001, 1001, wcs)
67 self.exposure.setPhotoCalib(makePhotoCalibFromCalibZeroPoint(10 ** (0.4 * 30), 1.0))
68 ny, nx = self.exposure.image.array.shape
69 grid_y, grid_x = np.mgrid[(-ny + 1) // 2 : ny // 2 + 1, (-nx + 1) // 2 : nx // 2 + 1]
70 self.exposure.image.array[:] += rng.normal(scale=sigma, size=self.exposure.image.array.shape)
71 self.exposure.image.array += pedestal
72 self.exposure.image.array += coef_x * grid_x
73 self.exposure.image.array += coef_y * grid_y
74 self.exposure.image.array += coef_x2 * grid_x * grid_x
75 self.exposure.image.array += coef_xy * grid_x * grid_y
76 self.exposure.image.array += coef_y2 * grid_y * grid_y
77 self.exposure.info.setVisitInfo(VisitInfo(id=12345))
78 camera = CameraWrapper().camera
79 detector = camera[10]
80 self.exposure.setDetector(detector)
81 for mask_plane in [
82 "BAD",
83 "CR",
84 "CROSSTALK",
85 "EDGE",
86 "NO_DATA",
87 "SAT",
88 "SUSPECT",
89 "UNMASKEDNAN",
90 "NEIGHBOR",
91 ]:
92 _ = self.exposure.mask.addMaskPlane(mask_plane)
93 self.exposure.variance.array.fill(1.0)
95 # Make a table of extended PSF candidate stars
96 corners = self.exposure.wcs.pixelToSky([Point2D(x) for x in self.exposure.getBBox().getCorners()])
97 corner_ras = [corner.getRa().asDegrees() for corner in corners]
98 corner_decs = [corner.getDec().asDegrees() for corner in corners]
99 num_stars = 10
100 ras = rng.uniform(np.min(corner_ras), np.max(corner_ras), num_stars)
101 decs = rng.uniform(np.min(corner_decs), np.max(corner_decs), num_stars)
102 sky_coords = [SpherePoint(ra, dec, degrees) for ra, dec in zip(ras, decs)]
103 pixel_coords = self.exposure.wcs.skyToPixel(sky_coords)
104 pixel_x = [coord.getX() for coord in pixel_coords]
105 pixel_y = [coord.getY() for coord in pixel_coords]
106 mags = rng.uniform(10, 20, num_stars)
107 fluxes = [self.exposure.photoCalib.magnitudeToInstFlux(mag) for mag in mags]
108 mm_coords = detector.transform(pixel_coords, PIXELS, FOCAL_PLANE)
109 mm_coords_x = np.array([mm_coord.x for mm_coord in mm_coords])
110 mm_coords_y = np.array([mm_coord.y for mm_coord in mm_coords])
111 radius_mm = np.sqrt(mm_coords_x**2 + mm_coords_y**2)
112 theta_radians = np.arctan2(mm_coords_y, mm_coords_x)
113 self.extended_psf_candidate_table = Table(
114 {
115 "id": np.arange(num_stars),
116 "coord_ra": ras,
117 "coord_dec": decs,
118 "phot_g_mean_flux": fluxes,
119 "mag": mags,
120 "pixel_x": pixel_x,
121 "pixel_y": pixel_y,
122 "radius_mm": radius_mm,
123 "angle_radians": theta_radians,
124 }
125 )
127 # Make a synthetic star
128 cutout_radius = 25
129 grid_y, grid_x = np.mgrid[-cutout_radius : cutout_radius + 1, -cutout_radius : cutout_radius + 1]
130 dist_from_center = np.sqrt(grid_x**2 + grid_y**2)
131 sigma = 1.5
132 psf_array = np.exp(-(dist_from_center**2) / (2 * sigma**2))
133 psf_array /= np.sum(psf_array)
134 fixed_kernel = FixedKernel(ImageD(psf_array))
135 kernel_psf = KernelPsf(fixed_kernel)
136 self.exposure.setPsf(kernel_psf)
137 psf = kernel_psf.computeKernelImage(kernel_psf.getAveragePosition())
139 # Add synthetic stars to the exposure
140 footprints = ImageF(self.exposure.getBBox())
141 for candidate_id, candidate in enumerate(self.extended_psf_candidate_table):
142 bbox_star = Box2I(Point2I(candidate["pixel_x"], candidate["pixel_y"]), Extent2I(1, 1)).dilatedBy(
143 cutout_radius
144 )
145 bbox_star_clipped = bbox_star.clippedTo(self.exposure.getBBox())
146 candidate_im = MaskedImageF(bbox_star)
147 candidate_im.image.array = candidate["phot_g_mean_flux"] * psf.getArray()
148 candidate_im = candidate_im[bbox_star_clipped]
149 detection_threshold = self.exposure.getPhotoCalib().magnitudeToInstFlux(25)
150 detected = candidate_im.image.array > detection_threshold
151 footprints[bbox_star_clipped].array = detected * (candidate_id + 1)
152 _ = candidate_im.mask.addMaskPlane("DETECTED")
153 candidate_im.mask.array[detected] |= candidate_im.mask.getPlaneBitMask("DETECTED")
154 candidate_im.variance.array.fill(1.0)
155 self.exposure.maskedImage[bbox_star_clipped] += candidate_im
156 self.footprints = footprints.array
158 # Run the cutout task
159 extendedPsfCutoutConfig = ExtendedPsfCutoutConfig()
160 extendedPsfCutoutTask = ExtendedPsfCutoutTask(config=extendedPsfCutoutConfig)
161 self.extended_psf_candidates = extendedPsfCutoutTask._get_extended_psf_candidates(
162 input_exposure=self.exposure,
163 input_background=None,
164 footprints=self.footprints,
165 extended_psf_candidate_table=self.extended_psf_candidate_table,
166 )
168 def test_extendedPsfCutout(self):
169 """Test ExtendedPsfCutoutTask."""
170 assert self.extended_psf_candidates is not None
171 self.assertAlmostEqual(
172 float(self.extended_psf_candidates.metadata["FOCAL_PLANE_RADIUS_MM_MIN"]), 5.22408977, 7
173 )
174 self.assertAlmostEqual(
175 float(self.extended_psf_candidates.metadata["FOCAL_PLANE_RADIUS_MM_MAX"]), 14.6045832, 7
176 )
177 self.assertEqual(len(self.extended_psf_candidates), len(self.extended_psf_candidate_table))
178 self.assertEqual(self.extended_psf_candidates[0].star_info.visit, 12345)
179 self.assertEqual(self.extended_psf_candidates[0].star_info.detector, 10)
181 for candidate, candidate_entry in zip(
182 self.extended_psf_candidates, self.extended_psf_candidate_table
183 ):
184 self.assertEqual(candidate.star_info.ref_id, candidate_entry["id"])
185 self.assertEqual(candidate.star_info.ref_mag, candidate_entry["mag"])
186 self.assertEqual(candidate.star_info.position_x, candidate_entry["pixel_x"])
187 self.assertEqual(candidate.star_info.position_y, candidate_entry["pixel_y"])
188 self.assertIsInstance(candidate.psf_kernel_image, Image)
189 self.assertEqual(candidate.psf_kernel_image.array.ndim, 2)
190 self.assertGreater(candidate.psf_kernel_image.array.shape[0], 0)
191 self.assertGreater(candidate.psf_kernel_image.array.shape[1], 0)
192 self.assertTrue(np.isfinite(candidate.psf_kernel_image.array).all())
193 self.assertGreater(candidate.psf_kernel_image.array.sum(), 0.0)
194 focal_plane_radius = cast(u.Quantity, candidate.star_info.focal_plane_radius)
195 focal_plane_angle = cast(u.Quantity, candidate.star_info.focal_plane_angle)
196 self.assertEqual(focal_plane_radius.to_value(u.mm), candidate_entry["radius_mm"])
197 self.assertEqual(focal_plane_angle.to_value(u.rad), candidate_entry["angle_radians"])
200def setup_module(module):
201 lsst.utils.tests.init()
204class MemoryTestCase(lsst.utils.tests.MemoryTestCase):
205 pass
208if __name__ == "__main__": 208 ↛ 209line 208 didn't jump to line 209 because the condition on line 208 was never true
209 lsst.utils.tests.init()
210 unittest.main()