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

Classes

class  PixelFinder
 

Functions

template<typename VertexIterator >
UnitVector3d centroid (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
Circle boundingCircle (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
Box boundingBox (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
Box3d boundingBox3d (VertexIterator const begin, VertexIterator const end)
 
template<typename VertexIterator >
bool contains (VertexIterator const begin, VertexIterator const end, UnitVector3d const &v)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, Box const &b)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, Circle const &c)
 
template<typename VertexIterator1 , typename VertexIterator2 >
Relationship relate (VertexIterator1 const begin1, VertexIterator1 const end1, VertexIterator2 const begin2, VertexIterator2 const end2)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, ConvexPolygon const &p)
 
template<typename VertexIterator >
Relationship relate (VertexIterator const begin, VertexIterator const end, Ellipse const &e)
 
template<template< typename, bool > class Finder, bool InteriorOnly>
RangeSet findPixels (Region const &r, size_t maxRanges, int level)
 

Function Documentation

◆ boundingBox()

template<typename VertexIterator >
Box lsst::sphgeom::detail::boundingBox ( VertexIterator const begin,
VertexIterator const end )

Definition at line 92 of file ConvexPolygonImpl.h.

92 {
93 Angle const eps(5.0e-10); // ~ 0.1 milli-arcseconds
94 Box bbox;
95 VertexIterator i = std::prev(end);
96 VertexIterator j = begin;
97 bool haveCW = false;
98 bool haveCCW = false;
99 // Compute the bounding box for each vertex. When converting a Vector3d
100 // to a LonLat, the relative error on the longitude is about 4*2^-53,
101 // and the relative error on the latitude is about twice that (assuming
102 // std::atan2 and std::sqrt accurate to within 1 ulp). We convert each
103 // vertex to a conservative bounding box for its spherical coordinates,
104 // and compute a bounding box for the union of all these boxes.
105 //
106 // Furthermore, the latitude range of an edge can be greater than the
107 // latitude range of its endpoints - this occurs when the minimum or
108 // maximum latitude point on the great circle defined by the edge vertices
109 // lies in the edge interior.
110 for (; j != end; i = j, ++j) {
111 LonLat p(*j);
112 bbox.expandTo(Box(p, eps, eps));
113 if (!haveCW || !haveCCW) {
114 int o = orientationZ(*i, *j);
115 haveCCW = haveCCW || (o > 0);
116 haveCW = haveCW || (o < 0);
117 }
118 // Compute the plane normal for edge i, j.
119 Vector3d n = (*i).robustCross(*j);
120 // Compute a vector v with positive z component that lies on both the
121 // edge plane and on the plane defined by the z axis and the edge plane
122 // normal. This is the direction of maximum latitude for the great
123 // circle containing the edge, and -v is the direction of minimum
124 // latitude.
125 //
126 // TODO(smm): Do a proper error analysis.
127 Vector3d v(-n.x() * n.z(),
128 -n.y() * n.z(),
129 n.x() * n.x() + n.y() * n.y());
130 if (v != Vector3d()) {
131 // The plane defined by the z axis and n has normal
132 // (-n.y(), n.x(), 0.0). Compute the dot product of this plane
133 // normal with vertices i and j.
134 double zni = i->y() * n.x() - i->x() * n.y();
135 double znj = j->y() * n.x() - j->x() * n.y();
136 // Check if v or -v is in the edge interior.
137 if (zni > 0.0 && znj < 0.0) {
138 bbox = Box(bbox.getLon(), bbox.getLat().expandedTo(
139 LonLat::latitudeOf(v) + eps));
140 } else if (zni < 0.0 && znj > 0.0) {
141 bbox = Box(bbox.getLon(), bbox.getLat().expandedTo(
142 LonLat::latitudeOf(-v) - eps));
143 }
144 }
145 }
146 // If this polygon contains a pole, its bounding box must contain all
147 // longitudes.
148 if (!haveCW) {
149 Box northPole(Box::allLongitudes(), AngleInterval(Angle(0.5 * PI)));
150 bbox.expandTo(northPole);
151 } else if (!haveCCW) {
152 Box southPole(Box::allLongitudes(), AngleInterval(Angle(-0.5 * PI)));
153 bbox.expandTo(southPole);
154 }
155 return bbox;
156}
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
int end
Angle represents an angle in radians.
Definition Angle.h:50
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition Box.h:61
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition LonLat.h:55
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition Vector3d.h:51
double x() const
Definition Vector3d.h:73
double y() const
Definition Vector3d.h:75
double z() const
Definition Vector3d.h:77
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
T prev(T... args)

