Coverage for tests/test_psfs.py: 59%

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 RoundtripFits, RoundtripJson, RoundtripNdf, compare_psf_to_legacy 

26 

27try: 

28 import h5py # noqa: F401 

29 

30 HAVE_H5PY = True 

31except ImportError: 

32 HAVE_H5PY = False 

33 

34try: 

35 from lsst.afw.detection import Psf as LegacyPsf 

36except ImportError: 

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

38 

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

40 

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

42 

43 

44@pytest.fixture(scope="session") 

45def legacy_piff_psf_and_bbox() -> tuple[LegacyPsf, Box]: 

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

47 

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

49 unavailable. 

50 """ 

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

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

53 try: 

54 import piff # noqa: F401 

55 

56 from lsst.afw.image import ExposureFitsReader 

57 except ImportError: 

58 pytest.skip("'piff' or 'lsst.afw.image' could not be imported.") 

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

60 reader = ExposureFitsReader(filename) 

61 legacy_psf = reader.readPsf() 

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

63 return legacy_psf, bounds 

64 

65 

66@pytest.fixture(scope="session") 

67def legacy_psfex_psf_and_bbox() -> tuple[LegacyPsf, Box]: 

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

69 

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

71 unavailable. 

72 """ 

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

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

75 try: 

76 from lsst.afw.image import ExposureFitsReader 

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

78 except ImportError: 

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

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

81 reader = ExposureFitsReader(filename) 

82 legacy_psf = reader.readPsf() 

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

84 return legacy_psf, bounds 

85 

86 

87def test_gaussian_psf_repr_pinned() -> None: 

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

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

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

91 

92 

93def test_gaussian_psf_describe() -> None: 

94 """GaussianPointSpreadFunction._describe returns a Report.""" 

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

96 assert isinstance(psf, DescribableMixin) 

97 report = psf._describe() 

98 assert isinstance(report, Report) 

99 assert report.type_name == "GaussianPointSpreadFunction" 

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

101 assert "sigma" in labels 

102 assert "stamp_size" in labels 

103 assert "bounds" in labels 

104 

105 

106def test_psf_base_describe() -> None: 

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

108 

109 GaussianPointSpreadFunction overrides _describe, so a minimal concrete 

110 subclass is used to exercise the actual base implementation. 

111 """ 

112 

113 class _MinimalPSF(PointSpreadFunction): 

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

115 

116 @property 

117 def bounds(self) -> Box: 

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

119 

120 @property 

121 def kernel_bbox(self) -> Box: 

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

123 

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

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

126 arr[2, 2] = 1.0 

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

128 

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

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

131 

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

133 return self.kernel_bbox 

134 

135 psf = _MinimalPSF() 

136 report = psf._describe() 

137 assert report.type_name == "_MinimalPSF" 

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

139 assert "bounds" in labels 

140 assert "kernel_bbox" in labels 

141 

142 

143def test_gaussian() -> None: 

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

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

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

147 assert psf.bounds == bounds 

148 

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

150 assert kernel.bbox == psf.kernel_bbox 

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

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

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

154 

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

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

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

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

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

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

161 

162 with RoundtripFits(psf) as roundtrip: 

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

164 

165 with pytest.raises(ValueError): 

166 # Even stamp size. 

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

168 

169 with pytest.raises(ValueError): 

170 # Negative stamp size. 

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

172 

173 with pytest.raises(ValueError): 

174 # Negative sigma. 

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

176 

177 

178def test_piff_writer_normalizes_tuple_metadata(): # intentionally untyped 

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

180 writer = _ArchivePiffWriter() 

181 writer.write_struct( 

182 "interp", 

183 { 

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

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

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

187 }, 

188 ) 

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

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

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

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

193 with warnings.catch_warnings(): 

194 warnings.simplefilter("error", UserWarning) 

195 model.model_dump_json() 

196 

197 

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

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

200 and converting it back to a legacy PSF. 

201 """ 

202 from piff import PSF 

203 

204 legacy_psf, bounds = legacy_piff_psf_and_bbox 

205 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

206 assert isinstance(psf, PiffWrapper) 

207 assert psf.bounds == bounds 

208 assert isinstance(psf.piff_psf, PSF) 

209 compare_psf_to_legacy(psf, legacy_psf) 

210 with RoundtripFits(psf) as roundtrip1: 

211 pass 

212 compare_psf_to_legacy(roundtrip1.result, legacy_psf) 

213 with RoundtripJson(psf) as roundtrip2: 

214 pass 

215 compare_psf_to_legacy(roundtrip2.result, legacy_psf) 

216 legacy_psf_2 = roundtrip1.result.to_legacy() 

217 compare_psf_to_legacy(psf, legacy_psf_2) 

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

219 

220 

221@skip_no_h5py 

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

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

224 legacy_psf, bounds = legacy_piff_psf_and_bbox 

225 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

226 with RoundtripNdf(psf) as roundtrip: 

227 pass 

228 compare_psf_to_legacy(roundtrip.result, legacy_psf) 

229 

230 

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

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

233 JSON. 

234 """ 

235 from lsst.meas.extensions.psfex import PsfexPsf 

236 

237 legacy_psf, bounds = legacy_psfex_psf_and_bbox 

238 psf = PointSpreadFunction.from_legacy(legacy_psf, bounds) 

239 assert isinstance(psf, PSFExWrapper) 

240 assert psf.bounds == bounds 

241 assert isinstance(psf.legacy_psf, PsfexPsf) 

242 compare_psf_to_legacy(psf, legacy_psf) 

243 with RoundtripFits(psf) as roundtrip1: 

244 pass 

245 compare_psf_to_legacy(roundtrip1.result, legacy_psf) 

246 with RoundtripJson(psf) as roundtrip2: 

247 pass 

248 compare_psf_to_legacy(roundtrip2.result, legacy_psf)