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
Ellipse.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2014-2015 AURA/LSST.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <https://www.lsstcorp.org/LegalNotices/>.
21  */
22 
25 
26 #include "lsst/sphgeom/Ellipse.h"
27 
28 #include <cmath>
29 #include <ostream>
30 #include <stdexcept>
31 
32 #include "lsst/sphgeom/Box.h"
33 #include "lsst/sphgeom/Box3d.h"
34 #include "lsst/sphgeom/Circle.h"
36 #include "lsst/sphgeom/codec.h"
37 
38 
39 namespace lsst {
40 namespace sphgeom {
41 
42 Ellipse::Ellipse(UnitVector3d const & f1, UnitVector3d const & f2, Angle alpha) :
43  _a(alpha.asRadians() - 0.5 * PI)
44 {
45  if (alpha.isNan()) {
46  throw std::invalid_argument("Invalid ellipse semi-axis angle");
47  }
48  if (f1 == f2) {
49  _gamma = Angle(0.0);
50  } else if (f1 == -f2) {
51  _gamma = Angle(0.5 * PI);
52  } else {
53  _gamma = 0.5 * NormalizedAngle(f1, f2);
54  }
55  if (isEmpty()) {
56  *this = empty();
57  return;
58  } else if (isFull()) {
59  *this = full();
60  return;
61  }
62  if (_gamma.asRadians() == 0.0) {
63  // The foci are identical, so this ellipse is a circle centered at
64  // the common focal point. Pick an orthonormal basis containing f1
65  // and use it to construct an orthogonal matrix that maps f1 to
66  // (0, 0, 1).
68  UnitVector3d b1 = UnitVector3d(f1.cross(b0));
69  _S = Matrix3d(b0.x(), b0.y(), b0.z(),
70  b1.x(), b1.y(), b1.z(),
71  f1.x(), f1.y(), f1.z());
72  _b = _a;
73  _tana = std::fabs(tan(_a));
74  _tanb = _tana;
75  return;
76  }
77  // _gamma != 0 implies that f1 - f2 != 0. Also, if f1 = -f2 then
78  // _gamma = PI/2, and the ellipse must either empty or full. So
79  // at this stage f1 + f2 != 0.
80  Vector3d b0 = f1 - f2;
81  Vector3d b2 = f1 + f2;
82  Vector3d b1 = b0.cross(b2);
83  b0.normalize();
84  b1.normalize();
85  b2.normalize();
86  _S = Matrix3d(b0.x(), b0.y(), b0.z(),
87  b1.x(), b1.y(), b1.z(),
88  b2.x(), b2.y(), b2.z());
89  // Compute _b.
90  double r = std::min(1.0, std::max(-1.0, cos(alpha) / cos(_gamma)));
91  _b = Angle(std::acos(r) - 0.5 * PI);
92  if (_a.asRadians() <= 0.0 && _b > _a) {
93  _b = _a;
94  } else if (_a.asRadians() > 0.0 && _b < _a) {
95  _b = _a;
96  }
97  _tana = std::fabs(tan(_a));
98  _tanb = std::fabs(tan(_b));
99  return;
100 }
101 
103  Angle alpha,
104  Angle beta,
106 {
107  if (!std::isfinite(orientation.asRadians())) {
108  throw std::invalid_argument("Invalid ellipse orientation");
109  }
110  if (alpha.isNan() ||
111  beta.isNan() ||
112  (alpha.asRadians() < 0.5 * PI && beta.asRadians() >= 0.5 * PI) ||
113  (alpha.asRadians() > 0.5 * PI && beta.asRadians() <= 0.5 * PI) ||
114  (alpha.asRadians() == 0.5 * PI && beta.asRadians() != 0.5 * PI)) {
115  throw std::invalid_argument("Invalid ellipse semi-axis angle(s)");
116  }
117  if (alpha.asRadians() < 0.0 || beta.asRadians() < 0.0) {
118  *this = empty();
119  return;
120  } else if (alpha.asRadians() > PI || beta.asRadians() > PI ||
121  (alpha.asRadians() == PI && beta.asRadians() == PI)) {
122  *this = full();
123  return;
124  }
125  if (alpha == beta) {
126  // The ellipse is a circle. Pick an orthonormal basis containing the
127  // center and use it to construct some orthogonal matrix that maps the
128  // center to (0, 0, 1).
130  UnitVector3d b1 = UnitVector3d(center.cross(b0));
131  _S = Matrix3d(b0.x(), b0.y(), b0.z(),
132  b1.x(), b1.y(), b1.z(),
133  center.x(), center.y(), center.z());
134  _a = alpha - Angle(0.5 * PI);
135  _b = _a;
136  _gamma = Angle(0.0);
137  _tana = std::fabs(tan(_a));
138  _tanb = _tana;
139  return;
140  }
141  if ((alpha.asRadians() < 0.5 * PI && alpha < beta) ||
142  (alpha.asRadians() > 0.5 * PI && alpha > beta)) {
143  std::swap(alpha, beta);
144  orientation = orientation + Angle(0.5 * PI);
145  }
146  UnitVector3d b0 =
148  UnitVector3d b1 = UnitVector3d(b0.cross(center));
149  _S = Matrix3d(b0.x(), b0.y(), b0.z(),
150  b1.x(), b1.y(), b1.z(),
151  center.x(), center.y(), center.z());
152  _a = alpha - Angle(0.5 * PI);
153  _b = beta - Angle(0.5 * PI);
154  double d = std::min(1.0, std::max(-1.0, cos(alpha) / cos(beta)));
155  _gamma = Angle(std::acos(d));
156  _tana = std::fabs(tan(_a));
157  _tanb = std::fabs(tan(_b));
158 }
159 
160 bool Ellipse::contains(UnitVector3d const & v) const {
161  UnitVector3d const c = getCenter();
162  double vdotc = v.dot(c);
163  Vector3d u;
164  double scz;
165  // To maintain high accuracy for very small and very large ellipses,
166  // decompose v as v = u ± c near ±c. Then S v = S u ± S c, and
167  // S c = (0, 0, 1).
168  if (vdotc > 0.5) {
169  u = v - c;
170  scz = 1.0;
171  } else if (vdotc < -0.5) {
172  u = v + c;
173  scz = -1.0;
174  } else {
175  u = v;
176  scz = 0.0;
177  }
178  u = _S * u;
179  double x = u.x() * _tana;
180  double y = u.y() * _tanb;
181  double z = u.z() + scz;
182  double d = (x * x + y * y) - z * z;
183  if (_a.asRadians() > 0.0) {
184  return z >= 0.0 || d >= 0.0;
185  } else {
186  return z >= 0.0 && d <= 0.0;
187  }
188 }
189 
191  // For now, simply return the bounding box of the ellipse bounding circle.
192  //
193  // Improving on this seems difficult, mainly because error bounds must be
194  // computed to guarantee that the resulting box tightly bounds the ellipse.
195  // In case this ends up being important and someone wants to go down
196  // this route in the future, here are some thoughts on how to proceed.
197  //
198  // First of all, if the ellipse contains a pole, its bounding box must
199  // contain all longitudes. Otherwise, consider the plane spanned by the
200  // basis vectors u = (0, 0, 1) and v = (cos θ, sin θ, 0). To obtain
201  // longitude bounds for the ellipse, we must find values of θ for which
202  // this plane is tangent to the elliptical cone that defines the spherical
203  // ellipse boundary. This is the case when the plane intersects the cone
204  // in a single line.
205  //
206  // Let μ v + λ u be the direction of the line of intersection, where
207  // μ, λ ∈ ℝ. We know that μ ≠ 0 because u is not on the ellipse boundary,
208  // so fix μ = 1. Let Q be the symmetric matrix representation of the
209  // ellipse. Then:
210  //
211  // (v + λ u)ᵀ Q (v + λ u) = 0
212  //
213  // Expanding gives:
214  //
215  // (vᵀ Q v) + λ (uᵀ Q v) + λ (vᵀ Q u) + λ² (uᵀ Q u) = 0
216  //
217  // By the symmetry of Q:
218  //
219  // λ² (uᵀ Q u) + 2λ (uᵀ Q v) + (vᵀ Q v) = 0
220  //
221  // This is a quadratic equation which has one solution exactly when its
222  // discriminant is 0, i.e. when:
223  //
224  // (uᵀ Q v)² - (uᵀ Q u) (vᵀ Q v) = 0
225  //
226  // Substituting for u and v and simplifying yields:
227  //
228  // (Q₁₂² - Q₁₁Q₂₂)tan²θ + 2(Q₀₂Q₁₂ - Q₀₁Q₂₂)tan θ + (Q₀₂² - Q₀₀Q₂₂) = 0
229  //
230  // Finding the latitude bounds can be done by parameterizing the ellipse
231  // boundary and then solving for the zeros of the derivative of z with
232  // respect to the parameter. This looks to be more involved than the
233  // longitude bound calculation, and I haven't worked through the details.
235 }
236 
239 }
240 
242  Angle r = std::max(getAlpha(), getBeta()) + 2.0 * Angle(MAX_ASIN_ERROR);
243  return Circle(getCenter(), r);
244 }
245 
247  return getBoundingCircle().relate(b) & (DISJOINT | WITHIN);
248 }
249 
250 // For now, implement ellipse-circle and ellipse-ellipse relation
251 // computation by approximating ellipses via their bounding circles.
252 //
253 // It should be possible to improve on this using the following algorithm to
254 // compute ellipse-ellipse intersection points.
255 //
256 // Ellipses that are neither empty nor full have quadratic forms with
257 // symmetric matrix representations P, Q of full rank. Consider the matrix
258 // pencil μ P + λ Q, where neither μ, λ ∈ ℝ is zero. This is a family of
259 // quadratic forms, where every member passes through the intersection points
260 // of P and Q. The scalars μ, λ are homogeneous, meaning that the quadratic
261 // forms for (μ, λ) and (k μ, k λ) are identical for k ≠ 0, so we can fix μ = 1.
262 //
263 // If we can find λ such that R = P - λ Q is rank deficient, then the resulting
264 // quadratic form corresponds to a line or to a pair of (possibly identical)
265 // planes. The intersection of a plane or a line with either P or Q is
266 // then easily computed (see below). Finding λ is a generalized eigenvalue
267 // problem, and one way to solve it is to note that:
268 //
269 // det(P - λ Q) = 0 → det(PQ⁻¹ - λ I) = 0
270 //
271 // so that the values of λ we are interested in are the eigenvalues of PQ⁻¹.
272 // Use one of these eigenvalues to construct R. Then:
273 //
274 // 1) If R has rank 0, P and Q are equivalent up to scale and the corresponding
275 // quadratic forms are identical.
276 //
277 // 2) If R has rank 1, then the quadratic form R can be factorized as
278 // (ax + by + cz)² = 0. Why? There is an eigen-decomposition of R,
279 // R = M D Mᵀ, where D is diagonal, M is orthogonal, D₀₀ is the only
280 // non-zero eigenvalue and the first column of M is the unit eigenvector
281 // n for D₀₀. Now
282 //
283 // vᵀ R v = (Mᵀ v)ᵀ D (Mᵀ v) = (n·v)² = 0
284 //
285 // This just says that there is some basis in which the quadratic form
286 // is x² = 0, which is the plane defined by x = 0.
287 //
288 // 3) Otherwise, R has rank 2, and corresponds to a line or to a pair of
289 // planes. To see why, again consider its eigen-decomposition R = M D Mᵀ,
290 // where M is orthogonal, D is diagonal with non-zero eigenvalues D₀₀ and
291 // D₁₁, and the first two columns of M are the corresponding eigenvectors.
292 //
293 // If D₀₀ and D₁₁ have the same sign, then there is a basis in which
294 // the quadratic form is a² x² + b² y² = 0, i.e. the z-axis. In the
295 // original basis this is a line perpendicular to the two eigenvectors of
296 // R where P and Q are tangent.
297 //
298 // If D₀₀ and D₁₁ have opposite signs, then there is some basis in which
299 // the quadratic form reduces to a² x² - b² y² = (a x - b y)(a x + b y) = 0;
300 // that is, solutions lie on two planes.
301 //
302 // To compute the intersection of an ellipse with quadratic form Q and a
303 // plane with unit normal n, we must solve vᵀ Q v = 0 subject to v·n = 0.
304 // Pick two unit vectors b₀, b₁ orthogonal to n, write v = μ b₀ + λ b₁, and
305 // substitute into the equation of the quadratic form to obtain:
306 //
307 // μ² (b₀ᵀ Q b₀) + 2μλ (b₁ᵀ Q b₀) + λ² (b₁ᵀ Q b₁) = 0
308 //
309 // This only has solutions when:
310 //
311 // (b₁ᵀ Q b₀)² ≥ (b₀ᵀ Q b₀) (b₁ᵀ Q b₁)
312 //
313 // As in the case of the ellipse bounding box computation, actually using the
314 // above in an implementation of relate() involves a non-trivial error analysis.
315 //
316 // TODO(smm): investigate the theory of matrix pencils and generalized
317 // eigenvalue problems. A couple unanswered questions this could shed some
318 // light on are:
319 //
320 // - What algorithm should be used to solve for the generalized eigenvalues?
321 // Note that very small and very large ellipses will have matrix
322 // representations with very large condition numbers.
323 // - Which of the generalized eigenvalues should be chosen / leads to the most
324 // accurate computation? Is there some usefully exploitable relationship
325 // between them and the degenerate quadratic forms they engender?
326 
328  return getBoundingCircle().relate(c) & (DISJOINT | WITHIN);
329 }
330 
332  return getBoundingCircle().relate(p) & (DISJOINT | WITHIN);
333 }
334 
336  return getBoundingCircle().relate(e.getBoundingCircle()) & DISJOINT;
337 }
338 
340  std::vector<uint8_t> buffer;
341  uint8_t tc = TYPE_CODE;
342  buffer.reserve(ENCODED_SIZE);
343  buffer.push_back(tc);
344  for (int r = 0; r < 3; ++r) {
345  for (int c = 0; c < 3; ++c) {
346  encodeDouble(_S(r, c), buffer);
347  }
348  }
349  encodeDouble(_a.asRadians(), buffer);
350  encodeDouble(_b.asRadians(), buffer);
351  encodeDouble(_gamma.asRadians(), buffer);
352  encodeDouble(_tana, buffer);
353  encodeDouble(_tanb, buffer);
354  return buffer;
355 }
356 
357 std::unique_ptr<Ellipse> Ellipse::decode(uint8_t const * buffer, size_t n) {
358  if (buffer == nullptr || n != ENCODED_SIZE || buffer[0] != TYPE_CODE) {
359  throw std::runtime_error("Byte-string is not an encoded Ellipse");
360  }
361  std::unique_ptr<Ellipse> ellipse(new Ellipse);
362  ++buffer;
363  double m00 = decodeDouble(buffer); buffer += 8;
364  double m01 = decodeDouble(buffer); buffer += 8;
365  double m02 = decodeDouble(buffer); buffer += 8;
366  double m10 = decodeDouble(buffer); buffer += 8;
367  double m11 = decodeDouble(buffer); buffer += 8;
368  double m12 = decodeDouble(buffer); buffer += 8;
369  double m20 = decodeDouble(buffer); buffer += 8;
370  double m21 = decodeDouble(buffer); buffer += 8;
371  double m22 = decodeDouble(buffer); buffer += 8;
372  ellipse->_S = Matrix3d(m00, m01, m02,
373  m10, m11, m12,
374  m20, m21, m22);
375  double a = decodeDouble(buffer); buffer += 8;
376  double b = decodeDouble(buffer); buffer += 8;
377  double gamma = decodeDouble(buffer); buffer += 8;
378  ellipse->_a = Angle(a);
379  ellipse->_b = Angle(b);
380  ellipse->_gamma = Angle(gamma);
381  double tana = decodeDouble(buffer); buffer += 8;
382  double tanb = decodeDouble(buffer); buffer += 8;
383  ellipse->_tana = tana;
384  ellipse->_tanb = tanb;
385  return ellipse;
386 }
387 
389  os << "{\"Ellipse\": [" << e.getTransformMatrix() << ", "
390  << e.getAlpha() << ", " << e.getBeta() << "]}";
391  return os;
392 }
393 
394 }} // namespace lsst::sphgeom
This file declares a class for representing axis-aligned bounding boxes in ℝ³.
double x
This file declares a class for representing circular regions on the unit sphere.
This file declares a class for representing convex polygons with great circle edges on the unit spher...
double z
Definition: Match.cc:44
std::ostream * os
Definition: Schema.cc:557
int y
Definition: SpanSet.cc:48
table::Key< int > b
table::Key< int > a
T acos(T... args)
Angle represents an angle in radians.
Definition: Angle.h:43
bool isNan() const
isNan returns true if the angle value is NaN.
Definition: Angle.h:91
double asRadians() const
asRadians returns the value of this angle in units of radians.
Definition: Angle.h:85
Box3d represents a box in ℝ³.
Definition: Box3d.h:42
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition: Box.h:54
Circle is a circular region on the unit sphere that contains its boundary.
Definition: Circle.h:46
Box getBoundingBox() const override
getBoundingBox returns a bounding-box for this region.
Definition: Circle.cc:211
Relationship relate(UnitVector3d const &v) const
Definition: Circle.cc:268
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
Definition: Circle.cc:219
ConvexPolygon is a closed convex polygon on the unit sphere.
Definition: ConvexPolygon.h:57
Ellipse is an elliptical region on the sphere.
Definition: Ellipse.h:170
Matrix3d const & getTransformMatrix() const
getTransformMatrix returns the orthogonal matrix that maps vectors to the basis in which the quadrati...
Definition: Ellipse.h:230
static constexpr uint8_t TYPE_CODE
Definition: Ellipse.h:172
Angle getAlpha() const
getAlpha returns α, the first semi-axis length of the ellipse.
Definition: Ellipse.h:252
static Ellipse full()
Definition: Ellipse.h:176
Angle getBeta() const
getBeta returns β, the second semi-axis length of the ellipse.
Definition: Ellipse.h:257
static std::unique_ptr< Ellipse > decode(std::vector< uint8_t > const &s)
Definition: Ellipse.h:303
Relationship relate(Region const &r) const override
Definition: Ellipse.h:289
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Definition: Ellipse.cc:241
Box getBoundingBox() const override
getBoundingBox returns a bounding-box for this region.
Definition: Ellipse.cc:190
std::vector< uint8_t > encode() const override
encode serializes this region into an opaque byte string.
Definition: Ellipse.cc:339
Ellipse()
This constructor creates an empty ellipse.
Definition: Ellipse.h:179
virtual bool contains(UnitVector3d const &) const=0
contains tests whether the given unit vector is inside this region.
bool isFull() const
Definition: Ellipse.h:221
static Ellipse empty()
Definition: Ellipse.h:174
bool isEmpty() const
Definition: Ellipse.h:219
UnitVector3d getCenter() const
getCenter returns the center of the ellipse as a unit vector.
Definition: Ellipse.h:233
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
Definition: Ellipse.cc:237
A 3x3 matrix with real entries stored in double precision.
Definition: Matrix3d.h:38
NormalizedAngle is an angle that lies in the range [0, 2π), with one exception - a NormalizedAngle ca...
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
double dot(Vector3d const &v) const
dot returns the inner product of this unit vector and v.
Definition: UnitVector3d.h:152
UnitVector3d rotatedAround(UnitVector3d const &k, Angle a) const
rotatedAround returns a copy of this unit vector, rotated around the unit vector k by angle a accordi...
Definition: UnitVector3d.h:193
static UnitVector3d orthogonalTo(Vector3d const &v)
orthogonalTo returns an arbitrary unit vector that is orthogonal to v.
Definition: UnitVector3d.cc:34
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this unit vector and v.
Definition: UnitVector3d.h:155
static UnitVector3d northFrom(Vector3d const &v)
northFrom returns the unit vector orthogonal to v that points "north" from v.
Definition: UnitVector3d.cc:51
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
double x() const
Definition: Vector3d.h:66
double y() const
Definition: Vector3d.h:68
double normalize()
normalize scales this vector to have unit norm and returns its norm prior to scaling.
Definition: Vector3d.cc:41
double z() const
Definition: Vector3d.h:70
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this vector and v.
Definition: Vector3d.h:101
This file contains simple helper functions for encoding and decoding primitive types to/from byte str...
T fabs(T... args)
T isfinite(T... args)
T max(T... args)
T min(T... args)
lsst::geom::Angle Angle
Definition: misc.h:33
constexpr double MAX_ASIN_ERROR
Definition: constants.h:45
void encodeDouble(double item, std::vector< uint8_t > &buffer)
encode appends an IEEE double in little-endian byte order to the end of buffer.
Definition: codec.h:46
std::ostream & operator<<(std::ostream &, Angle const &)
Definition: Angle.cc:34
int orientation(UnitVector3d const &a, UnitVector3d const &b, UnitVector3d const &c)
orientation computes and returns the orientations of 3 unit vectors a, b and c.
Definition: orientation.cc:135
double tan(Angle const &a)
Definition: Angle.h:104
double decodeDouble(uint8_t const *buffer)
decode extracts an IEEE double from the 8 byte little-endian byte sequence in buffer.
Definition: codec.h:66
double cos(Angle const &a)
Definition: Angle.h:103
constexpr double PI
Definition: constants.h:36
A base class for image defects.
T push_back(T... args)
T reserve(T... args)
This file declares a class for representing longitude/latitude angle boxes on the unit sphere.
This file declares a class for representing elliptical regions on the unit sphere.
T swap(T... args)