logician's updates

r1/interval.go

@@ -34,6 +34,23 @@ func EmptyInterval() Interval { return Interval{1, 0} }
 // IntervalFromPoint returns an interval representing a single point.
 func IntervalFromPoint(p float64) Interval { return Interval{p, p} }
 
+// Convenience method to construct an interval from two points that are
+// already ordered from lo to hi
+func IntervalFromEndpoints(p1, p2 float64) Interval {
+	return Interval{p1, p2}
+}
+
+// Convenience method to construct the minimal interval containing the two
+// given points. This is equivalent to starting with an empty interval and
+// calling AddPoint() twice, but it is more efficient.
+func IntervalFromPointPair(p1, p2 float64) Interval {
+	if p1 <= p2 {
+		return Interval{p1, p2}
+	} else {
+		return Interval{p2, p1}
+	}
+}
+
 // IsEmpty reports whether the interval is empty.
 func (i Interval) IsEmpty() bool { return i.Lo > i.Hi }
 

r2/vector.go

@@ -0,0 +1,104 @@
+/*
+ * Copyright 2005 Google Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package r2
+
+import (
+	"fmt"
+	"math"
+)
+
+/**
+ * r2.Vector represents a vector in the two-dimensional space. It defines the
+ * basic geometrical operations for 2D vectors, e.g. cross product, addition,
+ * norm, comparison etc.
+ *
+ */
+type Vector struct {
+	X, Y float64
+}
+
+func (v Vector) String() string { return fmt.Sprintf("(%v, %v)", v.X, v.Y) }
+
+// Norm returns the vector's norm.
+func (v Vector) Norm() float64 { return math.Sqrt(v.Dot(v)) }
+
+// Norm2 returns the square of the norm.
+func (v Vector) Norm2() float64 { return v.Dot(v) }
+
+// Normalize returns a unit vector in the same direction as v.
+func (v Vector) Normalize() Vector {
+	if v == (Vector{0, 0}) {
+		return v
+	}
+	return v.Mul(1 / v.Norm())
+}
+
+// Abs returns the vector with nonnegative components.
+func (v Vector) Abs() Vector { return Vector{math.Abs(v.X), math.Abs(v.Y)} }
+
+// Neg returns the negated vector
+func (v Vector) Neg() Vector { return Vector{-v.X, -v.Y} }
+
+// Add returns the standard vector sum of v and ov.
+func (v Vector) Add(ov Vector) Vector { return Vector{v.X + ov.X, v.Y + ov.Y} }
+
+// Sub returns the standard vector difference of v and ov.
+func (v Vector) Sub(ov Vector) Vector { return Vector{v.X - ov.X, v.Y - ov.Y} }
+
+// Mul returns the standard scalar product of v and m.
+func (v Vector) Mul(m float64) Vector { return Vector{v.X * m, v.Y * m} }
+
+// Mul returns the standard scalar product of v and m.
+func (v Vector) Div(m float64) Vector { return Vector{v.X / m, v.Y / m} }
+
+// Dot returns the standard dot product of v and ov.
+func (v Vector) Dot(ov Vector) float64 { return v.X*ov.X + v.Y*ov.Y }
+
+// Cross returns the standard cross product of v and ov.
+func (v Vector) Cross(ov Vector) float64 {
+	return v.X*ov.Y - v.Y*ov.X
+}
+
+func (v Vector) Equals(other Vector) bool {
+	return v.X == other.X && v.Y == other.Y
+}
+
+func (v Vector) LessThan(vb Vector) bool {
+	if v.X < vb.X {
+		return true
+	}
+	if vb.X < v.X {
+		return false
+	}
+	if v.Y < vb.Y {
+		return true
+	}
+	if vb.Y < v.Y {
+		return false
+	}
+	return false
+}
+
+func (v Vector) CompareTo(other Vector) int {
+	if v.LessThan(other) {
+		return -1
+	} else {
+		if v.Equals(other) {
+			return 0
+		}
+		return 1
+	}
+}

r3/vector.go

@@ -28,6 +28,37 @@ type Vector struct {
 	X, Y, Z float64
 }
 
