s2: Fix issues with Polygon.Invert.

The counters were not being properly cleared during the invert, so this
led to index out of bounds panics.

Also pulled in a few small changes to stay in sync with the C++ code for
Polycon.ContainsPoint.

Fixes #68.
Fixes #69.

Signed-off-by: David Symonds <dsymonds@golang.org>

Author: rsned

s2/loop.go

@@ -72,6 +72,7 @@ type Loop struct {
 func LoopFromPoints(pts []Point) *Loop {
 	l := &Loop{
 		vertices: pts,
+		index:    NewShapeIndex(),
 	}
 
 	l.initOriginAndBound()
@@ -96,6 +97,7 @@ func LoopFromCell(c Cell) *Loop {
 			c.Vertex(2),
 			c.Vertex(3),
 		},
+		index: NewShapeIndex(),
 	}
 
 	l.initOriginAndBound()
@@ -587,16 +589,24 @@ func (l *Loop) bruteForceContainsPoint(p Point) bool {
 
 // ContainsPoint returns true if the loop contains the point.
 func (l *Loop) ContainsPoint(p Point) bool {
-	// Empty and full loops don't need a special case, but invalid loops with
-	// zero vertices do, so we might as well handle them all at once.
-	if len(l.vertices) < 3 {
-		return l.originInside
+	if !l.index.IsFresh() && !l.bound.ContainsPoint(p) {
+		return false
 	}
 
-	// For small loops, and during initial construction, it is faster to just
-	// check all the crossing.
+	// For small loops it is faster to just check all the crossings.  We also
+	// use this method during loop initialization because InitOriginAndBound()
+	// calls Contains() before InitIndex().  Otherwise, we keep track of the
+	// number of calls to Contains() and only build the index when enough calls
+	// have been made so that we think it is worth the effort.  Note that the
+	// code below is structured so that if many calls are made in parallel only
+	// one thread builds the index, while the rest continue using brute force
+	// until the index is actually available.
+
 	const maxBruteForceVertices = 32
-	if len(l.vertices) < maxBruteForceVertices || l.index == nil {
+	// TODO(roberts): add unindexed contains calls tracking
+
+	if len(l.index.shapes) == 0 || // Index has not been initialized yet.
+		len(l.vertices) <= maxBruteForceVertices {
 		return l.bruteForceContainsPoint(p)
 	}
 
@@ -802,20 +812,23 @@ func (l *Loop) TurningAngle() float64 {
 		compensation = (oldSum - sum) + angle
 		n--
 	}
-	return float64(dir) * float64(sum+compensation)
+
+	const maxCurvature = 2*math.Pi - 4*dblEpsilon
+
+	return math.Max(-maxCurvature, math.Min(maxCurvature, float64(dir)*float64(sum+compensation)))
 }
 
 // turningAngleMaxError return the maximum error in TurningAngle. The value is not
 // constant; it depends on the loop.
 func (l *Loop) turningAngleMaxError() float64 {
 	// The maximum error can be bounded as follows:
-	//   2.24 * dblEpsilon    for RobustCrossProd(b, a)
-	//   2.24 * dblEpsilon    for RobustCrossProd(c, b)
+	//   3.00 * dblEpsilon    for RobustCrossProd(b, a)
+	//   3.00 * dblEpsilon    for RobustCrossProd(c, b)
 	//   3.25 * dblEpsilon    for Angle()
 	//   2.00 * dblEpsilon    for each addition in the Kahan summation
 	//   ------------------
-	//   9.73 * dblEpsilon
-	maxErrorPerVertex := 9.73 * dblEpsilon
+	//  11.25 * dblEpsilon
+	maxErrorPerVertex := 11.25 * dblEpsilon
 	return maxErrorPerVertex * float64(len(l.vertices))
 }
 
@@ -1275,12 +1288,12 @@ func (l *Loop) decode(d *decoder) {
 		l.vertices[i].Y = d.readFloat64()
 		l.vertices[i].Z = d.readFloat64()
 	}
+	l.index = NewShapeIndex()
 	l.originInside = d.readBool()
 	l.depth = int(d.readUint32())
 	l.bound.decode(d)
 	l.subregionBound = ExpandForSubregions(l.bound)
 
-	l.index = NewShapeIndex()
 	l.index.Add(l)
 }
 
@@ -1359,6 +1372,7 @@ func (l *Loop) decodeCompressed(d *decoder, snapLevel int) {
 		return
 	}
 
+	l.index = NewShapeIndex()
 	l.originInside = (properties & originInside) != 0
 
 	l.depth = int(d.readUvarint())
@@ -1373,7 +1387,6 @@ func (l *Loop) decodeCompressed(d *decoder, snapLevel int) {
 		l.initBound()
 	}
 
-	l.index = NewShapeIndex()
 	l.index.Add(l)
 }
 

s2/polygon.go

@@ -194,10 +194,12 @@ func (p *Polygon) Invert() {
 	// The empty and full polygons are handled specially.
 	if p.IsEmpty() {
 		*p = *FullPolygon()
+		p.initLoopProperties()
 		return
 	}
 	if p.IsFull() {
 		*p = Polygon{}
+		p.initLoopProperties()
 		return
 	}
 
@@ -244,6 +246,7 @@ func (p *Polygon) Invert() {
 			newLoops = append(newLoops, l)
 		}
 	}
+
 	p.loops = newLoops
 	p.initLoopProperties()
 }
@@ -383,6 +386,7 @@ func (p *Polygon) initOneLoop() {
 
 // initLoopProperties sets the properties for polygons with multiple loops.
 func (p *Polygon) initLoopProperties() {
+	p.numVertices = 0
 	// the loops depths are set by initNested/initOriented prior to this.
 	p.bound = EmptyRect()
 	p.hasHoles = false
@@ -402,6 +406,8 @@ func (p *Polygon) initLoopProperties() {
 // initEdgesAndIndex performs the shape related initializations and adds the final
 // polygon to the index.
 func (p *Polygon) initEdgesAndIndex() {
+	p.numEdges = 0
+	p.cumulativeEdges = nil
 	if p.IsFull() {
 		return
 	}
@@ -600,13 +606,8 @@ func (p *Polygon) ContainsPoint(point Point) bool {
 		return inside
 	}
 
-	// Otherwise, look up the ShapeIndex cell containing this point.
-	it := p.index.Iterator()
-	if !it.LocatePoint(point) {
-		return false
-	}
-
-	return p.iteratorContainsPoint(it, point)
+	// Otherwise we look up the ShapeIndex cell containing this point.
+	return NewContainsPointQuery(p.index, VertexModelSemiOpen).Contains(point)
 }
 
 // ContainsCell reports whether the polygon contains the given cell.

s2/polygon_test.go

@@ -1107,6 +1107,23 @@ func TestPolygonArea(t *testing.T) {
 	}
 }
 
+func TestPolygonInvert(t *testing.T) {
+	origin := PointFromLatLng(LatLngFromDegrees(0, 0))
+	pt := PointFromLatLng(LatLngFromDegrees(30, 30))
+	p := PolygonFromLoops([]*Loop{
+		RegularLoop(origin, 1000/earthRadiusKm, 100),
+	})
+
+	if p.ContainsPoint(pt) {
+		t.Errorf("polygon contains point outside itself")
+	}
+
+	p.Invert()
+	if !p.ContainsPoint(pt) {
+		t.Errorf("inverted polygon does not contain point that is inside itself")
+	}
+}
+
 // TODO(roberts): Remaining Tests
 // TestInit
 // TestMultipleInit