LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Static Public Attributes | List of all members
lsst::sphgeom::Ellipse Class Reference

Ellipse is an elliptical region on the sphere. More...

#include <Ellipse.h>

Inheritance diagram for lsst::sphgeom::Ellipse:
lsst::sphgeom::Region

Public Member Functions

 Ellipse ()
 This constructor creates an empty ellipse.
 
 Ellipse (Circle const &c)
 This constructor creates an ellipse corresponding to the given circle.
 
 Ellipse (UnitVector3d const &v, Angle alpha=Angle(0.0))
 This constructor creates an ellipse corresponding to the circle with the given center and opening angle.
 
 Ellipse (UnitVector3d const &f1, UnitVector3d const &f2, Angle alpha)
 This constructor creates an ellipse with the given foci and semi-axis angle.
 
 Ellipse (UnitVector3d const &center, Angle alpha, Angle beta, Angle orientation)
 This constructor creates an ellipse with the given center, semi-axis angles, and orientation.
 
bool operator== (Ellipse const &e) const
 
bool operator!= (Ellipse const &e) const
 
bool isEmpty () const
 
bool isFull () const
 
bool isGreatCircle () const
 
bool isCircle () const
 
Matrix3d const & getTransformMatrix () const
 getTransformMatrix returns the orthogonal matrix that maps vectors to the basis in which the quadratic form corresponding to this ellipse is diagonal.
 
UnitVector3d getCenter () const
 getCenter returns the center of the ellipse as a unit vector.
 
UnitVector3d getF1 () const
 getF1 returns the first focal point of the ellipse.
 
UnitVector3d getF2 () const
 getF2 returns the second focal point of the ellipse.
 
Angle getAlpha () const
 getAlpha returns α, the first semi-axis length of the ellipse.
 
Angle getBeta () const
 getBeta returns β, the second semi-axis length of the ellipse.
 
Angle getGamma () const
 getGamma returns ɣ ∈ [0, π/2], half of the angle between the foci.
 
Ellipsecomplement ()
 complement sets this ellipse to the closure of its complement.
 
Ellipse complemented () const
 complemented returns the closure of the complement of this ellipse.
 
std::unique_ptr< Regionclone () const override
 clone returns a deep copy of this region.
 
Box getBoundingBox () const override
 getBoundingBox returns a bounding-box for this region.
 
Box3d getBoundingBox3d () const override
 getBoundingBox3d returns a 3-dimensional bounding-box for this region.
 
Circle getBoundingCircle () const override
 getBoundingCircle returns a bounding-circle for this region.
 
bool contains (UnitVector3d const &v) const override
 contains tests whether the given unit vector is inside this region.
 
Relationship relate (Region const &r) const override
 
Relationship relate (Box const &) const override
 
Relationship relate (Circle const &) const override
 
Relationship relate (ConvexPolygon const &) const override
 
Relationship relate (Ellipse const &) const override
 
std::vector< uint8_t > encode () const override
 encode serializes this region into an opaque byte string.
 
virtual bool contains (UnitVector3d const &) const=0
 contains tests whether the given unit vector is inside this region.
 
bool contains (double x, double y, double z) const
 contains tests whether the unit vector defined by the given (not necessarily normalized) coordinates is inside this region.
 
bool contains (double lon, double lat) const
 contains tests whether the unit vector defined by the given longitude and latitude coordinates (in radians) is inside this region.
 

Static Public Member Functions

static Ellipse empty ()
 
static Ellipse full ()
 
static std::unique_ptr< Ellipsedecode (std::vector< uint8_t > const &s)
 
static std::unique_ptr< Ellipsedecode (uint8_t const *buffer, size_t n)
 

Static Public Attributes

static constexpr uint8_t TYPE_CODE = 'e'
 

Detailed Description

Ellipse is an elliptical region on the sphere.

Mathematical Definition

A spherical ellipse is defined as the set of unit vectors v such that:

d(v,f₁) + d(v,f₂) ≤ 2α                           (Eq. 1)

where f₁ and f₂ are unit vectors corresponding to the foci of the ellipse, d is the function that returns the angle between its two input vectors, and α is a constant.

If 2α < d(f₁,f₂), no point in S² satisfies the inequality, and the ellipse is empty. If f₁ = f₂, the ellipse is a circle with opening angle α. The ellipse defined by foci -f₁ and -f₂, and angle π - α satisfies:

