LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
ConvexPolygon.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2014-2015 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 
25 
27 
28 #include <ostream>
29 #include <stdexcept>
30 
31 #include "lsst/sphgeom/codec.h"
33 
34 #include "ConvexPolygonImpl.h"
35 
36 
37 namespace lsst {
38 namespace sphgeom {
39 
40 namespace {
41 
42 char const * const FOUND_ANTIPODAL_POINT =
43  "The convex hull of the given point set is the "
44  "entire unit sphere";
45 
46 char const * const NOT_ENOUGH_POINTS =
47  "The convex hull of a point set containing less than "
48  "3 distinct, non-coplanar points is not a convex polygon";
49 
50 // `findPlane` rearranges the entries of `points` such that the first two
51 // entries are distinct. If it is unable to do so, or if it encounters any
52 // antipodal points in the process, an exception is thrown. It returns an
53 // iterator to the first point in the input that was not consumed during
54 // the search.
57 {
58  if (points.empty()) {
59  throw std::invalid_argument(NOT_ENOUGH_POINTS);
60  }
61  // Starting with the first point v0, find a distinct second point v1.
62  UnitVector3d & v0 = points[0];
63  UnitVector3d & v1 = points[1];
66  for (; v != end; ++v) {
67  if (v0 == -*v) {
68  throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
69  }
70  if (v0 != *v) {
71  break;
72  }
73  }
74  if (v == end) {
75  throw std::invalid_argument(NOT_ENOUGH_POINTS);
76  }
77  v1 = *v;
78  return ++v;
79 }
80 
81 // `findTriangle` rearranges the entries of `points` such that the first
82 // three entries have counter-clockwise orientation. If it is unable to do so,
83 // or if it encounters any antipodal points in the process, an exception is
84 // thrown. It returns an iterator to the first point in the input that was not
85 // consumed during the search.
88 {
89  std::vector<UnitVector3d>::iterator v = findPlane(points);
90  std::vector<UnitVector3d>::iterator const end = points.end();
91  UnitVector3d & v0 = points[0];
92  UnitVector3d & v1 = points[1];
93  // Note that robustCross() gives a non-zero result for distinct,
94  // non-antipodal inputs, and that normalization never maps a non-zero
95  // vector to the zero vector.
96  UnitVector3d n(v0.robustCross(v1));
97  for (; v != end; ++v) {
98  int ccw = orientation(v0, v1, *v);
99  if (ccw > 0) {
100  // We found a counter-clockwise triangle.
101  break;
102  } else if (ccw < 0) {
103  // We found a clockwise triangle. Swap the first two vertices to
104  // flip its orientation.
105  std::swap(v0, v1);
106  break;
107  }
108  // v, v0 and v1 are coplanar.
109  if (*v == v0 || *v == v1) {
110  continue;
111  }
112  if (*v == -v0 || *v == -v1) {
113  throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
114  }
115  // At this point, v, v0 and v1 are distinct and non-antipodal.
116  // If v is in the interior of the great circle segment (v0, v1),
117  // discard v and continue.
118  int v0v = orientation(n, v0, *v);
119  int vv1 = orientation(n, *v, v1);
120  if (v0v == vv1) {
121  continue;
122  }
123  int v0v1 = orientation(n, v0, v1);
124  // If v1 is in the interior of (v0, v), replace v1 with v and continue.
125  if (v0v1 == -vv1) {
126  v1 = *v; continue;
127  }
128  // If v0 is in the interior of (v, v1), replace v0 with v and continue.
129  if (-v0v == v0v1) {
130  v0 = *v; continue;
131  }
132  // Otherwise, (v0, v1) ∪ (v1, v) and (v, v0) ∪ (v0, v1) both span
133  // more than π radians of the great circle defined by the v0 and v1,
134  // so there is a pair of antipodal points in the corresponding great
135  // circle segment.
136  throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
137  }
138  if (v == end) {
139  throw std::invalid_argument(NOT_ENOUGH_POINTS);
140  }
141  points[2] = *v;
142  return ++v;
143 }
144 
145 void computeHull(std::vector<UnitVector3d> & points) {
146  typedef std::vector<UnitVector3d>::iterator VertexIterator;
147  VertexIterator hullEnd = points.begin() + 3;
148  VertexIterator const end = points.end();
149  // Start with a triangular hull.
150  for (VertexIterator v = findTriangle(points); v != end; ++v) {
151  // Compute the hull of the current hull and v.
152  //
153  // Notice that if v is in the current hull, v can be ignored. If
154  // -v is in the current hull, then the hull of v and the current hull
155  // is not a convex polygon.
156  //
157  // Otherwise, let i and j be the first and second end-points of an edge
158  // in the current hull. We define the orientation of vertex j with
159  // respect to vertex v as orientation(v, i, j). When neither v or -v
160  // is in the current hull, there must be a sequence of consecutive
161  // hull vertices that do not have counter-clockwise orientation with
162  // respect to v. Insert v before the first vertex in this sequence
163  // and remove all vertices in the sequence except the last to obtain
164  // a new, larger convex hull.
165  VertexIterator i = hullEnd - 1;
166  VertexIterator j = points.begin();
167  // toCCW is the vertex before the transition to counter-clockwise
168  // orientation with respect to v.
169  VertexIterator toCCW = hullEnd;
170  // fromCCW is the vertex at which orientation with respect to v
171  // changes from counter-clockwise to clockwise or coplanar. It may
172  // equal toCCW.
173  VertexIterator fromCCW = hullEnd;
174  // Compute the orientation of the first point in the current hull
175  // with respect to v.
176  bool const firstCCW = orientation(*v, *i, *j) > 0;
177  bool prevCCW = firstCCW;
178  // Compute the orientation of points in the current hull with respect
179  // to v, starting with the second point. Update toCCW / fromCCW when
180  // we transition to / from counter-clockwise orientation.
181  for (i = j, ++j; j != hullEnd; i = j, ++j) {
182  if (orientation(*v, *i, *j) > 0) {
183  if (!prevCCW) {
184  toCCW = i;
185  prevCCW = true;
186  }
187  } else if (prevCCW) {
188  fromCCW = j;
189  prevCCW = false;
190  }
191  }
192  // Now that we know the orientation of the last point in the current
193  // hull with respect to v, consider the first point in the current hull.
194  if (firstCCW) {
195  if (!prevCCW) {
196  toCCW = i;
197  }
198  } else if (prevCCW) {
199  fromCCW = points.begin();
200  }
201  // Handle the case where there is never a transition to / from
202  // counter-clockwise orientation.
203  if (toCCW == hullEnd) {
204  // If all vertices are counter-clockwise with respect to v,
205  // v is inside the current hull and can be ignored.
206  if (firstCCW) {
207  continue;
208  }
209  // Otherwise, no vertex is counter-clockwise with respect to v,
210  // so that -v is inside the current hull.
211  throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
212  }
213  // Insert v into the current hull at fromCCW, and remove
214  // all vertices between fromCCW and toCCW.
215  if (toCCW < fromCCW) {
216  if (toCCW != points.begin()) {
217  fromCCW = std::copy(toCCW, fromCCW, points.begin());
218  }
219  *fromCCW++ = *v;
220  hullEnd = fromCCW;
221  } else if (toCCW > fromCCW) {
222  *fromCCW++ = *v;
223  if (toCCW != fromCCW) {
224  hullEnd = std::copy(toCCW, hullEnd, fromCCW);
225  }
226  } else {
227  if (fromCCW == points.begin()) {
228  *hullEnd = *v;
229  } else {
230  UnitVector3d u = *v;
231  std::copy_backward(fromCCW, hullEnd, hullEnd + 1);
232  *fromCCW = u;
233  }
234  ++hullEnd;
235  }
236  }
237  points.erase(hullEnd, end);
238 }
239 
240 // TODO(smm): for all of this to be fully rigorous, we must prove that no two
241 // UnitVector3d objects u and v are exactly colinear unless u == v or u == -v.
242 // It's not clear that this is true. For example, (1, 0, 0) and (1 + ε, 0, 0)
243 // are colinear. This means that UnitVector3d should probably always normalize
244 // on construction. Currently, it does not normalize when created from a LonLat,
245 // and also contains some escape-hatches for performance. The normalize()
246 // function implementation may also need to be revisited.
247 
248 // TODO(smm): This implementation is quadratic. It would be nice to implement
249 // a fast hull merging algorithm, which could then be used to implement Chan's
250 // algorithm.
251 
252 } // unnamed namespace
253 
254 
255 ConvexPolygon::ConvexPolygon(std::vector<UnitVector3d> const & points) :
256  _vertices(points)
257 {
258  computeHull(_vertices);
259 }
260 
262  if (this == &p) {
263  return true;
264  }
265  if (_vertices.size() != p._vertices.size()) {
266  return false;
267  }
268  VertexIterator i = _vertices.begin();
269  VertexIterator f = p._vertices.begin();
270  VertexIterator const ep = p._vertices.end();
271  // Find vertex f of p equal to the first vertex of this polygon.
272  for (; f != ep; ++f) {
273  if (*i == *f) {
274  break;
275  }
276  }
277  if (f == ep) {
278  // No vertex of p is equal to the first vertex of this polygon.
279  return false;
280  }
281  // Now, compare all vertices.
282  ++i;
283  for (VertexIterator j = f + 1; j != ep; ++i, ++j) {
284  if (*i != *j) {
285  return false;
286  }
287  }
288  for (VertexIterator j = p._vertices.begin(); j != f; ++i, ++j) {
289  if (*i != *j) {
290  return false;
291  }
292  }
293  return true;
294 }
295 
297  return detail::centroid(_vertices.begin(), _vertices.end());
298 }
299 
301  return detail::boundingCircle(_vertices.begin(), _vertices.end());
302 }
303 
305  return detail::boundingBox(_vertices.begin(), _vertices.end());
306 }
307 
309  return detail::boundingBox3d(_vertices.begin(), _vertices.end());
310 }
311 
312 bool ConvexPolygon::contains(UnitVector3d const & v) const {
313  return detail::contains(_vertices.begin(), _vertices.end(), v);
314 }
315 
316 bool ConvexPolygon::contains(Region const & r) const {
317  return (relate(r) & CONTAINS) != 0;
318 }
319 
320 bool ConvexPolygon::isDisjointFrom(Region const & r) const {
321  return (relate(r) & DISJOINT) != 0;
322 }
323 
324 bool ConvexPolygon::intersects(Region const & r) const {
325  return !isDisjointFrom(r);
326 }
327 
328 bool ConvexPolygon::isWithin(Region const & r) const {
329  return (relate(r) & WITHIN) != 0;
330 }
331 
333  return detail::relate(_vertices.begin(), _vertices.end(), b);
334 }
335 
337  return detail::relate(_vertices.begin(), _vertices.end(), c);
338 }
339 
341  return detail::relate(_vertices.begin(), _vertices.end(), p);
342 }
343 
345  return detail::relate(_vertices.begin(), _vertices.end(), e);
346 }
347 
349  std::vector<uint8_t> buffer;
350  uint8_t tc = TYPE_CODE;
351  buffer.reserve(1 + 24 * _vertices.size());
352  buffer.push_back(tc);
353  for (UnitVector3d const & v: _vertices) {
354  encodeDouble(v.x(), buffer);
355  encodeDouble(v.y(), buffer);
356  encodeDouble(v.z(), buffer);
357  }
358  return buffer;
359 }
360 
362  size_t n)
363 {
364  if (buffer == nullptr || *buffer != TYPE_CODE ||
365  n < 1 + 24*3 || (n - 1) % 24 != 0) {
366  throw std::runtime_error("Byte-string is not an encoded ConvexPolygon");
367  }
369  ++buffer;
370  size_t nv = (n - 1) / 24;
371  poly->_vertices.reserve(nv);
372  for (size_t i = 0; i < nv; ++i, buffer += 24) {
373  poly->_vertices.push_back(UnitVector3d::fromNormalized(
374  decodeDouble(buffer),
375  decodeDouble(buffer + 8),
376  decodeDouble(buffer + 16)
377  ));
378  }
379  return poly;
380 }
381 
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 }
391 
392 }} // namespace lsst::sphgeom
bool isDisjointFrom(Region const &r) const
T empty(T... args)
T copy(T... args)
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
Low-level polynomials (including special polynomials) in C++.
Definition: Basis1d.h:26
Box boundingBox(VertexIterator const begin, VertexIterator const end)
T swap(T... args)
bool operator==(ConvexPolygon const &p) const
Two convex polygons are equal iff they contain the same points.
table::Key< int > b
This file contains simple helper functions for encoding and decoding primitive types to/from byte str...
std::ostream & operator<<(std::ostream &, Angle const &)
Definition: Angle.cc:34
T end(T... args)
bool isWithin(Region const &r) const
Box represents a rectangle in spherical coordinate space that contains its boundary.
Definition: Box.h:54
UnitVector3d getCentroid() const
The centroid of a polygon is its center of mass projected onto S², assuming a uniform mass distributi...
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
Ellipse is an elliptical region on the sphere.
Definition: Ellipse.h:169
T push_back(T... args)
bool intersects(Region const &r) const
Relationship relate(Region const &r) const override
A base class for image defects.
std::vector< UnitVector3d > const & getVertices() const
Definition: ConvexPolygon.h:99
UnitVector3d centroid(VertexIterator const begin, VertexIterator const end)
Region is a minimal interface for 2-dimensional regions on the unit sphere.
Definition: Region.h:79
Relationship relate(VertexIterator const begin, VertexIterator const end, Box const &b)
T erase(T... args)
This file contains the meat of the ConvexPolygon implementation.
Circle is a circular region on the unit sphere that contains its boundary.
Definition: Circle.h:46
ConvexPolygon is a closed convex polygon on the unit sphere.
Definition: ConvexPolygon.h:57
double decodeDouble(uint8_t const *buffer)
decode extracts an IEEE double from the 8 byte little-endian byte sequence in buffer.
Definition: codec.h:59
T copy_backward(T... args)
bool contains(UnitVector3d const &v) const override
STL class.
STL class.
Box getBoundingBox() const override
getBoundingBox returns a bounding-box for this region.
T begin(T... args)
This file declares functions for orienting points on the sphere.
static constexpr uint8_t TYPE_CODE
Definition: ConvexPolygon.h:59
Box3d represents a box in ℝ³.
Definition: Box3d.h:42
std::vector< uint8_t > encode() const override
encode serializes this region into an opaque byte string.
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
void encodeDouble(double item, std::vector< uint8_t > &buffer)
encode appends an IEEE double in little-endian byte order to the end of buffer.
Definition: codec.h:38
std::ostream * os
Definition: Schema.cc:746
This file declares a class for representing convex polygons with great circle edges on the unit spher...
STL class.
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
Definition: UnitVector3d.h:82
STL class.
Box3d boundingBox3d(VertexIterator const begin, VertexIterator const end)
int end
Circle boundingCircle(VertexIterator const begin, VertexIterator const end)
T reserve(T... args)
static std::unique_ptr< ConvexPolygon > decode(std::vector< uint8_t > const &s)
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.