LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
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. More...
 
 Ellipse (Circle const &c)
 This constructor creates an ellipse corresponding to the given circle. More...
 
 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. More...
 
 Ellipse (UnitVector3d const &f1, UnitVector3d const &f2, Angle alpha)
 This constructor creates an ellipse with the given foci and semi-axis angle. More...
 
 Ellipse (UnitVector3d const &center, Angle alpha, Angle beta, Angle orientation)
 This constructor creates an ellipse with the given center, semi-axis angles, and orientation. More...
 
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. More...
 
UnitVector3d getCenter () const
 getCenter returns the center of the ellipse as a unit vector. More...
 
UnitVector3d getF1 () const
 getF1 returns the first focal point of the ellipse. More...
 
UnitVector3d getF2 () const
 getF2 returns the second focal point of the ellipse. More...
 
Angle getAlpha () const
 getAlpha returns α, the first semi-axis length of the ellipse. More...
 
Angle getBeta () const
 getBeta returns β, the second semi-axis length of the ellipse. More...
 
Angle getGamma () const
 getGamma returns ɣ ∈ [0, π/2], half of the angle between the foci. More...
 
Ellipsecomplement ()
 complement sets this ellipse to the closure of its complement. More...
 
Ellipse complemented () const
 complemented returns the closure of the complement of this ellipse. More...
 
std::unique_ptr< Regionclone () const override
 clone returns a deep copy of this region. More...
 
Box getBoundingBox () const override
 getBoundingBox returns a bounding-box for this region. More...
 
Box3d getBoundingBox3d () const override
 getBoundingBox3d returns a 3-dimensional bounding-box for this region. More...
 
Circle getBoundingCircle () const override
 getBoundingCircle returns a bounding-circle for this region. More...
 
bool contains (UnitVector3d const &v) const override
 contains tests whether the given unit vector is inside this region. More...
 
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. More...
 
virtual bool contains (UnitVector3d const &) const=0
 contains tests whether the given unit vector is inside this region. More...
 
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. More...
 
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. More...
 
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. More...
 
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. More...
 

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 169 of file Ellipse.h.

Constructor & Destructor Documentation

◆ Ellipse() [1/5]

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

This constructor creates an empty ellipse.

Definition at line 178 of file Ellipse.h.

178  :
179  _S(1.0),
180  _a(-2.0),
181  _b(-2.0),
182  _gamma(0.0),
185  {}

◆ Ellipse() [2/5]

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

This constructor creates an ellipse corresponding to the given circle.

Definition at line 188 of file Ellipse.h.

188  {
189  *this = Ellipse(c.getCenter(), c.getCenter(), c.getOpeningAngle());
190  }
Ellipse()
This constructor creates an empty ellipse.
Definition: Ellipse.h:178

◆ 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 194 of file Ellipse.h.

194  {
195  *this = Ellipse(v, v, alpha);
196  }

◆ 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 42 of file Ellipse.cc.

42  :
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).
67  UnitVector3d b0 = UnitVector3d::orthogonalTo(f1);
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 }
T acos(T... args)
double asRadians() const
asRadians returns the value of this angle in units of radians.
Definition: Angle.h:85
static Ellipse full()
Definition: Ellipse.h:175
bool isFull() const
Definition: Ellipse.h:220
static Ellipse empty()
Definition: Ellipse.h:173
bool isEmpty() const
Definition: Ellipse.h:218
static UnitVector3d orthogonalTo(Vector3d const &v)
orthogonalTo returns an arbitrary unit vector that is orthogonal to v.
Definition: UnitVector3d.cc:34
T fabs(T... args)
T max(T... args)
T min(T... args)
lsst::geom::Angle Angle
Definition: misc.h:33
double tan(Angle const &a)
Definition: Angle.h:104
double cos(Angle const &a)
Definition: Angle.h:103
constexpr double PI
Definition: constants.h:36

◆ 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 102 of file Ellipse.cc.

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).
129  UnitVector3d b0 = UnitVector3d::orthogonalTo(center);
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 }
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 northFrom(Vector3d const &v)
northFrom returns the unit vector orthogonal to v that points "north" from v.
Definition: UnitVector3d.cc:51
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.
Definition: orientation.cc:135
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 276 of file Ellipse.h.

276  {
277  return std::unique_ptr<Ellipse>(new Ellipse(*this));
278  }

◆ complement()

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

complement sets this ellipse to the closure of its complement.

Definition at line 263 of file Ellipse.h.

263  {
264  _S = Matrix3d(-_S(0,0), -_S(0,1), -_S(0,2),
265  _S(1,0), _S(1,1), _S(1,2),
266  -_S(2,0), -_S(2,1), -_S(2,2));
267  _a = -_a;
268  _b = -_b;
269  return *this;
270  }

◆ complemented()

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

complemented returns the closure of the complement of this ellipse.

Definition at line 273 of file Ellipse.h.

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

◆ contains() [1/6]

bool lsst::sphgeom::Region::contains

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

Definition at line 104 of file Region.cc.