d(v,-f₁) + d(v,-f₂) ≤ 2(π - α)                 →
π - d(v,f₁) + π - d(v,f₂) ≤ 2π - 2α            →
d(v,f₁) + d(v,f₂) ≥ 2α

In other words, it is the closure of the complement of the ellipse defined by f₁, f₂ and α. Therefore if 2π - 2α ≤ d(f₁,f₂), all points in S² satisfy Eq 1. and we say that the ellipse is full.

Consider now the equation d(v,f₁) + d(v,f₂) = 2α for v ∈ ℝ³. We know that

cos(d(v,fᵢ)) = (v·fᵢ)/(‖v‖‖fᵢ‖)
             = (v·fᵢ)/‖v‖            (since ‖fᵢ‖ = 1)

and, because sin²θ + cos²θ = 1 and ‖v‖² = v·v,

sin(d(v,fᵢ)) = √(v·v - (v·fᵢ)²)/‖v‖

Starting with:

d(v,f₁) + d(v,f₂) = 2α

we take the cosine of both sides, apply the angle sum identity for cosine, and substitute the expressions above to obtain:

cos(d(v,f₁) + d(v,f₂)) = cos 2α                                   →
cos(d(v,f₁)) cos(d(v,f₂)) - sin(d(v,f₁)) sin(d(v,f₂)) = cos 2α    →
(v·f₁) (v·f₂) - √(v·v - (v·f₁)²) √(v·v - (v·f₂)²) = cos 2α (v·v)

Rearranging to place the square roots on the RHS, squaring both sides, and simplifying:

((v·f₁) (v·f₂) - cos 2α (v·v))² = (v·v - (v·f₁)²) (v·v - (v·f₂)²) →
cos²2α (v·v) - 2 cos 2α (v·f₁) (v·f₂) = (v·v) - (v·f₁)² - (v·f₂)² →

sin²2α (v·v) + 2 cos 2α (v·f₁) (v·f₂) - (v·f₁)² - (v·f₂)² = 0   (Eq. 2)

Note in particular that if α = π/2, the above simplifies to:

 (v·f₁ + v·f₂)² = 0    ↔    v·(f₁ + f₂) = 0

That is, the equation describes the great circle obtained by intersecting S² with the plane having normal vector f₁ + f₂.

Writing v = (x, y, z) and substituting into Eq. 2, we see that the LHS is a homogeneous polynomial of degree two in 3 variables, or a ternary quadratic form. The matrix representation of this quadratic form is the symmetric 3 by 3 matrix Q such that:

vᵀ Q v = 0

is equivalent to Eq. 2. Consider now the orthonormal basis vectors:

b₀ = (f₁ - f₂)/‖f₁ - f₂‖
b₁ = (f₁ × f₂)/‖f₁ × f₂‖
b₂ = (f₁ + f₂)/‖f₁ + f₂‖

where x denotes the vector cross product. Let S be the matrix with these basis vectors as rows. Given coordinates u in this basis, we have v = Sᵀ u, and:

(Sᵀ u)ᵀ Q (Sᵀ u) = 0    ↔    uᵀ (S Q Sᵀ) u = 0

We now show that D = S Q Sᵀ is diagonal. Let d(f₁,f₂) = 2ɣ. The coordinates of f₁ and f₂ in this new basis are f₁ = (sin ɣ, 0, cos ɣ) and f₂ = (-sin ɣ, 0, cos ɣ). Writing u = (x, y, z) and substituting into Eq. 2:

 sin²2α (u·u) + 2 cos 2α (u·f₁) (u·f₂) - (u·f₁)² - (u·f₂)² = 0

we obtain:

 (sin²2α - 2 cos 2α sin²ɣ - 2 sin²ɣ) x² + (sin²2α) y² +
 (sin²2α + 2 cos 2α cos²ɣ - 2 cos²ɣ) z² = 0

Now sin²2α = 4 sin²α cos²α, cos 2α = cos²α - sin²α, so that:

 (cos²α (sin²α - sin²ɣ)) x² + (sin²α cos²α) y² +
 (sin²α (cos²α - cos²ɣ)) z² = 0

Dividing by sin²α (cos²ɣ - cos²α), and letting cos β = cos α / cos ɣ:

  x² cot²α + y² cot²β - z² = 0              (Eq. 3)

