Coverage for python/lsst/images/psfs/_gaussian.py: 52%

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 

14__all__ = ( 

15 "GaussianPSFSerializationModel", 

16 "GaussianPointSpreadFunction", 

17) 

18 

19from functools import cached_property 

20from typing import Any, ClassVar 

21 

22import numpy as np 

23import pydantic 

24 

25from lsst.images._image import Image 

26 

27from .. import serialization 

28from .._concrete_bounds import BoundsSerializationModel 

29from .._geom import Bounds, Box 

30from ..describe import DescribeOptions, Report, ReportField 

31from ..utils import round_half_up 

32from ._base import PointSpreadFunction 

33 

34 

35class GaussianPointSpreadFunction(PointSpreadFunction): 

36 """A PSF with a spatially-invariant circular Gaussian profile. 

37 

38 Parameters 

39 ---------- 

40 sigma 

41 Standard deviation of the Gaussian profile in pixels. 

42 bounds 

43 The pixel-coordinate region where the model can safely be 

44 evaluated. 

45 stamp_size 

46 Side length in pixels of the PSF image stamps; must be a positive 

47 odd number. 

48 """ 

49 

50 def __init__(self, sigma: float, bounds: Bounds, stamp_size: int) -> None: 

51 if sigma <= 0: 

52 raise ValueError(f"sigma must be positive; got {sigma}.") 

53 if stamp_size <= 0: 

54 raise ValueError(f"stamp_size must be positive; got {stamp_size}.") 

55 if stamp_size % 2 != 1: 

56 raise ValueError(f"stamp_size must be odd number, got {stamp_size}") 

57 self.sigma = float(sigma) 

58 self._stamp_size = stamp_size 

59 self._bounds = bounds 

60 self._sigma2 = self.sigma * self.sigma 

61 

62 def __eq__(self, other: Any) -> bool: 

63 if not isinstance(other, GaussianPointSpreadFunction): 63 ↛ 64line 63 didn't jump to line 64 because the condition on line 63 was never true

64 return NotImplemented 

65 if self.sigma != other.sigma: 65 ↛ 66line 65 didn't jump to line 66 because the condition on line 65 was never true

66 return False 

67 if self._stamp_size != other._stamp_size: 67 ↛ 68line 67 didn't jump to line 68 because the condition on line 67 was never true

68 return False 

69 if self._bounds != other._bounds: 69 ↛ 70line 69 didn't jump to line 70 because the condition on line 69 was never true

70 return False 

71 return True 

72 

73 def _describe(self, options: DescribeOptions = DescribeOptions(), /) -> Report: 

74 """Return a `Report` describing this Gaussian PSF. 

75 

76 Parameters 

77 ---------- 

78 options : `DescribeOptions`, optional 

79 Unused; accepted for interface compatibility. 

80 """ 

81 return Report( 

82 type_name="GaussianPointSpreadFunction", 

83 summary=f"Gaussian PSF sigma={self.sigma}", 

84 fields=[ 

85 ReportField(label="sigma", value=self.sigma, positional=True), 

86 ReportField(label="stamp_size", value=self._stamp_size), 

87 ReportField(label="bounds", value=self._bounds, repr_value=repr(self._bounds)), 

88 ], 

89 ) 

90 

91 @property 

92 def bounds(self) -> Bounds: 

93 return self._bounds 

94 

95 @cached_property 

96 def kernel_bbox(self) -> Box: 

97 r = self._stamp_size // 2 

98 return Box.factory[-r : r + 1, -r : r + 1] 

99 

100 @cached_property 

101 def _centered_coordinates(self) -> np.ndarray: 

102 r = self._stamp_size // 2 

103 return np.arange(-r, r + 1, dtype=np.float64) 

104 

105 @cached_property 

106 def _kernel_array(self) -> np.ndarray: 

107 profile = np.exp(-0.5 * np.square(self._centered_coordinates) / self._sigma2) 

108 kernel = np.multiply.outer(profile, profile) 

109 kernel /= kernel.sum() 

110 return kernel 

111 

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