+func (v Vector) GetAxis(axis int) float64 {
+	switch axis {
+	case 0:
+		return v.X
+	case 1:
+		return v.Y
+	case 2:
+		return v.Z
+	default:
+		return v.GetAxis((axis + 3) % 3)
+	}
+}
+
+/** Return the index of the largest component fabs */
+func (v Vector) LargestAbsComponent() int {
+	temp := v.Abs()
+	if temp.X > temp.Y {
+		if temp.X > temp.Z {
+			return 0
+		} else {
+			return 2
+		}
+	} else {
+		if temp.Y > temp.Z {
+			return 1
+		} else {
+			return 2
+		}
+	}
+}
+
 // ApproxEqual reports whether v and ov are equal within a small epsilon.
 func (v Vector) ApproxEqual(ov Vector) bool {
 	const epsilon = 1e-14
@@ -59,14 +90,20 @@ func (v Vector) IsUnit() bool {
 // Abs returns the vector with nonnegative components.
 func (v Vector) Abs() Vector { return Vector{math.Abs(v.X), math.Abs(v.Y), math.Abs(v.Z)} }
 
+// Neg returns the negated vector
+func (v Vector) Neg() Vector { return Vector{-v.X, -v.Y, -v.Z} }
+
 // Add returns the standard vector sum of v and ov.
 func (v Vector) Add(ov Vector) Vector { return Vector{v.X + ov.X, v.Y + ov.Y, v.Z + ov.Z} }
 
 // Sub returns the standard vector difference of v and ov.
 func (v Vector) Sub(ov Vector) Vector { return Vector{v.X - ov.X, v.Y - ov.Y, v.Z - ov.Z} }
 
 // Mul returns the standard scalar product of v and m.
-func (v Vector) Mul(m float64) Vector { return Vector{m * v.X, m * v.Y, m * v.Z} }
+func (v Vector) Mul(m float64) Vector { return Vector{v.X * m, v.Y * m, v.Z * m} }
+
+// Mul returns the standard scalar product of v and m.
+func (v Vector) Div(m float64) Vector { return Vector{v.X / m, v.Y / m, v.Z / m} }
 
 // Dot returns the standard dot product of v and ov.
 func (v Vector) Dot(ov Vector) float64 { return v.X*ov.X + v.Y*ov.Y + v.Z*ov.Z }

s1/interval.go

@@ -47,6 +47,24 @@ func IntervalFromEndpoints(lo, hi float64) Interval {
 	return i
 }
 
+// Convenience method to construct the minimal interval containing the two
+// given points. This is equivalent to starting with an empty interval and
+// calling AddPoint() twice, but it is more efficient.
+func IntervalFromPointPair(p1, p2 float64) Interval {
+	// assert (Math.abs(p1) <= S2.M_PI && Math.abs(p2) <= S2.M_PI);
+	if p1 == -math.Pi {
+		p1 = math.Pi
+	}
+	if p2 == -math.Pi {
+		p2 = math.Pi
+	}
+	if positiveDistance(p1, p2) <= math.Pi {
+		return Interval{p1, p2}
+	} else {
+		return Interval{p2, p1}
+	}
+}
+
 // EmptyInterval returns an empty interval.
 func EmptyInterval() Interval { return Interval{math.Pi, -math.Pi} }
 

s2/areacentroid.go

@@ -0,0 +1,33 @@
+package s2
+
+/**
+ * The area of an interior, i.e. the region on the left side of an odd
+ * number of loops and optionally a centroid.
+ * The area is between 0 and 4*Pi. If it has a centroid, it is
+ * the true centroid of the interiord multiplied by the area of the shape.
+ * Note that the centroid may not be contained by the shape.
+ *
+ * @author dbentley@google.com (Daniel Bentley)
+ */
+type AreaCentroid struct {
+	area     float64
+	centroid *Point
+}
+
+func NewAreaCentroid(area float64, centroid *Point) AreaCentroid {
+	return AreaCentroid{
+		area:     area,
+		centroid: centroid,
+	}
+}
+
+func (ac AreaCentroid) GetArea() float64 {
+	return ac.area
+}
+
+func (ac AreaCentroid) GetCentroid() Point {
+	if ac.centroid == nil {
+		panic("no centroid")
+	}
+	return *ac.centroid
+}

s2/cap.go

@@ -118,6 +118,11 @@ func (c Cap) IsFull() bool {
 	return c.height == fullHeight
 }
 
+// Center returns the cap's center or axis.
+func (c Cap) Center() Point {
+	return c.center
+}
+
 // Radius returns the cap's radius.
 func (c Cap) Radius() s1.Angle {
 	if c.IsEmpty() {
@@ -245,7 +250,7 @@ func (c Cap) RectBound() Rect {
 // a reasonable epsilon from the other cap.
 func (c Cap) ApproxEqual(other Cap) bool {
 	const epsilon = 1e-14
-	return c.center.ApproxEqual(other.center) &&
+	return c.center.ApproxEquals(other.center, epsilon) &&
 		math.Abs(c.height-other.height) <= epsilon ||
 		c.IsEmpty() && other.height <= epsilon ||
 		other.IsEmpty() && c.height <= epsilon ||
@@ -312,6 +317,88 @@ func radiusToHeight(r s1.Angle) float64 {
 
 }
 
+// ContainsCell reports whether the region completely contains the given region.
+// It returns false if containment could not be determined.
+func (c Cap) ContainsCell(cell Cell) bool {
+	vertices := make([]Point, 4)
+	for k := 0; k < 4; k++ {
+		vertices[k] = cell.Vertex(k)
+		if !c.ContainsPoint(vertices[k]) {
+			return false
+		}
+	}
+	return !c.Complement().intersects(cell, vertices)
+}
+
+// IntersectsCell reports whether the region intersects the given cell or
+// if intersection could not be determined. It returns false if the region
+// does not intersect.
+func (c Cap) IntersectsCell(cell Cell) bool {
+	vertices := make([]Point, 4)
+	for k := 0; k < 4; k++ {
+		vertices[k] = cell.Vertex(k)
+		if c.ContainsPoint(vertices[k]) {
+			return true
+		}
+	}
+	return c.intersects(cell, vertices)
+}
+
+/**
+ * Return true if the cap intersects 'cell', given that the cap vertices have
+ * alrady been checked.
+ */
+func (c Cap) intersects(cell Cell, vertices []Point) bool {
+	// Return true if this cap intersects any point of 'cell' excluding its
+	// vertices (which are assumed to already have been checked).
+
+	// If the cap is a hemisphere or larger, the cell and the complement of the
+	// cap are both convex. Therefore since no vertex of the cell is contained,
+	// no other interior point of the cell is contained either.
+	if c.height >= 1 {
+		return false
+	}
+
+	// We need to check for empty caps due to the axis check just below.
+	if c.IsEmpty() {
+		return false
+	}
+
+	// Optimization: return true if the cell contains the cap axis. (This
+	// allows half of the edge checks below to be skipped.)
+	if cell.ContainsPoint(c.center) {
+		return true
+	}
+
+	// At this point we know that the cell does not contain the cap axis,
+	// and the cap does not contain any cell vertex. The only way that they
+	// can intersect is if the cap intersects the interior of some edge.
+
+	sin2Angle := c.height * (2 - c.height) // sin^2(capAngle)
+	for k := 0; k < 4; k++ {
+		edge := cell.EdgeRaw(k)
+		dot := c.center.Dot(edge.Vector)
+		if dot > 0 {
+			// The axis is in the interior half-space defined by the edge. We don't
+			// need to consider these edges, since if the cap intersects this edge
+			// then it also intersects the edge on the opposite side of the cell
+			// (because we know the axis is not contained with the cell).
+			continue
+		}
+		// The Norm2() factor is necessary because "edge" is not normalized.
+		if dot*dot > sin2Angle*edge.Norm2() {
+			return false // Entire cap is on the exterior side of this edge.
+		}
+		// Otherwise, the great circle containing this edge intersects
+		// the interior of the cap. We just need to check whether the point
+		// of closest approach occurs between the two edge endpoints.
+		dir := edge.Cross(c.center.Vector)
+		if dir.Dot(vertices[k].Vector) < 0 && dir.Dot(vertices[(k+1)&3].Vector) > 0 {
+			return true
+		}
+	}
+	return false
+}
+
 // TODO(roberts): Differences from C++
-// Intersects(S2Cell), Contains(S2Cell), MayIntersect(S2Cell)
 // Centroid, Union