This says that the non-zero elements of S Q Sᵀ are on the diagonal and equal to (cot²α, cot²β, -1) up to scale. In other words, the boundary of a spherical ellipse is given by the intersection of S² and an elliptical cone in ℝ³ passing through the origin. Because z = 0 → x,y = 0 it is evident that the boundary of a spherical ellipse is hemispherical.

If 0 < α < π/2, then β ≤ α, and α is the semi-major axis angle of the ellipse while β is the semi-minor axis angle.

If α = π/2, then the spherical ellipse corresponds to a hemisphere.

If π/2 < α < π, then β ≥ α, and α is the semi-minor axis angle of the ellipse, while β is the semi-major axis angle.

Implementation

Internal state consists of the orthogonal transformation matrix S that maps the ellipse center to (0, 0, 1), as well as |cot α| and |cot β| (enough to reconstruct D, and hence Q), and α, β, ɣ.

In fact, a = α - π/2, b = β - π/2 are stored instead of α and β. This is for two reasons. The first is that when taking the complement of an ellipse, α is mapped to π - α but a is mapped to -a (and b → -b). As a result, taking the complement can be implemented using only changes of sign, and is therefore exact. The other reason is that |cot(α)| = |tan(a)|, and tan is more convenient numerically. In particular, cot(0) is undefined, but tan is finite since a is rational and cannot be exactly equal to ±π/2.

Definition at line 177 of file Ellipse.h.

Constructor & Destructor Documentation

◆ Ellipse() [1/5]

lsst::sphgeom::Ellipse::Ellipse ( )
inline

This constructor creates an empty ellipse.

Definition at line 186 of file Ellipse.h.

186 :
187 _S(1.0),
188 _a(-2.0),
189 _b(-2.0),
190 _gamma(0.0),
193 {}

◆ Ellipse() [2/5]

lsst::sphgeom::Ellipse::Ellipse ( Circle const & c)
inlineexplicit

This constructor creates an ellipse corresponding to the given circle.

Definition at line 196 of file Ellipse.h.

196 {
197 *this = Ellipse(c.getCenter(), c.getCenter(), c.getOpeningAngle());
198 }
Ellipse()
This constructor creates an empty ellipse.
Definition Ellipse.h:186

◆ Ellipse() [3/5]

lsst::sphgeom::Ellipse::Ellipse ( UnitVector3d const & v,
Angle alpha = Angle(0.0) )
inlineexplicit

This constructor creates an ellipse corresponding to the circle with the given center and opening angle.

Definition at line 202 of file Ellipse.h.

202 {
203 *this = Ellipse(v, v, alpha);
204 }

◆ Ellipse() [4/5]

lsst::sphgeom::Ellipse::Ellipse ( UnitVector3d const & f1,
UnitVector3d const & f2,
Angle alpha )

This constructor creates an ellipse with the given foci and semi-axis angle.

Definition at line 49 of file Ellipse.cc.

49 :
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).
74 UnitVector3d b0 = UnitVector3d::orthogonalTo(f1);
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}
T acos(T... args)
double asRadians() const
asRadians returns the value of this angle in units of radians.
Definition Angle.h:92
static Ellipse full()
Definition Ellipse.h:183
bool isFull() const
Definition Ellipse.h:228
static Ellipse empty()
Definition Ellipse.h:181
bool isEmpty() const
Definition Ellipse.h:226
static UnitVector3d orthogonalTo(Vector3d const &v)
orthogonalTo returns an arbitrary unit vector that is orthogonal to v.
T fabs(T... args)
T max(T... args)
T min(T... args)
double tan(Angle const &a)
Definition Angle.h:111
double cos(Angle const &a)
Definition Angle.h:110
constexpr double PI
Definition constants.h:43

◆ Ellipse() [5/5]

lsst::sphgeom::Ellipse::Ellipse ( UnitVector3d const & center,
Angle alpha,
Angle beta,
Angle orientation )

This constructor creates an ellipse with the given center, semi-axis angles, and orientation.

The orientation is defined as the position angle (east of north) of the first axis with respect to the north pole. Note that both alpha and beta must be less than, greater than, or equal to PI/2.

Definition at line 109 of file Ellipse.cc.

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).
136 UnitVector3d b0 = UnitVector3d::orthogonalTo(center);
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}
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 northFrom(Vector3d const &v)
northFrom returns the unit vector orthogonal to v that points "north" from v.
T isfinite(T... args)
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.
T swap(T... args)

