LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
ConvexPolygonImpl.h
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2014-2016 AURA/LSST.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <https://www.lsstcorp.org/LegalNotices/>.
21  */
22 
23 #ifndef LSST_SPHGEOM_CONVEXPOLYGONIMPL_H_
24 #define LSST_SPHGEOM_CONVEXPOLYGONIMPL_H_
25 
34 
35 #include "lsst/sphgeom/Box.h"
36 #include "lsst/sphgeom/Box3d.h"
37 #include "lsst/sphgeom/Circle.h"
38 #include "lsst/sphgeom/Ellipse.h"
40 #include "lsst/sphgeom/utils.h"
41 
42 
43 namespace lsst {
44 namespace sphgeom {
45 namespace detail {
46 
47 template <typename VertexIterator>
48 UnitVector3d centroid(VertexIterator const begin, VertexIterator const end) {
49  // The center of mass is obtained via trivial generalization of
50  // the formula for spherical triangles from:
51  //
52  // The centroid and inertia tensor for a spherical triangle
53  // John E. Brock
54  // 1974, Naval Postgraduate School, Monterey Calif.
55  Vector3d cm;
56  VertexIterator i = std::prev(end);
57  VertexIterator j = begin;
58  for (; j != end; i = j, ++j) {
59  Vector3d v = (*i).robustCross(*j);
60  double s = 0.5 * v.normalize();
61  double c = (*i).dot(*j);
62  double a = (s == 0.0 && c == 0.0) ? 0.0 : std::atan2(s, c);
63  cm += v * a;
64  }
65  return UnitVector3d(cm);
66 }
67 
68 template <typename VertexIterator>
69 Circle boundingCircle(VertexIterator const begin, VertexIterator const end) {
70  UnitVector3d c = centroid(begin, end);
71  // Compute the maximum squared chord length between the centroid and
72  // all vertices.
73  VertexIterator i = begin;
74  double cl2 = 0.0;
75  for (; i != end; ++i) {
76  cl2 = std::max(cl2, (*i - c).getSquaredNorm());
77  }
78  // Add double the maximum squared-chord-length error, so that the
79  // bounding circle we return also reliably CONTAINS this polygon.
80  return Circle(c, cl2 + 2.0 * MAX_SQUARED_CHORD_LENGTH_ERROR);
81 }
82 
83 
84 template <typename VertexIterator>
85 Box boundingBox(VertexIterator const begin, VertexIterator const end) {
86  Angle const eps(5.0e-10); // ~ 0.1 milli-arcseconds
87  Box bbox;
88  VertexIterator i = std::prev(end);
89  VertexIterator j = begin;
90  bool haveCW = false;
91  bool haveCCW = false;
92  // Compute the bounding box for each vertex. When converting a Vector3d
93  // to a LonLat, the relative error on the longitude is about 4*2^-53,
94  // and the relative error on the latitude is about twice that (assuming
95  // std::atan2 and std::sqrt accurate to within 1 ulp). We convert each
96  // vertex to a conservative bounding box for its spherical coordinates,
97  // and compute a bounding box for the union of all these boxes.
98  //
99  // Furthermore, the latitude range of an edge can be greater than the
100  // latitude range of its endpoints - this occurs when the minimum or
101  // maximum latitude point on the great circle defined by the edge vertices
102  // lies in the edge interior.
103  for (; j != end; i = j, ++j) {
104  LonLat p(*j);
105  bbox.expandTo(Box(p, eps, eps));
106  if (!haveCW || !haveCCW) {
107  int o = orientationZ(*i, *j);
108  haveCCW = haveCCW || (o > 0);
109  haveCW = haveCW || (o < 0);
110  }
111  // Compute the plane normal for edge i, j.
112  Vector3d n = (*i).robustCross(*j);
113  // Compute a vector v with positive z component that lies on both the
114  // edge plane and on the plane defined by the z axis and the edge plane
115  // normal. This is the direction of maximum latitude for the great
116  // circle containing the edge, and -v is the direction of minimum
117  // latitude.
118  //
119  // TODO(smm): Do a proper error analysis.
120  Vector3d v(-n.x() * n.z(),
121  -n.y() * n.z(),
122  n.x() * n.x() + n.y() * n.y());
123  if (v != Vector3d()) {
124  // The plane defined by the z axis and n has normal
125  // (-n.y(), n.x(), 0.0). Compute the dot product of this plane
126  // normal with vertices i and j.
127  double zni = i->y() * n.x() - i->x() * n.y();
128  double znj = j->y() * n.x() - j->x() * n.y();
129  // Check if v or -v is in the edge interior.
130  if (zni > 0.0 && znj < 0.0) {
131  bbox = Box(bbox.getLon(), bbox.getLat().expandedTo(
132  LonLat::latitudeOf(v) + eps));
133  } else if (zni < 0.0 && znj > 0.0) {
134  bbox = Box(bbox.getLon(), bbox.getLat().expandedTo(
135  LonLat::latitudeOf(-v) - eps));
136  }
137  }
138  }
139  // If this polygon contains a pole, its bounding box must contain all
140  // longitudes.
141  if (!haveCW) {
142  Box northPole(Box::allLongitudes(), AngleInterval(Angle(0.5 * PI)));
143  bbox.expandTo(northPole);
144  } else if (!haveCCW) {
145  Box southPole(Box::allLongitudes(), AngleInterval(Angle(-0.5 * PI)));
146  bbox.expandTo(southPole);
147  }
148  return bbox;
149 }
150 
151 template <typename VertexIterator>
152 Box3d boundingBox3d(VertexIterator const begin, VertexIterator const end) {
153  static double const maxError = 1.0e-14;
154  // Compute the extrema of all vertex coordinates.
155  VertexIterator j = begin;
156  double emin[3] = { j->x(), j->y(), j->z() };
157  double emax[3] = { j->x(), j->y(), j->z() };
158  for (++j; j != end; ++j) {
159  for (int i = 0; i < 3; ++i) {
160  double v = j->operator()(i);
161  emin[i] = std::min(emin[i], v);
162  emax[i] = std::max(emax[i], v);
163  }
164  }
165  // Compute the extrema of all edges.
166  //
167  // It can be shown that the great circle with unit normal vector
168  // n = (n₀, n₁, n₂) has extrema in x at:
169  //
170  // (∓√(1 - n₀²), ±n₁n₀/√(1 - n₀²), ±n₂n₀/√(1 - n₀²))
171  //
172  // in y at:
173  //
174  // (±n₀n₁/√(1 - n₁²), ∓√(1 - n₁²), ±n₂n₁/√(1 - n₁²))
175  //
176  // and in z at
177  //
178  // (±n₀n₂/√(1 - n₂²), ±n₁n₂/√(1 - n₂²), ∓√(1 - n₂²))
179  //
180  // Compute these vectors for each edge, determine whether they lie in
181  // the edge, and update the extrema if so. Rounding errors in these
182  // computations are compensated for by expanding the bounding box
183  // prior to returning it.
184  j = std::prev(end);
185  VertexIterator k = begin;
186  for (; k != end; j = k, ++k) {
187  UnitVector3d n(j->robustCross(*k));
188  for (int i = 0; i < 3; ++i) {
189  double ni = n(i);
190  double d = std::fabs(1.0 - ni * ni);
191  if (d > 0.0) {
192  Vector3d e(i == 0 ? -d : n.x() * ni,
193  i == 1 ? -d : n.y() * ni,
194  i == 2 ? -d : n.z() * ni);
195  // If e or -e lies in the lune defined by the half great
196  // circle passing through n and a and the half great circle
197  // passing through n and b, the edge contains an extremum.
198  Vector3d v = e.cross(n);
199  double vdj = v.dot(*j);
200  double vdk = v.dot(*k);
201  if (vdj >= 0.0 && vdk <= 0.0) {
202  emin[i] = std::min(emin[i], -std::sqrt(d));
203  }
204  if (vdj <= 0.0 && vdk >= 0.0) {
205  emax[i] = std::max(emax[i], std::sqrt(d));
206  }
207  }
208  }
209  }
210  // Check whether the standard basis vectors and their antipodes
211  // are inside this polygon.
212  bool a[3] = { true, true, true };
213  bool b[3] = { true, true, true };
214  j = std::prev(end);
215  k = begin;
216  for (; k != end; j = k, ++k) {
217  // Test the standard basis vectors against the plane defined by
218  // vertices (j, k). Note that orientation(-x, *j, *k) =
219  // -orientation(x, *j, *k).
220  int ox = orientationX(*j, *k);
221  a[0] = a[0] && (ox <= 0);
222  b[0] = b[0] && (ox >= 0);
223  int oy = orientationY(*j, *k);
224  a[1] = a[1] && (oy <= 0);
225  b[1] = b[1] && (oy >= 0);
226  int oz = orientationZ(*j, *k);
227  a[2] = a[2] && (oz <= 0);
228  b[2] = b[2] && (oz >= 0);
229  }
230  // At this point, b[i] is true iff the standard basis vector eᵢ
231  // is inside all the half spaces defined by the polygon edges.
232  // Similarly, a[i] is true iff -eᵢ is inside the same half spaces.
233  for (int i = 0; i < 3; ++i) {
234  emin[i] = a[i] ? -1.0 : std::max(-1.0, emin[i] - maxError);
235  emax[i] = b[i] ? 1.0 : std::min(1.0, emax[i] + maxError);
236  }
237  return Box3d(Interval1d(emin[0], emax[0]),
238  Interval1d(emin[1], emax[1]),
239  Interval1d(emin[2], emax[2]));
240 }
241 
242 template <typename VertexIterator>
243 bool contains(VertexIterator const begin,
244  VertexIterator const end,
245  UnitVector3d const & v)
246 {
247  VertexIterator i = std::prev(end);
248  VertexIterator j = begin;
249  for (; j != end; i = j, ++j) {
250  if (orientation(v, *i, *j) < 0) {
251  return false;
252  }
253  }
254  return true;
255 }
256 
257 template <typename VertexIterator>
258 Relationship relate(VertexIterator const begin,
259  VertexIterator const end,
260  Box const & b)
261 {
262  // TODO(smm): be more accurate when computing box relations.
263  return boundingBox(begin, end).relate(b) & (DISJOINT | WITHIN);
264 }
265 
266 template <typename VertexIterator>
267 Relationship relate(VertexIterator const begin,
268  VertexIterator const end,
269  Circle const & c)
270 {
271  if (c.isEmpty()) {
272  return CONTAINS | DISJOINT;
273  }
274  if (c.isFull()) {
275  return WITHIN;
276  }
277  // Determine whether or not the circle and polygon boundaries intersect.
278  // If the polygon vertices are not all inside or all outside of c, then the
279  // boundaries cross.
280  bool inside = false;
281  for (VertexIterator v = begin; v != end; ++v) {
282  double d = (*v - c.getCenter()).getSquaredNorm();
283  if (std::fabs(d - c.getSquaredChordLength()) <
285  // A polygon vertex is close to the circle boundary.
286  return INTERSECTS;
287  }
288  bool b = d < c.getSquaredChordLength();
289  if (v == begin) {
290  inside = b;
291  } else if (inside != b) {
292  // There are box vertices both inside and outside of c.
293  return INTERSECTS;
294  }
295  }
296  if (inside) {
297  // All polygon vertices are inside c. Look for points in the polygon
298  // edge interiors that are outside c.
299  for (VertexIterator a = std::prev(end), b = begin; b != end; a = b, ++b) {
300  Vector3d n = a->robustCross(*b);
301  double d = getMaxSquaredChordLength(c.getCenter(), *a, *b, n);
302  if (d > c.getSquaredChordLength() -
304  return INTERSECTS;
305  }
306  }
307  // The polygon boundary is conclusively inside c. It may still be the
308  // case that the circle punches a hole in the polygon. We check that
309  // the polygon does not contain the complement of c by testing whether
310  // or not it contains the anti-center of c.
311  if (contains(begin, end, -c.getCenter())) {
312  return INTERSECTS;
313  }
314  return WITHIN;
315  }
316  // All polygon vertices are outside c. Look for points in the polygon edge
317  // interiors that are inside c.
318  for (VertexIterator a = std::prev(end), b = begin; b != end; a = b, ++b) {
319  Vector3d n = a->robustCross(*b);
320  double d = getMinSquaredChordLength(c.getCenter(), *a, *b, n);
322  return INTERSECTS;
323  }
324  }
325  // The polygon boundary is conclusively outside of c. If the polygon
326  // contains the circle center, then the polygon contains c. Otherwise, the
327  // polygon and circle are disjoint.
328  if (contains(begin, end, c.getCenter())) {
329  return CONTAINS;
330  }
331  return DISJOINT;
332 }
333 
334 template <typename VertexIterator1,
335  typename VertexIterator2>
336 Relationship relate(VertexIterator1 const begin1,
337  VertexIterator1 const end1,
338  VertexIterator2 const begin2,
339  VertexIterator2 const end2)
340 {
341  // TODO(smm): Make this more performant. Instead of the current quadratic
342  // implementation, it should be possible to determine whether the boundaries
343  // intersect by adapting the following method to the sphere:
344  //
345  // A new linear algorithm for intersecting convex polygons
346  // Computer Graphics and Image Processing, Volume 19, Issue 1, May 1982, Page 92
347  // Joseph O'Rourke, Chi-Bin Chien, Thomas Olson, David Naddor
348  //
349  // http://www.sciencedirect.com/science/article/pii/0146664X82900235
350  bool all1 = true;
351  bool any1 = false;
352  bool all2 = true;
353  bool any2 = false;
354  for (VertexIterator1 i = begin1; i != end1; ++i) {
355  bool b = contains(begin2, end2, *i);
356  all1 = b && all1;
357  any1 = b || any1;
358  }
359  for (VertexIterator2 j = begin2; j != end2; ++j) {
360  bool b = contains(begin1, end1, *j);
361  all2 = b && all2;
362  any2 = b || any2;
363  }
364  if (all1 || all2) {
365  // All vertices of one or both polygons are inside the other
366  return (all1 ? WITHIN : INTERSECTS) | (all2 ? CONTAINS : INTERSECTS);
367  }
368  if (any1 || any2) {
369  // The polygons have at least one point in common.
370  return INTERSECTS;
371  }
372  // No vertex of either polygon is inside the other. Consider all
373  // possible edge pairs and look for a crossing.
374  for (VertexIterator1 a = std::prev(end1), b = begin1;
375  b != end1; a = b, ++b) {
376  for (VertexIterator2 c = std::prev(end2), d = begin2;
377  d != end2; c = d, ++d) {
378  int acd = orientation(*a, *c, *d);
379  int bdc = orientation(*b, *d, *c);
380  if (acd == bdc && acd != 0) {
381  int cba = orientation(*c, *b, *a);
382  int dab = orientation(*d, *a, *b);
383  if (cba == dab && cba == acd) {
384  // Found a non-degenerate edge crossing
385  return INTERSECTS;
386  }
387  }
388  }
389  }
390  return DISJOINT;
391 }
392 
393 template <typename VertexIterator>
394 Relationship relate(VertexIterator const begin,
395  VertexIterator const end,
396  ConvexPolygon const & p)
397 {
398  return relate(begin, end, p.getVertices().begin(), p.getVertices().end());
399 }
400 
401 template <typename VertexIterator>
402 Relationship relate(VertexIterator const begin,
403  VertexIterator const end,
404  Ellipse const & e)
405 {
406  return relate(begin, end, e.getBoundingCircle()) & (CONTAINS | DISJOINT);
407 }
408 
409 }}} // namespace lsst::sphgeom::detail
410 
411 #endif // LSST_SPHGEOM_CONVEXPOLYGONIMPL_H_
double normalize()
normalize scales this vector to have unit norm and returns its norm prior to scaling.
Definition: Vector3d.cc:41
This file declares a class for representing circular regions on the unit sphere.
Box & expandTo(LonLat const &x)
Definition: Box.h:246
static NormalizedAngleInterval allLongitudes()
allLongitudes returns a normalized angle interval containing all valid longitude angles.
Definition: Box.h:81
This file declares a class for representing longitude/latitude angle boxes on the unit sphere...
bool contains(VertexIterator const begin, VertexIterator const end, UnitVector3d const &v)
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
Box boundingBox(VertexIterator const begin, VertexIterator const end)
double getSquaredChordLength() const
getSquaredChordLength returns the squared length of chords between the circle center and points on th...
Definition: Circle.h:123
AngleInterval represents closed intervals of arbitrary angles.
Definition: AngleInterval.h:39
constexpr double MAX_SQUARED_CHORD_LENGTH_ERROR
Definition: constants.h:50
table::Key< int > b
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
Definition: orientation.cc:237
Interval1d represents closed intervals of ℝ.
Definition: Interval1d.h:39
table::Key< int > a
This file declares a class for representing axis-aligned bounding boxes in ℝ³.
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition: Box.h:54
double y() const
Definition: Vector3d.h:68
T atan2(T... args)
Relationship relate(LonLat const &p) const
Definition: Box.h:297
T prev(T... args)
T min(T... args)
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
Ellipse is an elliptical region on the sphere.
Definition: Ellipse.h:169
double x() const
Definition: Vector3d.h:66
A base class for image defects.
std::vector< UnitVector3d > const & getVertices() const
Definition: ConvexPolygon.h:99
UnitVector3d centroid(VertexIterator const begin, VertexIterator const end)
NormalizedAngleInterval const & getLon() const
getLon returns the longitude interval of this box.
Definition: Box.h:145
Relationship relate(VertexIterator const begin, VertexIterator const end, Box const &b)
constexpr double PI
Definition: constants.h:36
This file declares functions for orienting points on the sphere.
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:37
Derived expandedTo(Scalar x) const
Definition: Interval.h:218
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
This file declares miscellaneous utility functions.
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
Definition: orientation.cc:232
table::Box2IKey bbox
Definition: Detector.cc:169
T fabs(T... args)
This file declares a class for representing elliptical regions on the unit sphere.
T max(T... args)
solver_t * s
Circle is a circular region on the unit sphere that contains its boundary.
Definition: Circle.h:46
bool isFull() const
Definition: Circle.h:114
ConvexPolygon is a closed convex polygon on the unit sphere.
Definition: ConvexPolygon.h:57
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Definition: Ellipse.cc:241
Angle represents an angle in radians.
Definition: Angle.h:43
AngleInterval const & getLat() const
getLat returns the latitude interval of this box.
Definition: Box.h:148
UnitVector3d const & getCenter() const
getCenter returns the center of this circle as a unit vector.
Definition: Circle.h:118
LonLat represents a spherical coordinate (longitude/latitude angle) pair.
Definition: LonLat.h:48
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this vector and v.
Definition: Vector3d.h:101
Box3d represents a box in ℝ³.
Definition: Box3d.h:42
double z() const
Definition: Vector3d.h:70
lsst::geom::Angle Angle
Definition: misc.h:33
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
T sqrt(T... args)
Interval1d const & x() const
Definition: Box3d.h:123
bool isEmpty() const
Definition: Circle.h:109
STL class.
Box3d boundingBox3d(VertexIterator const begin, VertexIterator const end)
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
Definition: orientation.cc:227
int end
Circle boundingCircle(VertexIterator const begin, VertexIterator const end)
double dot(Vector3d const &v) const
dot returns the inner product of this vector and v.
Definition: Vector3d.h:73