Coverage for tests/test_psfs.py: 65%

1# This file is part of lsst-images. 

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# Use of this source code is governed by a 3-clause BSD-style 

10# license that can be found in the LICENSE file. 

11 

12from __future__ import annotations 

13 

14import os 

15import warnings 

16from typing import Any 

17 

18import numpy as np 

19import pytest 

20 

21from lsst.images import Box, Image 

22from lsst.images.describe import DescribableMixin, Report 

23from lsst.images.psfs import GaussianPointSpreadFunction, PiffWrapper, PointSpreadFunction, PSFExWrapper 

24from lsst.images.psfs._piff import _ArchivePiffWriter 

25from lsst.images.tests import ( 

26 RoundtripFits, 

27 RoundtripJson, 

28 RoundtripNdf, 

29 compare_psf_to_legacy, 

30 reset_afw_mask_planes, # noqa: F401 

31) 

32 

33try: 

34 import h5py # noqa: F401 

35 

36 HAVE_H5PY = True 

37except ImportError: 

38 HAVE_H5PY = False 

39 

40try: 

41 from lsst.afw.detection import Psf as LegacyPsf 

42except ImportError: 

43 type LegacyPsf = Any # type: ignore[no-redef] 

44 

45EXTERNAL_DATA_DIR = os.environ.get("TESTDATA_IMAGES_DIR", None) 

46 

47skip_no_h5py = pytest.mark.skipif(not HAVE_H5PY, reason="h5py is not installed") 

48 

49 

50@pytest.fixture 

51def legacy_piff_psf_and_bbox(reset_afw_mask_planes: None) -> tuple[LegacyPsf, Box]: # noqa: F811 

52 """Return a legacy-wrapped Piff PSF and its bounding box. 

53 

54 Skips if TESTDATA_IMAGES_DIR is unset, piff is unavailable, or afw is 

55 unavailable. 

56 """ 

57 # reset_afw_mask_planes will have already skipped if afw is not available. 

58 from lsst.afw.image import ExposureFitsReader 

59 

60 if EXTERNAL_DATA_DIR is None: 60 ↛ 62line 60 didn't jump to line 62 because the condition on line 60 was always true

61 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.") 

62 try: 

63 import piff # noqa: F401 

64 except ImportError: 

65 pytest.skip("'piff' could not be imported.") 

66 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "visit_image.fits") 

67 reader = ExposureFitsReader(filename) 

68 legacy_psf = reader.readPsf() 

69 bounds = Box.from_legacy(reader.readBBox()) 

70 return legacy_psf, bounds 

71 

72 

73@pytest.fixture 

74def legacy_psfex_psf_and_bbox(reset_afw_mask_planes: None) -> tuple[LegacyPsf, Box]: # noqa: F811 

75 """Return a legacy PSFEx PSF and its bounding box 

76 

77 Skips if TESTDATA_IMAGES_DIR is unset, afw is unavailable, or psfex is 

78 unavailable. 

79 """ 

80 if EXTERNAL_DATA_DIR is None: 80 ↛ 82line 80 didn't jump to line 82 because the condition on line 80 was always true

81 pytest.skip("TESTDATA_IMAGES_DIR is not in the environment.") 

82 try: 

83 from lsst.afw.image import ExposureFitsReader 

84 from lsst.meas.extensions.psfex import PsfexPsf # noqa: F401 

85 except ImportError: 

86 pytest.skip("'lsst.afw.image' or 'lsst.meas.extensions.psfex' could not be imported.") 

87 filename = os.path.join(EXTERNAL_DATA_DIR, "dp2", "legacy", "preliminary_visit_image.fits") 

88 reader = ExposureFitsReader(filename) 

89 legacy_psf = reader.readPsf() 

90 bounds = Box.from_legacy(reader.readBBox()) 

91 return legacy_psf, bounds 

92 

93 

94def test_gaussian_psf_repr_pinned() -> None: 

95 """GaussianPointSpreadFunction repr matches its documented form.""" 

96 psf = GaussianPointSpreadFunction(2.5, bounds=Box.factory[0:10, 0:10], stamp_size=21) 

97 assert repr(psf) == f"GaussianPointSpreadFunction(2.5, stamp_size=21, bounds={psf.bounds!r})" 

98 

99 

100def test_gaussian_psf_describe() -> None: 

101 """GaussianPointSpreadFunction._describe returns a Report.""" 

102 psf = GaussianPointSpreadFunction(2.5, bounds=Box.factory[0:10, 0:10], stamp_size=21) 