Member Function Documentation

◆ clone()

std::unique_ptr< Region > lsst::sphgeom::Ellipse::clone ( ) const
inlineoverridevirtual

clone returns a deep copy of this region.

Implements lsst::sphgeom::Region.

Definition at line 284 of file Ellipse.h.

284 {
285 return std::unique_ptr<Ellipse>(new Ellipse(*this));
286 }

◆ complement()

Ellipse & lsst::sphgeom::Ellipse::complement ( )
inline

complement sets this ellipse to the closure of its complement.

Definition at line 271 of file Ellipse.h.

271 {
272 _S = Matrix3d(-_S(0,0), -_S(0,1), -_S(0,2),
273 _S(1,0), _S(1,1), _S(1,2),
274 -_S(2,0), -_S(2,1), -_S(2,2));
275 _a = -_a;
276 _b = -_b;
277 return *this;
278 }

◆ complemented()

Ellipse lsst::sphgeom::Ellipse::complemented ( ) const
inline

complemented returns the closure of the complement of this ellipse.

Definition at line 281 of file Ellipse.h.

281{ return Ellipse(*this).complement(); }

◆ contains() [1/4]

bool lsst::sphgeom::Region::contains ( double lon,
double lat ) const

contains tests whether the unit vector defined by the given longitude and latitude coordinates (in radians) is inside this region.

Definition at line 111 of file Region.cc.

51 {
52 return contains(UnitVector3d(LonLat::fromRadians(lon, lat)));
53}
bool contains(UnitVector3d const &v) const override
contains tests whether the given unit vector is inside this region.
Definition Ellipse.cc:167
static LonLat fromRadians(double lon, double lat)
Definition LonLat.h:62

◆ contains() [2/4]

bool lsst::sphgeom::Region::contains ( double x,
double y,
double z ) const

contains tests whether the unit vector defined by the given (not necessarily normalized) coordinates is inside this region.

Definition at line 107 of file Region.cc.

47 {
48 return contains(UnitVector3d(x, y, z));
49}
double z
Definition Match.cc:44
int y
Definition SpanSet.cc:48

◆ contains() [3/4]

virtual bool lsst::sphgeom::Region::contains ( UnitVector3d const & ) const
virtual

contains tests whether the given unit vector is inside this region.

Implements lsst::sphgeom::Region.

◆ contains() [4/4]

bool lsst::sphgeom::Ellipse::contains ( UnitVector3d const & ) const
overridevirtual

contains tests whether the given unit vector is inside this region.

Implements lsst::sphgeom::Region.

Definition at line 167 of file Ellipse.cc.

167 {
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}
UnitVector3d getCenter() const
getCenter returns the center of the ellipse as a unit vector.
Definition Ellipse.h:240

◆ decode() [1/2]

static std::unique_ptr< Ellipse > lsst::sphgeom::Ellipse::decode ( std::vector< uint8_t > const & s)
inlinestatic

decode deserializes an Ellipse from a byte string produced by encode.

Definition at line 310 of file Ellipse.h.

310 {
311 return decode(s.data(), s.size());
312 }
static std::unique_ptr< Ellipse > decode(std::vector< uint8_t > const &s)
Definition Ellipse.h:310

◆ decode() [2/2]

std::unique_ptr< Ellipse > lsst::sphgeom::Ellipse::decode ( uint8_t const * buffer,
size_t n )
static

decode deserializes an Ellipse from a byte string produced by encode.

Definition at line 364 of file Ellipse.cc.

364 {
365 if (buffer == nullptr || n != ENCODED_SIZE || buffer[0] != TYPE_CODE) {
366 throw std::runtime_error("Byte-string is not an encoded Ellipse");
367 }
368 std::unique_ptr<Ellipse> ellipse(new Ellipse);
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}
table::Key< int > b
static constexpr uint8_t TYPE_CODE
Definition Ellipse.h:179
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

◆ empty()

static Ellipse lsst::sphgeom::Ellipse::empty ( )
inlinestatic

Definition at line 181 of file Ellipse.h.

181{ return Ellipse(); }

◆ encode()

std::vector< uint8_t > lsst::sphgeom::Ellipse::encode ( ) const
overridevirtual

encode serializes this region into an opaque byte string.

Byte strings emitted by encode can be deserialized with decode.