◆ boundingBox3d()

template<typename VertexIterator >
Box3d lsst::sphgeom::detail::boundingBox3d ( VertexIterator const begin,
VertexIterator const end )

Definition at line 159 of file ConvexPolygonImpl.h.

159 {
160 static double const maxError = 1.0e-14;
161 // Compute the extrema of all vertex coordinates.
162 VertexIterator j = begin;
163 double emin[3] = { j->x(), j->y(), j->z() };
164 double emax[3] = { j->x(), j->y(), j->z() };
165 for (++j; j != end; ++j) {
166 for (int i = 0; i < 3; ++i) {
167 double v = j->operator()(i);
168 emin[i] = std::min(emin[i], v);
169 emax[i] = std::max(emax[i], v);
170 }
171 }
172 // Compute the extrema of all edges.
173 //
174 // It can be shown that the great circle with unit normal vector
175 // n = (n₀, n₁, n₂) has extrema in x at:
176 //
177 // (∓√(1 - n₀²), ±n₁n₀/√(1 - n₀²), ±n₂n₀/√(1 - n₀²))
178 //
179 // in y at:
180 //
181 // (±n₀n₁/√(1 - n₁²), ∓√(1 - n₁²), ±n₂n₁/√(1 - n₁²))
182 //
183 // and in z at
184 //
185 // (±n₀n₂/√(1 - n₂²), ±n₁n₂/√(1 - n₂²), ∓√(1 - n₂²))
186 //
187 // Compute these vectors for each edge, determine whether they lie in
188 // the edge, and update the extrema if so. Rounding errors in these
189 // computations are compensated for by expanding the bounding box
190 // prior to returning it.
191 j = std::prev(end);
192 VertexIterator k = begin;
193 for (; k != end; j = k, ++k) {
194 UnitVector3d n(j->robustCross(*k));
195 for (int i = 0; i < 3; ++i) {
196 double ni = n(i);
197 double d = std::fabs(1.0 - ni * ni);
198 if (d > 0.0) {
199 Vector3d e(i == 0 ? -d : n.x() * ni,
200 i == 1 ? -d : n.y() * ni,
201 i == 2 ? -d : n.z() * ni);
202 // If e or -e lies in the lune defined by the half great
203 // circle passing through n and a and the half great circle
204 // passing through n and b, the edge contains an extremum.
205 Vector3d v = e.cross(n);
206 double vdj = v.dot(*j);
207 double vdk = v.dot(*k);
208 if (vdj >= 0.0 && vdk <= 0.0) {
209 emin[i] = std::min(emin[i], -std::sqrt(d));
210 }
211 if (vdj <= 0.0 && vdk >= 0.0) {
212 emax[i] = std::max(emax[i], std::sqrt(d));
213 }
214 }
215 }
216 }
217 // Check whether the standard basis vectors and their antipodes
218 // are inside this polygon.
219 bool a[3] = { true, true, true };
220 bool b[3] = { true, true, true };
221 j = std::prev(end);
222 k = begin;
223 for (; k != end; j = k, ++k) {
224 // Test the standard basis vectors against the plane defined by
225 // vertices (j, k). Note that orientation(-x, *j, *k) =
226 // -orientation(x, *j, *k).
227 int ox = orientationX(*j, *k);
228 a[0] = a[0] && (ox <= 0);
229 b[0] = b[0] && (ox >= 0);
230 int oy = orientationY(*j, *k);
231 a[1] = a[1] && (oy <= 0);
232 b[1] = b[1] && (oy >= 0);
233 int oz = orientationZ(*j, *k);
234 a[2] = a[2] && (oz <= 0);
235 b[2] = b[2] && (oz >= 0);
236 }
237 // At this point, b[i] is true iff the standard basis vector eᵢ
238 // is inside all the half spaces defined by the polygon edges.
239 // Similarly, a[i] is true iff -eᵢ is inside the same half spaces.
240 for (int i = 0; i < 3; ++i) {
241 emin[i] = a[i] ? -1.0 : std::max(-1.0, emin[i] - maxError);
242 emax[i] = b[i] ? 1.0 : std::min(1.0, emax[i] + maxError);
243 }
244 return Box3d(Interval1d(emin[0], emax[0]),
245 Interval1d(emin[1], emax[1]),
246 Interval1d(emin[2], emax[2]));
247}
double z
Definition Match.cc:44
int y
Definition SpanSet.cc:48
table::Key< int > b
T begin(T... args)
T fabs(T... args)
T max(T... args)
T min(T... args)
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
T sqrt(T... args)

