Loading [MathJax]/extensions/tex2jax.js
LSST Applications g0fba68d861+83433b07ee,g16d25e1f1b+23bc9e47ac,g1ec0fe41b4+3ea9d11450,g1fd858c14a+9be2b0f3b9,g2440f9efcc+8c5ae1fdc5,g35bb328faa+8c5ae1fdc5,g4a4af6cd76+d25431c27e,g4d2262a081+c74e83464e,g53246c7159+8c5ae1fdc5,g55585698de+1e04e59700,g56a49b3a55+92a7603e7a,g60b5630c4e+1e04e59700,g67b6fd64d1+3fc8cb0b9e,g78460c75b0+7e33a9eb6d,g786e29fd12+668abc6043,g8352419a5c+8c5ae1fdc5,g8852436030+60e38ee5ff,g89139ef638+3fc8cb0b9e,g94187f82dc+1e04e59700,g989de1cb63+3fc8cb0b9e,g9d31334357+1e04e59700,g9f33ca652e+0a83e03614,gabe3b4be73+8856018cbb,gabf8522325+977d9fabaf,gb1101e3267+8b4b9c8ed7,gb89ab40317+3fc8cb0b9e,gc0af124501+57ccba3ad1,gcf25f946ba+60e38ee5ff,gd6cbbdb0b4+1cc2750d2e,gd794735e4e+7be992507c,gdb1c4ca869+be65c9c1d7,gde0f65d7ad+c7f52e58fe,ge278dab8ac+6b863515ed,ge410e46f29+3fc8cb0b9e,gf35d7ec915+97dd712d81,gf5e32f922b+8c5ae1fdc5,gf618743f1b+747388abfa,gf67bdafdda+3fc8cb0b9e,w.2025.18
LSST Data Management Base Package
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
radii.cc
Go to the documentation of this file.
1// -*- lsst-c++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008, 2009, 2010 LSST Corporation.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24#include "lsst/afw/geom/ellipses/radii.h"
25#include <complex>
26
27namespace lsst {
28namespace afw {
29namespace geom {
30namespace ellipses {
31
32void DeterminantRadius::assignFromQuadrupole(double ixx, double iyy, double ixy, Distortion& distortion) {
33 _value = std::pow(ixx * iyy - ixy * ixy, 0.25);
34 distortion.setE1((ixx - iyy) / (ixx + iyy));
35 distortion.setE2(2.0 * ixy / (ixx + iyy));
36}
37
38BaseCore::Jacobian DeterminantRadius::dAssignFromQuadrupole(double ixx, double iyy, double ixy,
39 Distortion& distortion) {
40 double xx_yy = ixx + iyy;
41 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
42 _value = std::pow(ixx * iyy - ixy * ixy, 0.25);
43 distortion.setE1((ixx - iyy) / xx_yy);
44 distortion.setE2(2.0 * ixy / xx_yy);
45 result.block<2, 2>(0, 0).setConstant(2.0 / (xx_yy * xx_yy));
46 result(0, 0) *= iyy;
47 result(0, 1) *= -ixx;
48 result(1, 0) *= -ixy;
49 result(1, 1) *= -ixy;
50 result(1, 2) = 2.0 / xx_yy;
51 result.row(2).setConstant(1.0 / (4.0 * _value * _value * _value));
52 result(2, 0) *= iyy;
53 result(2, 1) *= ixx;
54 result(2, 2) *= -2.0 * ixy;
55 return result;
56}
57
58void DeterminantRadius::assignToQuadrupole(Distortion const& distortion, double& ixx, double& iyy,
59 double& ixy) const {
60 double den = std::sqrt(1.0 - std::norm(distortion.getComplex()));
61 double r2 = _value * _value;
62 ixx = r2 * (1.0 + distortion.getE1()) / den;
63 iyy = r2 * (1.0 - distortion.getE1()) / den;
64 ixy = r2 * distortion.getE2() / den;
65}
66
67BaseCore::Jacobian DeterminantRadius::dAssignToQuadrupole(Distortion const& distortion, double& ixx,
68 double& iyy, double& ixy) const {
69 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
70 double den = std::sqrt(1.0 - std::norm(distortion.getComplex()));
71 result.col(2).setConstant(2.0 * _value / den);
72 double r2 = _value * _value;
73 ixx = r2 * (1.0 + distortion.getE1()) / den;
74 iyy = r2 * (1.0 - distortion.getE1()) / den;
75 ixy = r2 * distortion.getE2() / den;
76 result.block<3, 2>(0, 0).setConstant(r2 / (den * den * den));
77 result(0, 0) *= (distortion.getE1() + 1.0 - distortion.getE2() * distortion.getE2());
78 result(1, 0) *= (distortion.getE1() - 1.0 + distortion.getE2() * distortion.getE2());
79 result(2, 0) *= distortion.getE1() * distortion.getE2();
80 result(0, 1) *= distortion.getE2() * (1.0 + distortion.getE1());
81 result(1, 1) *= distortion.getE2() * (1.0 - distortion.getE1());
82 result(2, 1) *= (1.0 - distortion.getE1() * distortion.getE1());
83 result(0, 2) *= (1.0 + distortion.getE1());
84 result(1, 2) *= (1.0 - distortion.getE1());
85 result(2, 2) *= distortion.getE2();
86 return result;
87}
88
89void TraceRadius::assignFromQuadrupole(double ixx, double iyy, double ixy, Distortion& distortion) {
90 _value = std::sqrt(0.5 * (ixx + iyy));
91 distortion.setE1((ixx - iyy) / (ixx + iyy));
92 distortion.setE2(2.0 * ixy / (ixx + iyy));
93}
94
95BaseCore::Jacobian TraceRadius::dAssignFromQuadrupole(double ixx, double iyy, double ixy,
96 Distortion& distortion) {
97 double xx_yy = ixx + iyy;
98 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
99 _value = std::sqrt(0.5 * xx_yy);
100 distortion.setE1((ixx - iyy) / xx_yy);
101 distortion.setE2(2.0 * ixy / xx_yy);
102 result.block<2, 2>(0, 0).setConstant(2.0 / (xx_yy * xx_yy));
103 result(0, 0) *= iyy;
104 result(0, 1) *= -ixx;
105 result(1, 0) *= -ixy;
106 result(1, 1) *= -ixy;
107 result(1, 2) = 2.0 / xx_yy;
108 result(2, 0) = 0.25 / _value;
109 result(2, 1) = 0.25 / _value;
110 return result;
111}
112
113void TraceRadius::assignToQuadrupole(Distortion const& distortion, double& ixx, double& iyy,
114 double& ixy) const {
115 double r2 = _value * _value;
116 ixx = r2 * (1.0 + distortion.getE1());
117 iyy = r2 * (1.0 - distortion.getE1());
118 ixy = r2 * distortion.getE2();
119}
120
121BaseCore::Jacobian TraceRadius::dAssignToQuadrupole(Distortion const& distortion, double& ixx, double& iyy,
122 double& ixy) const {
123 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
124 result.col(2).setConstant(2.0 * _value);
125 result(0, 2) *= (1.0 + distortion.getE1());
126 result(1, 2) *= (1.0 - distortion.getE1());
127 result(2, 2) *= distortion.getE2();
128 double r2 = _value * _value;
129 ixx = r2 * (1.0 + distortion.getE1());
130 iyy = r2 * (1.0 - distortion.getE1());
131 ixy = r2 * distortion.getE2();
132 result(0, 0) = r2;
133 result(1, 0) = -r2;
134 result(2, 1) = r2;
135 return result;
136}
137
138void LogDeterminantRadius::assignFromQuadrupole(double ixx, double iyy, double ixy, Distortion& distortion) {
139 _value = 0.25 * std::log(ixx * iyy - ixy * ixy);
140 distortion.setE1((ixx - iyy) / (ixx + iyy));
141 distortion.setE2(2.0 * ixy / (ixx + iyy));
142}
143
144BaseCore::Jacobian LogDeterminantRadius::dAssignFromQuadrupole(double ixx, double iyy, double ixy,
145 Distortion& distortion) {
146 double xx_yy = ixx + iyy;
147 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
148 double det = ixx * iyy - ixy * ixy;
149 _value = 0.25 * std::log(det);
150 distortion.setE1((ixx - iyy) / xx_yy);
151 distortion.setE2(2.0 * ixy / xx_yy);
152 result.block<2, 2>(0, 0).setConstant(2.0 / (xx_yy * xx_yy));
153 result(0, 0) *= iyy;
154 result(0, 1) *= -ixx;
155 result(1, 0) *= -ixy;
156 result(1, 1) *= -ixy;
157 result(1, 2) = 2.0 / xx_yy;
158 result(2, 0) = 0.25 * iyy / det;
159 result(2, 1) = 0.25 * ixx / det;
160 result(2, 2) = -0.5 * ixy / det;
161 return result;
162}
163
164void LogDeterminantRadius::assignToQuadrupole(Distortion const& distortion, double& ixx, double& iyy,
165 double& ixy) const {
166 double den = std::sqrt(1.0 - std::norm(distortion.getComplex()));
167 double r2 = std::exp(2.0 * _value);
168 ixx = r2 * (1.0 + distortion.getE1()) / den;
169 iyy = r2 * (1.0 - distortion.getE1()) / den;
170 ixy = r2 * distortion.getE2() / den;
171}
172
173BaseCore::Jacobian LogDeterminantRadius::dAssignToQuadrupole(Distortion const& distortion, double& ixx,
174 double& iyy, double& ixy) const {
175 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
176 double den = std::sqrt(1.0 - std::norm(distortion.getComplex()));
177 result.col(2).setConstant(2.0 * _value / den);
178 double r2 = std::exp(2.0 * _value);
179 ixx = r2 * (1.0 + distortion.getE1()) / den;
180 iyy = r2 * (1.0 - distortion.getE1()) / den;
181 ixy = r2 * distortion.getE2() / den;
182 result.block<3, 2>(0, 0).setConstant(r2 / (den * den * den));
183 result(0, 0) *= distortion.getE1() + 1.0 - distortion.getE2() * distortion.getE2();
184 result(1, 0) *= distortion.getE1() - 1.0 + distortion.getE2() * distortion.getE2();
185 result(2, 0) *= distortion.getE1() * distortion.getE2();
186 result(0, 1) *= distortion.getE2() * (1.0 + distortion.getE1());
187 result(1, 1) *= distortion.getE2() * (1.0 - distortion.getE1());
188 result(2, 1) *= 1.0 - distortion.getE1() * distortion.getE1();
189 result(0, 2) = 2.0 * ixx;
190 result(1, 2) = 2.0 * iyy;
191 result(2, 2) = 2.0 * ixy;
192 return result;
193}
194
195void LogTraceRadius::assignFromQuadrupole(double ixx, double iyy, double ixy, Distortion& distortion) {
196 _value = 0.5 * std::log(0.5 * (ixx + iyy));
197 distortion.setE1((ixx - iyy) / (ixx + iyy));
198 distortion.setE2(2.0 * ixy / (ixx + iyy));
199}
200
201BaseCore::Jacobian LogTraceRadius::dAssignFromQuadrupole(double ixx, double iyy, double ixy,
202 Distortion& distortion) {
203 double xx_yy = ixx + iyy;
204 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
205 _value = 0.5 * std::log(0.5 * xx_yy);
206 distortion.setE1((ixx - iyy) / xx_yy);
207 distortion.setE2(2.0 * ixy / xx_yy);
208 result.block<2, 2>(0, 0).setConstant(2.0 / (xx_yy * xx_yy));
209 result(0, 0) *= iyy;
210 result(0, 1) *= -ixx;
211 result(1, 0) *= -ixy;
212 result(1, 1) *= -ixy;
213 result(1, 2) = 2.0 / xx_yy;
214 result(2, 0) = 0.5 / xx_yy;
215 result(2, 1) = 0.5 / xx_yy;
216 return result;
217}
218
219void LogTraceRadius::assignToQuadrupole(Distortion const& distortion, double& ixx, double& iyy,
220 double& ixy) const {
221 double r2 = std::exp(2.0 * _value);
222 ixx = r2 * (1.0 + distortion.getE1());
223 iyy = r2 * (1.0 - distortion.getE1());
224 ixy = r2 * distortion.getE2();
225}
226
227BaseCore::Jacobian LogTraceRadius::dAssignToQuadrupole(Distortion const& distortion, double& ixx, double& iyy,
228 double& ixy) const {
229 BaseCore::Jacobian result = BaseCore::Jacobian::Zero();
230 double r2 = std::exp(2.0 * _value);
231 ixx = r2 * (1.0 + distortion.getE1());
232 iyy = r2 * (1.0 - distortion.getE1());
233 ixy = r2 * distortion.getE2();
234 result(0, 2) = 2.0 * ixx;
235 result(1, 2) = 2.0 * iyy;
236 result(2, 2) = 2.0 * ixy;
237 result(0, 0) = r2;
238 result(1, 0) = -r2;
239 result(2, 1) = r2;
240 return result;
241}
242} // namespace ellipses
243} // namespace geom
244} // namespace afw
245} // namespace lsst
lsst::afw::geom::ellipses::BaseCore::Jacobian
Eigen::Matrix3d Jacobian
Parameter Jacobian matrix type.
Definition BaseCore.h:64
lsst::afw::geom::ellipses::Distortion
A complex ellipticity with magnitude .
Definition Distortion.h:44
std::exp
T exp(T... args)
std::log
T log(T... args)
std::pow
T pow(T... args)
std::sqrt
T sqrt(T... args)
Generated on Sat May 3 2025 07:30:01 for LSST Applications by doxygen 1.13.2