103 assert isinstance(psf, DescribableMixin) 

104 report = psf._describe() 

105 assert isinstance(report, Report) 

106 assert report.type_name == "GaussianPointSpreadFunction" 

107 labels = {f.label for f in report.fields} 

108 assert "sigma" in labels 

109 assert "stamp_size" in labels 

110 assert "bounds" in labels 

111 

112 

113def test_psf_base_describe() -> None: 

114 """Base _describe uses type(self).__name__ and includes bounds/kernel_bbox. 

115 

116 GaussianPointSpreadFunction overrides _describe, so a minimal concrete 

117 subclass is used to exercise the actual base implementation. 

118 """ 

119 

120 class _MinimalPSF(PointSpreadFunction): 

121 """Minimal concrete subclass that exercises the base _describe path.""" 

122 

123 @property 

124 def bounds(self) -> Box: 

125 return Box.factory[-5:5, -5:5] 

126 

127 @property 

128 def kernel_bbox(self) -> Box: 

129 return Box.factory[-2:3, -2:3] 

130 

131 def compute_kernel_image(self, *, x: float, y: float) -> Image: 

132 arr = np.zeros((5, 5)) 

133 arr[2, 2] = 1.0 

134 return Image(arr, bbox=self.kernel_bbox) 

135 

136 def compute_stellar_image(self, *, x: float, y: float) -> Image: 

137 return self.compute_kernel_image(x=x, y=y) 

138 

139 def compute_stellar_bbox(self, *, x: float, y: float) -> Box: 

140 return self.kernel_bbox 

141 

142 psf = _MinimalPSF() 

143 report = psf._describe() 

144 assert report.type_name == "_MinimalPSF" 

145 labels = {f.label for f in report.fields} 

146 assert "bounds" in labels 

147 assert "kernel_bbox" in labels 

148 

149 

150def test_gaussian() -> None: 

151 """Test the built-in Gaussian PSF implementation.""" 

152 bounds = Box.factory[-1024:1024, -2048:2048] 

153 psf = GaussianPointSpreadFunction(2.5, bounds=bounds, stamp_size=33) 

154 assert psf.bounds == bounds 

155 

156 kernel = psf.compute_kernel_image(x=5.0, y=3.0) 

157 assert kernel.bbox == psf.kernel_bbox 

158 assert abs(float(kernel.array.sum()) - 1.0) < 1e-6 

159 center = kernel.array.shape[0] // 2 

160 assert np.unravel_index(np.argmax(kernel.array), kernel.array.shape) == (center, center) 

161 

162 stellar = psf.compute_stellar_image(x=5.25, y=3.75) 

163 assert stellar.bbox == psf.compute_stellar_bbox(x=5.25, y=3.75) 

164 assert abs(float(stellar.array.sum()) - 1.0) < 1e-6 

165 assert stellar.array[center - 1, center] > stellar.array[center + 1, center] 

166 assert stellar.array[center, center] > stellar.array[center, center - 1] 

167 assert stellar.array[center, center] > stellar.array[center - 1, center] 

168 

169 with RoundtripFits(psf) as roundtrip: 

170 assert roundtrip.result == psf, f"{roundtrip.result} != {psf}" 

171 

172 with pytest.raises(ValueError, match="stamp_size must be odd"): 

173 # Even stamp size. 

174 GaussianPointSpreadFunction(2.5, bounds=bounds, stamp_size=32) 

175 

176 with pytest.raises(ValueError, match="stamp_size must be positive"): 

177 # Negative stamp size. 

178 GaussianPointSpreadFunction(2.5, bounds=bounds, stamp_size=-33) 

179 

180 with pytest.raises(ValueError, match="sigma must be positive"): 

181 # Negative sigma. 

182 GaussianPointSpreadFunction(-2.5, bounds=bounds, stamp_size=33) 

183 

184 

185def test_gaussian_from_legacy(reset_afw_mask_planes: None) -> None: # noqa: F811 

186 """Test converting a legacy GaussianPsf to and from 

187 GaussianPointSpreadFunction. 

188 """ 

189 # reset_afw_mask_planes will have already skipped if afw is not available. 

190 from lsst.afw.detection import GaussianPsf 

191 from lsst.geom import Point2D 

192 

193 bounds = Box.factory[0:64, 0:48] 

194 legacy_psf = GaussianPsf(33, 33, 2.5) 

195 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

196 assert isinstance(psf, GaussianPointSpreadFunction) 

197 assert psf == GaussianPointSpreadFunction(2.5, bounds=bounds, stamp_size=33) 