Implements lsst::sphgeom::Region.

Definition at line 346 of file Ellipse.cc.

346 {
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}
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
T push_back(T... args)
T reserve(T... args)

◆ full()

static Ellipse lsst::sphgeom::Ellipse::full ( )
inlinestatic

Definition at line 183 of file Ellipse.h.

183{ return Ellipse().complement(); }

◆ getAlpha()

Angle lsst::sphgeom::Ellipse::getAlpha ( ) const
inline

getAlpha returns α, the first semi-axis length of the ellipse.

It is negative for empty ellipses, ≥ π for full ellipses and in [0, π) otherwise.

Definition at line 259 of file Ellipse.h.

259{ return Angle(0.5 * PI) + _a; }

◆ getBeta()

Angle lsst::sphgeom::Ellipse::getBeta ( ) const
inline

getBeta returns β, the second semi-axis length of the ellipse.

It is negative for empty ellipses, ≥ π for full ellipses and in [0, π) otherwise.

Definition at line 264 of file Ellipse.h.

264{ return Angle(0.5 * PI) + _b; }

◆ getBoundingBox()

Box lsst::sphgeom::Ellipse::getBoundingBox ( ) const
overridevirtual

getBoundingBox returns a bounding-box for this region.

Implements lsst::sphgeom::Region.

Definition at line 197 of file Ellipse.cc.

197 {
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}
Box getBoundingBox() const override
getBoundingBox returns a bounding-box for this region.
Definition Circle.cc:221
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Definition Ellipse.cc:248

◆ getBoundingBox3d()

Box3d lsst::sphgeom::Ellipse::getBoundingBox3d ( ) const
overridevirtual

getBoundingBox3d returns a 3-dimensional bounding-box for this region.

Implements lsst::sphgeom::Region.

Definition at line 244 of file Ellipse.cc.

244 {
246}
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
Definition Circle.cc:229

◆ getBoundingCircle()

Circle lsst::sphgeom::Ellipse::getBoundingCircle ( ) const
overridevirtual

getBoundingCircle returns a bounding-circle for this region.

Implements lsst::sphgeom::Region.

Definition at line 248 of file Ellipse.cc.

248 {
249 Angle r = std::max(getAlpha(), getBeta()) + 2.0 * Angle(MAX_ASIN_ERROR);
250 return Circle(getCenter(), r);
251}
Angle getAlpha() const
getAlpha returns α, the first semi-axis length of the ellipse.
Definition Ellipse.h:259
Angle getBeta() const
getBeta returns β, the second semi-axis length of the ellipse.
Definition Ellipse.h:264
constexpr double MAX_ASIN_ERROR
Definition constants.h:52

◆ getCenter()

UnitVector3d lsst::sphgeom::Ellipse::getCenter ( ) const
inline

getCenter returns the center of the ellipse as a unit vector.

Definition at line 240 of file Ellipse.h.

240 {
241 return UnitVector3d::fromNormalized(_S(2,0), _S(2,1), _S(2,2));
242 }
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.

◆ getF1()

UnitVector3d lsst::sphgeom::Ellipse::getF1 ( ) const
inline

getF1 returns the first focal point of the ellipse.

Definition at line 245 of file Ellipse.h.

245 {
246 UnitVector3d n = UnitVector3d::fromNormalized(_S(1,0), _S(1,1), _S(1,2));
247 return getCenter().rotatedAround(n, -_gamma);
248 }

◆ getF2()

UnitVector3d lsst::sphgeom::Ellipse::getF2 ( ) const
inline

getF2 returns the second focal point of the ellipse.

Definition at line 251 of file Ellipse.h.

251 {
252 UnitVector3d n = UnitVector3d::fromNormalized(_S(1,0), _S(1,1), _S(1,2));
253 return getCenter().rotatedAround(n, _gamma);
254 }

◆ getGamma()

Angle lsst::sphgeom::Ellipse::getGamma ( ) const
inline

getGamma returns ɣ ∈ [0, π/2], half of the angle between the foci.

The return value is arbitrary for empty and full ellipses.

Definition at line 268 of file Ellipse.h.

268{ return _gamma; }

◆ getTransformMatrix()

Matrix3d const & lsst::sphgeom::Ellipse::getTransformMatrix ( ) const
inline

getTransformMatrix returns the orthogonal matrix that maps vectors to the basis in which the quadratic form corresponding to this ellipse is diagonal.

