1 // Copyright 2018 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s2 16 17 import ( 18 "sort" 19 20 "github.com/golang/geo/r3" 21 ) 22 23 // ConvexHullQuery builds the convex hull of any collection of points, 24 // polylines, loops, and polygons. It returns a single convex loop. 25 // 26 // The convex hull is defined as the smallest convex region on the sphere that 27 // contains all of your input geometry. Recall that a region is "convex" if 28 // for every pair of points inside the region, the straight edge between them 29 // is also inside the region. In our case, a "straight" edge is a geodesic, 30 // i.e. the shortest path on the sphere between two points. 31 // 32 // Containment of input geometry is defined as follows: 33 // 34 // - Each input loop and polygon is contained by the convex hull exactly 35 // (i.e., according to Polygon's Contains(Polygon)). 36 // 37 // - Each input point is either contained by the convex hull or is a vertex 38 // of the convex hull. (Recall that S2Loops do not necessarily contain their 39 // vertices.) 40 // 41 // - For each input polyline, the convex hull contains all of its vertices 42 // according to the rule for points above. (The definition of convexity 43 // then ensures that the convex hull also contains the polyline edges.) 44 // 45 // To use this type, call the various Add... methods to add your input geometry, and 46 // then call ConvexHull. Note that ConvexHull does *not* reset the 47 // state; you can continue adding geometry if desired and compute the convex 48 // hull again. If you want to start from scratch, simply create a new 49 // ConvexHullQuery value. 50 // 51 // This implement Andrew's monotone chain algorithm, which is a variant of the 52 // Graham scan (see https://en.wikipedia.org/wiki/Graham_scan). The time 53 // complexity is O(n log n), and the space required is O(n). In fact only the 54 // call to "sort" takes O(n log n) time; the rest of the algorithm is linear. 55 // 56 // Demonstration of the algorithm and code: 57 // en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain 58 // 59 // This type is not safe for concurrent use. 60 type ConvexHullQuery struct { 61 bound Rect 62 points []Point 63 } 64 65 // NewConvexHullQuery creates a new ConvexHullQuery. 66 func NewConvexHullQuery() *ConvexHullQuery { 67 return &ConvexHullQuery{ 68 bound: EmptyRect(), 69 } 70 } 71 72 // AddPoint adds the given point to the input geometry. 73 func (q *ConvexHullQuery) AddPoint(p Point) { 74 q.bound = q.bound.AddPoint(LatLngFromPoint(p)) 75 q.points = append(q.points, p) 76 } 77 78 // AddPolyline adds the given polyline to the input geometry. 79 func (q *ConvexHullQuery) AddPolyline(p *Polyline) { 80 q.bound = q.bound.Union(p.RectBound()) 81 q.points = append(q.points, (*p)...) 82 } 83 84 // AddLoop adds the given loop to the input geometry. 85 func (q *ConvexHullQuery) AddLoop(l *Loop) { 86 q.bound = q.bound.Union(l.RectBound()) 87 if l.isEmptyOrFull() { 88 return 89 } 90 q.points = append(q.points, l.vertices...) 91 } 92 93 // AddPolygon adds the given polygon to the input geometry. 94 func (q *ConvexHullQuery) AddPolygon(p *Polygon) { 95 q.bound = q.bound.Union(p.RectBound()) 96 for _, l := range p.loops { 97 // Only loops at depth 0 can contribute to the convex hull. 98 if l.depth == 0 { 99 q.AddLoop(l) 100 } 101 } 102 } 103 104 // CapBound returns a bounding cap for the input geometry provided. 105 // 106 // Note that this method does not clear the geometry; you can continue 107 // adding to it and call this method again if desired. 108 func (q *ConvexHullQuery) CapBound() Cap { 109 // We keep track of a rectangular bound rather than a spherical cap because 110 // it is easy to compute a tight bound for a union of rectangles, whereas it 111 // is quite difficult to compute a tight bound around a union of caps. 112 // Also, polygons and polylines implement CapBound() in terms of 113 // RectBound() for this same reason, so it is much better to keep track 114 // of a rectangular bound as we go along and convert it at the end. 115 // 116 // TODO(roberts): We could compute an optimal bound by implementing Welzl's 117 // algorithm. However we would still need to have special handling of loops 118 // and polygons, since if a loop spans more than 180 degrees in any 119 // direction (i.e., if it contains two antipodal points), then it is not 120 // enough just to bound its vertices. In this case the only convex bounding 121 // cap is FullCap(), and the only convex bounding loop is the full loop. 122 return q.bound.CapBound() 123 } 124 125 // ConvexHull returns a Loop representing the convex hull of the input geometry provided. 126 // 127 // If there is no geometry, this method returns an empty loop containing no 128 // points. 129 // 130 // If the geometry spans more than half of the sphere, this method returns a 131 // full loop containing the entire sphere. 132 // 133 // If the geometry contains 1 or 2 points, or a single edge, this method 134 // returns a very small loop consisting of three vertices (which are a 135 // superset of the input vertices). 136 // 137 // Note that this method does not clear the geometry; you can continue 138 // adding to the query and call this method again. 139 func (q *ConvexHullQuery) ConvexHull() *Loop { 140 c := q.CapBound() 141 if c.Height() >= 1 { 142 // The bounding cap is not convex. The current bounding cap 143 // implementation is not optimal, but nevertheless it is likely that the 144 // input geometry itself is not contained by any convex polygon. In any 145 // case, we need a convex bounding cap to proceed with the algorithm below 146 // (in order to construct a point "origin" that is definitely outside the 147 // convex hull). 148 return FullLoop() 149 } 150 151 // Remove duplicates. We need to do this before checking whether there are 152 // fewer than 3 points. 153 x := make(map[Point]bool) 154 r, w := 0, 0 // read/write indexes 155 for ; r < len(q.points); r++ { 156 if x[q.points[r]] { 157 continue 158 } 159 q.points[w] = q.points[r] 160 x[q.points[r]] = true 161 w++ 162 } 163 q.points = q.points[:w] 164 165 // This code implements Andrew's monotone chain algorithm, which is a simple 166 // variant of the Graham scan. Rather than sorting by x-coordinate, instead 167 // we sort the points in CCW order around an origin O such that all points 168 // are guaranteed to be on one side of some geodesic through O. This 169 // ensures that as we scan through the points, each new point can only 170 // belong at the end of the chain (i.e., the chain is monotone in terms of 171 // the angle around O from the starting point). 172 origin := Point{c.Center().Ortho()} 173 sort.Slice(q.points, func(i, j int) bool { 174 return RobustSign(origin, q.points[i], q.points[j]) == CounterClockwise 175 }) 176 177 // Special cases for fewer than 3 points. 178 switch len(q.points) { 179 case 0: 180 return EmptyLoop() 181 case 1: 182 return singlePointLoop(q.points[0]) 183 case 2: 184 return singleEdgeLoop(q.points[0], q.points[1]) 185 } 186 187 // Generate the lower and upper halves of the convex hull. Each half 188 // consists of the maximal subset of vertices such that the edge chain 189 // makes only left (CCW) turns. 190 lower := q.monotoneChain() 191 192 // reverse the points 193 for left, right := 0, len(q.points)-1; left < right; left, right = left+1, right-1 { 194 q.points[left], q.points[right] = q.points[right], q.points[left] 195 } 196 upper := q.monotoneChain() 197 198 // Remove the duplicate vertices and combine the chains. 199 lower = lower[:len(lower)-1] 200 upper = upper[:len(upper)-1] 201 lower = append(lower, upper...) 202 203 return LoopFromPoints(lower) 204 } 205 206 // monotoneChain iterates through the points, selecting the maximal subset of points 207 // such that the edge chain makes only left (CCW) turns. 208 func (q *ConvexHullQuery) monotoneChain() []Point { 209 var output []Point 210 for _, p := range q.points { 211 // Remove any points that would cause the chain to make a clockwise turn. 212 for len(output) >= 2 && RobustSign(output[len(output)-2], output[len(output)-1], p) != CounterClockwise { 213 output = output[:len(output)-1] 214 } 215 output = append(output, p) 216 } 217 return output 218 } 219 220 // singlePointLoop constructs a 3-vertex polygon consisting of "p" and two nearby 221 // vertices. Note that ContainsPoint(p) may be false for the resulting loop. 222 func singlePointLoop(p Point) *Loop { 223 const offset = 1e-15 224 d0 := p.Ortho() 225 d1 := p.Cross(d0) 226 vertices := []Point{ 227 p, 228 {p.Add(d0.Mul(offset)).Normalize()}, 229 {p.Add(d1.Mul(offset)).Normalize()}, 230 } 231 return LoopFromPoints(vertices) 232 } 233 234 // singleEdgeLoop constructs a loop consisting of the two vertices and their midpoint. 235 func singleEdgeLoop(a, b Point) *Loop { 236 // If the points are exactly antipodal we return the full loop. 237 // 238 // Note that we could use the code below even in this case (which would 239 // return a zero-area loop that follows the edge AB), except that (1) the 240 // direction of AB is defined using symbolic perturbations and therefore is 241 // not predictable by ordinary users, and (2) Loop disallows anitpodal 242 // adjacent vertices and so we would need to use 4 vertices to define the 243 // degenerate loop. (Note that the Loop antipodal vertex restriction is 244 // historical and now could easily be removed, however it would still have 245 // the problem that the edge direction is not easily predictable.) 246 if a.Add(b.Vector) == (r3.Vector{}) { 247 return FullLoop() 248 } 249 250 // Construct a loop consisting of the two vertices and their midpoint. We 251 // use Interpolate() to ensure that the midpoint is very close to 252 // the edge even when its endpoints nearly antipodal. 253 vertices := []Point{a, b, Interpolate(0.5, a, b)} 254 loop := LoopFromPoints(vertices) 255 // The resulting loop may be clockwise, so invert it if necessary. 256 loop.Normalize() 257 return loop 258 } 259