LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+b203dec576,g18429d2f64+358861cd2c,g199a45376c+0ba108daf9,g1fd858c14a+dd066899e3,g262e1987ae+ebfced1d55,g29ae962dfc+72fd90588e,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+b668f15bc5,g4595892280+3897dae354,g47891489e3+abcf9c3559,g4d44eb3520+fb4ddce128,g53246c7159+8c5ae1fdc5,g67b6fd64d1+abcf9c3559,g67fd3c3899+1f72b5a9f7,g74acd417e5+cb6b47f07b,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+abcf9c3559,g8d7436a09f+bcf525d20c,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+6054cc57f1,g97be763408+06f794da49,g9dd6db0277+1f72b5a9f7,ga681d05dcb+7e36ad54cd,gabf8522325+735880ea63,gac2eed3f23+abcf9c3559,gb89ab40317+abcf9c3559,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+1f72b5a9f7,gdab6d2f7ff+cb6b47f07b,gdc713202bf+1f72b5a9f7,gdfd2d52018+8225f2b331,ge365c994fd+375fc21c71,ge410e46f29+abcf9c3559,geaed405ab2+562b3308c0,gf9a733ac38+8c5ae1fdc5,w.2025.35
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
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 Wed Sep 3 2025 07:01:18 for LSST Applications by doxygen 1.13.2