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
ConvexPolygonImpl.h
Go to the documentation of this file.
1/*
2 * This file is part of sphgeom.
3 *
4 * Developed for the LSST Data Management System.
5 * This product includes software developed by the LSST Project
6 * (http://www.lsst.org).
7 * See the COPYRIGHT file at the top-level directory of this distribution
8 * for details of code ownership.
9 *
10 * This software is dual licensed under the GNU General Public License and also
11 * under a 3-clause BSD license. Recipients may choose which of these licenses
12 * to use; please see the files gpl-3.0.txt and/or bsd_license.txt,
13 * respectively. If you choose the GPL option then the following text applies
14 * (but note that there is still no warranty even if you opt for BSD instead):
15 *
16 * This program is free software: you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation, either version 3 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program. If not, see <http://www.gnu.org/licenses/>.
28 */
29
30#ifndef LSST_SPHGEOM_CONVEXPOLYGONIMPL_H_
31#define LSST_SPHGEOM_CONVEXPOLYGONIMPL_H_
32
41
42#include "lsst/sphgeom/Box.h"
43#include "lsst/sphgeom/Box3d.h"
44#include "lsst/sphgeom/Circle.h"
47#include "lsst/sphgeom/utils.h"
48
49
50namespace lsst {
51namespace sphgeom {
52namespace detail {
53
54template <typename VertexIterator>
55UnitVector3d centroid(VertexIterator const begin, VertexIterator const end) {
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}
74
75template <typename VertexIterator>
76Circle boundingCircle(VertexIterator const begin, VertexIterator const end) {
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}
89
90
91template <typename VertexIterator>
92Box boundingBox(VertexIterator const begin, VertexIterator const end) {
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}
157
158template <typename VertexIterator>
159Box3d boundingBox3d(VertexIterator const begin, VertexIterator const end) {
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}
248
249template <typename VertexIterator>
250bool contains(VertexIterator const begin,
251 VertexIterator const end,
252 UnitVector3d const & v)
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}
263
264template <typename VertexIterator>
265Relationship relate(VertexIterator const begin,
266 VertexIterator const end,
267 Box const & b)
268{
269 // TODO(smm): be more accurate when computing box relations.
270 return boundingBox(begin, end).relate(b) & (DISJOINT | WITHIN);
271}
272
273template <typename VertexIterator>
274Relationship relate(VertexIterator const begin,
275 VertexIterator const end,
276 Circle const & c)
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()) <
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() -
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);
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}
340
341template <typename VertexIterator1,
342 typename VertexIterator2>
343Relationship relate(VertexIterator1 const begin1,
344 VertexIterator1 const end1,
345 VertexIterator2 const begin2,
346 VertexIterator2 const end2)
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}
399
400template <typename VertexIterator>
401Relationship relate(VertexIterator const begin,
402 VertexIterator const end,
403 ConvexPolygon const & p)
404{
405 return relate(begin, end, p.getVertices().begin(), p.getVertices().end());
406}
407
408template <typename VertexIterator>
409Relationship relate(VertexIterator const begin,
410 VertexIterator const end,
411 Ellipse const & e)
412{
413 return relate(begin, end, e.getBoundingCircle()) & (CONTAINS | DISJOINT);
414}
415
416}}} // namespace lsst::sphgeom::detail
417
418#endif // LSST_SPHGEOM_CONVEXPOLYGONIMPL_H_
AmpInfoBoxKey bbox
Definition Amplifier.cc:117
int end
This file declares a class for representing axis-aligned bounding boxes in ℝ³.
This file declares a class for representing circular regions on the unit sphere.
table::Key< int > b
table::Key< int > a
T atan2(T... args)
Angle represents an angle in radians.
Definition Angle.h:50
AngleInterval represents closed intervals of arbitrary angles.
Box3d represents a box in ℝ³.
Definition Box3d.h:49
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition Box.h:61
Relationship relate(LonLat const &p) const
Definition Box.h:304
static NormalizedAngleInterval allLongitudes()
allLongitudes returns a normalized angle interval containing all valid longitude angles.
Definition Box.h:88
Circle is a circular region on the unit sphere that contains its boundary.
Definition Circle.h:53
bool isEmpty() const
Definition Circle.h:116
bool isFull() const
Definition Circle.h:121
double getSquaredChordLength() const
getSquaredChordLength returns the squared length of chords between the circle center and points on th...
Definition Circle.h:130
UnitVector3d const & getCenter() const
getCenter returns the center of this circle as a unit vector.
Definition Circle.h:125
ConvexPolygon is a closed convex polygon on the unit sphere.
std::vector< UnitVector3d > const & getVertices() const
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
Interval1d represents closed intervals of ℝ.
Definition Interval1d.h:47
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition LonLat.h:55
static Angle latitudeOf(Vector3d const &v)
latitudeOf returns the latitude of the point on the unit sphere corresponding to the direction of v.
Definition LonLat.cc:44
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition Vector3d.h:51
double dot(Vector3d const &v) const
dot returns the inner product of this vector and v.
Definition Vector3d.h:80
double x() const
Definition Vector3d.h:73
double y() const
Definition Vector3d.h:75
double normalize()
normalize scales this vector to have unit norm and returns its norm prior to scaling.
Definition Vector3d.cc:48
double z() const
Definition Vector3d.h:77
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this vector and v.
Definition Vector3d.h:108
T fabs(T... args)
T max(T... args)
T min(T... args)
bool contains(VertexIterator const begin, VertexIterator const end, UnitVector3d const &v)
Circle boundingCircle(VertexIterator const begin, VertexIterator const end)
Box3d boundingBox3d(VertexIterator const begin, VertexIterator const end)
UnitVector3d centroid(VertexIterator const begin, VertexIterator const end)
Relationship relate(VertexIterator const begin, VertexIterator const end, Box const &b)
Box boundingBox(VertexIterator const begin, VertexIterator const end)
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
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
constexpr double MAX_SQUARED_CHORD_LENGTH_ERROR
Definition constants.h:57
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.
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
constexpr double PI
Definition constants.h:43
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
T prev(T... args)
This file declares a class for representing longitude/latitude angle boxes on the unit sphere.
This file declares a class for representing elliptical regions on the unit sphere.
This file declares functions for orienting points on the sphere.
This file declares miscellaneous utility functions.
T sqrt(T... args)