Coverage for python / lsst / analysis / tools / actions / vector / calcMomentSize.py: 34%

1# This file is part of analysis_tools. 

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# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21from __future__ import annotations 

22 

23__all__ = ("CalcMomentSize",) 

24 

25import numpy as np 

26 

27from lsst.pex.config import Field, FieldValidationError 

28from lsst.pex.config.choiceField import ChoiceField 

29 

30from ...interfaces import KeyedData, KeyedDataSchema, Vector, VectorAction 

31from ...math import power, sqrt 

32 

33 

34class CalcMomentSize(VectorAction): 

35 r"""Calculate a size based on 2D moments. 

36 

37 Given a 2x2 matrix of moments (i.e. moment of inertia), two sizes can be 

38 defined as follows: 

39 

40 Determinant radius: :math:`(I_{xx}I_{yy}-I_{xy}^2)^{\frac{1}{4}}` 

41 Trace radius: :math:`\sqrt{(I_{xx}+I_{yy})/2}` 

42 

43 The square of size measure is typically expressed either as the arithmetic 

44 mean of the eigenvalues of the moment matrix (trace radius) or as the 

45 geometric mean of the eigenvalues (determinant radius), which can be 

46 specified using the `sizeType` parameter. Both of these measures 

47 correspond to the :math:`\sigma^2` parameter for a 2D Gaussian. 

48 

49 Notes 

50 ----- 

51 Since lensing preserves surface brightness, the determinant radius relates 

52 the magnification cleanly as it is derived from the area of isophotes, but 

53 have a slightly higher chance of being NaNs for noisy moment estimates. 

54 """ 

55 

56 colXx = Field[str]( 

57 doc="The column name to get the xx shape component from.", 

58 default="{band}_ixx", 

59 ) 

60 

61 colYy = Field[str]( 

62 doc="The column name to get the yy shape component from.", 

63 default="{band}_iyy", 

64 ) 

65 

66 colXy = Field[str]( 

67 doc="The column name to get the xy shape component from.", 

68 default="{band}_ixy", 

69 optional=True, 

70 ) 

71 

72 is_covariance = Field[bool]( 

73 doc="Whether the fields are for a covariance matrix. If False, the XX/YY/XY terms are instead" 

74 " assumed to map to sigma_x/sigma_y/rho.", 

75 default=True, 

76 ) 

77 

78 sizeType = ChoiceField[str]( 

79 doc="The type of size to calculate", 

80 default="determinant", 

81 optional=False, 

82 allowed={ 

83 "trace": r"Trace radius :math:`\sqrt{(I_{xx}+I_{yy})/2}`", 

84 "determinant": r"Determinant radius :math:`(I_{xx}I_{yy}-I_{xy}^2)^{\frac{1}{4}}`", 

85 }, 

86 ) 

87 

88 def getInputSchema(self) -> KeyedDataSchema: 

89 if self.sizeType == "trace": 

90 return ( 

91 (self.colXx, Vector), 

92 (self.colYy, Vector), 

93 ) 

94 else: 

95 return ( 

96 (self.colXx, Vector), 

97 (self.colYy, Vector), 

98 (self.colXy, Vector), 

99 ) # type: ignore 

100 

101 def __call__(self, data: KeyedData, **kwargs) -> Vector: 

102 xx = np.array(data[self.colXx.format(**kwargs)]) 

103 yy = np.array(data[self.colYy.format(**kwargs)]) 

104 if self.sizeType == "trace": 

105 if not self.is_covariance: 

106 xx *= xx 

107 yy *= yy 

108 size = sqrt(0.5 * (xx + yy)) 

109 else: 

110 xy_sq = np.array(data[self.colXy.format(**kwargs)]) ** 2 

111 if not self.is_covariance: 

112 xx *= xx 

113 yy *= yy 

114 # cov = rho * sigma_x * sigma_y 

115 # this term needs to be cov^2 

116 xy_sq *= xx * yy 

117 size = power(xx * yy - xy_sq, 0.25) 

118 

119 return size 

120 

121 def validate(self): 

122 super().validate() 

123 if self.sizeType == "determinant" and self.colXy is None: 

124 msg = "colXy is required for determinant-type size" 

125 raise FieldValidationError(self.__class__.colXy, self, msg)