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