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

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 SerializableBounds 

29from .._geom import Bounds, Box 

30from ..utils import round_half_up 

31from ._base import PointSpreadFunction 

32 

33 

34class GaussianPointSpreadFunction(PointSpreadFunction): 

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

36 

37 Parameters 

38 ---------- 

39 sigma 

40 Standard deviation of the Gaussian profile in pixels. 

41 bounds 

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

43 evaluated. 

44 stamp_size 

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

46 odd number. 

47 """ 

48 

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

50 if sigma <= 0: 

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

52 if stamp_size <= 0: 

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

54 if stamp_size % 2 != 1: 

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

56 self.sigma = float(sigma) 

57 self._stamp_size = stamp_size 

58 self._bounds = bounds 

59 self._sigma2 = self.sigma * self.sigma 

60 

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

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

63 return NotImplemented 

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

65 return False 

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

67 return False 

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

69 return False 

70 return True 

71 

72 def __repr__(self) -> str: 

73 return ( 

74 f"GaussianPointSpreadFunction({self.sigma}, " 

75 f"stamp_size={self._stamp_size}, bounds={self._bounds!r})" 

76 ) 

77 

78 @property 

79 def bounds(self) -> Bounds: 

80 return self._bounds 

81 

82 @cached_property 

83 def kernel_bbox(self) -> Box: 

84 r = self._stamp_size // 2 

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

86 

87 @cached_property 

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

89 r = self._stamp_size // 2 

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

91 

92 @cached_property 

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

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

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

96 kernel /= kernel.sum() 

97 return kernel 

98 

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

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

101 

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

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

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

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

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

107 image /= image.sum() 

108 return Image(image, bbox=bbox) 

109 

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

111 r = self._stamp_size // 2 

112 xi = round_half_up(x) 

113 yi = round_half_up(y) 

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

115 

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

117 return GaussianPSFSerializationModel( 

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

119 ) 

120 

121 @staticmethod 

122 def _get_archive_tree_type( 

123 pointer_type: type[pydantic.BaseModel], 

124 ) -> type[GaussianPSFSerializationModel]: 

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

126 type that uses the given pointer type. 

127 """ 

128 return GaussianPSFSerializationModel 

129 

130 

131class GaussianPSFSerializationModel(serialization.ArchiveTree): 

132 SCHEMA_NAME: ClassVar[str] = "gaussian_psf" 

133 SCHEMA_VERSION: ClassVar[str] = "1.0.0" 

134 MIN_READ_VERSION: ClassVar[int] = 1 

135 PUBLIC_TYPE: ClassVar[type] = GaussianPointSpreadFunction 

136 

137 sigma: float = pydantic.Field( 

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

139 ) 

140 stamp_size: int = pydantic.Field( 

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

142 ) 

143 bounds: SerializableBounds = pydantic.Field( 

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

145 ) 

146 

147 def deserialize( 

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

149 ) -> GaussianPointSpreadFunction: 

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

151 raise serialization.InvalidParameterError( 

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

153 ) 

154 return GaussianPointSpreadFunction( 

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

156 )