LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
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 
30 namespace lsst {
31 namespace afw {
32 namespace geom {
33 namespace ellipses {
34 
35 namespace {
36 
38 
39 RegistryMap& getRegistry() {
40  static RegistryMap instance;
41  return instance;
42 }
43 
44 std::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 
67 std::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 
95 void 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 
103 void 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 
111 double 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 
123 double BaseCore::getTraceRadius() const {
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 
144  ParameterVector r;
145  writeParameters(r.data());
146  return r;
147 }
148 
150 
151 bool 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
166 BaseCore& 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 
183 void 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 
193 BaseCore::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 
218 void 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 
231 BaseCore::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
table::Key< std::string > name
Definition: Amplifier.cc:116
int end
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:48
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)
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
T sin(T... args)
T sqrt(T... args)