LSST Applications 26.0.0,g0265f82a02+6660c170cc,g07994bdeae+30b05a742e,g0a0026dc87+17526d298f,g0a60f58ba1+17526d298f,g0e4bf8285c+96dd2c2ea9,g0ecae5effc+c266a536c8,g1e7d6db67d+6f7cb1f4bb,g26482f50c6+6346c0633c,g2bbee38e9b+6660c170cc,g2cc88a2952+0a4e78cd49,g3273194fdb+f6908454ef,g337abbeb29+6660c170cc,g337c41fc51+9a8f8f0815,g37c6e7c3d5+7bbafe9d37,g44018dc512+6660c170cc,g4a941329ef+4f7594a38e,g4c90b7bd52+5145c320d2,g58be5f913a+bea990ba40,g635b316a6c+8d6b3a3e56,g67924a670a+bfead8c487,g6ae5381d9b+81bc2a20b4,g93c4d6e787+26b17396bd,g98cecbdb62+ed2cb6d659,g98ffbb4407+81bc2a20b4,g9ddcbc5298+7f7571301f,ga1e77700b3+99e9273977,gae46bcf261+6660c170cc,gb2715bf1a1+17526d298f,gc86a011abf+17526d298f,gcf0d15dbbd+96dd2c2ea9,gdaeeff99f8+0d8dbea60f,gdb4ec4c597+6660c170cc,ge23793e450+96dd2c2ea9,gf041782ebf+171108ac67
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
result
py::object result
Definition _schema.cc:429
lsst::afw::geom::ellipses::BaseCore::Jacobian
Eigen::Matrix3d Jacobian
Parameter Jacobian matrix type.
Definition BaseCore.h:64
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 Thu Dec 14 2023 04:18:41 for LSST Applications by doxygen 1.9.7