lsst::sphgeom Namespace Reference

Namespaces
	detail
	python
	version

Classes
class	Angle
	Angle represents an angle in radians. More...
class	AngleInterval
	AngleInterval represents closed intervals of arbitrary angles. More...
class	BigInteger
	BigInteger is an arbitrary precision signed integer class. More...
class	Box
	Box represents a rectangle in spherical coordinate space that contains its boundary. More...
class	Box3d
	Box3d represents a box in ℝ³. More...
class	Chunker
	Chunker subdivides the unit sphere into longitude-latitude boxes. More...
class	Circle
	Circle is a circular region on the unit sphere that contains its boundary. More...
class	ConvexPolygon
	ConvexPolygon is a closed convex polygon on the unit sphere. More...
class	Ellipse
	Ellipse is an elliptical region on the sphere. More...
class	HtmPixelization
	HtmPixelization provides HTM indexing of points and regions. More...
class	Interval
	Interval represents a closed interval of the real numbers by its upper and lower bounds. More...
class	Interval1d
	Interval1d represents closed intervals of ℝ. More...
class	LonLat
	LonLat represents a spherical coordinate (longitude/latitude angle) pair. More...
class	Matrix3d
	A 3x3 matrix with real entries stored in double precision. More...
class	Mq3cPixelization
	Mq3cPixelization provides modified Q3C indexing of points and regions. More...
class	NormalizedAngle
	NormalizedAngle is an angle that lies in the range [0, 2π), with one exception - a NormalizedAngle can be NaN. More...
class	NormalizedAngleInterval
	NormalizedAngleInterval represents closed intervals of normalized angles, i.e. More...
class	Pixelization
	A Pixelization (or partitioning) of the sphere is a mapping between points on the sphere and a set of pixels (a.k.a. More...
class	Q3cPixelization
	Q3cPixelization provides Q3C indexing of points and regions. More...
class	RangeSet
	A RangeSet is a set of unsigned 64 bit integers. More...
class	Region
	Region is a minimal interface for 2-dimensional regions on the unit sphere. More...
struct	SubChunks
	SubChunks represents a set of sub-chunks of a particular chunk. More...
class	UnitVector3d
	UnitVector3d is a unit vector in ℝ³ with components stored in double precision. More...
class	Vector3d
	Vector3d is a vector in ℝ³ with components stored in double precision. More...

Typedefs
using	Relationship = std::bitset< 3 >
	Relationship describes how two sets are related. More...

Functions
Angle	operator* (double a, Angle const &b)
std::ostream &	operator<< (std::ostream &, Angle const &)
double	sin (Angle const &a)
double	cos (Angle const &a)
double	tan (Angle const &a)
Angle	abs (Angle const &a)
std::ostream &	operator<< (std::ostream &, AngleInterval const &)
std::ostream &	operator<< (std::ostream &, Box const &)
std::ostream &	operator<< (std::ostream &, Box3d const &)
std::ostream &	operator<< (std::ostream &, Circle const &)
void	encodeDouble (double item, std::vector< uint8_t > &buffer)
	encode appends an IEEE double in little-endian byte order to the end of buffer. More...
double	decodeDouble (uint8_t const *buffer)
	decode extracts an IEEE double from the 8 byte little-endian byte sequence in buffer. More...
std::ostream &	operator<< (std::ostream &, ConvexPolygon const &)
uint64_t	mortonIndex (uint32_t x, uint32_t y)
	mortonIndex interleaves the bits of x and y. More...
std::tuple< uint32_t, uint32_t >	mortonIndexInverse (uint64_t z)
	mortonIndexInverse separates the even and odd bits of z. More...
uint64_t	mortonToHilbert (uint64_t z, int m)
	mortonToHilbert converts the 2m-bit Morton index z to the corresponding Hilbert index. More...
uint64_t	hilbertToMorton (uint64_t h, int m)
	hilbertToMorton converts the 2m-bit Hilbert index h to the corresponding Morton index. More...
uint64_t	hilbertIndex (uint32_t x, uint32_t y, int m)
	hilbertIndex returns the index of (x, y) in a 2-D Hilbert curve. More...
std::tuple< uint32_t, uint32_t >	hilbertIndexInverse (uint64_t h, int m)
	hilbertIndexInverse returns the point (x, y) with Hilbert index h, where x and y are m bit integers. More...
std::ostream &	operator<< (std::ostream &, Ellipse const &)
std::ostream &	operator<< (std::ostream &, Interval1d const &)
std::ostream &	operator<< (std::ostream &, LonLat const &)
std::ostream &	operator<< (std::ostream &, Matrix3d const &)
NormalizedAngle const &	abs (NormalizedAngle const &a)
std::ostream &	operator<< (std::ostream &, NormalizedAngleInterval const &)
int	orientationExact (Vector3d const &a, Vector3d const &b, Vector3d const &c)
	orientationExact computes and returns the orientations of 3 vectors a, b and c, which need not be normalized but are assumed to have finite components. More...
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. More...
int	orientationX (UnitVector3d const &b, UnitVector3d const &c)
	orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c). More...
int	orientationY (UnitVector3d const &b, UnitVector3d const &c)
	orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c). More...
int	orientationZ (UnitVector3d const &b, UnitVector3d const &c)
	orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c). More...