◆ boundingCircle()

template<typename VertexIterator >
Circle lsst::sphgeom::detail::boundingCircle ( VertexIterator const begin,
VertexIterator const end )

Definition at line 76 of file ConvexPolygonImpl.h.

76 {
77 UnitVector3d c = centroid(begin, end);
78 // Compute the maximum squared chord length between the centroid and
79 // all vertices.
80 VertexIterator i = begin;
81 double cl2 = 0.0;
82 for (; i != end; ++i) {
83 cl2 = std::max(cl2, (*i - c).getSquaredNorm());
84 }
85 // Add double the maximum squared-chord-length error, so that the
86 // bounding circle we return also reliably CONTAINS this polygon.
87 return Circle(c, cl2 + 2.0 * MAX_SQUARED_CHORD_LENGTH_ERROR);
88}
Circle is a circular region on the unit sphere that contains its boundary.
Definition Circle.h:53
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
UnitVector3d centroid(VertexIterator const begin, VertexIterator const end)

◆ centroid()

template<typename VertexIterator >
UnitVector3d lsst::sphgeom::detail::centroid ( VertexIterator const begin,
VertexIterator const end )

Definition at line 55 of file ConvexPolygonImpl.h.

55 {
56 // The center of mass is obtained via trivial generalization of
57 // the formula for spherical triangles from:
58 //
59 // The centroid and inertia tensor for a spherical triangle
60 // John E. Brock
61 // 1974, Naval Postgraduate School, Monterey Calif.
62 Vector3d cm;
63 VertexIterator i = std::prev(end);
64 VertexIterator j = begin;
65 for (; j != end; i = j, ++j) {
66 Vector3d v = (*i).robustCross(*j);
67 double s = 0.5 * v.normalize();
68 double c = (*i).dot(*j);
69 double a = (s == 0.0 && c == 0.0) ? 0.0 : std::atan2(s, c);
70 cm += v * a;
71 }
72 return UnitVector3d(cm);
73}
table::Key< int > a
double normalize()
normalize scales this vector to have unit norm and returns its norm prior to scaling.
Definition Vector3d.cc:48
STL namespace.

◆ contains()

template<typename VertexIterator >
bool lsst::sphgeom::detail::contains ( VertexIterator const begin,
VertexIterator const end,
UnitVector3d const & v )

Definition at line 250 of file ConvexPolygonImpl.h.

253{
254 VertexIterator i = std::prev(end);
255 VertexIterator j = begin;
256 for (; j != end; i = j, ++j) {
257 if (orientation(v, *i, *j) < 0) {
258 return false;
259 }
260 }
261 return true;
262}
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.

