LSST Applications g070148d5b3+33e5256705,g0d53e28543+25c8b88941,g0da5cf3356+2dd1178308,g1081da9e2a+62d12e78cb,g17e5ecfddb+7e422d6136,g1c76d35bf8+ede3a706f7,g295839609d+225697d880,g2e2c1a68ba+cc1f6f037e,g2ffcdf413f+853cd4dcde,g38293774b4+62d12e78cb,g3b44f30a73+d953f1ac34,g48ccf36440+885b902d19,g4b2f1765b6+7dedbde6d2,g5320a0a9f6+0c5d6105b6,g56b687f8c9+ede3a706f7,g5c4744a4d9+ef6ac23297,g5ffd174ac0+0c5d6105b6,g6075d09f38+66af417445,g667d525e37+2ced63db88,g670421136f+2ced63db88,g71f27ac40c+2ced63db88,g774830318a+463cbe8d1f,g7876bc68e5+1d137996f1,g7985c39107+62d12e78cb,g7fdac2220c+0fd8241c05,g96f01af41f+368e6903a7,g9ca82378b8+2ced63db88,g9d27549199+ef6ac23297,gabe93b2c52+e3573e3735,gb065e2a02a+3dfbe639da,gbc3249ced9+0c5d6105b6,gbec6a3398f+0c5d6105b6,gc9534b9d65+35b9f25267,gd01420fc67+0c5d6105b6,geee7ff78d7+a14128c129,gf63283c776+ede3a706f7,gfed783d017+0c5d6105b6,w.2022.47
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 Nov 24 2022 08:01:18 for LSST Applications by doxygen 1.9.5