Definition at line 237 of file Ellipse.h.

237{ return _S; }

◆ isCircle()

bool lsst::sphgeom::Ellipse::isCircle ( ) const
inline

Definition at line 232 of file Ellipse.h.

232{ return _a == _b; }

◆ isEmpty()

bool lsst::sphgeom::Ellipse::isEmpty ( ) const
inline

Definition at line 226 of file Ellipse.h.

226{ return Angle(0.5 * PI) + _a < _gamma; }

◆ isFull()

bool lsst::sphgeom::Ellipse::isFull ( ) const
inline

Definition at line 228 of file Ellipse.h.

228{ return Angle(0.5 * PI) - _a <= _gamma; }

◆ isGreatCircle()

bool lsst::sphgeom::Ellipse::isGreatCircle ( ) const
inline

Definition at line 230 of file Ellipse.h.

230{ return _a.asRadians() == 0.0; }

◆ operator!=()

bool lsst::sphgeom::Ellipse::operator!= ( Ellipse const & e) const
inline

Definition at line 224 of file Ellipse.h.

224{ return !(*this == e); }

◆ operator==()

bool lsst::sphgeom::Ellipse::operator== ( Ellipse const & e) const
inline

Definition at line 220 of file Ellipse.h.

220 {
221 return _S == e._S && _a == e._a && _b == e._b;
222 }

◆ relate() [1/5]

Relationship lsst::sphgeom::Ellipse::relate ( Box const & ) const
overridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

  • Bit S & DISJOINT is set only if A and B do not have any points in common.
  • Bit S & CONTAINS is set only if A contains all points in B.
  • Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 253 of file Ellipse.cc.

253 {
254 return getBoundingCircle().relate(b) & (DISJOINT | WITHIN);
255}
Relationship relate(UnitVector3d const &v) const
Definition Circle.cc:278

◆ relate() [2/5]

Relationship lsst::sphgeom::Ellipse::relate ( Circle const & ) const
overridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

  • Bit S & DISJOINT is set only if A and B do not have any points in common.
  • Bit S & CONTAINS is set only if A contains all points in B.
  • Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 334 of file Ellipse.cc.

334 {
335 return getBoundingCircle().relate(c) & (DISJOINT | WITHIN);
336}

◆ relate() [3/5]

Relationship lsst::sphgeom::Ellipse::relate ( ConvexPolygon const & ) const
overridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

  • Bit S & DISJOINT is set only if A and B do not have any points in common.
  • Bit S & CONTAINS is set only if A contains all points in B.
  • Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 338 of file Ellipse.cc.

338 {
339 return getBoundingCircle().relate(p) & (DISJOINT | WITHIN);
340}

◆ relate() [4/5]

Relationship lsst::sphgeom::Ellipse::relate ( Ellipse const & ) const
overridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

  • Bit S & DISJOINT is set only if A and B do not have any points in common.
  • Bit S & CONTAINS is set only if A contains all points in B.
  • Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 342 of file Ellipse.cc.

342 {
343 return getBoundingCircle().relate(e.getBoundingCircle()) & DISJOINT;
344}

◆ relate() [5/5]

Relationship lsst::sphgeom::Ellipse::relate ( Region const & ) const
inlineoverridevirtual

relate computes the spatial relationships between this region A and another region B. The return value S is a bitset with the following properties:

  • Bit S & DISJOINT is set only if A and B do not have any points in common.
  • Bit S & CONTAINS is set only if A contains all points in B.
  • Bit S & WITHIN is set only if B contains all points in A.

Said another way: if the CONTAINS, WITHIN or DISJOINT bit is set, then the corresponding spatial relationship between the two regions holds conclusively. If it is not set, the relationship may or may not hold.

These semantics allow for conservative relationship computations. In particular, a Region may choose to implement relate by replacing itself and/or the argument with a simplified bounding region.

Implements lsst::sphgeom::Region.

Definition at line 296 of file Ellipse.h.

296 {
297 // Dispatch on the type of r.
298 return invert(r.relate(*this));
299 }
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.

Member Data Documentation

◆ TYPE_CODE

constexpr uint8_t lsst::sphgeom::Ellipse::TYPE_CODE = 'e'
staticconstexpr

Definition at line 179 of file Ellipse.h.


The documentation for this class was generated from the following files: