LSST Applications g2079a07aa2+86d27d4dc4,g2305ad1205+a659bff248,g2bbee38e9b+3c60f8fe34,g337abbeb29+3c60f8fe34,g33d1c0ed96+3c60f8fe34,g3502564af9+d77d6d1350,g3a166c0a6a+3c60f8fe34,g487adcacf7+25d9892218,g4be5004598+d77d6d1350,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+4d81263f9a,g5cd07815a0+980d2b1c3b,g607f77f49a+d77d6d1350,g858d7b2824+d77d6d1350,g88963caddf+83e433e629,g99cad8db69+a4d3c48eeb,g9ddcbc5298+9a081db1e4,ga1e77700b3+bcf1af89ad,ga57fefb910+9a39d7b2d7,gae0086650b+585e252eca,gb065fddaf9+4f9fd82a2c,gb0e22166c9+60f28cb32d,gb363559e06+d84b1d3d07,gb3b7280ab2+4563d032e1,gb4b16eec92+babe958938,gba4ed39666+c2a2e4ac27,gbb8dafda3b+ed6854b564,gc120e1dc64+b72d212f87,gc28159a63d+3c60f8fe34,gc3e9b769f7+921dbcd359,gcf0d15dbbd+9a39d7b2d7,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+3c60f8fe34,w.2024.21
LSST Data Management Base Package
Loading...
Searching...
No Matches
ConvexPolygon.cc
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
32
34
35#include <ostream>
36#include <stdexcept>
37
38#include "lsst/sphgeom/codec.h"
40
41#include "ConvexPolygonImpl.h"
42
43
44namespace lsst {
45namespace sphgeom {
46
47namespace {
48
49char const * const FOUND_ANTIPODAL_POINT =
50 "The convex hull of the given point set is the "
51 "entire unit sphere";
52
53char const * const NOT_ENOUGH_POINTS =
54 "The convex hull of a point set containing less than "
55 "3 distinct, non-coplanar points is not a convex polygon";
56
57// `findPlane` rearranges the entries of `points` such that the first two
58// entries are distinct. If it is unable to do so, or if it encounters any
59// antipodal points in the process, an exception is thrown. It returns an
60// iterator to the first point in the input that was not consumed during
61// the search.
62std::vector<UnitVector3d>::iterator findPlane(
64{
65 if (points.empty()) {
66 throw std::invalid_argument(NOT_ENOUGH_POINTS);
67 }
68 // Starting with the first point v0, find a distinct second point v1.
69 UnitVector3d & v0 = points[0];
70 UnitVector3d & v1 = points[1];
71 std::vector<UnitVector3d>::iterator const end = points.end();
72 std::vector<UnitVector3d>::iterator v = points.begin() + 1;
73 for (; v != end; ++v) {
74 if (v0 == -*v) {
75 throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
76 }
77 if (v0 != *v) {
78 break;
79 }
80 }
81 if (v == end) {
82 throw std::invalid_argument(NOT_ENOUGH_POINTS);
83 }
84 v1 = *v;
85 return ++v;
86}
87
88// `findTriangle` rearranges the entries of `points` such that the first
89// three entries have counter-clockwise orientation. If it is unable to do so,
90// or if it encounters any antipodal points in the process, an exception is
91// thrown. It returns an iterator to the first point in the input that was not
92// consumed during the search.
93std::vector<UnitVector3d>::iterator findTriangle(
95{
96 std::vector<UnitVector3d>::iterator v = findPlane(points);
97 std::vector<UnitVector3d>::iterator const end = points.end();
98 UnitVector3d & v0 = points[0];
99 UnitVector3d & v1 = points[1];
100 // Note that robustCross() gives a non-zero result for distinct,
101 // non-antipodal inputs, and that normalization never maps a non-zero
102 // vector to the zero vector.
103 UnitVector3d n(v0.robustCross(v1));
104 for (; v != end; ++v) {
105 int ccw = orientation(v0, v1, *v);
106 if (ccw > 0) {
107 // We found a counter-clockwise triangle.
108 break;
109 } else if (ccw < 0) {
110 // We found a clockwise triangle. Swap the first two vertices to
111 // flip its orientation.
112 std::swap(v0, v1);
113 break;
114 }
115 // v, v0 and v1 are coplanar.
116 if (*v == v0 || *v == v1) {
117 continue;
118 }
119 if (*v == -v0 || *v == -v1) {
120 throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
121 }
122 // At this point, v, v0 and v1 are distinct and non-antipodal.
123 // If v is in the interior of the great circle segment (v0, v1),
124 // discard v and continue.
125 int v0v = orientation(n, v0, *v);
126 int vv1 = orientation(n, *v, v1);
127 if (v0v == vv1) {
128 continue;
129 }
130 int v0v1 = orientation(n, v0, v1);
131 // If v1 is in the interior of (v0, v), replace v1 with v and continue.
132 if (v0v1 == -vv1) {
133 v1 = *v; continue;
134 }
135 // If v0 is in the interior of (v, v1), replace v0 with v and continue.
136 if (-v0v == v0v1) {
137 v0 = *v; continue;
138 }
139 // Otherwise, (v0, v1) ∪ (v1, v) and (v, v0) ∪ (v0, v1) both span
140 // more than π radians of the great circle defined by the v0 and v1,
141 // so there is a pair of antipodal points in the corresponding great
142 // circle segment.
143 throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
144 }
145 if (v == end) {
146 throw std::invalid_argument(NOT_ENOUGH_POINTS);
147 }
148 points[2] = *v;
149 return ++v;
150}
151
152void computeHull(std::vector<UnitVector3d> & points) {
153 typedef std::vector<UnitVector3d>::iterator VertexIterator;
154 VertexIterator hullEnd = points.begin() + 3;
155 VertexIterator const end = points.end();
156 // Start with a triangular hull.
157 for (VertexIterator v = findTriangle(points); v != end; ++v) {
158 // Compute the hull of the current hull and v.
159 //
160 // Notice that if v is in the current hull, v can be ignored. If
161 // -v is in the current hull, then the hull of v and the current hull
162 // is not a convex polygon.
163 //
164 // Otherwise, let i and j be the first and second end-points of an edge
165 // in the current hull. We define the orientation of vertex j with
166 // respect to vertex v as orientation(v, i, j). When neither v or -v
167 // is in the current hull, there must be a sequence of consecutive
168 // hull vertices that do not have counter-clockwise orientation with
169 // respect to v. Insert v before the first vertex in this sequence
170 // and remove all vertices in the sequence except the last to obtain
171 // a new, larger convex hull.
172 VertexIterator i = hullEnd - 1;
173 VertexIterator j = points.begin();
174 // toCCW is the vertex before the transition to counter-clockwise
175 // orientation with respect to v.
176 VertexIterator toCCW = hullEnd;
177 // fromCCW is the vertex at which orientation with respect to v
178 // changes from counter-clockwise to clockwise or coplanar. It may
179 // equal toCCW.
180 VertexIterator fromCCW = hullEnd;
181 // Compute the orientation of the first point in the current hull
182 // with respect to v.
183 bool const firstCCW = orientation(*v, *i, *j) > 0;
184 bool prevCCW = firstCCW;
185 // Compute the orientation of points in the current hull with respect
186 // to v, starting with the second point. Update toCCW / fromCCW when
187 // we transition to / from counter-clockwise orientation.
188 for (i = j, ++j; j != hullEnd; i = j, ++j) {
189 if (orientation(*v, *i, *j) > 0) {
190 if (!prevCCW) {
191 toCCW = i;
192 prevCCW = true;
193 }
194 } else if (prevCCW) {
195 fromCCW = j;
196 prevCCW = false;
197 }
198 }
199 // Now that we know the orientation of the last point in the current
200 // hull with respect to v, consider the first point in the current hull.
201 if (firstCCW) {
202 if (!prevCCW) {
203 toCCW = i;
204 }
205 } else if (prevCCW) {
206 fromCCW = points.begin();
207 }
208 // Handle the case where there is never a transition to / from
209 // counter-clockwise orientation.
210 if (toCCW == hullEnd) {
211 // If all vertices are counter-clockwise with respect to v,
212 // v is inside the current hull and can be ignored.
213 if (firstCCW) {
214 continue;
215 }
216 // Otherwise, no vertex is counter-clockwise with respect to v,
217 // so that -v is inside the current hull.
218 throw std::invalid_argument(FOUND_ANTIPODAL_POINT);
219 }
220 // Insert v into the current hull at fromCCW, and remove
221 // all vertices between fromCCW and toCCW.
222 if (toCCW < fromCCW) {
223 if (toCCW != points.begin()) {
224 fromCCW = std::copy(toCCW, fromCCW, points.begin());
225 }
226 *fromCCW++ = *v;
227 hullEnd = fromCCW;
228 } else if (toCCW > fromCCW) {
229 *fromCCW++ = *v;
230 if (toCCW != fromCCW) {
231 hullEnd = std::copy(toCCW, hullEnd, fromCCW);
232 }
233 } else {
234 if (fromCCW == points.begin()) {
235 *hullEnd = *v;
236 } else {
237 UnitVector3d u = *v;
238 std::copy_backward(fromCCW, hullEnd, hullEnd + 1);
239 *fromCCW = u;
240 }
241 ++hullEnd;
242 }
243 }
244 points.erase(hullEnd, end);
245}
246
247// TODO(smm): for all of this to be fully rigorous, we must prove that no two
248// UnitVector3d objects u and v are exactly colinear unless u == v or u == -v.
249// It's not clear that this is true. For example, (1, 0, 0) and (1 + ε, 0, 0)
250// are colinear. This means that UnitVector3d should probably always normalize
251// on construction. Currently, it does not normalize when created from a LonLat,
252// and also contains some escape-hatches for performance. The normalize()
253// function implementation may also need to be revisited.
254
255// TODO(smm): This implementation is quadratic. It would be nice to implement
256// a fast hull merging algorithm, which could then be used to implement Chan's
257// algorithm.
258
259} // unnamed namespace
260
261
262ConvexPolygon::ConvexPolygon(std::vector<UnitVector3d> const & points) :
263 _vertices(points)
264{
265 computeHull(_vertices);
266}
267
269 if (this == &p) {
270 return true;
271 }
272 if (_vertices.size() != p._vertices.size()) {
273 return false;
274 }
275 VertexIterator i = _vertices.begin();
276 VertexIterator f = p._vertices.begin();
277 VertexIterator const ep = p._vertices.end();
278 // Find vertex f of p equal to the first vertex of this polygon.
279 for (; f != ep; ++f) {
280 if (*i == *f) {
281 break;
282 }
283 }
284 if (f == ep) {
285 // No vertex of p is equal to the first vertex of this polygon.
286 return false;
287 }
288 // Now, compare all vertices.
289 ++i;
290 for (VertexIterator j = f + 1; j != ep; ++i, ++j) {
291 if (*i != *j) {
292 return false;
293 }
294 }
295 for (VertexIterator j = p._vertices.begin(); j != f; ++i, ++j) {
296 if (*i != *j) {
297 return false;
298 }
299 }
300 return true;
301}
302
304 return detail::centroid(_vertices.begin(), _vertices.end());
305}
306
308 return detail::boundingCircle(_vertices.begin(), _vertices.end());
309}
310
312 return detail::boundingBox(_vertices.begin(), _vertices.end());
313}
314
316 return detail::boundingBox3d(_vertices.begin(), _vertices.end());
317}
318
320 return detail::contains(_vertices.begin(), _vertices.end(), v);
321}
322
323bool ConvexPolygon::contains(Region const & r) const {
324 return (relate(r) & CONTAINS) != 0;
325}
326
328 return (relate(r) & DISJOINT) != 0;
329}
330
331bool ConvexPolygon::intersects(Region const & r) const {
332 return !isDisjointFrom(r);
333}
334
335bool ConvexPolygon::isWithin(Region const & r) const {
336 return (relate(r) & WITHIN) != 0;
337}
338
340 return detail::relate(_vertices.begin(), _vertices.end(), b);
341}
342
344 return detail::relate(_vertices.begin(), _vertices.end(), c);
345}
346
348 return detail::relate(_vertices.begin(), _vertices.end(), p);
349}
350
352 return detail::relate(_vertices.begin(), _vertices.end(), e);
353}
354
357 uint8_t tc = TYPE_CODE;
358 buffer.reserve(1 + 24 * _vertices.size());
359 buffer.push_back(tc);
360 for (UnitVector3d const & v: _vertices) {
361 encodeDouble(v.x(), buffer);
362 encodeDouble(v.y(), buffer);
363 encodeDouble(v.z(), buffer);
364 }
365 return buffer;
366}
367
369 size_t n)
370{
371 if (buffer == nullptr || *buffer != TYPE_CODE ||
372 n < 1 + 24*3 || (n - 1) % 24 != 0) {
373 throw std::runtime_error("Byte-string is not an encoded ConvexPolygon");
374 }
376 ++buffer;
377 size_t nv = (n - 1) / 24;
378 poly->_vertices.reserve(nv);
379 for (size_t i = 0; i < nv; ++i, buffer += 24) {
380 poly->_vertices.push_back(UnitVector3d::fromNormalized(
381 decodeDouble(buffer),
382 decodeDouble(buffer + 8),
383 decodeDouble(buffer + 16)
384 ));
385 }
386 return poly;
387}
388
390 typedef std::vector<UnitVector3d>::const_iterator VertexIterator;
391 VertexIterator v = p.getVertices().begin();
392 VertexIterator const end = p.getVertices().end();
393 os << "{\"ConvexPolygon\": [" << *v;
394 for (++v; v != end; ++v) { os << ", " << *v; }
395 os << "]}";
396 return os;
397}
398
399}} // namespace lsst::sphgeom
int end
This file declares a class for representing convex polygons with great circle edges on the unit spher...
This file contains the meat of the ConvexPolygon implementation.
std::ostream * os
Definition Schema.cc:557
table::Key< int > b
T begin(T... args)
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
Circle is a circular region on the unit sphere that contains its boundary.
Definition Circle.h:53
ConvexPolygon is a closed convex polygon on the unit sphere.
Box3d getBoundingBox3d() const override
getBoundingBox3d returns a 3-dimensional bounding-box for this region.
bool intersects(Region const &r) const
UnitVector3d getCentroid() const
The centroid of a polygon is its center of mass projected onto S², assuming a uniform mass distributi...
bool contains(UnitVector3d const &v) const override
bool isDisjointFrom(Region const &r) const
Relationship relate(Region const &r) const override
Circle getBoundingCircle() const override
getBoundingCircle returns a bounding-circle for this region.
std::vector< uint8_t > encode() const override
encode serializes this region into an opaque byte string.
static constexpr uint8_t TYPE_CODE
Box getBoundingBox() const override
getBoundingBox returns a bounding-box for this region.
static std::unique_ptr< ConvexPolygon > decode(std::vector< uint8_t > const &s)
bool operator==(ConvexPolygon const &p) const
Two convex polygons are equal iff they contain the same points.
std::vector< UnitVector3d > const & getVertices() const
bool isWithin(Region const &r) const
Ellipse is an elliptical region on the sphere.
Definition Ellipse.h:177
Region is a minimal interface for 2-dimensional regions on the unit sphere.
Definition Region.h:86
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
static UnitVector3d fromNormalized(Vector3d const &v)
fromNormalized returns the unit vector equal to v, which is assumed to be normalized.
This file contains simple helper functions for encoding and decoding primitive types to/from byte str...
T copy_backward(T... args)
T copy(T... args)
T empty(T... args)
T end(T... args)
T erase(T... args)
Low-level polynomials (including special polynomials) in C++.
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)
void encodeDouble(double item, std::vector< uint8_t > &buffer)
encodeDouble appends an IEEE double in little-endian byte order to the end of buffer.
Definition codec.h:56
std::ostream & operator<<(std::ostream &, Angle const &)
Definition Angle.cc:41
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.
double decodeDouble(uint8_t const *buffer)
decodeDouble extracts an IEEE double from the 8 byte little-endian byte sequence in buffer.
Definition codec.h:76
T push_back(T... args)
T reserve(T... args)
This file declares functions for orienting points on the sphere.
T swap(T... args)