◆ findPixels()

template<template< typename, bool > class Finder, bool InteriorOnly>
RangeSet lsst::sphgeom::detail::findPixels ( Region const & r,
size_t maxRanges,
int level )

Definition at line 157 of file PixelFinder.h.

157 {
158 RangeSet s;
159 Circle const * c = nullptr;
160 Ellipse const * e = nullptr;
161 Box const * b = nullptr;
162 if ((c = dynamic_cast<Circle const *>(&r))) {
163 Finder<Circle, InteriorOnly> find(s, *c, level, maxRanges);
164 find();
165 } else if ((e = dynamic_cast<Ellipse const *>(&r))) {
166 Finder<Circle, InteriorOnly> find(
167 s, e->getBoundingCircle(), level, maxRanges);
168 find();
169 } else if ((b = dynamic_cast<Box const *>(&r))) {
170 Finder<Box, InteriorOnly> find(s, *b, level, maxRanges);
171 find();
172 } else {
173 Finder<ConvexPolygon, InteriorOnly> find(
174 s, dynamic_cast<ConvexPolygon const &>(r), level, maxRanges);
175 find();
176 }
177 return s;
178}
Ellipse is an elliptical region on the sphere.
Definition Ellipse.h:177
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Definition Ellipse.cc:248
A RangeSet is a set of unsigned 64 bit integers.
Definition RangeSet.h:106

◆ relate() [1/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const begin,
VertexIterator const end,
Box const & b )

Definition at line 265 of file ConvexPolygonImpl.h.

268{
269 // TODO(smm): be more accurate when computing box relations.
270 return boundingBox(begin, end).relate(b) & (DISJOINT | WITHIN);
271}
Relationship relate(LonLat const &p) const
Definition Box.h:304
Box boundingBox(VertexIterator const begin, VertexIterator const end)

◆ relate() [2/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const begin,
VertexIterator const end,
Circle const & c )

Definition at line 274 of file ConvexPolygonImpl.h.

277{
278 if (c.isEmpty()) {
279 return CONTAINS | DISJOINT;
280 }
281 if (c.isFull()) {
282 return WITHIN;
283 }
284 // Determine whether or not the circle and polygon boundaries intersect.
285 // If the polygon vertices are not all inside or all outside of c, then the
286 // boundaries cross.
287 bool inside = false;
288 for (VertexIterator v = begin; v != end; ++v) {
289 double d = (*v - c.getCenter()).getSquaredNorm();
290 if (std::fabs(d - c.getSquaredChordLength()) <
291 MAX_SQUARED_CHORD_LENGTH_ERROR) {
292 // A polygon vertex is close to the circle boundary.
293 return INTERSECTS;
294 }
295 bool b = d < c.getSquaredChordLength();
296 if (v == begin) {
297 inside = b;
298 } else if (inside != b) {
299 // There are box vertices both inside and outside of c.
300 return INTERSECTS;
301 }
302 }
303 if (inside) {
304 // All polygon vertices are inside c. Look for points in the polygon
305 // edge interiors that are outside c.
306 for (VertexIterator a = std::prev(end), b = begin; b != end; a = b, ++b) {
307 Vector3d n = a->robustCross(*b);
308 double d = getMaxSquaredChordLength(c.getCenter(), *a, *b, n);
309 if (d > c.getSquaredChordLength() -
310 MAX_SQUARED_CHORD_LENGTH_ERROR) {
311 return INTERSECTS;
312 }
313 }
314 // The polygon boundary is conclusively inside c. It may still be the
315 // case that the circle punches a hole in the polygon. We check that
316 // the polygon does not contain the complement of c by testing whether
317 // or not it contains the anti-center of c.
318 if (contains(begin, end, -c.getCenter())) {
319 return INTERSECTS;
320 }
321 return WITHIN;
322 }
323 // All polygon vertices are outside c. Look for points in the polygon edge
324 // interiors that are inside c.
325 for (VertexIterator a = std::prev(end), b = begin; b != end; a = b, ++b) {
326 Vector3d n = a->robustCross(*b);
327 double d = getMinSquaredChordLength(c.getCenter(), *a, *b, n);
328 if (d < c.getSquaredChordLength() + MAX_SQUARED_CHORD_LENGTH_ERROR) {
329 return INTERSECTS;
330 }
331 }
332 // The polygon boundary is conclusively outside of c. If the polygon
333 // contains the circle center, then the polygon contains c. Otherwise, the
334 // polygon and circle are disjoint.
335 if (contains(begin, end, c.getCenter())) {
336 return CONTAINS;
337 }
338 return DISJOINT;
339}
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:65
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:43

