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

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 

21 

22import numpy as np 

23import pydantic 

24 

25from lsst.images._image import Image 

26 

27from .. import serialization 

28from .._geom import Bounds, Box, SerializableBounds 

29from ._base import PointSpreadFunction 

30 

31 

32class GaussianPointSpreadFunction(PointSpreadFunction): 

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

34 

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

36 if sigma <= 0: 

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

38 if stamp_size <= 0: 

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

40 if stamp_size % 2 != 1: 

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

42 self.sigma = float(sigma) 

43 self._stamp_size = stamp_size 

44 self._bounds = bounds 

45 self._sigma2 = self.sigma * self.sigma 

46 

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

48 if not isinstance(other, GaussianPointSpreadFunction): 

49 return NotImplemented 

50 if self.sigma != other.sigma: 

51 return False 

52 if self._stamp_size != other._stamp_size: 

53 return False 

54 if self._bounds != other._bounds: 

55 return False 

56 return True 

57 

58 def __repr__(self) -> str: 

59 return ( 

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

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

62 ) 

63 

64 @property 

65 def bounds(self) -> Bounds: 

66 return self._bounds 

67 

68 @cached_property 

69 def kernel_bbox(self) -> Box: 

70 r = self._stamp_size // 2 

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

72 

73 @cached_property 

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

75 r = self._stamp_size // 2 

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

77 

78 @cached_property 

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

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

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

82 kernel /= kernel.sum() 

83 return kernel 

84 

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

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

87 

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

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

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

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

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

93 image /= image.sum() 

94 return Image(image, bbox=bbox) 

95 

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

97 r = self._stamp_size // 2 

98 xi = round(x) 

99 yi = round(y) 

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

101 

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

103 return GaussianPSFSerializationModel( 

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

105 ) 

106 

107 @classmethod 

108 def deserialize( 

109 cls, model: GaussianPSFSerializationModel, archive: serialization.InputArchive[Any] 

110 ) -> GaussianPointSpreadFunction: 

111 return cls(sigma=model.sigma, bounds=Bounds.deserialize(model.bounds), stamp_size=model.stamp_size) 

112 

113 @staticmethod 

114 def _get_archive_tree_type( 

115 pointer_type: type[pydantic.BaseModel], 

116 ) -> type[GaussianPSFSerializationModel]: 

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

118 type that uses the given pointer type. 

119 """ 

120 return GaussianPSFSerializationModel 

121 

122 

123class GaussianPSFSerializationModel(serialization.ArchiveTree): 

124 sigma: float = pydantic.Field( 

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

126 ) 

127 stamp_size: int = pydantic.Field( 

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

129 ) 

130 bounds: SerializableBounds = pydantic.Field( 

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

132 )