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