void	swap (RangeSet &a, RangeSet &b)
std::ostream &	operator<< (std::ostream &, RangeSet const &)
Relationship	invert (Relationship r)
	Given the relationship between two sets A and B (i.e. More...
std::ostream &	operator<< (std::ostream &, UnitVector3d const &)
double	getMinSquaredChordLength (Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
	Let p be the unit vector closest to v that lies on the plane with normal n in the direction of the cross product of a and b. More...
double	getMaxSquaredChordLength (Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
	Let p be the unit vector furthest from v that lies on the plane with normal n in the direction of the cross product of a and b. More...
Angle	getMinAngleToCircle (Angle x, Angle c)
	getMinAngleToCircle returns the minimum angular separation between a point at latitude x and the points on the circle of constant latitude c. More...
Angle	getMaxAngleToCircle (Angle x, Angle c)
	getMaxAngleToCircle returns the maximum angular separation between a point at latitude x and the points on the circle of constant latitude c. More...
Vector3d	getWeightedCentroid (UnitVector3d const &v0, UnitVector3d const &v1, UnitVector3d const &v2)
	getWeightedCentroid returns the center of mass of the given spherical triangle (assuming a uniform mass distribution over the triangle surface), weighted by the triangle area. More...
Vector3d	operator* (double s, Vector3d const &v)
std::ostream &	operator<< (std::ostream &, Vector3d const &)

Variables
constexpr double	PI = 3.1415926535897932384626433832795
constexpr double	ONE_OVER_PI = 0.318309886183790671537767526745
constexpr double	RAD_PER_DEG = 0.0174532925199432957692369076849
constexpr double	DEG_PER_RAD = 57.2957795130823208767981548141
constexpr double	MAX_ASIN_ERROR = 1.5e-8
constexpr double	MAX_SQUARED_CHORD_LENGTH_ERROR = 2.5e-15
constexpr double	EPSILON = 1.1102230246251565e-16

Typedef Documentation

Relationship

using lsst::sphgeom::Relationship = typedef std::bitset<3>

Relationship describes how two sets are related.

Definition at line 35 of file Relationship.h.

Function Documentation

abs() [1/2]

Angle lsst::sphgeom::abs ( Angle const &  a )

inline

Definition at line 106 of file Angle.h.

{ return Angle(std::fabs(a.asRadians())); }

abs() [2/2]

NormalizedAngle const& lsst::sphgeom::abs ( NormalizedAngle const &  a )

inline

Definition at line 150 of file NormalizedAngle.h.

{ return a; }

cos()

double lsst::sphgeom::cos ( Angle const &  a )

inline

Definition at line 103 of file Angle.h.

{ return std::cos(a.asRadians()); }

decodeDouble()

double lsst::sphgeom::decodeDouble ( uint8_t const *  buffer )

inline

decode extracts an IEEE double from the 8 byte little-endian byte sequence in buffer.

Definition at line 59 of file codec.h.

                                                   {
#if defined(__x86_64__)
    // x86-64 is little endian and supports unaligned loads.
    return *reinterpret_cast<double const *>(buffer);
#else
    union { uint64_t u; double d };
    u = static_cast<uint64_t>(buffer[0]) +
        (static_cast<uint64_t>(buffer[1]) << 8) +
        (static_cast<uint64_t>(buffer[2]) << 16) +
        (static_cast<uint64_t>(buffer[3]) << 24) +
        (static_cast<uint64_t>(buffer[4]) << 32) +
        (static_cast<uint64_t>(buffer[5]) << 40) +
        (static_cast<uint64_t>(buffer[6]) << 48) +
        (static_cast<uint64_t>(buffer[7]) << 56);
    return d;
#endif
}

encodeDouble()

void lsst::sphgeom::encodeDouble ( double  item, std::vector< uint8_t > &  buffer  )

inline

encode appends an IEEE double in little-endian byte order to the end of buffer.

Definition at line 38 of file codec.h.

                                                                   {
#if defined(__x86_64__)
    // x86-64 is little endian.
    auto ptr = reinterpret_cast<uint8_t const *>(&item);
    buffer.insert(buffer.end(), ptr, ptr + 8);
#else
    union { uint64_t u; double d };
    d = item;
    buffer.push_back(static_cast<uint8_t>(value));
    buffer.push_back(static_cast<uint8_t>(value >> 8));
    buffer.push_back(static_cast<uint8_t>(value >> 16));
    buffer.push_back(static_cast<uint8_t>(value >> 24));
    buffer.push_back(static_cast<uint8_t>(value >> 32));
    buffer.push_back(static_cast<uint8_t>(value >> 40));
    buffer.push_back(static_cast<uint8_t>(value >> 48));
    buffer.push_back(static_cast<uint8_t>(value >> 56));
#endif
}

getMaxAngleToCircle()

Angle lsst::sphgeom::getMaxAngleToCircle ( Angle  x, Angle  c  )

inline

getMaxAngleToCircle returns the maximum angular separation between a point at latitude x and the points on the circle of constant latitude c.

Definition at line 68 of file utils.h.

                                                   {
    Angle a = getMinAngleToCircle(x, c);
    if (abs(x) <= abs(c)) {
        return a + Angle(PI) - 2.0 * abs(c);
    }
    if (a < abs(x)) {
        return Angle(PI) - 2.0 * abs(c) - a;
    }
    return Angle(PI) + 2.0 * abs(c) - a;
}

getMaxSquaredChordLength()

double lsst::sphgeom::getMaxSquaredChordLength	(	Vector3d const &	v,
		Vector3d const &	a,
		Vector3d const &	b,
		Vector3d const &	n
	)

Let p be the unit vector furthest from v that lies on the plane with normal n in the direction of the cross product of a and b.

If p is in the interior of the great circle segment from a to b, then this helper function returns the squared chord length between p and v. Otherwise it returns 0 - the minimum squared chord length between any pair of points on the sphere.

Definition at line 58 of file utils.cc.

{
    Vector3d vxn = v.cross(n);
    if (vxn.dot(a) < 0.0 && vxn.dot(b) > 0.0) {
        // v is in the lune defined by the half great circle passing through
        // n and -a and the half great circle passing through n and -b, so p
        // is in the interior of the great circle segment from a to b. The
        // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
        // and ‖v‖ ‖n‖ cos θ = -‖v × n‖. The desired squared chord length is
        // 4 sin²(θ/2).
        double s = std::fabs(v.dot(n));
        double c = - vxn.getNorm();
        double d = std::sin(0.5 * std::atan2(s, c));
        return 4.0 * d * d;
    }
    return 0.0;
}

getMinAngleToCircle()

Angle lsst::sphgeom::getMinAngleToCircle ( Angle  x, Angle  c  )

inline

getMinAngleToCircle returns the minimum angular separation between a point at latitude x and the points on the circle of constant latitude c.

Definition at line 62 of file utils.h.

                                                   {
    return abs(x - c);
}

getMinSquaredChordLength()

double lsst::sphgeom::getMinSquaredChordLength	(	Vector3d const &	v,
		Vector3d const &	a,
		Vector3d const &	b,
		Vector3d const &	n
	)

Let p be the unit vector closest to v that lies on the plane with normal n in the direction of the cross product of a and b.

If p is in the interior of the great circle segment from a to b, then this function returns the squared chord length between p and v. Otherwise it returns 4 - the maximum squared chord length between any pair of points on the unit sphere.

Definition at line 36 of file utils.cc.

{
    Vector3d vxn = v.cross(n);
    if (vxn.dot(a) > 0.0 && vxn.dot(b) < 0.0) {
        // v is in the lune defined by the half great circle passing through
        // n and a and the half great circle passing through n and b, so p
        // is in the interior of the great circle segment from a to b. The
        // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
        // and ‖v‖ ‖n‖ cos θ = ‖v × n‖. The desired squared chord length is
        // 4 sin²(θ/2).
        double s = std::fabs(v.dot(n));
        double c = vxn.getNorm();
        double theta = (c == 0.0) ? 0.5 * PI : std::atan(s / c);
        double d = std::sin(0.5 * theta);
        return 4.0 * d * d;
    }
    return 4.0;
}

getWeightedCentroid()

Vector3d lsst::sphgeom::getWeightedCentroid	(	UnitVector3d const &	v0,
		UnitVector3d const &	v1,
		UnitVector3d const &	v2
	)

getWeightedCentroid returns the center of mass of the given spherical triangle (assuming a uniform mass distribution over the triangle surface), weighted by the triangle area.

Definition at line 79 of file utils.cc.

{
    // For the details, see:
    //
    // The centroid and inertia tensor for a spherical triangle
    // John E. Brock
    // 1974, Naval Postgraduate School, Monterey Calif.
    //
    // https://openlibrary.org/books/OL25493734M/The_centroid_and_inertia_tensor_for_a_spherical_triangle
 
    Vector3d x01 = v0.robustCross(v1); // twice the cross product of v0 and v1
    Vector3d x12 = v1.robustCross(v2);
    Vector3d x20 = v2.robustCross(v0);
    double s01 = 0.5 * x01.normalize(); // sine of the angle between v0 and v1
    double s12 = 0.5 * x12.normalize();
    double s20 = 0.5 * x20.normalize();
    double c01 = v0.dot(v1); // cosine of the angle between v0 and v1
    double c12 = v1.dot(v2);
    double c20 = v2.dot(v0);
    double a0 = (s12 == 0.0 && c12 == 0.0) ? 0.0 : std::atan2(s12, c12);
    double a1 = (s20 == 0.0 && c20 == 0.0) ? 0.0 : std::atan2(s20, c20);
    double a2 = (s01 == 0.0 && c01 == 0.0) ? 0.0 : std::atan2(s01, c01);
    return 0.5 * (x01 * a2 + x12 * a0 + x20 * a1);
}

hilbertIndex()

uint64_t lsst::sphgeom::hilbertIndex ( uint32_t  x, uint32_t  y, int  m  )

inline

hilbertIndex returns the index of (x, y) in a 2-D Hilbert curve.

Only the m least significant bits of x and y are used. Computing the Hilbert index of a point has been measured to take 4 to 15 times as long as computing its Morton index on an Intel Core i7-3820QM CPU. With Xcode 7.3 and -O3, latency is ~19ns per call at a CPU frequency of 3.5 GHz.

Definition at line 349 of file curve.h.

                                                            {
    return mortonToHilbert(mortonIndex(x, y), m);
}

hilbertIndexInverse()

std::tuple<uint32_t, uint32_t> lsst::sphgeom::hilbertIndexInverse ( uint64_t  h, int  m  )

inline

hilbertIndexInverse returns the point (x, y) with Hilbert index h, where x and y are m bit integers.

Definition at line 361 of file curve.h.

                                                                           {
    return mortonIndexInverse(hilbertToMorton(h, m));
}

hilbertToMorton()

uint64_t lsst::sphgeom::hilbertToMorton ( uint64_t  h, int  m  )

inline

hilbertToMorton converts the 2m-bit Hilbert index h to the corresponding Morton index.

Definition at line 290 of file curve.h.

                                                   {
    alignas(64) static uint8_t const HILBERT_INVERSE_LUT_3[256] = {
        0x40, 0x02, 0x03, 0xc1, 0x04, 0x45, 0x47, 0x86,
        0x0c, 0x4d, 0x4f, 0x8e, 0xcb, 0x89, 0x88, 0x4a,
        0x20, 0x61, 0x63, 0xa2, 0x68, 0x2a, 0x2b, 0xe9,
        0x6c, 0x2e, 0x2f, 0xed, 0xa7, 0xe6, 0xe4, 0x25,
        0x30, 0x71, 0x73, 0xb2, 0x78, 0x3a, 0x3b, 0xf9,
        0x7c, 0x3e, 0x3f, 0xfd, 0xb7, 0xf6, 0xf4, 0x35,
        0xdf, 0x9d, 0x9c, 0x5e, 0x9b, 0xda, 0xd8, 0x19,
        0x93, 0xd2, 0xd0, 0x11, 0x54, 0x16, 0x17, 0xd5,
        0x00, 0x41, 0x43, 0x82, 0x48, 0x0a, 0x0b, 0xc9,
        0x4c, 0x0e, 0x0f, 0xcd, 0x87, 0xc6, 0xc4, 0x05,
        0x50, 0x12, 0x13, 0xd1, 0x14, 0x55, 0x57, 0x96,
        0x1c, 0x5d, 0x5f, 0x9e, 0xdb, 0x99, 0x98, 0x5a,
        0x70, 0x32, 0x33, 0xf1, 0x34, 0x75, 0x77, 0xb6,
        0x3c, 0x7d, 0x7f, 0xbe, 0xfb, 0xb9, 0xb8, 0x7a,
        0xaf, 0xee, 0xec, 0x2d, 0xe7, 0xa5, 0xa4, 0x66,
        0xe3, 0xa1, 0xa0, 0x62, 0x28, 0x69, 0x6b, 0xaa,
        0xff, 0xbd, 0xbc, 0x7e, 0xbb, 0xfa, 0xf8, 0x39,
        0xb3, 0xf2, 0xf0, 0x31, 0x74, 0x36, 0x37, 0xf5,
        0x9f, 0xde, 0xdc, 0x1d, 0xd7, 0x95, 0x94, 0x56,
        0xd3, 0x91, 0x90, 0x52, 0x18, 0x59, 0x5b, 0x9a,
        0x8f, 0xce, 0xcc, 0x0d, 0xc7, 0x85, 0x84, 0x46,
        0xc3, 0x81, 0x80, 0x42, 0x08, 0x49, 0x4b, 0x8a,
        0x60, 0x22, 0x23, 0xe1, 0x24, 0x65, 0x67, 0xa6,
        0x2c, 0x6d, 0x6f, 0xae, 0xeb, 0xa9, 0xa8, 0x6a,
        0xbf, 0xfe, 0xfc, 0x3d, 0xf7, 0xb5, 0xb4, 0x76,
        0xf3, 0xb1, 0xb0, 0x72, 0x38, 0x79, 0x7b, 0xba,
        0xef, 0xad, 0xac, 0x6e, 0xab, 0xea, 0xe8, 0x29,
        0xa3, 0xe2, 0xe0, 0x21, 0x64, 0x26, 0x27, 0xe5,
        0xcf, 0x8d, 0x8c, 0x4e, 0x8b, 0xca, 0xc8, 0x09,
        0x83, 0xc2, 0xc0, 0x01, 0x44, 0x06, 0x07, 0xc5,
        0x10, 0x51, 0x53, 0x92, 0x58, 0x1a, 0x1b, 0xd9,
        0x5c, 0x1e, 0x1f, 0xdd, 0x97, 0xd6, 0xd4, 0x15
    };
    uint64_t z = 0;
    uint64_t i = 0;
    for (m = 2 * m; m >= 6;) {
        m -= 6;
        uint8_t j = HILBERT_INVERSE_LUT_3[i | ((h >> m) & 0x3f)];
        z = (z << 6) | (j & 0x3f);
        i = j & 0xc0;
    }
    if (m != 0) {
        // m = 2 or 4
        int r = 6 - m;
        uint8_t j = HILBERT_INVERSE_LUT_3[i | ((h << r) & 0x3f)];
        z = (z << m) | ((j & 0x3f) >> r);
    }
    return z;
}

invert()

Relationship lsst::sphgeom::invert ( Relationship  r )

inline

Given the relationship between two sets A and B (i.e.

the output of A.relate(B)), invert returns the relationship between B and A (B.relate(A)).

Definition at line 55 of file Relationship.h.

                                           {
    // If A is disjoint from B, then B is disjoint from A. But if A contains B
    // then B is within A, so the corresponding bits must be swapped.
    return (r & DISJOINT) | ((r & CONTAINS) << 1) | ((r & WITHIN) >> 1);
}

log2() [1/2]

uint8_t lsst::sphgeom::log2 ( uint32_t  x )

inline

log2 returns the index of the most significant 1 bit in x. If x is 0, the return value is 0.

A beautiful algorithm to find this index is presented in:

Using de Bruijn Sequences to Index a 1 in a Computer Word C. E. Leiserson, H. Prokop, and K. H. Randall. http://supertech.csail.mit.edu/papers/debruijn.pdf

Definition at line 125 of file curve.h.

                                {
    // See https://graphics.stanford.edu/~seander/bithacks.html#IntegerLogDeBruijn
    alignas(32) static uint8_t const PERFECT_HASH_TABLE[32] = {
        0,  9,  1, 10, 13, 21,  2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
        8, 12, 20, 28, 15, 17, 24,  7, 19, 27, 23,  6, 26,  5, 4, 31
    };
    uint32_t const DE_BRUIJN_SEQUENCE = UINT32_C(0x07c4acdd);
    x |= (x >> 1);
    x |= (x >> 2);
    x |= (x >> 4);
    x |= (x >> 8);
    x |= (x >> 16);
    return PERFECT_HASH_TABLE[(DE_BRUIJN_SEQUENCE * x) >> 27];
}

log2() [2/2]

uint8_t lsst::sphgeom::log2 ( uint64_t  x )

inline

log2 returns the index of the most significant 1 bit in x. If x is 0, the return value is 0.

A beautiful algorithm to find this index is presented in:

Using de Bruijn Sequences to Index a 1 in a Computer Word C. E. Leiserson, H. Prokop, and K. H. Randall. http://supertech.csail.mit.edu/papers/debruijn.pdf

Definition at line 98 of file curve.h.

                                {
    alignas(64) static uint8_t const PERFECT_HASH_TABLE[64] = {
         0,  1,  2,  7,  3, 13,  8, 19,  4, 25, 14, 28,  9, 34, 20, 40,
         5, 17, 26, 38, 15, 46, 29, 48, 10, 31, 35, 54, 21, 50, 41, 57,
        63,  6, 12, 18, 24, 27, 33, 39, 16, 37, 45, 47, 30, 53, 49, 56,
        62, 11, 23, 32, 36, 44, 52, 55, 61, 22, 43, 51, 60, 42, 59, 58
    };
    uint64_t const DE_BRUIJN_SEQUENCE = UINT64_C(0x0218a392cd3d5dbf);
    // First ensure that all bits below the MSB are set.
    x |= (x >> 1);
    x |= (x >> 2);
    x |= (x >> 4);
    x |= (x >> 8);
    x |= (x >> 16);
    x |= (x >> 32);
    // Then, subtract them away.
    x = x - (x >> 1);
    // Multiplication by x is now a shift by the index i of the MSB.
    //
    // By definition, the value of the upper 6 bits of a 64-bit De Bruijn
    // sequence left shifted by i is different for every value of i in
    // [0, 64). It can therefore be used as an an index into a lookup table
    // that recovers i. In other words, (DE_BRUIJN_SEQUENCE * x) >> 58 is a
    // minimal perfect hash function for 64 bit powers of 2.
    return PERFECT_HASH_TABLE[(DE_BRUIJN_SEQUENCE * x) >> 58];
}

mortonIndex()

uint64_t lsst::sphgeom::mortonIndex ( uint32_t  x, uint32_t  y  )

inline

mortonIndex interleaves the bits of x and y.

The 32 even bits of the return value will be the bits of x, and the 32 odd bits those of y. This is the z-value of (x,y) defined by the Morton order function. See https://en.wikipedia.org/wiki/Z-order_curve for more information.

Definition at line 148 of file curve.h.

                                                        {
        // This is just a 64-bit extension of:
        // http://graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN
        uint64_t b = y;
        uint64_t a = x;
        b = (b | (b << 16)) & UINT64_C(0x0000ffff0000ffff);
        a = (a | (a << 16)) & UINT64_C(0x0000ffff0000ffff);
        b = (b | (b << 8)) & UINT64_C(0x00ff00ff00ff00ff);
        a = (a | (a << 8)) & UINT64_C(0x00ff00ff00ff00ff);
        b = (b | (b << 4)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
        a = (a | (a << 4)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
        b = (b | (b << 2)) & UINT64_C(0x3333333333333333);
        a = (a | (a << 2)) & UINT64_C(0x3333333333333333);
        b = (b | (b << 1)) & UINT64_C(0x5555555555555555);
        a = (a | (a << 1)) & UINT64_C(0x5555555555555555);
        return a | (b << 1);
    }

mortonIndexInverse()

std::tuple<uint32_t, uint32_t> lsst::sphgeom::mortonIndexInverse ( uint64_t  z )

inline

mortonIndexInverse separates the even and odd bits of z.

The 32 even bits of z are returned in the first element of the result tuple, and the 32 odd bits in the second. This is the inverse of mortonIndex().

Definition at line 195 of file curve.h.

                                                                       {
        uint64_t x = z & UINT64_C(0x5555555555555555);
        uint64_t y = (z >> 1) & UINT64_C(0x5555555555555555);
        x = (x | (x >> 1)) & UINT64_C(0x3333333333333333);
        y = (y | (y >> 1)) & UINT64_C(0x3333333333333333);
        x = (x | (x >> 2)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
        y = (y | (y >> 2)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
        x = (x | (x >> 4)) & UINT64_C(0x00ff00ff00ff00ff);
        y = (y | (y >> 4)) & UINT64_C(0x00ff00ff00ff00ff);
        x = (x | (x >> 8)) & UINT64_C(0x0000ffff0000ffff);
        y = (y | (y >> 8)) & UINT64_C(0x0000ffff0000ffff);
        return std::make_tuple(static_cast<uint32_t>(x | (x >> 16)),
                               static_cast<uint32_t>(y | (y >> 16)));
    }

mortonToHilbert()

uint64_t lsst::sphgeom::mortonToHilbert ( uint64_t  z, int  m  )

inline

mortonToHilbert converts the 2m-bit Morton index z to the corresponding Hilbert index.

Definition at line 236 of file curve.h.

                                                   {
    alignas(64) static uint8_t const HILBERT_LUT_3[256] = {
        0x40, 0xc3, 0x01, 0x02, 0x04, 0x45, 0x87, 0x46,
        0x8e, 0x8d, 0x4f, 0xcc, 0x08, 0x49, 0x8b, 0x4a,
        0xfa, 0x3b, 0xf9, 0xb8, 0x7c, 0xff, 0x3d, 0x3e,
        0xf6, 0x37, 0xf5, 0xb4, 0xb2, 0xb1, 0x73, 0xf0,
        0x10, 0x51, 0x93, 0x52, 0xde, 0x1f, 0xdd, 0x9c,
        0x54, 0xd7, 0x15, 0x16, 0x58, 0xdb, 0x19, 0x1a,
        0x20, 0x61, 0xa3, 0x62, 0xee, 0x2f, 0xed, 0xac,
        0x64, 0xe7, 0x25, 0x26, 0x68, 0xeb, 0x29, 0x2a,
        0x00, 0x41, 0x83, 0x42, 0xce, 0x0f, 0xcd, 0x8c,
        0x44, 0xc7, 0x05, 0x06, 0x48, 0xcb, 0x09, 0x0a,
        0x50, 0xd3, 0x11, 0x12, 0x14, 0x55, 0x97, 0x56,
        0x9e, 0x9d, 0x5f, 0xdc, 0x18, 0x59, 0x9b, 0x5a,
        0xba, 0xb9, 0x7b, 0xf8, 0xb6, 0xb5, 0x77, 0xf4,
        0x3c, 0x7d, 0xbf, 0x7e, 0xf2, 0x33, 0xf1, 0xb0,
        0x60, 0xe3, 0x21, 0x22, 0x24, 0x65, 0xa7, 0x66,
        0xae, 0xad, 0x6f, 0xec, 0x28, 0x69, 0xab, 0x6a,
        0xaa, 0xa9, 0x6b, 0xe8, 0xa6, 0xa5, 0x67, 0xe4,
        0x2c, 0x6d, 0xaf, 0x6e, 0xe2, 0x23, 0xe1, 0xa0,
        0x9a, 0x99, 0x5b, 0xd8, 0x96, 0x95, 0x57, 0xd4,
        0x1c, 0x5d, 0x9f, 0x5e, 0xd2, 0x13, 0xd1, 0x90,
        0x70, 0xf3, 0x31, 0x32, 0x34, 0x75, 0xb7, 0x76,
        0xbe, 0xbd, 0x7f, 0xfc, 0x38, 0x79, 0xbb, 0x7a,
        0xca, 0x0b, 0xc9, 0x88, 0x4c, 0xcf, 0x0d, 0x0e,
        0xc6, 0x07, 0xc5, 0x84, 0x82, 0x81, 0x43, 0xc0,
        0xea, 0x2b, 0xe9, 0xa8, 0x6c, 0xef, 0x2d, 0x2e,
        0xe6, 0x27, 0xe5, 0xa4, 0xa2, 0xa1, 0x63, 0xe0,
        0x30, 0x71, 0xb3, 0x72, 0xfe, 0x3f, 0xfd, 0xbc,
        0x74, 0xf7, 0x35, 0x36, 0x78, 0xfb, 0x39, 0x3a,
        0xda, 0x1b, 0xd9, 0x98, 0x5c, 0xdf, 0x1d, 0x1e,
        0xd6, 0x17, 0xd5, 0x94, 0x92, 0x91, 0x53, 0xd0,
        0x8a, 0x89, 0x4b, 0xc8, 0x86, 0x85, 0x47, 0xc4,
        0x0c, 0x4d, 0x8f, 0x4e, 0xc2, 0x03, 0xc1, 0x80
    };
    uint64_t h = 0;
    uint64_t i = 0;
    for (m = 2 * m; m >= 6;) {
        m -= 6;
        uint8_t j = HILBERT_LUT_3[i | ((z >> m) & 0x3f)];
        h = (h << 6) | (j & 0x3f);
        i = j & 0xc0;
    }
    if (m != 0) {
        // m = 2 or 4
        int r = 6 - m;
        uint8_t j = HILBERT_LUT_3[i | ((z << r) & 0x3f)];
        h = (h << m) | ((j & 0x3f) >> r);
    }
    return h;
}

operator*() [1/2]

Angle lsst::sphgeom::operator* ( double  a, Angle const &  b  )

inline

Definition at line 98 of file Angle.h.

{ return b * a; }

operator*() [2/2]

Vector3d lsst::sphgeom::operator* ( double  s, Vector3d const &  v  )

inline

Definition at line 165 of file Vector3d.h.

{ return v * s; }

operator<<() [1/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Angle const &	a
	)

Definition at line 34 of file Angle.cc.

                                                          {
    char buf[32];
    std::snprintf(buf, sizeof(buf), "%.17g", a.asRadians());
    return os << buf;
}

operator<<() [2/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		AngleInterval const &	i
	)

Definition at line 34 of file AngleInterval.cc.

                                                                  {
    return os << '[' << i.getA() << ", " << i.getB() << ']';
}

operator<<() [3/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Box const &	b
	)

Definition at line 475 of file Box.cc.

                                                        {
    return os << "{\"Box\": [" << b.getLon() << ", " << b.getLat() << "]}";
}

operator<<() [4/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Box3d const &	b
	)

Definition at line 34 of file Box3d.cc.

                                                          {
    return os << "{\"Box3d\": [" << b.x() << ", " << b.y() << ", " << b.z() << "]}";
}

operator<<() [5/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Circle const &	c
	)

Definition at line 354 of file Circle.cc.

                                                           {
    char tail[32];
    std::snprintf(tail, sizeof(tail), ", %.17g]}", c.getSquaredChordLength());
    return os << "{\"Circle\": [" << c.getCenter() << tail;
}

operator<<() [6/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		ConvexPolygon const &	p
	)

Definition at line 382 of file ConvexPolygon.cc.

                                                                  {
    typedef std::vector<UnitVector3d>::const_iterator VertexIterator;
    VertexIterator v = p.getVertices().begin();
    VertexIterator const end = p.getVertices().end();
    os << "{\"ConvexPolygon\": [" << *v;
    for (++v; v != end; ++v) { os << ", " << *v; }
    os << "]}";
    return os;
}

operator<<() [7/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Ellipse const &	e
	)

Definition at line 388 of file Ellipse.cc.

                                                            {
    os << "{\"Ellipse\": [" << e.getTransformMatrix() << ", "
       << e.getAlpha() << ", " << e.getBeta() << "]}";
    return os;
}

operator<<() [8/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Interval1d const &	i
	)

Definition at line 35 of file Interval1d.cc.

                                                               {
    char buf[64];
    std::snprintf(buf, sizeof(buf), "[%.17g, %.17g]", i.getA(), i.getB());
    return os << buf;
}

operator<<() [9/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		LonLat const &	p
	)

Definition at line 82 of file LonLat.cc.

                                                           {
    return os << '[' << p.getLon() << ", " << p.getLat() << ']';
}

operator<<() [10/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Matrix3d const &	m
	)

Definition at line 35 of file Matrix3d.cc.

                                                             {
    return os << '[' << m.getRow(0) << ", " << m.getRow(1) << ", " << m.getRow(2) << ']';
}

operator<<() [11/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		NormalizedAngleInterval const &	i
	)

Definition at line 284 of file NormalizedAngleInterval.cc.

{
    return os << '[' << i.getA() << ", " << i.getB() << ']';
}

operator<<() [12/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		RangeSet const &	s
	)

Definition at line 616 of file RangeSet.cc.

                                                             {
    os << "{\"RangeSet\": [";
    bool first = true;
    for (auto const & t: s) {
        if (!first) {
            os << ", ";
        }
        first = false;
        os << '[' << std::get<0>(t) << ", " << std::get<1>(t) << ']';
    }
    os << "]}";
    return os;
}

operator<<() [13/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		UnitVector3d const &	v
	)

Definition at line 73 of file UnitVector3d.cc.

                                                                 {
    return os << static_cast<Vector3d const &>(v);
}

operator<<() [14/14]

std::ostream & lsst::sphgeom::operator<<	(	std::ostream &	os,
		Vector3d const &	v
	)

Definition at line 133 of file Vector3d.cc.

                                                             {
    char buf[128];
    std::snprintf(buf, sizeof(buf), "[%.17g, %.17g, %.17g]",
                  v.x(), v.y(), v.z());
    return os << buf;
}

orientation()

int lsst::sphgeom::orientation	(	UnitVector3d const &	a,
		UnitVector3d const &	b,
		UnitVector3d const &	c
	)

orientation computes and returns the orientations of 3 unit vectors a, b and c.

The return value is +1 if the vectors are in counter-clockwise orientation, 0 if they are coplanar, colinear or identical, and -1 if they are in clockwise orientation.

This is equivalent to computing the sign of the scalar triple product a · (b x c), which is the sign of the determinant of the 3x3 matrix with a, b and c as columns/rows.

The implementation proceeds by first computing a double precision approximation, and then falling back to arbitrary precision arithmetic when necessary. Consequently, the result is exact.

Definition at line 135 of file orientation.cc.

{
    // This constant is a little more than 5ε, where ε = 2^-53. When multiplied
    // by the permanent of |M|, it gives an error bound on the determinant of
    // M. Here, M is a 3x3 matrix and |M| denotes the matrix obtained by
    // taking the absolute value of each of its components. The derivation of
    // this proceeds in the same manner as the derivation of the error bounds
    // in section 4.3 of:
    //
    //     Adaptive Precision Floating-Point Arithmetic
    //     and Fast Robust Geometric Predicates,
    //     Jonathan Richard Shewchuk,
    //     Discrete & Computational Geometry 18(3):305–363, October 1997.
    //
    // available online at http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf
    static double const relativeError = 5.6e-16;
    // Because all 3 unit vectors are normalized, the maximum absolute value of
    // any vector component, cross product component or dot product term in
    // the calculation is very close to 1. The permanent of |M| must therefore
    // be below 3 + c, where c is some small multiple of ε. This constant, a
    // little larger than 3 * 5ε, is an upper bound on the absolute error in
    // the determinant calculation.
    static double const maxAbsoluteError = 1.7e-15;
    // This constant accounts for floating point underflow (assuming hardware
    // without gradual underflow, just to be conservative) in the computation
    // of det(M). It is a little more than 14 * 2^-1022.
    static double const minAbsoluteError = 4.0e-307;
 
    double bycz = b.y() * c.z();
    double bzcy = b.z() * c.y();
    double bzcx = b.z() * c.x();
    double bxcz = b.x() * c.z();
    double bxcy = b.x() * c.y();
    double bycx = b.y() * c.x();
    double determinant = a.x() * (bycz - bzcy) +
                         a.y() * (bzcx - bxcz) +
                         a.z() * (bxcy - bycx);
    if (determinant > maxAbsoluteError) {
        return 1;
    } else if (determinant < -maxAbsoluteError) {
        return -1;
    }
    // Expend some more effort on what is hopefully a tighter error bound
    // before falling back on arbitrary precision arithmetic.
    double permanent = std::fabs(a.x()) * (std::fabs(bycz) + std::fabs(bzcy)) +
                       std::fabs(a.y()) * (std::fabs(bzcx) + std::fabs(bxcz)) +
                       std::fabs(a.z()) * (std::fabs(bxcy) + std::fabs(bycx));
    double maxError = relativeError * permanent + minAbsoluteError;
    if (determinant > maxError) {
        return 1;
    } else if (determinant < -maxError) {
        return -1;
    }
    // Avoid the slow path when any two inputs are identical or antipodal.
    if (a == b || b == c || a == c || a == -b || b == -c || a == -c) {
        return 0;
    }
    return orientationExact(a, b, c);
}

orientationExact()

int lsst::sphgeom::orientationExact	(	Vector3d const &	a,
		Vector3d const &	b,
		Vector3d const &	c
	)

orientationExact computes and returns the orientations of 3 vectors a, b and c, which need not be normalized but are assumed to have finite components.

The return value is +1 if the vectors a, b, and c are in counter-clockwise orientation, 0 if they are coplanar, colinear, or identical, and -1 if they are in clockwise orientation. The implementation uses arbitrary precision arithmetic to avoid floating point rounding error, underflow and overflow.

Definition at line 76 of file orientation.cc.

{
    // Product mantissa storage buffers.
    uint32_t mantissaBuffers[6][6];
    // Product mantissas.
    BigInteger mantissas[6] = {
        BigInteger(mantissaBuffers[0], 6),
        BigInteger(mantissaBuffers[1], 6),
        BigInteger(mantissaBuffers[2], 6),
        BigInteger(mantissaBuffers[3], 6),
        BigInteger(mantissaBuffers[4], 6),
        BigInteger(mantissaBuffers[5], 6)
    };
    BigFloat products[6] = {
        BigFloat(&mantissas[0]),
        BigFloat(&mantissas[1]),
        BigFloat(&mantissas[2]),
        BigFloat(&mantissas[3]),
        BigFloat(&mantissas[4]),
        BigFloat(&mantissas[5])
    };
    // An accumulator and its storage.
    uint32_t accumulatorBuffer[512];
    BigInteger accumulator(accumulatorBuffer,
                           sizeof(accumulatorBuffer) / sizeof(uint32_t));
    // Compute the products in the determinant. Performing all multiplication
    // up front means that each product mantissa occupies at most 3*53 bits.
    computeProduct(products[0], a.x(), b.y(), c.z());
    computeProduct(products[1], a.x(), b.z(), c.y());
    computeProduct(products[2], a.y(), b.z(), c.x());
    computeProduct(products[3], a.y(), b.x(), c.z());
    computeProduct(products[4], a.z(), b.x(), c.y());
    computeProduct(products[5], a.z(), b.y(), c.x());
    mantissas[1].negate();
    mantissas[3].negate();
    mantissas[5].negate();
    // Sort the array of products in descending exponent order.
    std::sort(products, products + 6, [](BigFloat const & a, BigFloat const & b) {
        return a.exponent > b.exponent;
    });
    // First, initialize the accumulator to the product with the highest
    // exponent, then add the remaining products. Prior to each addition, we
    // must shift the accumulated value so that its radix point lines up with
    // the the radix point of the product to add.
    //
    // More precisely, at each step we have an accumulated value A·2ʲ and a
    // product P·2ᵏ, and we update the accumulator to equal (A·2ʲ⁻ᵏ + P)·2ᵏ.
    // Because the products were sorted beforehand, j ≥ k and 2ʲ⁻ᵏ is an
    // integer.
    accumulator = *products[0].mantissa;
    for (int i = 1; i < 6; ++i) {
        accumulator.multiplyPow2(products[i - 1].exponent - products[i].exponent);
        accumulator.add(*products[i].mantissa);
    }
    return accumulator.getSign();
}

orientationX()

int lsst::sphgeom::orientationX	(	UnitVector3d const &	b,
		UnitVector3d const &	c
	)

orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).

Definition at line 227 of file orientation.cc.

                                                                 {
    int o = _orientationXYZ(b.y() * c.z(), b.z() * c.y());
    return (o != 0) ? o : orientationExact(UnitVector3d::X(), b, c);
}

orientationY()

int lsst::sphgeom::orientationY	(	UnitVector3d const &	b,
		UnitVector3d const &	c
	)

orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).

Definition at line 232 of file orientation.cc.

                                                                 {
    int o = _orientationXYZ(b.z() * c.x(), b.x() * c.z());
    return (o != 0) ? o : orientationExact(UnitVector3d::Y(), b, c);
}

orientationZ()

int lsst::sphgeom::orientationZ	(	UnitVector3d const &	b,
		UnitVector3d const &	c
	)

orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).

Definition at line 237 of file orientation.cc.

                                                                 {
    int o = _orientationXYZ(b.x() * c.y(), b.y() * c.x());
    return (o != 0) ? o : orientationExact(UnitVector3d::Z(), b, c);
}

sin()

double lsst::sphgeom::sin ( Angle const &  a )

inline

Definition at line 102 of file Angle.h.

{ return std::sin(a.asRadians()); }

swap()

void lsst::sphgeom::swap ( RangeSet &  a, RangeSet &  b  )

inline

Definition at line 608 of file RangeSet.h.

                                             {
    a.swap(b);
}

tan()

double lsst::sphgeom::tan ( Angle const &  a )

inline

Definition at line 104 of file Angle.h.

{ return std::tan(a.asRadians()); }

Variable Documentation

DEG_PER_RAD

constexpr double lsst::sphgeom::DEG_PER_RAD = 57.2957795130823208767981548141

constexpr

Definition at line 39 of file constants.h.

EPSILON

constexpr double lsst::sphgeom::EPSILON = 1.1102230246251565e-16

constexpr

Definition at line 54 of file constants.h.

MAX_ASIN_ERROR

constexpr double lsst::sphgeom::MAX_ASIN_ERROR = 1.5e-8

constexpr

Definition at line 45 of file constants.h.

MAX_SQUARED_CHORD_LENGTH_ERROR

constexpr double lsst::sphgeom::MAX_SQUARED_CHORD_LENGTH_ERROR = 2.5e-15

constexpr

Definition at line 50 of file constants.h.

ONE_OVER_PI

constexpr double lsst::sphgeom::ONE_OVER_PI = 0.318309886183790671537767526745

constexpr

Definition at line 37 of file constants.h.

PI

constexpr double lsst::sphgeom::PI = 3.1415926535897932384626433832795

constexpr

Definition at line 36 of file constants.h.

RAD_PER_DEG

constexpr double lsst::sphgeom::RAD_PER_DEG = 0.0174532925199432957692369076849

constexpr

Definition at line 38 of file constants.h.