◆ relate() [3/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const begin,
VertexIterator const end,
ConvexPolygon const & p )

Definition at line 401 of file ConvexPolygonImpl.h.

404{
405 return relate(begin, end, p.getVertices().begin(), p.getVertices().end());
406}
Relationship relate(VertexIterator const begin, VertexIterator const end, Box const &b)

◆ relate() [4/5]

template<typename VertexIterator >
Relationship lsst::sphgeom::detail::relate ( VertexIterator const begin,
VertexIterator const end,
Ellipse const & e )

Definition at line 409 of file ConvexPolygonImpl.h.

412{
413 return relate(begin, end, e.getBoundingCircle()) & (CONTAINS | DISJOINT);
414}

◆ relate() [5/5]

template<typename VertexIterator1 , typename VertexIterator2 >
Relationship lsst::sphgeom::detail::relate ( VertexIterator1 const begin1,
VertexIterator1 const end1,
VertexIterator2 const begin2,
VertexIterator2 const end2 )

Definition at line 343 of file ConvexPolygonImpl.h.

347{
348 // TODO(smm): Make this more performant. Instead of the current quadratic
349 // implementation, it should be possible to determine whether the boundaries
350 // intersect by adapting the following method to the sphere:
351 //
352 // A new linear algorithm for intersecting convex polygons
353 // Computer Graphics and Image Processing, Volume 19, Issue 1, May 1982, Page 92
354 // Joseph O'Rourke, Chi-Bin Chien, Thomas Olson, David Naddor
355 //
356 // http://www.sciencedirect.com/science/article/pii/0146664X82900235
357 bool all1 = true;
358 bool any1 = false;
359 bool all2 = true;
360 bool any2 = false;
361 for (VertexIterator1 i = begin1; i != end1; ++i) {
362 bool b = contains(begin2, end2, *i);
363 all1 = b && all1;
364 any1 = b || any1;
365 }
366 for (VertexIterator2 j = begin2; j != end2; ++j) {
367 bool b = contains(begin1, end1, *j);
368 all2 = b && all2;
369 any2 = b || any2;
370 }
371 if (all1 || all2) {
372 // All vertices of one or both polygons are inside the other
373 return (all1 ? WITHIN : INTERSECTS) | (all2 ? CONTAINS : INTERSECTS);
374 }
375 if (any1 || any2) {
376 // The polygons have at least one point in common.
377 return INTERSECTS;
378 }
379 // No vertex of either polygon is inside the other. Consider all
380 // possible edge pairs and look for a crossing.
381 for (VertexIterator1 a = std::prev(end1), b = begin1;
382 b != end1; a = b, ++b) {
383 for (VertexIterator2 c = std::prev(end2), d = begin2;
384 d != end2; c = d, ++d) {
385 int acd = orientation(*a, *c, *d);
386 int bdc = orientation(*b, *d, *c);
387 if (acd == bdc && acd != 0) {
388 int cba = orientation(*c, *b, *a);
389 int dab = orientation(*d, *a, *b);
390 if (cba == dab && cba == acd) {
391 // Found a non-degenerate edge crossing
392 return INTERSECTS;
393 }
394 }
395 }
396 }
397 return DISJOINT;
398}