Coverage for tests/test_psfs.py: 60%

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_piff_writer_normalizes_tuple_metadata(): # intentionally untyped 

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

187 writer = _ArchivePiffWriter() 

188 writer.write_struct( 

189 "interp", 

190 { 

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

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

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

194 }, 

195 ) 

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

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

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

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

200 with warnings.catch_warnings(): 

201 warnings.simplefilter("error", UserWarning) 

202 model.model_dump_json() 

203 

204 

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

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

207 and converting it back to a legacy PSF. 

208 """ 

209 from piff import PSF 

210 

211 legacy_psf, bounds = legacy_piff_psf_and_bbox 

212 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

213 assert isinstance(psf, PiffWrapper) 

214 assert psf.bounds == bounds 

215 assert isinstance(psf.piff_psf, PSF) 

216 compare_psf_to_legacy(psf, legacy_psf) 

217 with RoundtripFits(psf) as roundtrip1: 

218 pass 

219 compare_psf_to_legacy(roundtrip1.result, legacy_psf) 

220 with RoundtripJson(psf) as roundtrip2: 

221 pass 

222 compare_psf_to_legacy(roundtrip2.result, legacy_psf) 

223 legacy_psf_2 = roundtrip1.result.to_legacy() 

224 compare_psf_to_legacy(psf, legacy_psf_2) 

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

226 

227 

228@skip_no_h5py 

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

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

231 legacy_psf, bounds = legacy_piff_psf_and_bbox 

232 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

233 with RoundtripNdf(psf) as roundtrip: 

234 pass 

235 compare_psf_to_legacy(roundtrip.result, legacy_psf) 

236 

237 

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

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

240 JSON. 

241 """ 

242 from lsst.meas.extensions.psfex import PsfexPsf 

243 

244 legacy_psf, bounds = legacy_psfex_psf_and_bbox 

245 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

246 assert isinstance(psf, PSFExWrapper) 

247 assert psf.bounds == bounds 

248 assert isinstance(psf.legacy_psf, PsfexPsf) 

249 compare_psf_to_legacy(psf, legacy_psf) 

250 with RoundtripFits(psf) as roundtrip1: 

251 pass 

252 compare_psf_to_legacy(roundtrip1.result, legacy_psf) 

253 with RoundtripJson(psf) as roundtrip2: 

254 pass 

255 compare_psf_to_legacy(roundtrip2.result, legacy_psf)