113 return Image(self._kernel_array.copy(), bbox=self.kernel_bbox) 

114 

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

116 bbox = self.compute_stellar_bbox(x=x, y=y) 

117 x_profile = np.exp(-0.5 * np.square(bbox.x.arange - x) / self._sigma2) 

118 y_profile = np.exp(-0.5 * np.square(bbox.y.arange - y) / self._sigma2) 

119 image = np.multiply.outer(y_profile, x_profile) 

120 image /= image.sum() 

121 return Image(image, bbox=bbox) 

122 

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

124 r = self._stamp_size // 2 

125 xi = round_half_up(x) 

126 yi = round_half_up(y) 

127 return Box.factory[yi - r : yi + r + 1, xi - r : xi + r + 1] 

128 

129 @classmethod 

130 def from_legacy(cls, legacy_psf: Any, bounds: Bounds) -> GaussianPointSpreadFunction: 

131 """Convert from a legacy `lsst.afw.detection.GaussianPsf`. 

132 

133 Parameters 

134 ---------- 

135 legacy_psf 

136 Legacy Gaussian PSF, with a square stamp. 

137 bounds 

138 The region where this PSF model is valid. 

139 

140 Returns 

141 ------- 

142 `GaussianPointSpreadFunction` 

143 The converted PSF. 

144 

145 Raises 

146 ------ 

147 ValueError 

148 Raised if the stamp of ``legacy_psf`` is not square. 

149 """ 

150 dimensions = legacy_psf.getDimensions() 

151 if dimensions.x != dimensions.y: 

152 raise ValueError(f"Legacy GaussianPsf stamp must be square; got {dimensions}.") 

153 return cls(legacy_psf.getSigma(), bounds=bounds, stamp_size=dimensions.x) 

154 

155 def to_legacy(self) -> Any: 

156 """Convert to a legacy `lsst.afw.detection.GaussianPsf`.""" 

157 from lsst.afw.detection import GaussianPsf 

158 

159 return GaussianPsf(self._stamp_size, self._stamp_size, self.sigma) 

160 

161 def serialize(self, archive: serialization.OutputArchive[Any]) -> GaussianPSFSerializationModel: 

162 return GaussianPSFSerializationModel( 

163 sigma=self.sigma, stamp_size=self._stamp_size, bounds=self._bounds.serialize() 

164 ) 

165 

166 @staticmethod 

167 def _get_archive_tree_type( 

168 pointer_type: type[pydantic.BaseModel], 

169 ) -> type[GaussianPSFSerializationModel]: 

170 """Return the serialization model type for this object for an archive 

171 type that uses the given pointer type. 

172 """ 

173 return GaussianPSFSerializationModel 

174 

175 

176class GaussianPSFSerializationModel(serialization.ArchiveTree): 

177 SCHEMA_NAME: ClassVar[str] = "gaussian_psf" 

178 SCHEMA_VERSION: ClassVar[str] = "1.0.0" 

179 MIN_READ_VERSION: ClassVar[int] = 1 

180 PUBLIC_TYPE: ClassVar[type] = GaussianPointSpreadFunction 

181 

182 sigma: float = pydantic.Field( 

183 description="Gaussian sigma for the PSF.", 

184 ) 

185 stamp_size: int = pydantic.Field( 

186 description="Width of the (square) images returned by this PSF's methods." 

187 ) 

188 bounds: BoundsSerializationModel = pydantic.Field( 

189 description="The bounds object that represents the PSF's validity region." 

190 ) 

191 

192 def deserialize( 

193 self, archive: serialization.InputArchive[Any], **kwargs: Any 

194 ) -> GaussianPointSpreadFunction: 

195 if kwargs: 195 ↛ 196line 195 didn't jump to line 196 because the condition on line 195 was never true

196 raise serialization.InvalidParameterError( 

197 f"Unrecognized parameters for GaussianPointSpreadFunction: {set(kwargs.keys())}." 

198 ) 

199 return GaussianPointSpreadFunction( 

200 sigma=self.sigma, bounds=self.bounds.deserialize(), stamp_size=self.stamp_size 

201 )