LSSTApplications  16.0-10-g0ee56ad+5,16.0-11-ga33d1f2+5,16.0-12-g3ef5c14+3,16.0-12-g71e5ef5+18,16.0-12-gbdf3636+3,16.0-13-g118c103+3,16.0-13-g8f68b0a+3,16.0-15-gbf5c1cb+4,16.0-16-gfd17674+3,16.0-17-g7c01f5c+3,16.0-18-g0a50484+1,16.0-20-ga20f992+8,16.0-21-g0e05fd4+6,16.0-21-g15e2d33+4,16.0-22-g62d8060+4,16.0-22-g847a80f+4,16.0-25-gf00d9b8+1,16.0-28-g3990c221+4,16.0-3-gf928089+3,16.0-32-g88a4f23+5,16.0-34-gd7987ad+3,16.0-37-gc7333cb+2,16.0-4-g10fc685+2,16.0-4-g18f3627+26,16.0-4-g5f3a788+26,16.0-5-gaf5c3d7+4,16.0-5-gcc1f4bb+1,16.0-6-g3b92700+4,16.0-6-g4412fcd+3,16.0-6-g7235603+4,16.0-69-g2562ce1b+2,16.0-8-g14ebd58+4,16.0-8-g2df868b+1,16.0-8-g4cec79c+6,16.0-8-gadf6c7a+1,16.0-8-gfc7ad86,16.0-82-g59ec2a54a+1,16.0-9-g5400cdc+2,16.0-9-ge6233d7+5,master-g2880f2d8cf+3,v17.0.rc1
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:284
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 Fri Feb 15 2019 09:17:47 for LSSTApplications by   doxygen 1.8.13