198 # The kernels match to round-off, so compare_psf_to_legacy (which checks 

199 # exact equality) cannot be used. 

200 for x, y in [(10.0, 10.0), (10.3, 20.7)]: 

201 legacy_kernel = Image.from_legacy(legacy_psf.computeKernelImage(Point2D(x, y))) 

202 kernel = psf.compute_kernel_image(x=x, y=y) 

203 assert kernel.bbox == legacy_kernel.bbox 

204 np.testing.assert_allclose(kernel.array, legacy_kernel.array, rtol=0, atol=1e-15) 

205 legacy_psf_2 = psf.to_legacy() 

206 assert isinstance(legacy_psf_2, GaussianPsf) 

207 assert legacy_psf_2.getSigma() == legacy_psf.getSigma() 

208 assert legacy_psf_2.getDimensions() == legacy_psf.getDimensions() 

209 

210 with pytest.raises(ValueError, match="must be square"): 

211 PointSpreadFunction.from_legacy(GaussianPsf(33, 31, 2.5), bounds) 

212 

213 

214def test_piff_writer_normalizes_tuple_metadata(): # intentionally untyped 

215 """Test that Piff metadata is normalized to JSON-like values.""" 

216 writer = _ArchivePiffWriter() 

217 writer.write_struct( 

218 "interp", 

219 { 

220 "keys": ("u", "v"), 

221 "scale": np.float64(1.5), 

222 "flags": [np.bool_(True), np.int64(3)], 

223 }, 

224 ) 

225 model = writer.serialize(None) # type: ignore[arg-type] 

226 assert model.structs["interp"]["keys"] == ["u", "v"] 

227 assert model.structs["interp"]["scale"] == 1.5 

228 assert model.structs["interp"]["flags"] == [True, 3] 

229 with warnings.catch_warnings(): 

230 warnings.simplefilter("error", UserWarning) 

231 model.model_dump_json() 

232 

233 

234def test_piff(legacy_piff_psf_and_bbox: tuple[LegacyPsf, Box]) -> None: 

235 """Test round-tripping a legacy Piff PSF through FITS and JSON archives, 

236 and converting it back to a legacy PSF. 

237 """ 

238 from piff import PSF 

239 

240 legacy_psf, bounds = legacy_piff_psf_and_bbox 

241 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

242 assert isinstance(psf, PiffWrapper) 

243 assert psf.bounds == bounds 

244 assert isinstance(psf.piff_psf, PSF) 

245 compare_psf_to_legacy(psf, legacy_psf) 

246 with RoundtripFits(psf) as roundtrip1: 

247 pass 

248 compare_psf_to_legacy(roundtrip1.result, legacy_psf) 

249 with RoundtripJson(psf) as roundtrip2: 

250 pass 

251 compare_psf_to_legacy(roundtrip2.result, legacy_psf) 

252 legacy_psf_2 = roundtrip1.result.to_legacy() 

253 compare_psf_to_legacy(psf, legacy_psf_2) 

254 assert legacy_psf.getAveragePosition() == legacy_psf_2.getAveragePosition() 

255 

256 

257@skip_no_h5py 

258def test_piff_ndf_roundtrip(legacy_piff_psf_and_bbox: tuple[LegacyPsf, Box]) -> None: 

259 """Test round-tripping a legacy Piff PSF through an NDF archive.""" 

260 legacy_psf, bounds = legacy_piff_psf_and_bbox 

261 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

262 with RoundtripNdf(psf) as roundtrip: 

263 pass 

264 compare_psf_to_legacy(roundtrip.result, legacy_psf) 

265 

266 

267def test_psfex(legacy_psfex_psf_and_bbox: tuple[LegacyPsf, Box]) -> None: 

268 """Test wrapping a legacy PSFEx PSF and round-tripping through FITS and 

269 JSON. 

270 """ 

271 from lsst.meas.extensions.psfex import PsfexPsf 

272 

273 legacy_psf, bounds = legacy_psfex_psf_and_bbox 

274 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

275 assert isinstance(psf, PSFExWrapper) 

276 assert psf.bounds == bounds 

277 assert isinstance(psf.legacy_psf, PsfexPsf) 

278 compare_psf_to_legacy(psf, legacy_psf) 

279 with RoundtripFits(psf) as roundtrip1: 

280 pass 

281 compare_psf_to_legacy(roundtrip1.result, legacy_psf) 

282 with RoundtripJson(psf) as roundtrip2: 

283 pass 

284 compare_psf_to_legacy(roundtrip2.result, legacy_psf)