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
BaseCore.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 */
26#include "lsst/geom/Angle.h"
27#include <boost/format.hpp>
28#include <map>
29
30namespace lsst {
31namespace afw {
32namespace geom {
33namespace ellipses {
34
35namespace {
36
38
39RegistryMap& getRegistry() {
40 static RegistryMap instance;
41 return instance;
42}
43
44std::shared_ptr<BaseCore> getRegistryCopy(std::string const& name) {
45 RegistryMap::iterator i = getRegistry().find(name);
46 if (i == getRegistry().end()) {
48 (boost::format("Ellipse core with name '%s' not found in registry.") % name).str());
49 }
50 return i->second->clone();
51}
52
53} // namespace
54
56 std::shared_ptr<BaseCore> result = getRegistryCopy(name);
57 *result = Quadrupole();
58 return result;
59}
60
62 std::shared_ptr<BaseCore> result = getRegistryCopy(name);
63 result->setParameterVector(parameters);
64 return result;
65}
66
67std::shared_ptr<BaseCore> BaseCore::make(std::string const& name, double v1, double v2, double v3) {
68 std::shared_ptr<BaseCore> result = getRegistryCopy(name);
69 result->setParameterVector(ParameterVector(v1, v2, v3));
70 return result;
71}
72
74 std::shared_ptr<BaseCore> result = getRegistryCopy(name);
75 *result = other;
76 return result;
77}
78
80 std::shared_ptr<BaseCore> result = getRegistryCopy(name);
81 other.apply(*result);
82 return result;
83}
84
86 std::shared_ptr<BaseCore> result = getRegistryCopy(name);
87 other.apply(*result);
88 return result;
89}
90
92 getRegistry()[example->getName()] = example;
93}
94
95void BaseCore::grow(double buffer) {
96 double a, b, theta;
97 _assignToAxes(a, b, theta);
98 a += buffer;
99 b += buffer;
100 _assignFromAxes(a, b, theta);
101}
102
103void BaseCore::scale(double factor) {
104 double a, b, theta;
105 _assignToAxes(a, b, theta);
106 a *= factor;
107 b *= factor;
108 _assignFromAxes(a, b, theta);
109}
110
111double BaseCore::getArea() const {
112 double a, b, theta;
113 _assignToAxes(a, b, theta);
114 return a * b * lsst::geom::PI;
115}
116
118 double a, b, theta;
119 _assignToAxes(a, b, theta);
120 return std::sqrt(a * b);
121}
122
124 double ixx, iyy, ixy;
125 _assignToQuadrupole(ixx, iyy, ixy);
126 return std::sqrt(0.5 * (ixx + iyy));
127}
128
130 double a, b, theta;
131 _assignToAxes(a, b, theta);
132 double c = std::cos(theta);
133 double s = std::sin(theta);
134 c *= c;
135 s *= s;
136 b *= b;
137 a *= a;
138 lsst::geom::Extent2D dimensions(std::sqrt(b * s + a * c), std::sqrt(a * s + b * c));
139 dimensions *= 2;
140 return dimensions;
141}
142
145 writeParameters(r.data());
146 return r;
147}
148
150
151bool BaseCore::operator==(BaseCore const& other) const {
152 return getParameterVector() == other.getParameterVector() && getName() == other.getName();
153}
154
156 if (&other != this) {
157 // We use Axes instead of Quadrupole here because it allows us to copy Axes without
158 // implicitly normalizing them.
159 double a, b, theta;
160 other._assignToAxes(a, b, theta);
161 _assignFromAxes(a, b, theta);
162 }
163 return *this;
164}
165// Delegate to copy-constructor for backward-compatibility
166BaseCore& BaseCore::operator=(BaseCore&& other) { return *this = other; }
167
169 if (getName() == other.getName()) {
170 this->operator=(other);
171 return Jacobian::Identity();
172 }
173 // We use Quadrupole instead of Axes here because the ambiguity of the position angle
174 // in the circular case causes some of the Jacobians to/from Axes to be undefined for
175 // exact circles. Quadrupoles don't have that problem, and the Axes-to-Axes case is
176 // handled by the above if block.
177 double ixx, iyy, ixy;
178 Jacobian rhs = other._dAssignToQuadrupole(ixx, iyy, ixy);
179 Jacobian lhs = _dAssignFromQuadrupole(ixx, iyy, ixy);
180 return lhs * rhs;
181}
182
183void BaseCore::_assignQuadrupoleToAxes(double ixx, double iyy, double ixy, double& a, double& b,
184 double& theta) {
185 double xx_p_yy = ixx + iyy;
186 double xx_m_yy = ixx - iyy;
187 double t = std::sqrt(xx_m_yy * xx_m_yy + 4 * ixy * ixy);
188 a = std::sqrt(0.5 * (xx_p_yy + t));
189 b = std::sqrt(0.5 * (xx_p_yy - t));
190 theta = 0.5 * std::atan2(2.0 * ixy, xx_m_yy);
191}
192
193BaseCore::Jacobian BaseCore::_dAssignQuadrupoleToAxes(double ixx, double iyy, double ixy, double& a,
194 double& b, double& theta) {
195 double xx_p_yy = ixx + iyy;
196 double xx_m_yy = ixx - iyy;
197 double t2 = xx_m_yy * xx_m_yy + 4.0 * ixy * ixy;
198 Eigen::Vector3d dt2(2.0 * xx_m_yy, -2.0 * xx_m_yy, 8.0 * ixy);
199 double t = std::sqrt(t2);
200 a = std::sqrt(0.5 * (xx_p_yy + t));
201 b = std::sqrt(0.5 * (xx_p_yy - t));
202 theta = 0.5 * std::atan2(2.0 * ixy, xx_m_yy);
203 Jacobian m = Jacobian::Zero();
204 m(0, 0) = 0.25 * (1.0 + 0.5 * dt2[0] / t) / a;
205 m(0, 1) = 0.25 * (1.0 + 0.5 * dt2[1] / t) / a;
206 m(0, 2) = 0.25 * (0.5 * dt2[2] / t) / a;
207 m(1, 0) = 0.25 * (1.0 - 0.5 * dt2[0] / t) / b;
208 m(1, 1) = 0.25 * (1.0 - 0.5 * dt2[1] / t) / b;
209 m(1, 2) = 0.25 * (-0.5 * dt2[2] / t) / b;
210
211 m.row(2).setConstant(1.0 / (t * t));
212 m(2, 0) *= -ixy;
213 m(2, 1) *= ixy;
214 m(2, 2) *= xx_m_yy;
215 return m;
216}
217
218void BaseCore::_assignAxesToQuadrupole(double a, double b, double theta, double& ixx, double& iyy,
219 double& ixy) {
220 a *= a;
221 b *= b;
222 double c = std::cos(theta);
223 double s = std::sin(theta);
224 ixy = (a - b) * c * s;
225 c *= c;
226 s *= s;
227 ixx = c * a + s * b;
228 iyy = s * a + c * b;
229}
230
231BaseCore::Jacobian BaseCore::_dAssignAxesToQuadrupole(double a, double b, double theta, double& ixx,
232 double& iyy, double& ixy) {
233 Jacobian m;
234 m.col(0).setConstant(2 * a);
235 m.col(1).setConstant(2 * b);
236 a *= a;
237 b *= b;
238 m.col(2).setConstant(a - b);
239 double c = std::cos(theta);
240 double s = std::sin(theta);
241 double cs = c * s;
242 ixy = (a - b) * c * s;
243 c *= c;
244 s *= s;
245 ixx = c * a + s * b;
246 iyy = s * a + c * b;
247 m(0, 0) *= c;
248 m(0, 1) *= s;
249 m(0, 2) *= -2.0 * cs;
250 m(1, 0) *= s;
251 m(1, 1) *= c;
252 m(1, 2) *= 2.0 * cs;
253 m(2, 0) *= cs;
254 m(2, 1) *= -cs;
255 m(2, 2) *= (c - s);
256 return m;
257}
258} // namespace ellipses
259} // namespace geom
260} // namespace afw
261} // namespace lsst
py::object result
Definition _schema.cc:429
int end
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::PointKey< int > dimensions
int m
Definition SpanSet.cc:48
table::Key< int > b
table::Key< int > a
T atan2(T... args)
A temporary-only expression object for ellipse core convolution.
Definition Convolution.h:44
A temporary-only expression object for ellipse core transformations.
Definition Transformer.h:49
A base class for parametrizations of the "core" of an ellipse - the ellipticity and size.
Definition BaseCore.h:55
virtual void writeParameters(double *iter) const =0
Return the size of the bounding box for the ellipse core.
static Jacobian _dAssignAxesToQuadrupole(double a, double b, double theta, double &ixx, double &iyy, double &ixy)
Return the size of the bounding box for the ellipse core.
Definition BaseCore.cc:231
Eigen::Vector3d ParameterVector
Parameter vector type.
Definition BaseCore.h:63
static void registerSubclass(std::shared_ptr< BaseCore > const &example)
Return the size of the bounding box for the ellipse core.
Definition BaseCore.cc:91
void scale(double factor)
Scale the size of the ellipse core by the given factor.
Definition BaseCore.cc:103
virtual void readParameters(double const *iter)=0
Return the size of the bounding box for the ellipse core.
double getArea() const
Return the area of the ellipse core.
Definition BaseCore.cc:111
virtual void _assignToQuadrupole(double &ixx, double &iyy, double &ixy) const =0
Return the size of the bounding box for the ellipse core.
void setParameterVector(ParameterVector const &vector)
Set the core parameters from a vector.
Definition BaseCore.cc:149
virtual Jacobian _dAssignToQuadrupole(double &ixx, double &iyy, double &ixy) const =0
Return the size of the bounding box for the ellipse core.
virtual std::string getName() const =0
Return a string that identifies this parametrization.
virtual void _assignToAxes(double &a, double &b, double &theta) const =0
Return the size of the bounding box for the ellipse core.
double getDeterminantRadius() const
Return the radius defined as the 4th root of the determinant of the quadrupole matrix.
Definition BaseCore.cc:117
bool operator==(BaseCore const &other) const
Compare two ellipse cores for equality.
Definition BaseCore.cc:151
virtual Jacobian _dAssignFromQuadrupole(double ixx, double iyy, double ixy)=0
Return the size of the bounding box for the ellipse core.
lsst::geom::Extent2D computeDimensions() const
Return the size of the bounding box for the ellipse core.
Definition BaseCore.cc:129
static void _assignAxesToQuadrupole(double a, double b, double theta, double &ixx, double &iyy, double &ixy)
Return the size of the bounding box for the ellipse core.
Definition BaseCore.cc:218
static Jacobian _dAssignQuadrupoleToAxes(double ixx, double iyy, double ixy, double &a, double &b, double &theta)
Return the size of the bounding box for the ellipse core.
Definition BaseCore.cc:193
ParameterVector const getParameterVector() const
Return the core parameters as a vector.
Definition BaseCore.cc:143
Jacobian dAssign(BaseCore const &other)
Assign other to this and return the derivative of the conversion, d(this)/d(other).
Definition BaseCore.cc:168
virtual void _assignFromAxes(double a, double b, double theta)=0
Return the size of the bounding box for the ellipse core.
static std::shared_ptr< BaseCore > make(std::string const &name)
Definition BaseCore.cc:55
void grow(double buffer)
Increase the major and minor radii of the ellipse core by the given buffer.
Definition BaseCore.cc:95
static void _assignQuadrupoleToAxes(double ixx, double iyy, double ixy, double &a, double &b, double &theta)
Return the size of the bounding box for the ellipse core.
Definition BaseCore.cc:183
double getTraceRadius() const
Return the radius defined as the square root of one half the trace of the quadrupole matrix.
Definition BaseCore.cc:123
BaseCore & operator=(BaseCore const &other)
Set the parameters of this ellipse core from another.
Definition BaseCore.cc:155
Eigen::Matrix3d Jacobian
Parameter Jacobian matrix type.
Definition BaseCore.h:64
An ellipse core with quadrupole moments as parameters.
Definition Quadrupole.h:47
Reports invalid arguments.
Definition Runtime.h:66
T cos(T... args)
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
T sin(T... args)
T sqrt(T... args)