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

#include <Ellipse.h>

Inheritance diagram for lsst::sphgeom::Ellipse:

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 override
	`isEmpty` returns true when a region does not contain any points.

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.

Ellipse &	complement ()
	`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< Region >	clone () 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

TriState	overlaps (Region const &other) const override

TriState	overlaps (Box const &) const override

TriState	overlaps (Circle const &) const override

TriState	overlaps (ConvexPolygon const &) const override

TriState	overlaps (Ellipse const &) const override

std::vector< std::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::vector< std::unique_ptr< Region > >	getRegions (Region const &region)
	`getRegions` returns a vector of Region.


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

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


static std::unique_ptr< Region >	decodeBase64 (std::string const &s)

static std::unique_ptr< Region >	decodeBase64 (std::string_view const &s)


static TriState	decodeOverlapsBase64 (std::string const &s)

static TriState	decodeOverlapsBase64 (std::string_view const &s)

Static Public Attributes
static constexpr std::uint8_t	TYPE_CODE = 'e'

Static Protected Member Functions
static TriState	_relationship_to_overlaps (Relationship r)

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.

              :
        _S(1.0),
        _a(-2.0),
        _b(-2.0),
        _gamma(0.0),
        _tana(std::numeric_limits<double>::infinity()),
        _tanb(std::numeric_limits<double>::infinity())
    {}

◆ 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.

                                       {
        *this = Ellipse(c.getCenter(), c.getCenter(), c.getOpeningAngle());
    }

◆ 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.

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

◆ 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.

                                                                              :
    _a(alpha.asRadians() - 0.5 * PI)
{
    if (alpha.isNan()) {
        throw std::invalid_argument("Invalid ellipse semi-axis angle");
    }
    if (f1 == f2) {
        _gamma = Angle(0.0);
    } else if (f1 == -f2) {
        _gamma = Angle(0.5 * PI);
    } else {
        _gamma = 0.5 * NormalizedAngle(f1, f2);
    }
    if (isEmpty()) {
        *this = empty();
        return;
    } else if (isFull()) {
        *this = full();
        return;
    }
    if (_gamma.asRadians() == 0.0) {
        // The foci are identical, so this ellipse is a circle centered at
        // the common focal point. Pick an orthonormal basis containing f1
        // and use it to construct an orthogonal matrix that maps f1 to
        // (0, 0, 1).
        UnitVector3d b0 = UnitVector3d::orthogonalTo(f1);
        UnitVector3d b1 = UnitVector3d(f1.cross(b0));
        _S = Matrix3d(b0.x(), b0.y(), b0.z(),
                      b1.x(), b1.y(), b1.z(),
                      f1.x(), f1.y(), f1.z());
        _b = _a;
        _tana = std::fabs(tan(_a));
        _tanb = _tana;
        return;
    }
    // _gamma != 0 implies that f1 - f2 != 0. Also, if f1 = -f2 then
    // _gamma = PI/2, and the ellipse must either empty or full. So
    // at this stage f1 + f2 != 0.
    Vector3d b0 = f1 - f2;
    Vector3d b2 = f1 + f2;
    Vector3d b1 = b0.cross(b2);
    b0.normalize();
    b1.normalize();
    b2.normalize();
    _S = Matrix3d(b0.x(), b0.y(), b0.z(),
                  b1.x(), b1.y(), b1.z(),
                  b2.x(), b2.y(), b2.z());
    // Compute _b.
    double r = std::min(1.0, std::max(-1.0, cos(alpha) / cos(_gamma)));
    _b = Angle(std::acos(r) - 0.5 * PI);
    if (_a.asRadians() <= 0.0 && _b > _a) {
        _b = _a;
    } else if (_a.asRadians() > 0.0 && _b < _a) {
        _b = _a;
    }
    _tana = std::fabs(tan(_a));
    _tanb = std::fabs(tan(_b));
    return;
}

◆ 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.

{
    if (!std::isfinite(orientation.asRadians())) {
        throw std::invalid_argument("Invalid ellipse orientation");
    }
    if (alpha.isNan() ||
        beta.isNan() ||
        (alpha.asRadians() <  0.5 * PI && beta.asRadians() >= 0.5 * PI) ||
        (alpha.asRadians() >  0.5 * PI && beta.asRadians() <= 0.5 * PI) ||
        (alpha.asRadians() == 0.5 * PI && beta.asRadians() != 0.5 * PI)) {
        throw std::invalid_argument("Invalid ellipse semi-axis angle(s)");
    }
    if (alpha.asRadians() < 0.0 || beta.asRadians() < 0.0) {
        *this = empty();
        return;
    } else if (alpha.asRadians() > PI || beta.asRadians() > PI ||
               (alpha.asRadians() == PI && beta.asRadians() == PI)) {
        *this = full();
        return;
    }
    if (alpha == beta) {
        // The ellipse is a circle. Pick an orthonormal basis containing the
        // center and use it to construct some orthogonal matrix that maps the
        // center to (0, 0, 1).
        UnitVector3d b0 = UnitVector3d::orthogonalTo(center);
        UnitVector3d b1 = UnitVector3d(center.cross(b0));
        _S = Matrix3d(b0.x(), b0.y(), b0.z(),
                      b1.x(), b1.y(), b1.z(),
                      center.x(), center.y(), center.z());
        _a = alpha - Angle(0.5 * PI);
        _b = _a;
        _gamma = Angle(0.0);
        _tana = std::fabs(tan(_a));
        _tanb = _tana;
        return;
    }
    if ((alpha.asRadians() < 0.5 * PI && alpha < beta) ||
        (alpha.asRadians() > 0.5 * PI && alpha > beta)) {
        std::swap(alpha, beta);
        orientation = orientation + Angle(0.5 * PI);
    }
    UnitVector3d b0 =
        UnitVector3d::northFrom(center).rotatedAround(center, -orientation);
    UnitVector3d b1 = UnitVector3d(b0.cross(center));
    _S = Matrix3d(b0.x(), b0.y(), b0.z(),
                  b1.x(), b1.y(), b1.z(),
                  center.x(), center.y(), center.z());
    _a = alpha - Angle(0.5 * PI);
    _b = beta - Angle(0.5 * PI);
    double d = std::min(1.0, std::max(-1.0, cos(alpha) / cos(beta)));
    _gamma = Angle(std::acos(d));
    _tana = std::fabs(tan(_a));
    _tanb = std::fabs(tan(_b));
}

Member Function Documentation

◆ _relationship_to_overlaps()

static TriState lsst::sphgeom::Region::_relationship_to_overlaps ( Relationship r )

inlinestaticprotectedinherited

Definition at line 206 of file Region.h.

                                                              {
        // `relate` returns exact relation when specific bit is set, if it is
        // not then relation may be true or not.
        if ((r & DISJOINT) == DISJOINT) {
            return TriState(false);
        }
        if ((r & (WITHIN | CONTAINS)).any()) {
            return TriState(true);
        }
        return TriState();
    }

◆ 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.

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

◆ 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.

                           {
        _S = Matrix3d(-_S(0,0), -_S(0,1), -_S(0,2),
                       _S(1,0),  _S(1,1),  _S(1,2),
                      -_S(2,0), -_S(2,1), -_S(2,2));
        _a = -_a;
        _b = -_b;
        return *this;
    }

◆ 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 117 of file Region.cc.

                                                  {
    return contains(UnitVector3d(LonLat::fromRadians(lon, lat)));
}

◆ 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 113 of file Region.cc.

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

◆ 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.

                                                   {
    UnitVector3d const c = getCenter();
    double vdotc = v.dot(c);
    Vector3d u;
    double scz;
    // To maintain high accuracy for very small and very large ellipses,
    // decompose v as v = u ± c near ±c. Then S v = S u ± S c, and
    // S c = (0, 0, 1).
    if (vdotc > 0.5) {
        u = v - c;
        scz = 1.0;
    } else if (vdotc < -0.5) {
        u = v + c;
        scz = -1.0;
    } else {
        u = v;
        scz = 0.0;
    }
    u = _S * u;
    double x = u.x() * _tana;
    double y = u.y() * _tanb;
    double z = u.z() + scz;
    double d = (x * x + y * y) - z * z;
    if (_a.asRadians() > 0.0) {
        return z >= 0.0 || d >= 0.0;
    } else {
        return z >= 0.0 && d <= 0.0;
    }
}

◆ decode() [1/2]

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

static

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

Definition at line 388 of file Ellipse.cc.

                                                                          {
    if (buffer == nullptr || n != ENCODED_SIZE || buffer[0] != TYPE_CODE) {
        throw std::runtime_error("Byte-string is not an encoded Ellipse");
    }
    std::unique_ptr<Ellipse> ellipse(new Ellipse);
    ++buffer;
    double m00 = decodeDouble(buffer); buffer += 8;
    double m01 = decodeDouble(buffer); buffer += 8;
    double m02 = decodeDouble(buffer); buffer += 8;
    double m10 = decodeDouble(buffer); buffer += 8;
    double m11 = decodeDouble(buffer); buffer += 8;
    double m12 = decodeDouble(buffer); buffer += 8;
    double m20 = decodeDouble(buffer); buffer += 8;
    double m21 = decodeDouble(buffer); buffer += 8;
    double m22 = decodeDouble(buffer); buffer += 8;
    ellipse->_S = Matrix3d(m00, m01, m02,
                           m10, m11, m12,
                           m20, m21, m22);
    double a = decodeDouble(buffer); buffer += 8;
    double b = decodeDouble(buffer); buffer += 8;
    double gamma = decodeDouble(buffer); buffer += 8;
    ellipse->_a = Angle(a);
    ellipse->_b = Angle(b);
    ellipse->_gamma = Angle(gamma);
    double tana = decodeDouble(buffer); buffer += 8;
    double tanb = decodeDouble(buffer); buffer += 8;
    ellipse->_tana = tana;
    ellipse->_tanb = tanb;
    return ellipse;
}

◆ decode() [2/2]

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

inlinestatic

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

Definition at line 318 of file Ellipse.h.

                                                                          {
        return decode(s.data(), s.size());
    }

◆ decodeBase64() [1/2]

static std::unique_ptr< Region > lsst::sphgeom::Region::decodeBase64 ( std::string const & s )

inlinestaticinherited

decodeBase64 deserializes a Region from an ASCII string produced by encode and then base64-encoding that result.

This method also interprets ':' as a delimiter for the elements of a UnionRegion, to support cases where a union of region is constructed server-side in a database as a concatenation with that delimiter.

Definition at line 176 of file Region.h.

                                                                   {
        return decodeBase64(s);
    }

◆ decodeBase64() [2/2]

std::unique_ptr< Region > lsst::sphgeom::Region::decodeBase64 ( std::string_view const & s )

staticinherited

decodeBase64 deserializes a Region from an ASCII string produced by encode and then base64-encoding that result.

This method also interprets ':' as a delimiter for the elements of a UnionRegion, to support cases where a union of region is constructed server-side in a database as a concatenation with that delimiter.

Definition at line 93 of file Region.cc.

                                                                   {
    if (s.empty()) {
        return std::unique_ptr<UnionRegion>(new UnionRegion({}));
    }
    auto region_begin = s.begin();
    auto region_end = std::find(s.begin(), s.end(), ':');
    if (region_end != s.end()) {
        std::vector<std::unique_ptr<Region>> union_args;
        while (region_end != s.end()) {
            auto bytes = base64::decode_into<std::vector<std::uint8_t>>(region_begin, region_end);
            union_args.push_back(decode(bytes));
            region_begin = region_end;
            ++region_begin;
            region_end = std::find(region_begin, s.end(), ':');
        }
        auto bytes = base64::decode_into<std::vector<std::uint8_t>>(region_begin, region_end);
        union_args.push_back(decode(bytes));
        return std::unique_ptr<UnionRegion>(new UnionRegion(std::move(union_args)));
    } else {
        auto bytes = base64::decode_into<std::vector<std::uint8_t>>(region_begin, region_end);
        return decode(bytes);
    }
}

◆ decodeOverlapsBase64() [1/2]

static TriState lsst::sphgeom::Region::decodeOverlapsBase64 ( std::string const & s )

inlinestaticinherited

decodeOverlapsBase64 evaluates an encoded overlap expression.

A single overlap expression is formed by concatenating a pair of base64-encoded regions (Region::encode then base64 encoding) with '&' as the delimiter. Multiple such pairwise overlap expressions can then be concatenated with '|' as the delimiter to form the logical OR.

Definition at line 190 of file Region.h.

                                                              {
        return decodeOverlapsBase64(s);
    }

◆ decodeOverlapsBase64() [2/2]

TriState lsst::sphgeom::Region::decodeOverlapsBase64 ( std::string_view const & s )

staticinherited

decodeOverlapsBase64 evaluates an encoded overlap expression.

A single overlap expression is formed by concatenating a pair of base64-encoded regions (Region::encode then base64 encoding) with '&' as the delimiter. Multiple such pairwise overlap expressions can then be concatenated with '|' as the delimiter to form the logical OR.

Definition at line 117 of file Region.cc.

                                                              {
    TriState result(false);
    if (s.empty()) {
        // False makes the most sense as the limit of a logical OR of zero
        // terms (e.g. `any([])` in Python).
        return result;
    }
    auto begin = s.begin();
    while (result != true) {  // if result is known to be true, we're done.
        auto mid = std::find(begin, s.end(), '&');
        if (mid == s.end()) {
            throw std::runtime_error("No '&' found in encoded overlap expression term.");
        }
        auto a = Region::decode(base64::decode_into<std::vector<std::uint8_t>>(begin, mid));
        ++mid;
        auto end = std::find(mid, s.end(), '|');
        auto b = Region::decode(base64::decode_into<std::vector<std::uint8_t>>(mid, end));
        result = result | a->overlaps(*b);
        if (end == s.end()) {
            break;
        } else {
            begin = end;
            ++begin;
        }
    }
    return result;
}

◆ empty()

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

inlinestatic

Definition at line 181 of file Ellipse.h.

181{ return Ellipse(); }

◆ encode()

std::vector< std::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 370 of file Ellipse.cc.

                                            {
    std::vector<std::uint8_t> buffer;
    std::uint8_t tc = TYPE_CODE;
    buffer.reserve(ENCODED_SIZE);
    buffer.push_back(tc);
    for (int r = 0; r < 3; ++r) {
        for (int c = 0; c < 3; ++c) {
            encodeDouble(_S(r, c), buffer);
        }
    }
    encodeDouble(_a.asRadians(), buffer);
    encodeDouble(_b.asRadians(), buffer);
    encodeDouble(_gamma.asRadians(), buffer);
    encodeDouble(_tana, buffer);
    encodeDouble(_tanb, buffer);
    return buffer;
}

◆ 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.

                                  {
    // For now, simply return the bounding box of the ellipse bounding circle.
    //
    // Improving on this seems difficult, mainly because error bounds must be
    // computed to guarantee that the resulting box tightly bounds the ellipse.
    // In case this ends up being important and someone wants to go down
    // this route in the future, here are some thoughts on how to proceed.
    //
    // First of all, if the ellipse contains a pole, its bounding box must
    // contain all longitudes. Otherwise, consider the plane spanned by the
    // basis vectors u = (0, 0, 1) and v = (cos θ, sin θ, 0). To obtain
    // longitude bounds for the ellipse, we must find values of θ for which
    // this plane is tangent to the elliptical cone that defines the spherical
    // ellipse boundary. This is the case when the plane intersects the cone
    // in a single line.
    //
    // Let μ v + λ u be the direction of the line of intersection, where
    // μ, λ ∈ ℝ. We know that μ ≠ 0 because u is not on the ellipse boundary,
    // so fix μ = 1. Let Q be the symmetric matrix representation of the
    // ellipse. Then:
    //
    //    (v + λ u)ᵀ Q (v + λ u) = 0
    //
    // Expanding gives:
    //
    //    (vᵀ Q v) + λ (uᵀ Q v) + λ (vᵀ Q u) + λ² (uᵀ Q u) = 0
    //
    // By the symmetry of Q:
    //
    //    λ² (uᵀ Q u) + 2λ (uᵀ Q v) + (vᵀ Q v) = 0
    //
    // This is a quadratic equation which has one solution exactly when its
    // discriminant is 0, i.e. when:
    //
    //    (uᵀ Q v)² - (uᵀ Q u) (vᵀ Q v) = 0
    //
    // Substituting for u and v and simplifying yields:
    //
    //    (Q₁₂² - Q₁₁Q₂₂)tan²θ + 2(Q₀₂Q₁₂ - Q₀₁Q₂₂)tan θ + (Q₀₂² - Q₀₀Q₂₂) = 0
    //
    // Finding the latitude bounds can be done by parameterizing the ellipse
    // boundary and then solving for the zeros of the derivative of z with
    // respect to the parameter. This looks to be more involved than the
    // longitude bound calculation, and I haven't worked through the details.
    return getBoundingCircle().getBoundingBox();
}

◆ 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.

                                      {
    return getBoundingCircle().getBoundingBox3d();
}

◆ 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.

                                        {
    Angle r = std::max(getAlpha(), getBeta()) + 2.0 * Angle(MAX_ASIN_ERROR);
    return Circle(getCenter(), r);
}

◆ 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.

                                   {
        return UnitVector3d::fromNormalized(_S(2,0), _S(2,1), _S(2,2));
    }

◆ 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.

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

◆ 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.

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

◆ 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; }

◆ getRegions()

std::vector< std::unique_ptr< Region > > lsst::sphgeom::Region::getRegions ( Region const & region )

staticinherited

getRegions returns a vector of Region.

Definition at line 145 of file Region.cc.

                                                                        {
    std::vector<std::unique_ptr<Region>> result;
    if (auto union_region = dynamic_cast<UnionRegion const *>(&region)) {
        for(unsigned i = 0; i < union_region->nOperands(); ++i) {
            result.emplace_back(union_region->getOperand(i).clone());
        }
    } else if(auto intersection_region = dynamic_cast<IntersectionRegion const *>(&region)) {
        for(unsigned i = 0; i < intersection_region->nOperands(); ++i) {
            result.emplace_back(intersection_region->getOperand(i).clone());
        }
    } else {
        result.emplace_back(region.clone());
    }
    return result;
}

◆ 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

inlineoverridevirtual

isEmpty returns true when a region does not contain any points.

Implements lsst::sphgeom::Region.

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.

                                             {
        return _S == e._S && _a == e._a && _b == e._b;
    }

◆ overlaps() [1/5]

TriState lsst::sphgeom::Ellipse::overlaps ( Box const & ) const

overridevirtual

overlaps tests whether two regions overlap. This method returns a TriState object, when the value is true it means that regions definitely overlap, false means they are definitely disjont, and unknown state means that they may or may not overlap.

Implements lsst::sphgeom::Region.

Definition at line 346 of file Ellipse.cc.

                                             {
    // `relate` uses bounding circle approximation which is not exact and can
    // only reliably return DISJOINT and WITHIN.
    return _relationship_to_overlaps(relate(b));
}

◆ overlaps() [2/5]

TriState lsst::sphgeom::Ellipse::overlaps ( Circle const & ) const

overridevirtual

overlaps tests whether two regions overlap. This method returns a TriState object, when the value is true it means that regions definitely overlap, false means they are definitely disjont, and unknown state means that they may or may not overlap.

Implements lsst::sphgeom::Region.

Definition at line 352 of file Ellipse.cc.

                                                {
    // `relate` uses bounding circle approximation which is not exact and can
    // only reliably return DISJOINT and WITHIN.
    return _relationship_to_overlaps(relate(c));
}

◆ overlaps() [3/5]

TriState lsst::sphgeom::Ellipse::overlaps ( ConvexPolygon const & ) const

overridevirtual

overlaps tests whether two regions overlap. This method returns a TriState object, when the value is true it means that regions definitely overlap, false means they are definitely disjont, and unknown state means that they may or may not overlap.

Implements lsst::sphgeom::Region.

Definition at line 358 of file Ellipse.cc.

                                                       {
    // `relate` uses bounding circle approximation which is not exact and can
    // only reliably return DISJOINT and WITHIN.
    return _relationship_to_overlaps(relate(p));
}

◆ overlaps() [4/5]

TriState lsst::sphgeom::Ellipse::overlaps ( Ellipse const & ) const

overridevirtual

overlaps tests whether two regions overlap. This method returns a TriState object, when the value is true it means that regions definitely overlap, false means they are definitely disjont, and unknown state means that they may or may not overlap.

Implements lsst::sphgeom::Region.

Definition at line 364 of file Ellipse.cc.

                                                 {
    // `relate` uses bounding circle approximation which is not exact and can
    // only reliably return DISJOINT and WITHIN.
    return _relationship_to_overlaps(relate(e));
}

◆ overlaps() [5/5]

TriState lsst::sphgeom::Ellipse::overlaps ( Region const & other ) const

inlineoverridevirtual

overlaps tests whether two regions overlap. This method returns a TriState object, when the value is true it means that regions definitely overlap, false means they are definitely disjont, and unknown state means that they may or may not overlap.

Implements lsst::sphgeom::Region.

Definition at line 306 of file Ellipse.h.

                                                          {
        return other.overlaps(*this);
    }

◆ 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.

                                                {
    return getBoundingCircle().relate(b) & (DISJOINT | WITHIN);
}

◆ 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.

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

◆ 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.

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

◆ 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.

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

◆ 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.

                                                         {
        // Dispatch on the type of r.
        return invert(r.relate(*this));
    }

Member Data Documentation

◆ TYPE_CODE

constexpr std::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:

/j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/sphgeom/g3e8969e208+a8ce1bb630/include/lsst/sphgeom/Ellipse.h
/j/snowflake/release/lsstsw/stack/lsst-scipipe-9.0.0/Linux64/sphgeom/g3e8969e208+a8ce1bb630/src/Ellipse.cc

Public Member Functions

Static Public Member Functions

Static Public Attributes

Static Protected Member Functions

Detailed Description

Mathematical Definition

Implementation

Constructor & Destructor Documentation

◆ Ellipse() [1/5]

◆ Ellipse() [2/5]

◆ Ellipse() [3/5]

◆ Ellipse() [4/5]

◆ Ellipse() [5/5]

Member Function Documentation

◆ _relationship_to_overlaps()

◆ clone()

◆ complement()

◆ complemented()

◆ contains() [1/4]

◆ contains() [2/4]

◆ contains() [3/4]

◆ contains() [4/4]

◆ decode() [1/2]

◆ decode() [2/2]

◆ decodeBase64() [1/2]

◆ decodeBase64() [2/2]

◆ decodeOverlapsBase64() [1/2]

◆ decodeOverlapsBase64() [2/2]

◆ empty()

◆ encode()

◆ full()

◆ getAlpha()

◆ getBeta()

◆ getBoundingBox()

◆ getBoundingBox3d()

◆ getBoundingCircle()

◆ getCenter()

◆ getF1()

◆ getF2()

◆ getGamma()

◆ getRegions()

◆ getTransformMatrix()

◆ isCircle()

◆ isEmpty()

◆ isFull()

◆ isGreatCircle()

◆ operator!=()

◆ operator==()

◆ overlaps() [1/5]

◆ overlaps() [2/5]

◆ overlaps() [3/5]

◆ overlaps() [4/5]

◆ overlaps() [5/5]

◆ relate() [1/5]

◆ relate() [2/5]

◆ relate() [3/5]

◆ relate() [4/5]

◆ relate() [5/5]

Member Data Documentation

◆ TYPE_CODE