LSST Applications  21.0.0+c4f5df5339,21.0.0+e70536a077,21.0.0-1-ga51b5d4+7c60f8a6ea,21.0.0-10-g560fb7b+411cd868f8,21.0.0-10-gcf60f90+8c49d86aa0,21.0.0-13-gc485e61d+38156233bf,21.0.0-16-g7a993c7b9+1041c3824f,21.0.0-2-g103fe59+d9ceee3e5a,21.0.0-2-g1367e85+0b2f7db15a,21.0.0-2-g45278ab+e70536a077,21.0.0-2-g5242d73+0b2f7db15a,21.0.0-2-g7f82c8f+feb9862f5e,21.0.0-2-g8f08a60+9c9a9cfcc8,21.0.0-2-ga326454+feb9862f5e,21.0.0-2-gde069b7+bedfc5e1fb,21.0.0-2-gecfae73+417509110f,21.0.0-2-gfc62afb+0b2f7db15a,21.0.0-3-g21c7a62+a91f7c0b59,21.0.0-3-g357aad2+062581ff1a,21.0.0-3-g4be5c26+0b2f7db15a,21.0.0-3-g65f322c+85aa0ead76,21.0.0-3-g7d9da8d+c4f5df5339,21.0.0-3-gaa929c8+411cd868f8,21.0.0-3-gc44e71e+fd4029fd48,21.0.0-3-ge02ed75+5d9b90b8aa,21.0.0-38-g070523fc+44fda2b515,21.0.0-4-g591bb35+5d9b90b8aa,21.0.0-4-g88306b8+3cdc83ea97,21.0.0-4-gc004bbf+d52368b591,21.0.0-4-gccdca77+a5c54364a0,21.0.0-5-g7ebb681+81e2098694,21.0.0-5-gdf36809+87b8d260e6,21.0.0-6-g2d4f3f3+e70536a077,21.0.0-6-g4e60332+5d9b90b8aa,21.0.0-6-g5ef7dad+3f4e29eeae,21.0.0-7-gc8ca178+0f5e56d48f,21.0.0-9-g9eb8d17+cc2c7a81aa,master-gac4afde19b+5d9b90b8aa,w.2021.07
LSST Data Management Base Package
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 
27 namespace lsst {
28 namespace afw {
29 namespace geom {
30 namespace ellipses {
31 
32 void 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 
38 BaseCore::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 
58 void 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 
67 BaseCore::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 
89 void 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 
95 BaseCore::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 
113 void 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 
121 BaseCore::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 
138 void 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 
144 BaseCore::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 
164 void 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 
173 BaseCore::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 
195 void 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 
201 BaseCore::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 
219 void 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 
227 BaseCore::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)
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:194
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
std::pow
T pow(T... args)
std::sqrt
T sqrt(T... args)
Generated on Wed Feb 17 2021 08:13:24 for LSST Applications by doxygen 1.9.1