LSST Applications g063fba187b+cac8b7c890,g0f08755f38+6aee506743,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+b4475c5878,g1dcb35cd9c+8f9bc1652e,g20f6ffc8e0+6aee506743,g217e2c1bcf+73dee94bd0,g28da252d5a+1f19c529b9,g2bbee38e9b+3f2625acfc,g2bc492864f+3f2625acfc,g3156d2b45e+6e55a43351,g32e5bea42b+1bb94961c2,g347aa1857d+3f2625acfc,g35bb328faa+a8ce1bb630,g3a166c0a6a+3f2625acfc,g3e281a1b8c+c5dd892a6c,g3e8969e208+a8ce1bb630,g414038480c+5927e1bc1e,g41af890bb2+8a9e676b2a,g7af13505b9+809c143d88,g80478fca09+6ef8b1810f,g82479be7b0+f568feb641,g858d7b2824+6aee506743,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+2903d499ea,gb58c049af0+d64f4d3760,gc28159a63d+3f2625acfc,gcab2d0539d+b12535109e,gcf0d15dbbd+46a3f46ba9,gda6a2b7d83+46a3f46ba9,gdaeeff99f8+1711a396fd,ge79ae78c31+3f2625acfc,gef2f8181fd+0a71e47438,gf0baf85859+c1f95f4921,gfa517265be+6aee506743,gfa999e8aa5+17cd334064,w.2024.51
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 Sat Dec 21 2024 08:53:13 for LSST Applications by doxygen 1.10.0