lsst::sphgeom Namespace Reference

_yaml

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::ostreamoperator<< (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::ostreamoperator<< (std::ostream &, AngleInterval const &)

std::ostreamoperator<< (std::ostream &, Box const &)

std::ostreamoperator<< (std::ostream &, Box3d const &)

std::ostreamoperator<< (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::ostreamoperator<< (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::ostreamoperator<< (std::ostream &, Ellipse const &)

std::ostreamoperator<< (std::ostream &, Interval1d const &)

std::ostreamoperator<< (std::ostream &, LonLat const &)

std::ostreamoperator<< (std::ostream &, Matrix3d const &)

NormalizedAngle const & abs (NormalizedAngle const &a)

std::ostreamoperator<< (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...

template<typename Pybind11Class >
void defineClass (Pybind11Class &cls)

void swap (RangeSet &a, RangeSet &b)

std::ostreamoperator<< (std::ostream &, RangeSet const &)

Relationship invert (Relationship r)
Given the relationship between two sets A and B (i.e. More...

std::ostreamoperator<< (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::ostreamoperator<< (std::ostream &, Vector3d const &)

template<>
void defineClass (py::class_< Angle > &cls)

template<>
void defineClass (py::class_< AngleInterval, std::shared_ptr< AngleInterval >> &cls)

template<>
void defineClass (py::class_< Box, std::unique_ptr< Box >, Region > &cls)

template<>
void defineClass (py::class_< Box3d, std::shared_ptr< Box3d >> &cls)

template<>
void defineClass (py::class_< Chunker, std::shared_ptr< Chunker >> &cls)

template<>
void defineClass (py::class_< Circle, std::unique_ptr< Circle >, Region > &cls)

template<>
void defineClass (py::class_< ConvexPolygon, std::unique_ptr< ConvexPolygon >, Region > &cls)

void defineCurve (py::module &mod)

template<>
void defineClass (py::class_< Ellipse, std::unique_ptr< Ellipse >, Region > &cls)

template<>
void defineClass (py::class_< HtmPixelization, Pixelization > &cls)

template<>
void defineClass (py::class_< Interval1d, std::shared_ptr< Interval1d >> &cls)

template<>
void defineClass (py::class_< LonLat, std::shared_ptr< LonLat >> &cls)

template<>
void defineClass (py::class_< Matrix3d, std::shared_ptr< Matrix3d >> &cls)

template<>
void defineClass (py::class_< Mq3cPixelization, Pixelization > &cls)

template<>
void defineClass (py::class_< NormalizedAngle > &cls)

template<>
void defineClass (py::class_< NormalizedAngleInterval, std::shared_ptr< NormalizedAngleInterval >> &cls)

void defineOrientation (py::module &mod)

template<>
void defineClass (py::class_< Pixelization > &cls)

template<>
void defineClass (py::class_< Q3cPixelization, Pixelization > &cls)

template<>
void defineClass (py::class_< RangeSet, std::shared_ptr< RangeSet >> &cls)

template<>
void defineClass (py::class_< Region, std::unique_ptr< Region >> &cls)

void defineRelationship (py::module &mod)

void defineUtils (py::module &)

template<>
void defineClass (py::class_< UnitVector3d, std::shared_ptr< UnitVector3d >> &cls)

template<>
void defineClass (py::class_< Vector3d, std::shared_ptr< Vector3d >> &cls)

## Variables

constexpr double PI = 3.1415926535897932384626433832795

constexpr double ONE_OVER_PI = 0.318309886183790671537767526745

constexpr double MAX_ASIN_ERROR = 1.5e-8

constexpr double MAX_SQUARED_CHORD_LENGTH_ERROR = 2.5e-15

constexpr double EPSILON = 1.1102230246251565e-16

## ◆ Relationship

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

Relationship describes how two sets are related.

Definition at line 35 of file Relationship.h.

## ◆ abs() [1/2]

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

Definition at line 106 of file Angle.h.

## ◆ abs() [2/2]

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

Definition at line 150 of file NormalizedAngle.h.

150 { return a; }

## ◆ cos()

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

Definition at line 103 of file Angle.h.

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

59  {
60 #if defined(__x86_64__)
61  // x86-64 is little endian and supports unaligned loads.
62  return *reinterpret_cast<double const *>(buffer);
63 #else
64  union { uint64_t u; double d };
65  u = static_cast<uint64_t>(buffer[0]) +
66  (static_cast<uint64_t>(buffer[1]) << 8) +
67  (static_cast<uint64_t>(buffer[2]) << 16) +
68  (static_cast<uint64_t>(buffer[3]) << 24) +
69  (static_cast<uint64_t>(buffer[4]) << 32) +
70  (static_cast<uint64_t>(buffer[5]) << 40) +
71  (static_cast<uint64_t>(buffer[6]) << 48) +
72  (static_cast<uint64_t>(buffer[7]) << 56);
73  return d;
74 #endif
75 }

## ◆ defineClass() [1/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Angle > & cls )

Definition at line 36 of file _angle.cc.

36  {
37  cls.def_static("nan", &Angle::nan);
38  cls.def_static("fromDegrees", &Angle::fromDegrees);
40
41  cls.def(py::init<>());
43  cls.def(py::init<Angle>(), "angle"_a);
44  // Construct an Angle from a NormalizedAngle, enabling implicit
45  // conversion from NormalizedAngle to Angle in python via
46  // py::implicitly_convertible
47  cls.def(py::init(
48  [](NormalizedAngle &a) {
50  }),
51  "normalizedAngle"_a);
52
53  cls.def("__eq__", &Angle::operator==, py::is_operator());
54  cls.def("__ne__", &Angle::operator!=, py::is_operator());
55  cls.def("__lt__", &Angle::operator<, py::is_operator());
56  cls.def("__gt__", &Angle::operator>, py::is_operator());
57  cls.def("__le__", &Angle::operator<=, py::is_operator());
58  cls.def("__ge__", &Angle::operator>=, py::is_operator());
59
60  cls.def("__neg__", (Angle(Angle::*)() const) & Angle::operator-);
62  cls.def("__sub__",
63  (Angle(Angle::*)(Angle const &) const) & Angle::operator-,
64  py::is_operator());
65  cls.def("__mul__", &Angle::operator*, py::is_operator());
66  cls.def("__rmul__", &Angle::operator*, py::is_operator());
67  cls.def("__truediv__", (Angle(Angle::*)(double) const) & Angle::operator/,
68  py::is_operator());
69  cls.def("__truediv__",
70  (double (Angle::*)(Angle const &) const) & Angle::operator/,
71  py::is_operator());
72
74  cls.def("__isub__", &Angle::operator-=);
75  cls.def("__imul__", &Angle::operator*=);
76  cls.def("__itruediv__", &Angle::operator/=);
77
78  cls.def("asDegrees", &Angle::asDegrees);
80  cls.def("isNormalized", &Angle::isNormalized);
81  cls.def("isNan", &Angle::isNan);
82
83  cls.def("__str__", [](Angle const &self) {
85  });
86  cls.def("__repr__", [](Angle const &self) {
88  });
89
90  cls.def("__reduce__", [cls](Angle const &self) {
92  });
93 }

## ◆ defineClass() [2/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< AngleInterval, std::shared_ptr< AngleInterval >> & cls )

Definition at line 36 of file _angleInterval.cc.

37  {
38  python::defineInterval<decltype(cls), AngleInterval, Angle>(cls);
39
40  cls.def_static("fromDegrees", &AngleInterval::fromDegrees, "x"_a, "y"_a);
42  cls.def_static("empty", &AngleInterval::empty);
43  cls.def_static("full", &AngleInterval::full);
44
45  cls.def(py::init<>());
46  cls.def(py::init<Angle>(), "x"_a);
47  cls.def(py::init<Angle, Angle>(), "x"_a, "y"_a);
48  cls.def(py::init<AngleInterval const &>(), "interval"_a);
49
50  cls.def("__str__", [](AngleInterval const &self) {
51  return py::str("[{!s}, {!s}]")
53  });
54  cls.def("__repr__", [](AngleInterval const &self) {
57  });
58 }

## ◆ defineClass() [3/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Box, std::unique_ptr< Box >, Region > & cls )

Definition at line 56 of file _box.cc.

56  {
57  cls.attr("TYPE_CODE") = py::int_(Box::TYPE_CODE);
58
59  cls.def_static("fromDegrees", &Box::fromDegrees, "lon1"_a, "lat1"_a,
60  "lon2"_a, "lat2"_a);
62  "lon2"_a, "lat2"_a);
63  cls.def_static("empty", &Box::empty);
64  cls.def_static("full", &Box::full);
66  "lat"_a);
67  cls.def_static("allLongitudes", &Box::allLongitudes);
68  cls.def_static("allLatitudes", &Box::allLatitudes);
69
70  cls.def(py::init<>());
71  cls.def(py::init<LonLat const &>(), "point"_a);
72  cls.def(py::init<LonLat const &, LonLat const &>(), "point1"_a, "point2"_a);
73  cls.def(py::init<LonLat const &, Angle, Angle>(), "center"_a, "width"_a,
74  "height"_a);
75  cls.def(py::init<NormalizedAngleInterval const &, AngleInterval const &>(),
76  "lon"_a, "lat"_a);
77  cls.def(py::init<Box const &>(), "box"_a);
78
79  cls.def("__eq__", (bool (Box::*)(Box const &) const) & Box::operator==,
80  py::is_operator());
81  cls.def("__eq__", (bool (Box::*)(LonLat const &) const) & Box::operator==,
82  py::is_operator());
83  cls.def("__ne__", (bool (Box::*)(Box const &) const) & Box::operator!=,
84  py::is_operator());
85  cls.def("__ne__", (bool (Box::*)(LonLat const &) const) & Box::operator!=,
86  py::is_operator());
87  cls.def("__contains__",
88  (bool (Box::*)(LonLat const &) const) & Box::contains,
89  py::is_operator());
90  cls.def("__contains__", (bool (Box::*)(Box const &) const) & Box::contains,
91  py::is_operator());
92  // Rewrap this base class method since there are overloads in this subclass
93  cls.def("__contains__",
94  (bool (Box::*)(UnitVector3d const &) const) & Box::contains,
95  py::is_operator());
96
97  cls.def("getLon", &Box::getLon);
98  cls.def("getLat", &Box::getLat);
99  cls.def("isEmpty", &Box::isEmpty);
100  cls.def("isFull", &Box::isFull);
101  cls.def("getCenter", &Box::getCenter);
102  cls.def("getWidth", &Box::getWidth);
103  cls.def("getHeight", &Box::getHeight);
104  cls.def("contains", (bool (Box::*)(LonLat const &) const) & Box::contains);
105  cls.def("contains", (bool (Box::*)(Box const &) const) & Box::contains);
106  // Rewrap this base class method since there are overloads in this subclass
107  cls.def("contains",
108  (bool (Box::*)(UnitVector3d const &) const) & Box::contains);
109  cls.def("isDisjointFrom",
110  (bool (Box::*)(LonLat const &) const) & Box::isDisjointFrom);
111  cls.def("isDisjointFrom",
112  (bool (Box::*)(Box const &) const) & Box::isDisjointFrom);
113  cls.def("intersects",
114  (bool (Box::*)(LonLat const &) const) & Box::intersects);
115  cls.def("intersects", (bool (Box::*)(Box const &) const) & Box::intersects);
116  cls.def("isWithin", (bool (Box::*)(LonLat const &) const) & Box::isWithin);
117  cls.def("isWithin", (bool (Box::*)(Box const &) const) & Box::isWithin);
118  cls.def("clipTo", (Box & (Box::*)(LonLat const &)) & Box::clipTo);
119  cls.def("clipTo", (Box & (Box::*)(Box const &)) & Box::clipTo);
120  cls.def("clippedTo", (Box(Box::*)(LonLat const &) const) & Box::clippedTo);
121  cls.def("clippedTo", (Box(Box::*)(Box const &) const) & Box::clippedTo);
122  cls.def("expandTo", (Box & (Box::*)(LonLat const &)) & Box::expandTo);
123  cls.def("expandTo", (Box & (Box::*)(Box const &)) & Box::expandTo);
124  cls.def("expandedTo",
125  (Box(Box::*)(LonLat const &) const) & Box::expandedTo);
126  cls.def("expandedTo", (Box(Box::*)(Box const &) const) & Box::expandedTo);
127  cls.def("dilateBy", (Box & (Box::*)(Angle)) & Box::dilateBy, "angle"_a);
128  cls.def("dilateBy", (Box & (Box::*)(Angle, Angle)) & Box::dilateBy,
129  "width"_a, "height"_a);
130  cls.def("dilatedBy", (Box(Box::*)(Angle) const) & Box::dilatedBy,
131  "angle"_a);
132  cls.def("dilatedBy", (Box(Box::*)(Angle, Angle) const) & Box::dilatedBy,
133  "width"_a, "height"_a);
134  cls.def("erodeBy", (Box & (Box::*)(Angle)) & Box::erodeBy, "angle"_a);
135  cls.def("erodeBy", (Box & (Box::*)(Angle, Angle)) & Box::erodeBy, "width"_a,
136  "height"_a);
137  cls.def("erodedBy", (Box(Box::*)(Angle) const) & Box::erodedBy, "angle"_a);
138  cls.def("erodedBy", (Box(Box::*)(Angle, Angle) const) & Box::erodedBy,
139  "width"_a, "height"_a);
140  cls.def("getArea", &Box::getArea);
141  cls.def("relate",
142  (Relationship(Box::*)(LonLat const &) const) & Box::relate,
143  "point"_a);
144  // Rewrap this base class method since there are overloads in this subclass
145  cls.def("relate",
146  (Relationship(Box::*)(Region const &) const) & Box::relate,
147  "region"_a);
148
149  // Note that the Region interface has already been wrapped.
150
151  // The lambda is necessary for now; returning the unique pointer
152  // directly leads to incorrect results and crashes.
153  cls.def_static("decode",
154  [](py::bytes bytes) { return decode(bytes).release(); },
155  "bytes"_a);
156
157  cls.def("__str__", [](Box const &self) {
158  return py::str("Box({!s}, {!s})").format(self.getLon(), self.getLat());
159  });
160  cls.def("__repr__", [](Box const &self) {
161  return py::str("Box({!r}, {!r})").format(self.getLon(), self.getLat());
162  });
163  cls.def(py::pickle(
164  [](const Box &self) { return python::encode(self); },
165  [](py::bytes bytes) { return decode(bytes).release(); }));
166 }

## ◆ defineClass() [4/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Box3d, std::shared_ptr< Box3d >> & cls )

Definition at line 37 of file _box3d.cc.

37  {
38  cls.def_static("empty", &Box3d::empty);
39  cls.def_static("full", &Box3d::full);
40  cls.def_static("aroundUnitSphere", &Box3d::aroundUnitSphere);
41
42  cls.def(py::init<>());
43  cls.def(py::init<Vector3d const &>(), "vector"_a);
44  cls.def(py::init<Vector3d const &, Vector3d const &>(), "vector1"_a,
45  "vector2"_a);
46  cls.def(py::init<Vector3d const &, double, double, double>(), "center"_a,
47  "halfWidth"_a, "halfHeight"_a, "halfDepth"_a);
48  cls.def(py::init<Interval1d const &, Interval1d const &,
49  Interval1d const &>(),
50  "x"_a, "y"_a, "z"_a);
51  cls.def(py::init<Box3d const &>(), "box3d"_a);
52
53  cls.def("__eq__",
54  (bool (Box3d::*)(Box3d const &) const) & Box3d::operator==,
55  py::is_operator());
56  cls.def("__eq__",
57  (bool (Box3d::*)(Vector3d const &) const) & Box3d::operator==,
58  py::is_operator());
59  cls.def("__ne__",
60  (bool (Box3d::*)(Box3d const &) const) & Box3d::operator!=,
61  py::is_operator());
62  cls.def("__ne__",
63  (bool (Box3d::*)(Vector3d const &) const) & Box3d::operator!=,
64  py::is_operator());
65  cls.def("__contains__",
66  (bool (Box3d::*)(Vector3d const &) const) & Box3d::contains,
67  py::is_operator());
68  cls.def("__contains__",
69  (bool (Box3d::*)(Box3d const &) const) & Box3d::contains,
70  py::is_operator());
71  cls.def("__len__", [](Box3d const &self) { return py::int_(3); });
72  cls.def("__getitem__", [](Box3d const &self, py::int_ row) {
73  return self(static_cast<int>(python::convertIndex(3, row)));
74  });
75
76  cls.def("x", &Box3d::x);
77  cls.def("y", &Box3d::y);
78  cls.def("z", &Box3d::z);
79  cls.def("isEmpty", &Box3d::isEmpty);
80  cls.def("isFull", &Box3d::isFull);
81  cls.def("getCenter", &Box3d::getCenter);
82  cls.def("getWidth", &Box3d::getWidth);
83  cls.def("getHeight", &Box3d::getHeight);
84  cls.def("getDepth", &Box3d::getDepth);
85
86  cls.def("contains",
87  (bool (Box3d::*)(Vector3d const &) const) & Box3d::contains);
88  cls.def("contains",
89  (bool (Box3d::*)(Box3d const &) const) & Box3d::contains);
90  cls.def("isDisjointFrom",
91  (bool (Box3d::*)(Vector3d const &) const) & Box3d::isDisjointFrom);
92  cls.def("isDisjointFrom",
93  (bool (Box3d::*)(Box3d const &) const) & Box3d::isDisjointFrom);
94  cls.def("intersects",
95  (bool (Box3d::*)(Vector3d const &) const) & Box3d::intersects);
96  cls.def("intersects",
97  (bool (Box3d::*)(Box3d const &) const) & Box3d::intersects);
98  cls.def("isWithin",
99  (bool (Box3d::*)(Vector3d const &) const) & Box3d::isWithin);
100  cls.def("isWithin",
101  (bool (Box3d::*)(Box3d const &) const) & Box3d::isWithin);
102
103  cls.def("clipTo", (Box3d & (Box3d::*)(Vector3d const &)) & Box3d::clipTo);
104  cls.def("clipTo", (Box3d & (Box3d::*)(Box3d const &)) & Box3d::clipTo);
105  cls.def("clippedTo",
106  (Box3d(Box3d::*)(Vector3d const &) const) & Box3d::clippedTo);
107  cls.def("clippedTo",
108  (Box3d(Box3d::*)(Box3d const &) const) & Box3d::clippedTo);
109  cls.def("expandTo",
110  (Box3d & (Box3d::*)(Vector3d const &)) & Box3d::expandTo);
111  cls.def("expandTo", (Box3d & (Box3d::*)(Box3d const &)) & Box3d::expandTo);
112  cls.def("expandedTo",
113  (Box3d(Box3d::*)(Vector3d const &) const) & Box3d::expandedTo);
114  cls.def("expandedTo",
115  (Box3d(Box3d::*)(Box3d const &) const) & Box3d::expandedTo);
116
117  cls.def("dilateBy", (Box3d & (Box3d::*)(double)) & Box3d::dilateBy,
119  cls.def("dilateBy",
120  (Box3d & (Box3d::*)(double, double, double)) & Box3d::dilateBy,
121  "width"_a, "height"_a, "depth"_a);
122  cls.def("dilatedBy", (Box3d(Box3d::*)(double) const) & Box3d::dilatedBy,
124  cls.def("dilatedBy",
125  (Box3d(Box3d::*)(double, double, double) const) & Box3d::dilatedBy,
126  "width"_a, "height"_a, "depth"_a);
127  cls.def("erodeBy", (Box3d & (Box3d::*)(double)) & Box3d::erodeBy,
129  cls.def("erodeBy",
130  (Box3d & (Box3d::*)(double, double, double)) & Box3d::erodeBy,
131  "width"_a, "height"_a, "depth"_a);
132  cls.def("erodedBy", (Box3d(Box3d::*)(double) const) & Box3d::erodedBy,
134  cls.def("erodedBy",
135  (Box3d(Box3d::*)(double, double, double) const) & Box3d::erodedBy,
136  "width"_a, "height"_a, "depth"_a);
137
138  cls.def("relate",
139  (Relationship(Box3d::*)(Vector3d const &) const) & Box3d::relate);
140  cls.def("relate",
141  (Relationship(Box3d::*)(Box3d const &) const) & Box3d::relate);
142
143  cls.def("__str__", [](Box3d const &self) {
144  return py::str("[{!s},\n"
145  " {!s},\n"
146  " {!s}]")
147  .format(self.x(), self.y(), self.z());
148  });
149  cls.def("__repr__", [](Box3d const &self) {
150  return py::str("Box3d({!r},\n"
151  " {!r},\n"
152  " {!r})")
153  .format(self.x(), self.y(), self.z());
154  });
155  cls.def("__reduce__", [cls](Box3d const &self) {
156  return py::make_tuple(cls,
157  py::make_tuple(self.x(), self.y(), self.z()));
158  });
159 }

## ◆ defineClass() [5/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Chunker, std::shared_ptr< Chunker >> & cls )

Definition at line 43 of file _chunker.cc.

43  {
44  cls.def(py::init<int32_t, int32_t>(), "numStripes"_a,
45  "numSubStripesPerStripe"_a);
46
47  cls.def("__eq__", &Chunker::operator==, py::is_operator());
48  cls.def("__ne__", &Chunker::operator!=, py::is_operator());
49
52  &Chunker::getNumSubStripesPerStripe);
53
54  cls.def("getChunksIntersecting", &Chunker::getChunksIntersecting,
55  "region"_a);
56  cls.def("getSubChunksIntersecting",
57  [](Chunker const &self, Region const &region) {
58  py::list results;
59  for (auto const &sc : self.getSubChunksIntersecting(region)) {
60  results.append(py::make_tuple(sc.chunkId, sc.subChunkIds));
61  }
62  return results;
63  },
64  "region"_a);
65  cls.def("getAllChunks", &Chunker::getAllChunks);
66  cls.def("getAllSubChunks", &Chunker::getAllSubChunks, "chunkId"_a);
67
68  cls.def("getChunkBoundingBox", &Chunker::getChunkBoundingBox, "stripe"_a, "chunk"_a);
69  cls.def("getSubChunkBoundingBox", &Chunker::getSubChunkBoundingBox, "subStripe"_a, "subChunk"_a);
70
71  cls.def("getStripe", &Chunker::getStripe, "chunkId"_a);
72  cls.def("getChunk", &Chunker::getChunk, "chunkId"_a, "stripe"_a);
73
74
75  cls.def("__str__", &toString);
76  cls.def("__repr__", &toString);
77
78  cls.def("__reduce__", [cls](Chunker const &self) {
79  return py::make_tuple(cls,
80  py::make_tuple(self.getNumStripes(),
81  self.getNumSubStripesPerStripe()));
82  });
83 }

## ◆ defineClass() [6/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Circle, std::unique_ptr< Circle >, Region > & cls )

Definition at line 52 of file _circle.cc.

52  {
53  cls.attr("TYPE_CODE") = py::int_(Circle::TYPE_CODE);
54
55  cls.def_static("empty", &Circle::empty);
56  cls.def_static("full", &Circle::full);
57  cls.def_static("squaredChordLengthFor", &Circle::squaredChordLengthFor,
58  "openingAngle"_a);
59  cls.def_static("openingAngleFor", &Circle::openingAngleFor,
60  "squaredChordLength"_a);
61
62  cls.def(py::init<>());
63  cls.def(py::init<UnitVector3d const &>(), "center"_a);
64  cls.def(py::init<UnitVector3d const &, Angle>(), "center"_a, "angle"_a);
65  cls.def(py::init<UnitVector3d const &, double>(), "center"_a,
66  "squaredChordLength"_a);
67  cls.def(py::init<Circle const &>(), "circle"_a);
68
69  cls.def("__eq__", &Circle::operator==, py::is_operator());
70  cls.def("__ne__", &Circle::operator!=, py::is_operator());
71  cls.def("__contains__",
72  (bool (Circle::*)(Circle const &) const) & Circle::contains,
73  py::is_operator());
74  // Rewrap this base class method since there are overloads in this subclass
75  cls.def("__contains__",
76  (bool (Circle::*)(UnitVector3d const &) const) & Circle::contains,
77  py::is_operator());
78
79  cls.def("isEmpty", &Circle::isEmpty);
80  cls.def("isFull", &Circle::isFull);
81  cls.def("getCenter", &Circle::getCenter);
82  cls.def("getSquaredChordLength", &Circle::getSquaredChordLength);
83  cls.def("getOpeningAngle", &Circle::getOpeningAngle);
84  cls.def("contains",
85  (bool (Circle::*)(Circle const &) const) & Circle::contains);
86  // Rewrap this base class method since there are overloads in this subclass
87  cls.def("contains",
88  (bool (Circle::*)(UnitVector3d const &) const) & Circle::contains);
89
90  cls.def("isDisjointFrom",
91  (bool (Circle::*)(UnitVector3d const &) const) &
92  Circle::isDisjointFrom);
93  cls.def("isDisjointFrom",
94  (bool (Circle::*)(Circle const &) const) & Circle::isDisjointFrom);
95  cls.def("intersects",
96  (bool (Circle::*)(UnitVector3d const &) const) &
97  Circle::intersects);
98  cls.def("intersects",
99  (bool (Circle::*)(Circle const &) const) & Circle::intersects);
100  cls.def("isWithin",
101  (bool (Circle::*)(UnitVector3d const &) const) & Circle::isWithin);
102  cls.def("isWithin",
103  (bool (Circle::*)(Circle const &) const) & Circle::isWithin);
104  cls.def("clipTo",
105  (Circle & (Circle::*)(UnitVector3d const &)) & Circle::clipTo);
106  cls.def("clipTo", (Circle & (Circle::*)(Circle const &)) & Circle::clipTo);
107  cls.def("clippedTo",
108  (Circle(Circle::*)(UnitVector3d const &) const) &
109  Circle::clippedTo);
110  cls.def("clippedTo",
111  (Circle(Circle::*)(Circle const &) const) & Circle::clippedTo);
112  cls.def("expandTo",
113  (Circle & (Circle::*)(UnitVector3d const &)) & Circle::expandTo);
114  cls.def("expandTo",
115  (Circle & (Circle::*)(Circle const &)) & Circle::expandTo);
116  cls.def("expandedTo",
117  (Circle(Circle::*)(UnitVector3d const &) const) &
118  Circle::expandedTo);
119  cls.def("expandedTo",
120  (Circle(Circle::*)(Circle const &) const) & Circle::expandedTo);
125  cls.def("getArea", &Circle::getArea);
126  cls.def("complement", &Circle::complement);
127  cls.def("complemented", &Circle::complemented);
128
129  // Note that the Region interface has already been wrapped.
130
131  // The lambda is necessary for now; returning the unique pointer
132  // directly leads to incorrect results and crashes.
133  cls.def_static("decode",
134  [](py::bytes bytes) { return decode(bytes).release(); },
135  "bytes"_a);
136
137  cls.def("__str__", [](Circle const &self) {
138  return py::str("Circle({!s}, {!s})")
139  .format(self.getCenter(), self.getOpeningAngle());
140  });
141  cls.def("__repr__", [](Circle const &self) {
142  return py::str("Circle({!r}, {!r})")
143  .format(self.getCenter(), self.getOpeningAngle());
144  });
145  cls.def(py::pickle(
146  [](const Circle &self) { return python::encode(self); },
147  [](py::bytes bytes) { return decode(bytes).release(); }));
148 }

## ◆ defineClass() [7/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< ConvexPolygon, std::unique_ptr< ConvexPolygon >, Region > & cls )

Definition at line 53 of file _convexPolygon.cc.

54  {
55  cls.attr("TYPE_CODE") = py::int_(ConvexPolygon::TYPE_CODE);
56
57  cls.def_static("convexHull", &ConvexPolygon::convexHull, "points"_a);
58
59  cls.def(py::init<std::vector<UnitVector3d> const &>(), "points"_a);
60  // Do not wrap the two unsafe (3 and 4 vertex) constructors
61  cls.def(py::init<ConvexPolygon const &>(), "convexPolygon"_a);
62
63  cls.def("__eq__", &ConvexPolygon::operator==, py::is_operator());
64  cls.def("__ne__", &ConvexPolygon::operator!=, py::is_operator());
65
66  cls.def("getVertices", &ConvexPolygon::getVertices);
67  cls.def("getCentroid", &ConvexPolygon::getCentroid);
68
69  // Note that much of the Region interface has already been wrapped. Here are bits that have not:
70  cls.def("contains", py::overload_cast<UnitVector3d const &>(&ConvexPolygon::contains, py::const_));
71  cls.def("contains", py::overload_cast<Region const &>(&ConvexPolygon::contains, py::const_));
72  cls.def("isDisjointFrom", &ConvexPolygon::isDisjointFrom);
73  cls.def("intersects", &ConvexPolygon::intersects);
74  cls.def("isWithin", &ConvexPolygon::isWithin);
75
76  // The lambda is necessary for now; returning the unique pointer
77  // directly leads to incorrect results and crashes.
78  cls.def_static("decode",
79  [](py::bytes bytes) { return decode(bytes).release(); },
80  "bytes"_a);
81
82  cls.def("__repr__", [](ConvexPolygon const &self) {
83  return py::str("ConvexPolygon({!r})").format(self.getVertices());
84  });
85  cls.def(py::pickle(
86  [](const ConvexPolygon &self) { return python::encode(self); },
87  [](py::bytes bytes) { return decode(bytes).release(); }));
88 }

## ◆ defineClass() [8/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Ellipse, std::unique_ptr< Ellipse >, Region > & cls )

Definition at line 53 of file _ellipse.cc.

53  {
54  cls.attr("TYPE_CODE") = py::int_(Ellipse::TYPE_CODE);
55
56  cls.def_static("empty", &Ellipse::empty);
57  cls.def_static("full", &Ellipse::full);
58
59  cls.def(py::init<>());
60  cls.def(py::init<Circle const &>(), "circle"_a);
61  cls.def(py::init<UnitVector3d const &, Angle>(), "center"_a,
62  "angle"_a = Angle(0.0));
63  cls.def(py::init<UnitVector3d const &, UnitVector3d const &, Angle>(),
64  "focus1"_a, "focus2"_a, "alpha"_a);
65  cls.def(py::init<UnitVector3d const &, Angle, Angle, Angle>(), "center"_a,
66  "alpha"_a, "beta"_a, "orientation"_a);
67  cls.def(py::init<Ellipse const &>(), "ellipse"_a);
68
69  cls.def("__eq__", &Ellipse::operator==, py::is_operator());
70  cls.def("__ne__", &Ellipse::operator!=, py::is_operator());
71
72  cls.def("isEmpty", &Ellipse::isEmpty);
73  cls.def("isFull", &Ellipse::isFull);
74  cls.def("isGreatCircle", &Ellipse::isGreatCircle);
75  cls.def("isCircle", &Ellipse::isCircle);
76  cls.def("getTransformMatrix", &Ellipse::getTransformMatrix);
77  cls.def("getCenter", &Ellipse::getCenter);
78  cls.def("getF1", &Ellipse::getF1);
79  cls.def("getF2", &Ellipse::getF2);
80  cls.def("getAlpha", &Ellipse::getAlpha);
81  cls.def("getBeta", &Ellipse::getBeta);
82  cls.def("getGamma", &Ellipse::getGamma);
83  cls.def("complement", &Ellipse::complement);
84  cls.def("complemented", &Ellipse::complemented);
85
86  // Note that the Region interface has already been wrapped.
87
88  // The lambda is necessary for now; returning the unique pointer
89  // directly leads to incorrect results and crashes.
90  cls.def_static("decode",
91  [](py::bytes bytes) { return decode(bytes).release(); },
92  "bytes"_a);
93
94  cls.def("__str__", [](Ellipse const &self) {
95  return py::str("Ellipse({!s}, {!s}, {!s})")
96  .format(self.getF1(), self.getF2(), self.getAlpha());
97  });
98  cls.def("__repr__", [](Ellipse const &self) {
99  return py::str("Ellipse({!r}, {!r}, {!r})")
100  .format(self.getF1(), self.getF2(), self.getAlpha());
101  });
102  cls.def(py::pickle(
103  [](const Ellipse &self) { return python::encode(self); },
104  [](py::bytes bytes) { return decode(bytes).release(); }));
105 }

## ◆ defineClass() [9/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< HtmPixelization, Pixelization > & cls )

Definition at line 35 of file _htmPixelization.cc.

35  {
36  cls.attr("MAX_LEVEL") = py::int_(HtmPixelization::MAX_LEVEL);
37
38  cls.def_static("level", &HtmPixelization::level, "i"_a);
39  cls.def_static("triangle", &HtmPixelization::triangle, "i"_a);
40  cls.def_static("asString", &HtmPixelization::asString, "i"_a);
41
42  cls.def(py::init<int>(), "level"_a);
43  cls.def(py::init<HtmPixelization const &>(), "htmPixelization"_a);
44
45  cls.def("getLevel", &HtmPixelization::getLevel);
46
47  cls.def("__eq__",
48  [](HtmPixelization const &self, HtmPixelization const &other) {
49  return self.getLevel() == other.getLevel();
50  });
51  cls.def("__ne__",
52  [](HtmPixelization const &self, HtmPixelization const &other) {
53  return self.getLevel() != other.getLevel();
54  });
55  cls.def("__repr__", [](HtmPixelization const &self) {
56  return py::str("HtmPixelization({!s})").format(self.getLevel());
57  });
58  cls.def("__reduce__", [cls](HtmPixelization const &self) {
59  return py::make_tuple(cls, py::make_tuple(self.getLevel()));
60  });
61 }

## ◆ defineClass() [10/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Interval1d, std::shared_ptr< Interval1d >> & cls )

Definition at line 36 of file _interval1d.cc.

36  {
37  python::defineInterval<decltype(cls), Interval1d, double>(cls);
38
39  cls.def_static("empty", &Interval1d::empty);
40  cls.def_static("full", &Interval1d::full);
41
42  cls.def(py::init<>());
43  cls.def(py::init<double>(), "x"_a);
44  cls.def(py::init<double, double>(), "x"_a, "y"_a);
45  cls.def(py::init<Interval1d const &>(), "interval"_a);
46
47  cls.def("isFull", &Interval1d::isFull);
48
49  cls.def("__str__", [](Interval1d const &self) {
50  return py::str("[{!s}, {!s}]").format(self.getA(), self.getB());
51  });
52  cls.def("__repr__", [](Interval1d const &self) {
53  return py::str("Interval1d({!r}, {!r})")
54  .format(self.getA(), self.getB());
55  });
56 }

## ◆ defineClass() [11/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< LonLat, std::shared_ptr< LonLat >> & cls )

Definition at line 36 of file _lonLat.cc.

36  {
37  cls.def_static("fromDegrees", &LonLat::fromDegrees);
39  cls.def_static("latitudeOf", &LonLat::latitudeOf);
40  cls.def_static("longitudeOf", &LonLat::longitudeOf);
41
42  cls.def(py::init<>());
43  cls.def(py::init<LonLat const &>());
44  cls.def(py::init<NormalizedAngle, Angle>(), "lon"_a, "lat"_a);
45  cls.def(py::init<Vector3d const &>(), "vector"_a);
46
47  cls.def("__eq__", &LonLat::operator==, py::is_operator());
48  cls.def("__nq__", &LonLat::operator!=, py::is_operator());
49
50  cls.def("getLon", &LonLat::getLon);
51  cls.def("getLat", &LonLat::getLat);
52
53  cls.def("__len__", [](LonLat const &self) { return py::int_(2); });
54  cls.def("__getitem__", [](LonLat const &self, py::object key) {
55  auto t = py::make_tuple(self.getLon(), self.getLat());
56  return t.attr("__getitem__")(key);
57  });
58  cls.def("__iter__", [](LonLat const &self) {
59  auto t = py::make_tuple(self.getLon(), self.getLat());
60  return t.attr("__iter__")();
61  });
62
63  cls.def("__str__", [](LonLat const &self) {
64  return py::str("[{!s}, {!s}]")
66  });
67  cls.def("__repr__", [](LonLat const &self) {
70  });
71  cls.def("__reduce__", [cls](LonLat const &self) {
72  return py::make_tuple(cls,
73  py::make_tuple(self.getLon(), self.getLat()));
74  });
75 }

## ◆ defineClass() [12/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Matrix3d, std::shared_ptr< Matrix3d >> & cls )

Definition at line 42 of file _matrix3d.cc.

## ◆ defineClass() [13/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Mq3cPixelization, Pixelization > & cls )

Definition at line 35 of file _mq3cPixelization.cc.

35  {
36  cls.attr("MAX_LEVEL") = py::int_(Mq3cPixelization::MAX_LEVEL);
37
38  cls.def_static("level", &Mq3cPixelization::level);
40  cls.def_static("neighborhood", &Mq3cPixelization::neighborhood);
41  cls.def_static("asString", &Mq3cPixelization::asString);
42
43  cls.def(py::init<int>(), "level"_a);
44  cls.def(py::init<Mq3cPixelization const &>(), "mq3cPixelization"_a);
45
46  cls.def("getLevel", &Mq3cPixelization::getLevel);
47
48  cls.def("__eq__",
49  [](Mq3cPixelization const &self, Mq3cPixelization const &other) {
50  return self.getLevel() == other.getLevel();
51  });
52  cls.def("__ne__",
53  [](Mq3cPixelization const &self, Mq3cPixelization const &other) {
54  return self.getLevel() != other.getLevel();
55  });
56  cls.def("__repr__", [](Mq3cPixelization const &self) {
57  return py::str("Mq3cPixelization({!s})").format(self.getLevel());
58  });
59  cls.def("__reduce__", [cls](Mq3cPixelization const &self) {
60  return py::make_tuple(cls, py::make_tuple(self.getLevel()));
61  });
62 }

## ◆ defineClass() [14/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< NormalizedAngle > & cls )

Definition at line 37 of file _normalizedAngle.cc.

37  {
38  // Provide the equivalent of the NormalizedAngle to Angle C++ cast
39  // operator in Python
40  py::implicitly_convertible<NormalizedAngle, Angle>();
41
42  cls.def_static("nan", &NormalizedAngle::nan);
43  cls.def_static("fromDegrees", &NormalizedAngle::fromDegrees);
45  cls.def_static("between", &NormalizedAngle::between, "a"_a, "b"_a);
46  cls.def_static("center", &NormalizedAngle::center, "a"_a, "b"_a);
47
48  cls.def(py::init<>());
49  cls.def(py::init<NormalizedAngle const &>());
50  cls.def(py::init<Angle const &>());
52  cls.def(py::init<LonLat const &, LonLat const &>(), "a"_a, "b"_a);
53  cls.def(py::init<Vector3d const &, Vector3d const &>(), "a"_a, "b"_a);
54
55  cls.def("__eq__", &NormalizedAngle::operator==, py::is_operator());
56  cls.def("__ne__", &NormalizedAngle::operator!=, py::is_operator());
57  cls.def("__lt__", &NormalizedAngle::operator<, py::is_operator());
58  cls.def("__gt__", &NormalizedAngle::operator>, py::is_operator());
59  cls.def("__le__", &NormalizedAngle::operator<=, py::is_operator());
60  cls.def("__ge__", &NormalizedAngle::operator>=, py::is_operator());
61
62  cls.def("__neg__",
63  (Angle(NormalizedAngle::*)() const) & NormalizedAngle::operator-);
65  cls.def("__sub__",
66  (Angle(NormalizedAngle::*)(Angle const &) const) &
67  NormalizedAngle::operator-,
68  py::is_operator());
69  cls.def("__mul__", &NormalizedAngle::operator*, py::is_operator());
70  cls.def("__rmul__", &NormalizedAngle::operator*, py::is_operator());
71  cls.def("__truediv__",
72  (Angle(NormalizedAngle::*)(double) const) &
73  NormalizedAngle::operator/,
74  py::is_operator());
75  cls.def("__truediv__",
76  (double (NormalizedAngle::*)(Angle const &) const) &
77  NormalizedAngle::operator/,
78  py::is_operator());
79
80  cls.def("asDegrees", &NormalizedAngle::asDegrees);
82  cls.def("isNan", &NormalizedAngle::isNan);
83  cls.def("getAngleTo", &NormalizedAngle::getAngleTo);
84
85  cls.def("__str__", [](NormalizedAngle const &self) {
87  });
88  cls.def("__repr__", [](NormalizedAngle const &self) {
90  });
91
92  cls.def("__reduce__", [cls](NormalizedAngle const &self) {
94  });
95 }

## ◆ defineClass() [15/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< NormalizedAngleInterval, std::shared_ptr< NormalizedAngleInterval >> & cls )

Definition at line 36 of file _normalizedAngleInterval.cc.

37  {
38  python::defineInterval<decltype(cls), NormalizedAngleInterval,
39  NormalizedAngle>(cls);
40
41  cls.def_static("fromDegrees", &NormalizedAngleInterval::fromDegrees, "x"_a,
42  "y"_a);
44  "y"_a);
45  cls.def_static("empty", &NormalizedAngleInterval::empty);
46  cls.def_static("full", &NormalizedAngleInterval::full);
47
48  cls.def(py::init<>());
49  cls.def(py::init<Angle>(), "x"_a);
50  cls.def(py::init<NormalizedAngle>(), "x"_a);
51  cls.def(py::init<Angle, Angle>(), "x"_a, "y"_a);
52  cls.def(py::init<NormalizedAngle, NormalizedAngle>(), "x"_a, "y"_a);
53  cls.def(py::init<NormalizedAngleInterval const &>(), "angleInterval"_a);
54
55  cls.def("isEmpty", &NormalizedAngleInterval::isEmpty);
56  cls.def("isFull", &NormalizedAngleInterval::isFull);
57  cls.def("wraps", &NormalizedAngleInterval::wraps);
58
59  cls.def("__str__", [](NormalizedAngleInterval const &self) {
60  return py::str("[{!s}, {!s}]")
62  });
63  cls.def("__repr__", [](NormalizedAngleInterval const &self) {
65  " {!r})")
67  });
68 }

## ◆ defineClass() [16/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Pixelization > & cls )

Definition at line 37 of file _pixelization.cc.

37  {
38  cls.def("universe", &Pixelization::universe);
39  cls.def("pixel", &Pixelization::pixel, "i"_a);
40  cls.def("index", &Pixelization::index, "i"_a);
41  cls.def("toString", &Pixelization::toString, "i"_a);
42  cls.def("envelope", &Pixelization::envelope, "region"_a, "maxRanges"_a = 0);
43  cls.def("interior", &Pixelization::interior, "region"_a, "maxRanges"_a = 0);
44 }

## ◆ defineClass() [17/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Q3cPixelization, Pixelization > & cls )

Definition at line 35 of file _q3cPixelization.cc.

35  {
36  cls.attr("MAX_LEVEL") = py::int_(Q3cPixelization::MAX_LEVEL);
37
38  cls.def(py::init<int>(), "level"_a);
39  cls.def(py::init<Q3cPixelization const &>(), "q3cPixelization"_a);
40
41  cls.def("getLevel", &Q3cPixelization::getLevel);
43  cls.def("neighborhood", &Q3cPixelization::neighborhood);
44
45  cls.def("__eq__",
46  [](Q3cPixelization const &self, Q3cPixelization const &other) {
47  return self.getLevel() == other.getLevel();
48  });
49  cls.def("__ne__",
50  [](Q3cPixelization const &self, Q3cPixelization const &other) {
51  return self.getLevel() != other.getLevel();
52  });
53  cls.def("__repr__", [](Q3cPixelization const &self) {
54  return py::str("Q3cPixelization({!s})").format(self.getLevel());
55  });
56  cls.def("__reduce__", [cls](Q3cPixelization const &self) {
57  return py::make_tuple(cls, py::make_tuple(self.getLevel()));
58  });
59 }

## ◆ defineClass() [18/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< RangeSet, std::shared_ptr< RangeSet >> & cls )

Definition at line 87 of file _rangeSet.cc.

87  {
88  cls.def(py::init<>());
89  cls.def(py::init<uint64_t>(), "integer"_a);
90  cls.def(py::init([](uint64_t a, uint64_t b) {
91  return new RangeSet(a, b);
92  }),
93  "first"_a, "last"_a);
94  cls.def(py::init<RangeSet const &>(), "rangeSet"_a);
95  cls.def(py::init(
96  [](py::iterable iterable) {
97  return new RangeSet(makeRangeSet(iterable));
98  }),
99  "iterable"_a);
100  cls.def("__eq__", &RangeSet::operator==, py::is_operator());
101  cls.def("__ne__", &RangeSet::operator!=, py::is_operator());
102
103  cls.def("insert", (void (RangeSet::*)(uint64_t)) & RangeSet::insert,
104  "integer"_a);
105  cls.def("insert",
106  (void (RangeSet::*)(uint64_t, uint64_t)) & RangeSet::insert,
107  "first"_a, "last"_a);
108  cls.def("erase", (void (RangeSet::*)(uint64_t)) & RangeSet::erase,
109  "integer"_a);
110  cls.def("erase", (void (RangeSet::*)(uint64_t, uint64_t)) & RangeSet::erase,
111  "first"_a, "last"_a);
112
113  cls.def("complement", &RangeSet::complement);
114  cls.def("complemented", &RangeSet::complemented);
115  cls.def("intersection", &RangeSet::intersection, "rangeSet"_a);
116  // In C++, the set union function is named join because union is a keyword.
117  // Python does not suffer from the same restriction.
118  cls.def("union", &RangeSet::join, "rangeSet"_a);
119  cls.def("difference", &RangeSet::difference, "rangeSet"_a);
120  cls.def("symmetricDifference", &RangeSet::symmetricDifference,
121  "rangeSet"_a);
122  cls.def("__invert__", &RangeSet::operator~, py::is_operator());
123  cls.def("__and__", &RangeSet::operator&, py::is_operator());
124  cls.def("__or__", &RangeSet::operator|, py::is_operator());
125  cls.def("__sub__", &RangeSet::operator-, py::is_operator());
126  cls.def("__xor__", &RangeSet::operator^, py::is_operator());
127  cls.def("__iand__", &RangeSet::operator&=);
128  cls.def("__ior__", &RangeSet::operator|=);
129  cls.def("__isub__", &RangeSet::operator-=);
130  cls.def("__ixor__", &RangeSet::operator^=);
131
132  cls.def("__len__", &RangeSet::size);
133  cls.def("__getitem__", [](RangeSet const &self, py::int_ i) {
134  auto j = python::convertIndex(static_cast<ptrdiff_t>(self.size()), i);
135  return py::cast(self.begin()[j]);
136  });
137
138  cls.def("intersects",
139  (bool (RangeSet::*)(uint64_t) const) & RangeSet::intersects,
140  "integer"_a);
141  cls.def("intersects",
142  (bool (RangeSet::*)(uint64_t, uint64_t) const) &
143  RangeSet::intersects,
144  "first"_a, "last"_a);
145  cls.def("intersects",
146  (bool (RangeSet::*)(RangeSet const &) const) & RangeSet::intersects,
147  "rangeSet"_a);
148
149  cls.def("contains",
150  (bool (RangeSet::*)(uint64_t) const) & RangeSet::contains,
151  "integer"_a);
152  cls.def("contains",
153  (bool (RangeSet::*)(uint64_t, uint64_t) const) & RangeSet::contains,
154  "first"_a, "last"_a);
155  cls.def("contains",
156  (bool (RangeSet::*)(RangeSet const &) const) & RangeSet::contains,
157  "rangeSet"_a);
158  cls.def("__contains__",
159  (bool (RangeSet::*)(uint64_t) const) & RangeSet::contains,
160  "integer"_a, py::is_operator());
161  cls.def("__contains__",
162  (bool (RangeSet::*)(uint64_t, uint64_t) const) & RangeSet::contains,
163  "first"_a, "last"_a, py::is_operator());
164  cls.def("__contains__",
165  (bool (RangeSet::*)(RangeSet const &) const) & RangeSet::contains,
166  "rangeSet"_a, py::is_operator());
167
168  cls.def("isWithin",
169  (bool (RangeSet::*)(uint64_t) const) & RangeSet::isWithin,
170  "integer"_a);
171  cls.def("isWithin",
172  (bool (RangeSet::*)(uint64_t, uint64_t) const) & RangeSet::isWithin,
173  "first"_a, "last"_a);
174  cls.def("isWithin",
175  (bool (RangeSet::*)(RangeSet const &) const) & RangeSet::isWithin,
176  "rangeSet"_a);
177
178  cls.def("isDisjointFrom",
179  (bool (RangeSet::*)(uint64_t) const) & RangeSet::isDisjointFrom,
180  "integer"_a);
181  cls.def("isDisjointFrom",
182  (bool (RangeSet::*)(uint64_t, uint64_t) const) &
183  RangeSet::isDisjointFrom,
184  "first"_a, "last"_a);
185  cls.def("isDisjointFrom",
186  (bool (RangeSet::*)(RangeSet const &) const) &
187  RangeSet::isDisjointFrom,
188  "rangeSet"_a);
189
190  cls.def("simplify", &RangeSet::simplify, "n"_a);
191  cls.def("simplified", &RangeSet::simplified, "n"_a);
192  cls.def("scale", &RangeSet::scale, "factor"_a);
193  cls.def("scaled", &RangeSet::scaled, "factor"_a);
194  cls.def("fill", &RangeSet::fill);
195  cls.def("clear", &RangeSet::clear);
196  cls.def("empty", &RangeSet::empty);
197  cls.def("full", &RangeSet::full);
198  cls.def("size", &RangeSet::size);
199  cls.def("cardinality", &RangeSet::cardinality);
200  // max_size() and swap() are omitted. The former is a C++ container
201  // requirement, and the latter doesn't seem relevant to Python.
202  cls.def("isValid", &RangeSet::cardinality);
203  cls.def("ranges", &ranges);
204
205  cls.def("__str__",
206  [](RangeSet const &self) { return py::str(ranges(self)); });
207  cls.def("__repr__", [](RangeSet const &self) {
208  return py::str("RangeSet({!s})").format(ranges(self));
209  });
210
211  cls.def("__reduce__", [cls](RangeSet const &self) {
212  return py::make_tuple(cls, py::make_tuple(ranges(self)));
213  });
214 }

## ◆ defineClass() [19/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Region, std::unique_ptr< Region >> & cls )

Definition at line 52 of file _region.cc.

52  {
53  cls.def("clone", [](Region const &self) { return self.clone().release(); });
54  cls.def("getBoundingBox", &Region::getBoundingBox);
55  cls.def("getBoundingBox3d", &Region::getBoundingBox3d);
56  cls.def("getBoundingCircle", &Region::getBoundingCircle);
57  cls.def("contains", &Region::contains, "unitVector"_a);
58  cls.def("__contains__", &Region::contains, "unitVector"_a,
59  py::is_operator());
60  // The per-subclass relate() overloads are used to implement
61  // double-dispatch in C++, and are not needed in Python.
62  cls.def("relate",
63  (Relationship(Region::*)(Region const &) const) & Region::relate,
64  "region"_a);
65  cls.def("encode", &python::encode);
66
67  cls.def_static(
68  "decode",
69  [](py::bytes bytes) {
70  uint8_t const *buffer = reinterpret_cast<uint8_t const *>(
71  PYBIND11_BYTES_AS_STRING(bytes.ptr()));
72  size_t n =
73  static_cast<size_t>(PYBIND11_BYTES_SIZE(bytes.ptr()));
74  return Region::decode(buffer, n).release();
75  },
76  "bytes"_a);
77 }

## ◆ defineClass() [20/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< UnitVector3d, std::shared_ptr< UnitVector3d >> & cls )

Definition at line 39 of file _unitVector3d.cc.

39  {
40  // Provide the equivalent of the UnitVector3d to Vector3d C++ cast
41  // operator in Python
42  py::implicitly_convertible<UnitVector3d, Vector3d>();
43
44  cls.def_static(
45  "orthogonalTo",
46  (UnitVector3d(*)(Vector3d const &)) & UnitVector3d::orthogonalTo,
47  "vector"_a);
48  cls.def_static("orthogonalTo",
49  (UnitVector3d(*)(Vector3d const &, Vector3d const &)) &
50  UnitVector3d::orthogonalTo,
51  "vector1"_a, "vector2"_a);
52  cls.def_static("orthogonalTo",
53  (UnitVector3d(*)(NormalizedAngle const &)) &
54  UnitVector3d::orthogonalTo,
55  "meridian"_a);
56  cls.def_static("northFrom", &UnitVector3d::northFrom, "vector"_a);
57  cls.def_static("X", &UnitVector3d::X);
58  cls.def_static("Y", &UnitVector3d::Y);
59  cls.def_static("Z", &UnitVector3d::Z);
60  // The fromNormalized static factory functions are not exposed to
61  // Python, as they are easy to misuse and intended only for performance
62  // critical code (i.e. not Python).
63
64  cls.def(py::init<>());
65  cls.def(py::init<UnitVector3d const &>(), "unitVector"_a);
66  cls.def(py::init<Vector3d const &>(), "vector"_a);
67  cls.def(py::init<double, double, double>(), "x"_a, "y"_a, "z"_a);
68  cls.def(py::init<LonLat const &>(), "lonLat"_a);
69  cls.def(py::init<Angle, Angle>(), "lon"_a, "lat"_a);
70
71  cls.def("__eq__", &UnitVector3d::operator==, py::is_operator());
72  cls.def("__ne__", &UnitVector3d::operator!=, py::is_operator());
73  cls.def("__neg__",
74  (UnitVector3d(UnitVector3d::*)() const) & UnitVector3d::operator-);
76  cls.def("__sub__",
77  (Vector3d(UnitVector3d::*)(Vector3d const &) const) &
78  UnitVector3d::operator-,
79  py::is_operator());
80  cls.def("__mul__", &UnitVector3d::operator*, py::is_operator());
81  cls.def("__truediv__", &UnitVector3d::operator/, py::is_operator());
82
83  cls.def("x", &UnitVector3d::x);
84  cls.def("y", &UnitVector3d::y);
85  cls.def("z", &UnitVector3d::z);
86  cls.def("x", &UnitVector3d::dot);
87  cls.def("dot", &UnitVector3d::dot);
88  cls.def("cross", &UnitVector3d::cross);
89  cls.def("robustCross", &UnitVector3d::robustCross);
90  cls.def("cwiseProduct", &UnitVector3d::cwiseProduct);
91  cls.def("rotatedAround", &UnitVector3d::rotatedAround, "axis"_a, "angle"_a);
92
93  cls.def("__len__", [](UnitVector3d const &self) { return py::int_(3); });
94  cls.def("__getitem__", [](UnitVector3d const &self, py::int_ i) {
95  return self(python::convertIndex(3, i));
96  });
97
98  cls.def("__str__", [](UnitVector3d const &self) {
99  return py::str("[{!s}, {!s}, {!s}]")
100  .format(self.x(), self.y(), self.z());
101  });
102  cls.def("__repr__", [](UnitVector3d const &self) {
103  return py::str("UnitVector3d({!r}, {!r}, {!r})")
104  .format(self.x(), self.y(), self.z());
105  });
106
107  // Do not implement __reduce__ for pickling. Why? Given:
108  //
109  // u = UnitVector3d(x, y, z)
110  // v = UnitVector3d(u.x(), u.y(), u.z())
111  //
112  // u may not be identical to v, since the constructor normalizes its input
113  // components. Furthermore, UnitVector3d::fromNormalized is not visible to
114  // Python, and even if it were, pybind11 is currently incapable of returning
115  // a picklable reference to it.
116  cls.def(py::pickle([](UnitVector3d const &self) { return py::make_tuple(self.x(), self.y(), self.z()); },
117  [](py::tuple t) {
118  if (t.size() != 3) {
119  throw std::runtime_error("Tuple size = " + std::to_string(t.size()) +
120  "; must be 3 for a UnitVector3d");
121  }
122  return new UnitVector3d(UnitVector3d::fromNormalized(
123  t[0].cast<double>(), t[1].cast<double>(), t[2].cast<double>()));
124  }));
125 }

## ◆ defineClass() [21/22]

template<>
 void lsst::sphgeom::defineClass ( py::class_< Vector3d, std::shared_ptr< Vector3d >> & cls )

Definition at line 38 of file _vector3d.cc.

38  {
39  cls.def(py::init<>());
40  cls.def(py::init<double, double, double>(), "x"_a, "y"_a, "z"_a);
41  cls.def(py::init<Vector3d const &>(), "vector"_a);
42  // Construct a Vector3d from a UnitVector3d, enabling implicit
43  // conversion from UnitVector3d to Vector3d in python via
44  // py::implicitly_convertible
45  cls.def(py::init([](UnitVector3d const &u) {
46  return new Vector3d(u.x(), u.y(), u.z());
47  }));
48
49  cls.def("__eq__", &Vector3d::operator==, py::is_operator());
50  cls.def("__ne__", &Vector3d::operator!=, py::is_operator());
51  cls.def("__neg__", (Vector3d(Vector3d::*)() const) & Vector3d::operator-);
53  cls.def("__sub__",
54  (Vector3d(Vector3d::*)(Vector3d const &) const) &
55  Vector3d::operator-,
56  py::is_operator());
57  cls.def("__mul__", &Vector3d::operator*, py::is_operator());
58  cls.def("__truediv__", &Vector3d::operator/, py::is_operator());
59
61  cls.def("__isub__", &Vector3d::operator-=);
62  cls.def("__imul__", &Vector3d::operator*=);
63  cls.def("__itruediv__", &Vector3d::operator/=);
64
65  cls.def("x", &Vector3d::x);
66  cls.def("y", &Vector3d::y);
67  cls.def("z", &Vector3d::z);
68  cls.def("dot", &Vector3d::dot);
69  cls.def("getSquaredNorm", &Vector3d::getSquaredNorm);
70  cls.def("getNorm", &Vector3d::getNorm);
71  cls.def("isZero", &Vector3d::isZero);
72  cls.def("normalize", &Vector3d::normalize);
73  cls.def("isNormalized", &Vector3d::isNormalized);
74  cls.def("cross", &Vector3d::cross);
75  cls.def("cwiseProduct", &Vector3d::cwiseProduct);
76  cls.def("rotatedAround", &Vector3d::rotatedAround, "axis"_a, "angle"_a);
77
78  cls.def("__len__", [](Vector3d const &self) { return py::int_(3); });
79  cls.def("__getitem__", [](Vector3d const &self, py::int_ i) {
80  return self(python::convertIndex(3, i));
81  });
82
83  cls.def("__str__", [](Vector3d const &self) {
84  return py::str("[{!s}, {!s}, {!s}]")
85  .format(self.x(), self.y(), self.z());
86  });
87  cls.def("__repr__", [](Vector3d const &self) {
88  return py::str("Vector3d({!r}, {!r}, {!r})")
89  .format(self.x(), self.y(), self.z());
90  });
91
92  cls.def("__reduce__", [cls](Vector3d const &self) {
93  return py::make_tuple(cls,
94  py::make_tuple(self.x(), self.y(), self.z()));
95  });
96 }

## ◆ defineClass() [22/22]

template<typename Pybind11Class >
 void lsst::sphgeom::defineClass ( Pybind11Class & cls )

## ◆ defineCurve()

 void lsst::sphgeom::defineCurve ( py::module & mod )

Definition at line 32 of file _curve.cc.

32  {
33  mod.def("log2", (uint8_t(*)(uint64_t)) & log2);
34  mod.def("mortonIndex", (uint64_t(*)(uint32_t, uint32_t)) & mortonIndex,
35  "x"_a, "y"_a);
36  mod.def("mortonIndexInverse",
38  "z"_a);
39  mod.def("mortonToHilbert", &mortonToHilbert, "z"_a, "m"_a);
40  mod.def("hilbertToMorton", &hilbertToMorton, "h"_a, "m"_a);
41  mod.def("hilbertIndex",
42  (uint64_t(*)(uint32_t, uint32_t, int)) & hilbertIndex, "x"_a, "y"_a,
43  "m"_a);
44  mod.def("hilbertIndexInverse",
45  (std::tuple<uint32_t, uint32_t>(*)(uint64_t, int)) &
47  "h"_a, "m"_a);
48 }

## ◆ defineOrientation()

 void lsst::sphgeom::defineOrientation ( py::module & mod )

Definition at line 32 of file _orientation.cc.

32  {
33  mod.def("orientationExact", &orientationExact, "a"_a, "b"_a, "c"_a);
34  mod.def("orientation", &orientation, "a"_a, "b"_a, "c"_a);
35  mod.def("orientationX", &orientationX, "b"_a, "c"_a);
36  mod.def("orientationY", &orientationY, "b"_a, "c"_a);
37  mod.def("orientationZ", &orientationZ, "b"_a, "c"_a);
38 }

## ◆ defineRelationship()

 void lsst::sphgeom::defineRelationship ( py::module & mod )

Definition at line 32 of file _relationship.cc.

32  {
33  mod.attr("DISJOINT") = py::cast(DISJOINT.to_ulong());
34  mod.attr("INTERSECTS") = py::cast(INTERSECTS.to_ulong());
35  mod.attr("CONTAINS") = py::cast(CONTAINS.to_ulong());
36  mod.attr("WITHIN") = py::cast(WITHIN.to_ulong());
37
38  mod.def("invert", &invert, "relationship"_a);
39 }

## ◆ defineUtils()

 void lsst::sphgeom::defineUtils ( py::module & mod )

Definition at line 37 of file _utils.cc.

37  {
38  mod.def("getMinSquaredChordLength", &getMinSquaredChordLength, "v"_a, "a"_a,
39  "b"_a, "n"_a);
40  mod.def("getMaxSquaredChordLength", &getMaxSquaredChordLength, "v"_a, "a"_a,
41  "b"_a, "n"_a);
42  mod.def("getMinAngleToCircle", &getMinAngleToCircle, "x"_a, "c"_a);
43  mod.def("getMaxAngleToCircle", &getMaxAngleToCircle, "x"_a, "c"_a);
44  mod.def("getWeightedCentroid", &getWeightedCentroid, "vector0"_a,
45  "vector1"_a, "vector2"_a);
46 }

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

38  {
39 #if defined(__x86_64__)
40  // x86-64 is little endian.
41  auto ptr = reinterpret_cast<uint8_t const *>(&item);
42  buffer.insert(buffer.end(), ptr, ptr + 8);
43 #else
44  union { uint64_t u; double d };
45  d = item;
46  buffer.push_back(static_cast<uint8_t>(value));
47  buffer.push_back(static_cast<uint8_t>(value >> 8));
48  buffer.push_back(static_cast<uint8_t>(value >> 16));
49  buffer.push_back(static_cast<uint8_t>(value >> 24));
50  buffer.push_back(static_cast<uint8_t>(value >> 32));
51  buffer.push_back(static_cast<uint8_t>(value >> 40));
52  buffer.push_back(static_cast<uint8_t>(value >> 48));
53  buffer.push_back(static_cast<uint8_t>(value >> 56));
54 #endif
55 }

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

68  {
70  if (abs(x) <= abs(c)) {
71  return a + Angle(PI) - 2.0 * abs(c);
72  }
73  if (a < abs(x)) {
74  return Angle(PI) - 2.0 * abs(c) - a;
75  }
76  return Angle(PI) + 2.0 * abs(c) - a;
77 }

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

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

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

62  {
63  return abs(x - c);
64 }

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

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

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

82 {
83  // For the details, see:
84  //
85  // The centroid and inertia tensor for a spherical triangle
86  // John E. Brock
87  // 1974, Naval Postgraduate School, Monterey Calif.
88  //
89  // https://openlibrary.org/books/OL25493734M/The_centroid_and_inertia_tensor_for_a_spherical_triangle
90
91  Vector3d x01 = v0.robustCross(v1); // twice the cross product of v0 and v1
92  Vector3d x12 = v1.robustCross(v2);
93  Vector3d x20 = v2.robustCross(v0);
94  double s01 = 0.5 * x01.normalize(); // sine of the angle between v0 and v1
95  double s12 = 0.5 * x12.normalize();
96  double s20 = 0.5 * x20.normalize();
97  double c01 = v0.dot(v1); // cosine of the angle between v0 and v1
98  double c12 = v1.dot(v2);
99  double c20 = v2.dot(v0);
100  double a0 = (s12 == 0.0 && c12 == 0.0) ? 0.0 : std::atan2(s12, c12);
101  double a1 = (s20 == 0.0 && c20 == 0.0) ? 0.0 : std::atan2(s20, c20);
102  double a2 = (s01 == 0.0 && c01 == 0.0) ? 0.0 : std::atan2(s01, c01);
103  return 0.5 * (x01 * a2 + x12 * a0 + x20 * a1);
104 }

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

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

## ◆ hilbertIndexInverse()

 std::tuple 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.

361  {
363 }

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

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

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

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

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

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

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

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

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

148  {
149  // This is just a 64-bit extension of:
150  // http://graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN
151  uint64_t b = y;
152  uint64_t a = x;
153  b = (b | (b << 16)) & UINT64_C(0x0000ffff0000ffff);
154  a = (a | (a << 16)) & UINT64_C(0x0000ffff0000ffff);
155  b = (b | (b << 8)) & UINT64_C(0x00ff00ff00ff00ff);
156  a = (a | (a << 8)) & UINT64_C(0x00ff00ff00ff00ff);
157  b = (b | (b << 4)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
158  a = (a | (a << 4)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
159  b = (b | (b << 2)) & UINT64_C(0x3333333333333333);
160  a = (a | (a << 2)) & UINT64_C(0x3333333333333333);
161  b = (b | (b << 1)) & UINT64_C(0x5555555555555555);
162  a = (a | (a << 1)) & UINT64_C(0x5555555555555555);
163  return a | (b << 1);
164  }

## ◆ mortonIndexInverse()

 std::tuple 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.

195  {
196  uint64_t x = z & UINT64_C(0x5555555555555555);
197  uint64_t y = (z >> 1) & UINT64_C(0x5555555555555555);
198  x = (x | (x >> 1)) & UINT64_C(0x3333333333333333);
199  y = (y | (y >> 1)) & UINT64_C(0x3333333333333333);
200  x = (x | (x >> 2)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
201  y = (y | (y >> 2)) & UINT64_C(0x0f0f0f0f0f0f0f0f);
202  x = (x | (x >> 4)) & UINT64_C(0x00ff00ff00ff00ff);
203  y = (y | (y >> 4)) & UINT64_C(0x00ff00ff00ff00ff);
204  x = (x | (x >> 8)) & UINT64_C(0x0000ffff0000ffff);
205  y = (y | (y >> 8)) & UINT64_C(0x0000ffff0000ffff);
206  return std::make_tuple(static_cast<uint32_t>(x | (x >> 16)),
207  static_cast<uint32_t>(y | (y >> 16)));
208  }

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

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

## ◆ operator*() [1/2]

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

Definition at line 98 of file Angle.h.

98 { return b * a; }

## ◆ operator*() [2/2]

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

Definition at line 165 of file Vector3d.h.

165 { 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.

34  {
35  char buf[32];
37  return os << buf;
38 }

## ◆ operator<<() [2/14]

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

Definition at line 34 of file AngleInterval.cc.

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

## ◆ operator<<() [3/14]

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

Definition at line 475 of file Box.cc.

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

## ◆ operator<<() [4/14]

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

Definition at line 34 of file Box3d.cc.

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

## ◆ operator<<() [5/14]

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

Definition at line 354 of file Circle.cc.

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

## ◆ operator<<() [6/14]

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

Definition at line 382 of file ConvexPolygon.cc.

382  {
383  typedef std::vector<UnitVector3d>::const_iterator VertexIterator;
384  VertexIterator v = p.getVertices().begin();
385  VertexIterator const end = p.getVertices().end();
386  os << "{\"ConvexPolygon\": [" << *v;
387  for (++v; v != end; ++v) { os << ", " << *v; }
388  os << "]}";
389  return os;
390 }

## ◆ operator<<() [7/14]

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

Definition at line 388 of file Ellipse.cc.

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

## ◆ operator<<() [8/14]

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

Definition at line 35 of file Interval1d.cc.

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

## ◆ operator<<() [9/14]

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

Definition at line 82 of file LonLat.cc.

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

## ◆ operator<<() [10/14]

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

Definition at line 35 of file Matrix3d.cc.

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

## ◆ operator<<() [11/14]

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

Definition at line 284 of file NormalizedAngleInterval.cc.

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

## ◆ operator<<() [12/14]

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

Definition at line 616 of file RangeSet.cc.

616  {
617  os << "{\"RangeSet\": [";
618  bool first = true;
619  for (auto const & t: s) {
620  if (!first) {
621  os << ", ";
622  }
623  first = false;
624  os << '[' << std::get<0>(t) << ", " << std::get<1>(t) << ']';
625  }
626  os << "]}";
627  return os;
628 }

## ◆ operator<<() [13/14]

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

Definition at line 73 of file UnitVector3d.cc.

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

## ◆ operator<<() [14/14]

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

Definition at line 133 of file Vector3d.cc.

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

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

138 {
139  // This constant is a little more than 5ε, where ε = 2^-53. When multiplied
140  // by the permanent of |M|, it gives an error bound on the determinant of
141  // M. Here, M is a 3x3 matrix and |M| denotes the matrix obtained by
142  // taking the absolute value of each of its components. The derivation of
143  // this proceeds in the same manner as the derivation of the error bounds
144  // in section 4.3 of:
145  //
146  // Adaptive Precision Floating-Point Arithmetic
147  // and Fast Robust Geometric Predicates,
148  // Jonathan Richard Shewchuk,
149  // Discrete & Computational Geometry 18(3):305–363, October 1997.
150  //
151  // available online at http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf
152  static double const relativeError = 5.6e-16;
153  // Because all 3 unit vectors are normalized, the maximum absolute value of
154  // any vector component, cross product component or dot product term in
155  // the calculation is very close to 1. The permanent of |M| must therefore
156  // be below 3 + c, where c is some small multiple of ε. This constant, a
157  // little larger than 3 * 5ε, is an upper bound on the absolute error in
158  // the determinant calculation.
159  static double const maxAbsoluteError = 1.7e-15;
160  // This constant accounts for floating point underflow (assuming hardware
161  // without gradual underflow, just to be conservative) in the computation
162  // of det(M). It is a little more than 14 * 2^-1022.
163  static double const minAbsoluteError = 4.0e-307;
164
165  double bycz = b.y() * c.z();
166  double bzcy = b.z() * c.y();
167  double bzcx = b.z() * c.x();
168  double bxcz = b.x() * c.z();
169  double bxcy = b.x() * c.y();
170  double bycx = b.y() * c.x();
171  double determinant = a.x() * (bycz - bzcy) +
172  a.y() * (bzcx - bxcz) +
173  a.z() * (bxcy - bycx);
174  if (determinant > maxAbsoluteError) {
175  return 1;
176  } else if (determinant < -maxAbsoluteError) {
177  return -1;
178  }
179  // Expend some more effort on what is hopefully a tighter error bound
180  // before falling back on arbitrary precision arithmetic.
181  double permanent = std::fabs(a.x()) * (std::fabs(bycz) + std::fabs(bzcy)) +
182  std::fabs(a.y()) * (std::fabs(bzcx) + std::fabs(bxcz)) +
183  std::fabs(a.z()) * (std::fabs(bxcy) + std::fabs(bycx));
184  double maxError = relativeError * permanent + minAbsoluteError;
185  if (determinant > maxError) {
186  return 1;
187  } else if (determinant < -maxError) {
188  return -1;
189  }
190  // Avoid the slow path when any two inputs are identical or antipodal.
191  if (a == b || b == c || a == c || a == -b || b == -c || a == -c) {
192  return 0;
193  }
194  return orientationExact(a, b, c);
195 }

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

79 {
80  // Product mantissa storage buffers.
81  uint32_t mantissaBuffers[6][6];
82  // Product mantissas.
83  BigInteger mantissas[6] = {
84  BigInteger(mantissaBuffers[0], 6),
85  BigInteger(mantissaBuffers[1], 6),
86  BigInteger(mantissaBuffers[2], 6),
87  BigInteger(mantissaBuffers[3], 6),
88  BigInteger(mantissaBuffers[4], 6),
89  BigInteger(mantissaBuffers[5], 6)
90  };
91  BigFloat products[6] = {
92  BigFloat(&mantissas[0]),
93  BigFloat(&mantissas[1]),
94  BigFloat(&mantissas[2]),
95  BigFloat(&mantissas[3]),
96  BigFloat(&mantissas[4]),
97  BigFloat(&mantissas[5])
98  };
99  // An accumulator and its storage.
100  uint32_t accumulatorBuffer[512];
101  BigInteger accumulator(accumulatorBuffer,
102  sizeof(accumulatorBuffer) / sizeof(uint32_t));
103  // Compute the products in the determinant. Performing all multiplication
104  // up front means that each product mantissa occupies at most 3*53 bits.
105  computeProduct(products[0], a.x(), b.y(), c.z());
106  computeProduct(products[1], a.x(), b.z(), c.y());
107  computeProduct(products[2], a.y(), b.z(), c.x());
108  computeProduct(products[3], a.y(), b.x(), c.z());
109  computeProduct(products[4], a.z(), b.x(), c.y());
110  computeProduct(products[5], a.z(), b.y(), c.x());
111  mantissas[1].negate();
112  mantissas[3].negate();
113  mantissas[5].negate();
114  // Sort the array of products in descending exponent order.
115  std::sort(products, products + 6, [](BigFloat const & a, BigFloat const & b) {
116  return a.exponent > b.exponent;
117  });
118  // First, initialize the accumulator to the product with the highest
119  // exponent, then add the remaining products. Prior to each addition, we
120  // must shift the accumulated value so that its radix point lines up with
122  //
123  // More precisely, at each step we have an accumulated value A·2ʲ and a
124  // product P·2ᵏ, and we update the accumulator to equal (A·2ʲ⁻ᵏ + P)·2ᵏ.
125  // Because the products were sorted beforehand, j ≥ k and 2ʲ⁻ᵏ is an
126  // integer.
127  accumulator = *products[0].mantissa;
128  for (int i = 1; i < 6; ++i) {
129  accumulator.multiplyPow2(products[i - 1].exponent - products[i].exponent);
131  }
132  return accumulator.getSign();
133 }

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

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

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

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

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

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

## ◆ sin()

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

Definition at line 102 of file Angle.h.

## ◆ swap()

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

Definition at line 608 of file RangeSet.h.

608  {
609  a.swap(b);
610 }

## ◆ tan()

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

Definition at line 104 of file Angle.h.

## Variable Documentation

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.

constexpr

Definition at line 38 of file constants.h.

y
int y
Definition: SpanSet.cc:49
std::tan
T tan(T... args)
std::make_tuple
T make_tuple(T... args)
lsst::geom::PI
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:39
std::fabs
T fabs(T... args)
std::atan2
T atan2(T... args)
std::cos
T cos(T... args)
lsst::sphgeom::orientation
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
std::vector
STL class.
lsst::sphgeom::abs
Angle abs(Angle const &a)
Definition: Angle.h:106
lsst::afw::table._match.first
first
Definition: _match.py:76
std::tuple
lsst::log.log.logContinued.getLevel
def getLevel(loggername)
Definition: logContinued.py:178
lsst::sphgeom::log2
uint8_t log2(uint64_t x)
Definition: curve.h:98
lsst::afw::image::X
@ X
Definition: ImageUtils.h:36
std::sort
T sort(T... args)
end
int end
Definition: BoundedField.cc:105
std::vector::push_back
T push_back(T... args)
lsst::afw::table::Angle
lsst::geom::Angle Angle
Definition: misc.h:33
lsst::afw.display.ds9.dot
def dot(symb, c, r, frame=None, size=2, ctype=None, origin=afwImage.PARENT, *args, **kwargs)
Definition: ds9.py:101
lsst::geom::polynomials::simplified
PolynomialFunction1d simplified(ScaledPolynomialFunction1d const &f)
Calculate the standard polynomial function that is equivalent to a scaled standard polynomial functio...
Definition: PolynomialFunction1d.cc:32
std::snprintf
T snprintf(T... args)
astshim.fitsChanContinued.contains
def contains(self, name)
Definition: fitsChanContinued.py:127
z
double z
Definition: Match.cc:44
x
double x
Definition: ChebyshevBoundedField.cc:277
other
ItemVariant const * other
Definition: Schema.cc:56
std::to_string
T to_string(T... args)
lsst::sphgeom::orientationX
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
Definition: orientation.cc:227
std::atan
T atan(T... args)
lsst::sphgeom::getMinAngleToCircle
Angle getMinAngleToCircle(Angle x, Angle c)
getMinAngleToCircle returns the minimum angular separation between a point at latitude x and the poin...
Definition: utils.h:62
lsst::sphgeom::detail::relate
Relationship relate(VertexIterator const begin, VertexIterator const end, Box const &b)
Definition: ConvexPolygonImpl.h:258
lsst::afw::image::Y
@ Y
Definition: ImageUtils.h:36
std::runtime_error
STL class.
ptr
uint64_t * ptr
Definition: RangeSet.cc:88
lsst::sphgeom::hilbertIndexInverse
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.
Definition: curve.h:361
lsst::ip::isr::between
std::string between(std::string &s, char ldelim, char rdelim)
Definition: Isr.cc:100
pixel
table::PointKey< int > pixel
Definition: DeltaFunctionKernel.cc:101
lsst::sphgeom::getMaxSquaredChordLength
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...
Definition: utils.cc:58
lsst::sphgeom::hilbertIndex
uint64_t hilbertIndex(uint32_t x, uint32_t y, int m)
hilbertIndex returns the index of (x, y) in a 2-D Hilbert curve.
Definition: curve.h:349
lsst::sphgeom::Relationship
std::bitset< 3 > Relationship
Relationship describes how two sets are related.
Definition: Relationship.h:35
mantissa
BigInteger * mantissa
Definition: orientation.cc:40
b
table::Key< int > b
Definition: TransmissionCurve.cc:467
lsst::sphgeom::mortonIndexInverse
std::tuple< uint32_t, uint32_t > mortonIndexInverse(uint64_t z)
mortonIndexInverse separates the even and odd bits of z.
Definition: curve.h:195
lsst::sphgeom::mortonToHilbert
uint64_t mortonToHilbert(uint64_t z, int m)
mortonToHilbert converts the 2m-bit Morton index z to the corresponding Hilbert index.
Definition: curve.h:236
lsst::sphgeom::orientationY
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
Definition: orientation.cc:232
std::sin
T sin(T... args)
lsst::sphgeom::hilbertToMorton
uint64_t hilbertToMorton(uint64_t h, int m)
hilbertToMorton converts the 2m-bit Hilbert index h to the corresponding Morton index.
Definition: curve.h:290
exponent
int exponent
Definition: orientation.cc:41
lsst::sphgeom::getMaxAngleToCircle
Angle getMaxAngleToCircle(Angle x, Angle c)
getMaxAngleToCircle returns the maximum angular separation between a point at latitude x and the poin...
Definition: utils.h:68
os
std::ostream * os
Definition: Schema.cc:746
bytes
table::Key< table::Array< std::uint8_t > > bytes
Definition: Transform.cc:199
row
int row
Definition: CR.cc:145
list
daf::base::PropertyList * list
Definition: fits.cc:913
a
table::Key< int > a
Definition: TransmissionCurve.cc:466
std::begin
T begin(T... args)
lsst::sphgeom::orientationExact
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 nor...
Definition: orientation.cc:76
std::vector::insert
T insert(T... args)
key
Key< U > key
Definition: Schema.cc:281
lsst::afw.display.ds9.scale
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:109
lsst::sphgeom::orientationZ
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
Definition: orientation.cc:237
std::vector::end
T end(T... args)
lsst::sphgeom::getMinSquaredChordLength
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 cr...
Definition: utils.cc:36
lsst::utils.tests.init
def init()
Definition: tests.py:59
lsst::afw.display.ds9.erase
def erase(frame=None)
Definition: ds9.py:97
lsst::sphgeom::mortonIndex
uint64_t mortonIndex(uint32_t x, uint32_t y)
mortonIndex interleaves the bits of x and y.
Definition: curve.h:148
lsst::sphgeom::getWeightedCentroid
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 ma...
Definition: utils.cc:79
lsst::sphgeom::invert
Relationship invert(Relationship r)
Given the relationship between two sets A and B (i.e.
Definition: Relationship.h:55
m
int m
Definition: SpanSet.cc:49