43  {
44  return contains(UnitVector3d(LonLat::fromRadians(lon, lat)));
45 }
virtual bool contains(UnitVector3d const &) const=0
contains tests whether the given unit vector is inside this region.
static LonLat fromRadians(double lon, double lat)
Definition: LonLat.h:55

◆ contains() [2/6]

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

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

Definition at line 43 of file Region.cc.

43  {
44  return contains(UnitVector3d(LonLat::fromRadians(lon, lat)));
45 }
virtual bool contains(UnitVector3d const &) const =0
contains tests whether the given unit vector is inside this region.

◆ contains() [3/6]

bool lsst::sphgeom::Region::contains

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

Definition at line 100 of file Region.cc.

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

◆ contains() [4/6]

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

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

Definition at line 39 of file Region.cc.

39  {
40  return contains(UnitVector3d(x, y, z));
41 }

◆ contains() [5/6]

virtual bool lsst::sphgeom::Region::contains

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

◆ contains() [6/6]

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 160 of file Ellipse.cc.

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

◆ 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 302 of file Ellipse.h.

302  {
303  return decode(s.data(), s.size());
304  }
static std::unique_ptr< Ellipse > decode(std::vector< uint8_t > const &s)
Definition: Ellipse.h:302
T data(T... args)
T size(T... args)

◆ 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 357 of file Ellipse.cc.

357  {
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 }
table::Key< int > b
table::Key< int > a
static constexpr uint8_t TYPE_CODE
Definition: Ellipse.h:171
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

◆ empty()

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

Definition at line 173 of file Ellipse.h.

173 { 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 339 of file Ellipse.cc.

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

◆ full()

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

Definition at line 175 of file Ellipse.h.

175 { 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 251 of file Ellipse.h.

251 { 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 256 of file Ellipse.h.

256 { 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 190 of file Ellipse.cc.

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

◆ 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 237 of file Ellipse.cc.

237  {
239 }
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
Definition: Circle.cc:219

◆ getBoundingCircle()

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

getBoundingCircle returns a bounding-circle for this region.

Implements lsst::sphgeom::Region.

Definition at line 241 of file Ellipse.cc.

241  {
242  Angle r = std::max(getAlpha(), getBeta()) + 2.0 * Angle(MAX_ASIN_ERROR);
243  return Circle(getCenter(), r);
244 }
Angle getAlpha() const
getAlpha returns α, the first semi-axis length of the ellipse.
Definition: Ellipse.h:251
Angle getBeta() const
getBeta returns β, the second semi-axis length of the ellipse.
Definition: Ellipse.h:256
constexpr double MAX_ASIN_ERROR
Definition: constants.h:45

◆ getCenter()

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

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

Definition at line 232 of file Ellipse.h.

232  {
233  return UnitVector3d::fromNormalized(_S(2,0), _S(2,1), _S(2,2));
234  }
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
Definition: UnitVector3d.h:82

◆ getF1()

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

getF1 returns the first focal point of the ellipse.

Definition at line 237 of file Ellipse.h.

237  {
238  UnitVector3d n = UnitVector3d::fromNormalized(_S(1,0), _S(1,1), _S(1,2));
239  return getCenter().rotatedAround(n, -_gamma);
240  }

◆ getF2()

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

getF2 returns the second focal point of the ellipse.

Definition at line 243 of file Ellipse.h.

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

◆ 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 260 of file Ellipse.h.

260 { 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 229 of file Ellipse.h.

229 { return _S; }

◆ isCircle()

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

Definition at line 224 of file Ellipse.h.

224 { return _a == _b; }

◆ isEmpty()

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

Definition at line 218 of file Ellipse.h.

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

◆ isFull()

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

Definition at line 220 of file Ellipse.h.

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

◆ isGreatCircle()

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

Definition at line 222 of file Ellipse.h.

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

◆ operator!=()

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

Definition at line 216 of file Ellipse.h.

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

◆ operator==()

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

Definition at line 212 of file Ellipse.h.

212  {
213  return _S == e._S && _a == e._a && _b == e._b;
214  }

◆ 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 246 of file Ellipse.cc.

246  {
247  return getBoundingCircle().relate(b) & (DISJOINT | WITHIN);
248 }
Relationship relate(UnitVector3d const &v) const
Definition: Circle.cc:268

◆ 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 327 of file Ellipse.cc.

327  {
328  return getBoundingCircle().relate(c) & (DISJOINT | WITHIN);
329 }

◆ 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 331 of file Ellipse.cc.

331  {
332  return getBoundingCircle().relate(p) & (DISJOINT | WITHIN);
333 }

◆ 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 335 of file Ellipse.cc.

335  {
336  return getBoundingCircle().relate(e.getBoundingCircle()) & DISJOINT;
337 }

◆ 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 288 of file Ellipse.h.

288  {
289  // Dispatch on the type of r.
290  return invert(r.relate(*this));
291  }
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.
Definition: Relationship.h:55

Member Data Documentation

◆ TYPE_CODE

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

Definition at line 171 of file Ellipse.h.


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