Coverage for python / lsst / analysis / tools / actions / vector / ellipticity.py: 34%
90 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-30 09:27 +0000
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
23__all__ = (
24 "CalcE",
25 "CalcEDiff",
26 "CalcE1",
27 "CalcE2",
28)
30from typing import cast
32import numpy as np
34from lsst.pex.config import ChoiceField, Field, FieldValidationError
35from lsst.pex.config.configurableActions import ConfigurableActionField
37from ...interfaces import KeyedData, KeyedDataSchema, Vector, VectorAction
40class CalcE(VectorAction):
41 r"""Calculate a complex value representation of the ellipticity.
43 The complex ellipticity is typically defined as:
45 .. math::
46 e &= |e|\exp{(\mathrm{i}2\theta)} = e_1+\mathrm{i}e_2, \\
47 &= \frac{(I_{xx} - I_{yy}) + \mathrm{i}2I_{xy}}{I_{xx} + I_{yy}},
49 where :math:`\mathrm{i}` is the square root of -1 and :math:`I_{xx}`,
50 :math:`I_{yy}`, and :math:`I_{xy}` are second-order central moments.
51 This is sometimes referred to as distortion, and denoted in GalSim by
52 :math:`e=(e_1,e_2)` (see Eq. 4.4. of Bartelmann and Schneider, 2001 [1]_).
53 The other definition differs in normalization.
54 It is referred to as shear, and denoted by :math:`g=(g_{1},g_{2})`
55 in GalSim (see Eq. 4.10 of Bartelmann and Schneider, 2001 [1]_).
56 It is defined as
58 .. math::
60 g = \frac{(I_{xx} - I_{yy}) + \mathrm{i}2I_{xy}}
61 {I_{xx} + I_{yy} + 2\sqrt{(I_{xx}I_{yy}-I_{xy}^{2})}}.
63 The shear measure is unbiased in weak-lensing shear, but may exclude some
64 objects in the presence of noisy moment estimates. The distortion measure
65 is biased in weak-lensing distortion, but does not suffer from selection
66 artifacts.
68 References
69 ----------
70 .. [1] Bartelmann, M. and Schneider, P., “Weak gravitational lensing”,
71 Physics Reports, vol. 340, no. 4–5, pp. 291–472, 2001.
72 doi:10.1016/S0370-1573(00)00082-X;
73 https://arxiv.org/abs/astro-ph/9912508
75 Notes
76 -----
78 1. This is a shape measurement used for doing QA on the ellipticity
79 of the sources.
81 2. For plotting purposes we might want to plot quivers whose lengths
82 are proportional to :math:`|e|` and whose angles correspond to
83 :math:`\theta`.
84 If `halvePhaseAngle` config parameter is set to `True`, then the returned
85 quantity therefore corresponds to the complex quantity
86 :math:`|e|\exp{(\mathrm{i}\theta)}` or its real and imaginary parts
87 (depending on the `component`).
89 See Also
90 --------
91 CalcE1
92 CalcE2
93 """
95 colXx = Field[str](
96 doc="The column name to get the xx shape component from.",
97 default="{band}_ixx",
98 )
100 colYy = Field[str](
101 doc="The column name to get the yy shape component from.",
102 default="{band}_iyy",
103 )
105 colXy = Field[str](
106 doc="The column name to get the xy shape component from.",
107 default="{band}_ixy",
108 )
110 ellipticityType = ChoiceField[str](
111 doc="The type of ellipticity to calculate",
112 allowed={
113 "distortion": (
114 "Distortion, defined as " r":math:`(I_{xx}-I_{yy}+\mathrm{i}2I_{xy})/(I_{xx}+I_{yy})`"
115 ),
116 "shear": (
117 "Shear, defined as "
118 r":math:`(I_{xx}-I_{yy}+\mathrm{i}2I_{xy})/(I_{xx}+I_{yy}+2\sqrt{I_{xx}I_{yy}-I_{xy}^2})`"
119 ),
120 },
121 default="distortion",
122 optional=False,
123 )
125 halvePhaseAngle = Field[bool](
126 doc="Divide the phase angle by 2? Suitable for quiver plots.",
127 default=False,
128 )
130 component = ChoiceField[str](
131 doc="Which component of the ellipticity to return. If `None`, return complex ellipticity values.",
132 allowed={
133 "1": r":math:`e_1` or :math:`g_1` (depending on `ellipticityType`)",
134 "2": r":math:`e_2` or :math:`g_2` (depending on `ellipticityType`)",
135 None: r":math:`e_1 + \mathrm{i}e_2` or :math:`g_1 + \mathrm{i}g_2`"
136 " (depending on `ellipticityType`)",
137 },
138 )
140 def getInputSchema(self) -> KeyedDataSchema:
141 return ((self.colXx, Vector), (self.colXy, Vector), (self.colYy, Vector))
143 def __call__(self, data: KeyedData, **kwargs) -> Vector:
144 e = (data[self.colXx.format(**kwargs)] - data[self.colYy.format(**kwargs)]) + 1j * (
145 2 * data[self.colXy.format(**kwargs)]
146 )
147 denom = data[self.colXx.format(**kwargs)] + data[self.colYy.format(**kwargs)]
149 if self.ellipticityType == "shear":
150 denom += 2 * np.sqrt(
151 data[self.colXx.format(**kwargs)] * data[self.colYy.format(**kwargs)]
152 - data[self.colXy.format(**kwargs)] ** 2
153 )
154 e = cast(Vector, e)
155 denom = cast(Vector, denom)
157 e /= denom
159 if self.halvePhaseAngle:
160 # Ellipiticity is |e|*exp(i*2*theta), but we want to return
161 # |e|*exp(i*theta). So we multiply by |e| and take its square root
162 # instead of the more expensive trig calls.
163 e *= np.abs(e)
164 e = np.sqrt(e)
166 if self.component == "1":
167 return np.real(e)
168 elif self.component == "2":
169 return np.imag(e)
170 else:
171 return e
174class CalcEDiff(VectorAction):
175 r"""Calculate the difference of two ellipticities as a complex quantity.
177 The complex ellipticity (for both distortion-type and shear-type)
178 difference ( between :math:`e_A` and :math:`e_B` is defined as
179 :math:`e_{A}-e_{B}=\delta e=|\delta e|\exp{(\mathrm{i}2\theta_{\delta})}`
182 See Also
183 --------
184 CalcE
186 Notes
187 -----
189 1. This is a shape measurement used for doing QA on the ellipticity
190 of the sources.
192 2. The `ellipticityType` of `colA` and `colB` have to be the same.
194 3. For plotting purposes we might want to plot quivers whose lengths
195 are proportional to :math:`|\delta e|` and whose angles correspond to
196 :math:`\theta_\delta`.
197 If `halvePhaseAngle` config parameter is set to `True`, then the returned
198 quantity therefore corresponds to the complex quantity
199 :math:`|\delta e|\exp{(\mathrm{i}\theta_\delta)}`.
200 """
202 colA = ConfigurableActionField[VectorAction](
203 doc="Ellipticity to subtract from",
204 default=CalcE,
205 )
207 colB = ConfigurableActionField[VectorAction](
208 doc="Ellipticity to subtract",
209 dtype=VectorAction,
210 default=CalcE,
211 )
213 halvePhaseAngle = Field[bool](
214 doc="Divide the phase angle by 2? Suitable for quiver plots.",
215 default=False,
216 )
218 component = ChoiceField[str](
219 doc="Which component of the ellipticity to return. If `None`, return complex ellipticity values.",
220 allowed={
221 "1": r":math:`\delta e_1` or :math:`\delta g_1` (depending on the common `ellipiticityType`)",
222 "2": r":math:`\delta e_2` or :math:`\delta g_2` (depending on the common `ellipiticityType`)",
223 None: r":math:`\delta e_1+\mathrm{i}\delta e_2` or :math:`\delta g_1 \mathrm{i}\delta g_2`"
224 " (depending on the common `ellipticityType`)",
225 },
226 )
228 def getInputSchema(self) -> KeyedDataSchema:
229 yield from self.colA.getInputSchema()
230 yield from self.colB.getInputSchema()
232 def validate(self):
233 super().validate()
234 if self.colA.ellipticityType != self.colB.ellipticityType:
235 msg = "Both the ellipticities in CalcEDiff must have the same type."
236 raise FieldValidationError(self.colB.__class__.ellipticityType, self, msg)
238 def __call__(self, data: KeyedData, **kwargs) -> Vector:
239 eMeas = self.colA(data, **kwargs)
240 ePSF = self.colB(data, **kwargs)
241 eDiff = eMeas - ePSF
242 if self.halvePhaseAngle:
243 # Ellipiticity is |e|*exp(i*2*theta), but we want to return
244 # |e|*exp(j*theta). So we multiply by |e| and take its square root
245 # instead of the more expensive trig calls.
246 eDiff *= np.abs(eDiff)
247 eDiff = np.sqrt(eDiff)
249 if self.component == "1":
250 return np.real(eDiff)
251 elif self.component == "2":
252 return np.imag(eDiff)
253 else:
254 return eDiff
257class CalcE1(VectorAction):
258 r"""Calculate :math:`e_1` (distortion-type) or :math:`g_1` (shear-type).
260 The definitions are as follows:
262 .. math::
263 e_1&=(I_{xx}-I_{yy})/(I_{xx}+I_{yy}) \\
264 g_1&=(I_{xx}-I_{yy})/(I_{xx}+I_{yy}+2\sqrt{I_{xx}I_{yy}-I_{xy}^{2}}).
266 See Also
267 --------
268 CalcE
269 CalcE2
271 Notes
272 -----
273 This is a shape measurement used for doing QA on the ellipticity
274 of the sources.
275 """
277 colXx = Field[str](
278 doc="The column name to get the xx shape component from.",
279 default="{band}_ixx",
280 )
282 colYy = Field[str](
283 doc="The column name to get the yy shape component from.",
284 default="{band}_iyy",
285 )
287 colXy = Field[str](
288 doc="The column name to get the xy shape component from.",
289 default="{band}_ixy",
290 optional=True,
291 )
293 ellipticityType = ChoiceField[str](
294 doc="The type of ellipticity to calculate",
295 optional=False,
296 allowed={
297 "distortion": ("Distortion, measured as " r":math:`(I_{xx}-I_{yy})/(I_{xx}+I_{yy})`"),
298 "shear": (
299 "Shear, measured as " r":math:`(I_{xx}-I_{yy})/(I_{xx}+I_{yy}+2\sqrt{I_{xx}I_{yy}-I_{xy}^2})`"
300 ),
301 },
302 default="distortion",
303 )
305 def getInputSchema(self) -> KeyedDataSchema:
306 if self.ellipticityType == "distortion":
307 return (
308 (self.colXx, Vector),
309 (self.colYy, Vector),
310 )
311 else:
312 return (
313 (self.colXx, Vector),
314 (self.colYy, Vector),
315 (self.colXy, Vector),
316 )
318 def __call__(self, data: KeyedData, **kwargs) -> Vector:
319 denom = data[self.colXx.format(**kwargs)] + data[self.colYy.format(**kwargs)]
320 if self.ellipticityType == "shear":
321 denom += 2 * np.sqrt(
322 data[self.colXx.format(**kwargs)] * data[self.colYy.format(**kwargs)]
323 - data[self.colXy.format(**kwargs)] ** 2
324 )
325 e1 = (data[self.colXx.format(**kwargs)] - data[self.colYy.format(**kwargs)]) / denom
327 return cast(Vector, e1)
329 def validate(self):
330 super().validate()
331 if self.ellipticityType == "shear" and self.colXy is None:
332 msg = "colXy is required for shear-type shear ellipticity"
333 raise FieldValidationError(self.__class__.colXy, self, msg)
336class CalcE2(VectorAction):
337 r"""Calculate :math:`e_2` (distortion-type) or :math:`g_2` (shear-type).
339 The definitions are as follows:
341 .. math::
342 e_2 &= 2I_{xy}/(I_{xx}+I_{yy}) \\
343 g_2 &= 2I_{xy}/(I_{xx}+I_{yy}+2\sqrt{(I_{xx}I_{yy}-I_{xy}^{2})}).
345 See Also
346 --------
347 CalcE
348 CalcE1
350 Notes
351 -----
352 This is a shape measurement used for doing QA on the ellipticity
353 of the sources.
354 """
356 colXx = Field[str](
357 doc="The column name to get the xx shape component from.",
358 default="{band}_ixx",
359 )
361 colYy = Field[str](
362 doc="The column name to get the yy shape component from.",
363 default="{band}_iyy",
364 )
366 colXy = Field[str](
367 doc="The column name to get the xy shape component from.",
368 default="{band}_ixy",
369 )
371 ellipticityType = ChoiceField[str](
372 doc="The type of ellipticity to calculate",
373 optional=False,
374 allowed={
375 "distortion": ("Distortion, defined as " r":math:`2I_{xy}/(I_{xx}+I_{yy})`"),
376 "shear": r"Shear, defined as :math:`2I_{xy}/(I_{xx}+I_{yy}+2\sqrt{I_{xx}I_{yy}-I_{xy}^2})`",
377 },
378 default="distortion",
379 )
381 def getInputSchema(self) -> KeyedDataSchema:
382 return (
383 (self.colXx, Vector),
384 (self.colYy, Vector),
385 (self.colXy, Vector),
386 )
388 def __call__(self, data: KeyedData, **kwargs) -> Vector:
389 denom = data[self.colXx.format(**kwargs)] + data[self.colYy.format(**kwargs)]
390 if self.ellipticityType == "shear":
391 denom += 2 * np.sqrt(
392 data[self.colXx.format(**kwargs)] * data[self.colYy.format(**kwargs)]
393 - data[self.colXy.format(**kwargs)] ** 2
394 )
395 e2 = 2 * data[self.colXy.format(**kwargs)] / denom
396 return cast(Vector, e2)