LSSTApplications  18.1.0
LSSTDataManagementBasePackage
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
std::exp
T exp(T... args)
std::log
T log(T... args)
result
py::object result
Definition: schema.cc:418
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
lsst::afw::geom::ellipses::BaseCore::Jacobian
Eigen::Matrix3d Jacobian
Parameter Jacobian matrix type.
Definition: BaseCore.h:64
std::pow
T pow(T... args)
std::sqrt
T sqrt(T... args)
Generated on Thu Aug 8 2019 20:20:06 for LSSTApplications by   doxygen 1.8.13