More additions

Author: alex-alex

Sources/R1Interval.swift

@@ -26,20 +26,20 @@ public struct R1Interval: Equatable { //, Hashable {
 	public static let empty = R1Interval(lo: 1, hi: 0)
 
 	/// Convenience method to construct an interval containing a single point.
-	public static func from(point p: Double) -> R1Interval {
-		return R1Interval(lo: p, hi: p)
+	public init(point p: Double) {
+		self.init(lo: p, hi: p)
 	}
 
 	/**
 		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.
 	*/
-	public static func fromPointPair(p1: Double, p2: Double) -> R1Interval {
+	public init(p1: Double, p2: Double) {
 		if p1 <= p2 {
-			return R1Interval(lo: p1, hi: p2)
+			self.init(lo: p1, hi: p2)
 		} else {
-			return R1Interval(lo: p2, hi: p1)
+			self.init(lo: p2, hi: p1)
 		}
 	}
 
@@ -95,7 +95,7 @@ public struct R1Interval: Equatable { //, Hashable {
 	/// Expand the interval so that it contains the given point "p".
 	public func add(point p: Double) -> R1Interval {
 		if isEmpty {
-			return R1Interval.from(point: p)
+			return R1Interval(point: p)
 		} else if p < lo {
 			return R1Interval(lo: p, hi: hi)
 		} else if p > hi {

Sources/S1Interval.swift

@@ -223,7 +223,7 @@ public struct S1Interval {
 		[2,Pi] intersect, for example.
 	*/
 	public func intersects(with y: S1Interval) -> Bool {
-		if (isEmpty || y.isEmpty) {
+		if isEmpty || y.isEmpty {
 			return false
 		}
 		if isInverted {

Sources/S2.swift

@@ -20,10 +20,6 @@ internal extension Double {
 
 public struct S2 {
 
-	enum Error: ErrorProtocol {
-		case IllegalArgumentException
-	}
-	
 	// Together these flags define a cell orientation. If SWAP_MASK
 	// is true, then canonical traversal order is flipped around the
 	// diagonal (i.e. i and j are swapped with each other). If
@@ -68,8 +64,8 @@ public struct S2 {
 
 		- Returns: A bit mask containing some combination of {@link #SWAP_MASK} and {@link #INVERT_MASK}.
 	*/
-	public static func posToOrientation(position: Int) throws -> Int {
-		guard 0 <= position && position < 4 else { throw Error.IllegalArgumentException }
+	public static func posToOrientation(position: Int) -> Int {
+		precondition(0 <= position && position < 4)
 		return posToOrientation[position]
 	}
 
@@ -93,9 +89,9 @@ public struct S2 {
 
 		- Returns: The IJ-index where `0->(0,0), 1->(0,1), 2->(1,0), 3->(1,1)`.
 	*/
-	public static func posToIJ(orientation: Int, position: Int) throws -> Int {
-		guard 0 <= orientation && orientation < 4 else { throw Error.IllegalArgumentException }
-		guard 0 <= position && position < 4 else { throw Error.IllegalArgumentException }
+	public static func posToIJ(orientation: Int, position: Int) -> Int {
+		precondition(0 <= orientation && orientation < 4)
+		precondition(0 <= position && position < 4)
 		return posToIJ[orientation][position]
 	}
 
@@ -119,9 +115,9 @@ public struct S2 {
 
 		- Returns: The position of the subcell in the Hilbert traversal, in the range [0,3].
 	*/
-	public static func toPos(orientation: Int, ijIndex: Int) throws -> Int {
-		guard 0 <= orientation && orientation < 4 else { throw Error.IllegalArgumentException }
-		guard 0 <= ijIndex && ijIndex < 4 else { throw Error.IllegalArgumentException }
+	public static func toPos(orientation: Int, ijIndex: Int) -> Int {
+		precondition(0 <= orientation && orientation < 4)
+		precondition(0 <= ijIndex && ijIndex < 4)
 		return ijToPos[orientation][ijIndex]
 	}
 
@@ -310,7 +306,7 @@ public struct S2 {
 		let sa = b.angle(to: c)
 		let sb = c.angle(to: a)
 		let sc = a.angle(to: b)
-		let s = 0.5 * (sa + sb + sc);
+		let s = 0.5 * (sa + sb + sc)
 		if s >= 3e-4 {
 			// Consider whether Girard's formula might be more accurate.
 			let s2 = s * s
@@ -324,12 +320,7 @@ public struct S2 {
 			}
 		}
 		// Use l'Huilier's formula.
-		return 4
-		* atan(
-		sqrt(
-		max(0.0,
-		tan(0.5 * s) * tan(0.5 * (s - sa)) * tan(0.5 * (s - sb))
-		* tan(0.5 * (s - sc)))))
+		return 4 * atan(sqrt(max(0.0, tan(0.5 * s) * tan(0.5 * (s - sa)) * tan(0.5 * (s - sb)) * tan(0.5 * (s - sc)))))
 	}
 
 	/**
@@ -353,6 +344,77 @@ public struct S2 {
 		return area(a: a, b: b, c: c) * Double(robustCCW(a: a, b: b, c: c))
 	}
 
+	// About centroids:
+	// ----------------
+	//
+	// There are several notions of the "centroid" of a triangle. First, there
+	// // is the planar centroid, which is simply the centroid of the ordinary
+	// (non-spherical) triangle defined by the three vertices. Second, there is
+	// the surface centroid, which is defined as the intersection of the three
+	// medians of the spherical triangle. It is possible to show that this
+	// point is simply the planar centroid projected to the surface of the
+	// sphere. Finally, there is the true centroid (mass centroid), which is
+	// defined as the area integral over the spherical triangle of (x,y,z)
+	// divided by the triangle area. This is the point that the triangle would
+	// rotate around if it was spinning in empty space.
+	//
+	// The best centroid for most purposes is the true centroid. Unlike the
+	// planar and surface centroids, the true centroid behaves linearly as
+	// regions are added or subtracted. That is, if you split a triangle into
+	// pieces and compute the average of their centroids (weighted by triangle
+	// area), the result equals the centroid of the original triangle. This is
+	// not true of the other centroids.
+	//
+	// Also note that the surface centroid may be nowhere near the intuitive
+	// "center" of a spherical triangle. For example, consider the triangle
+	// with vertices A=(1,eps,0), B=(0,0,1), C=(-1,eps,0) (a quarter-sphere).
+	// The surface centroid of this triangle is at S=(0, 2*eps, 1), which is
+	// within a distance of 2*eps of the vertex B. Note that the median from A
+	// (the segment connecting A to the midpoint of BC) passes through S, since
+	// this is the shortest path connecting the two endpoints. On the other
+	// hand, the true centroid is at M=(0, 0.5, 0.5), which when projected onto
+	// the surface is a much more reasonable interpretation of the "center" of
+	// this triangle.
+	
+	/**
+		Return the centroid of the planar triangle ABC. This can be normalized to
+		unit length to obtain the "surface centroid" of the corresponding spherical
+		triangle, i.e. the intersection of the three medians. However, note that
+		for large spherical triangles the surface centroid may be nowhere near the
+		intuitive "center" (see example above).
+	*/
+	public static func planarCentroid(a: S2Point, b: S2Point, c: S2Point) -> S2Point {
+		return S2Point(x: (a.x + b.x + c.x) / 3.0, y: (a.y + b.y + c.y) / 3.0, z: (a.z + b.z + c.z) / 3.0)
+	}
+	
+	/**
+		Returns the true centroid of the spherical triangle ABC multiplied by the
+		signed area of spherical triangle ABC. The reasons for multiplying by the
+		signed area are (1) this is the quantity that needs to be summed to compute
+		the centroid of a union or difference of triangles, and (2) it's actually
+		easier to calculate this way.
+	*/
+	public static func trueCentroid(a: S2Point, b: S2Point, c: S2Point) -> S2Point {
+		// I couldn't find any references for computing the true centroid of a
+		// spherical triangle... I have a truly marvellous demonstration of this
+		// formula which this margin is too narrow to contain :)
+		
+		// assert (isUnitLength(a) && isUnitLength(b) && isUnitLength(c));
+		let sina = b.crossProd(c).norm
+		let sinb = c.crossProd(a).norm
+		let sinc = a.crossProd(b).norm
+		let ra = (sina == 0) ? 1 : (asin(sina) / sina)
+		let rb = (sinb == 0) ? 1 : (asin(sinb) / sinb)
+		let rc = (sinc == 0) ? 1 : (asin(sinc) / sinc)
+		
+		// Now compute a point M such that M.X = rX * det(ABC) / 2 for X in A,B,C.
+		let x = S2Point(x: a.x, y: b.x, z: c.x)
+		let y = S2Point(x: a.y, y: b.y, z: c.y)
+		let z = S2Point(x: a.z, y: b.z, z: c.z)
+		let r = S2Point(x: ra, y: rb, z: rc)
+		return S2Point(x: 0.5 * y.crossProd(z).dotProd(r), y: 0.5 * z.crossProd(x).dotProd(r), z: 0.5 * x.crossProd(y).dotProd(r))
+	}
+	
 	/**
 		Return true if the points A, B, C are strictly counterclockwise. Return
 		false if the points are clockwise or colinear (i.e. if they are all
@@ -632,6 +694,16 @@ public struct S2 {
 		return (robustCCW(a: a, b: b, c: c) > 0) ? outAngle : -outAngle
 	}
 
+	/// Return true if two points are within the given distance of each other (mainly useful for testing).
+	public static func approxEquals(a: S2Point, b: S2Point, maxError: Double = 1e-15) -> Bool {
+		return a.angle(to: b) <= maxError
+	}
+	
+	/// Return true if two points are within the given distance of each other (mainly useful for testing).
+	public static func approxEquals(a: Double, b: Double, maxError: Double = 1e-15) -> Bool {
+		return abs(a - b) <= maxError
+	}
+	
 	// Don't instantiate
 	private init() { }
 

Sources/S2AreaCentroid.swift

@@ -16,9 +16,9 @@
 public struct S2AreaCentroid {
 
 	public let area: Double
-	public let centroid: S2Point?
+	public let centroid: S2Point
 
-	public init(area: Double, centroid: S2Point? = nil) {
+	public init(area: Double, centroid: S2Point) {
 		self.area = area
 		self.centroid = centroid
 	}

Sources/S2Cell.swift

@@ -36,16 +36,16 @@ public struct S2Cell: S2Region, Equatable {
 	}
 
 	/// An S2Cell always corresponds to a particular S2CellId. The other constructors are just convenience methods.
-	public init(id: S2CellId) {
-		cellId = id
+	public init(cellId: S2CellId) {
+		self.cellId = cellId
 
 		var i = 0
 		var j = 0
 		var mOrientation: Int? = 0
 
-		face = UInt8(id.toFaceIJOrientation(i: &i, j: &j, orientation: &mOrientation))
+		face = UInt8(cellId.toFaceIJOrientation(i: &i, j: &j, orientation: &mOrientation))
 		orientation = UInt8(mOrientation!)
-		level = UInt8(id.level)
+		level = UInt8(cellId.level)
 
 		let cellSize = 1 << (S2CellId.maxLevel - Int(level))
 		var _uv: [[Double]] = [[0, 0], [0, 0]]
@@ -61,16 +61,16 @@ public struct S2Cell: S2Region, Equatable {
 
 	// This is a static method in order to provide named parameters.
 	public init(face: Int, pos: UInt8, level: Int) {
-		self.init(id: S2CellId(face: face, pos: Int64(pos), level: level))
+		self.init(cellId: S2CellId(face: face, pos: Int64(pos), level: level))
 	}
 
 	// Convenience methods.
 	public init(point: S2Point) {
-		self.init(id: S2CellId(point: point))
+		self.init(cellId: S2CellId(point: point))
 	}
 
 	public init(latlng: S2LatLng) {
-		self.init(id: S2CellId(latlng: latlng))
+		self.init(cellId: S2CellId(latlng: latlng))
 	}
 
 	public var isLeaf: Bool {
@@ -121,7 +121,7 @@ public struct S2Cell: S2Region, Equatable {
 
 		except that it is more than two times faster.
 	*/
-	public func subdivide() throws -> [S2Cell] {
+	public func subdivide() -> [S2Cell] {
 		// This function is equivalent to just iterating over the child cell ids
 		// and calling the S2Cell constructor, but it is about 2.5 times faster.
 
@@ -136,15 +136,15 @@ public struct S2Cell: S2Region, Equatable {
 		for pos in 0 ..< 4 {
 
 			var _uv: [[Double]] = [[0, 0], [0, 0]]
-			let ij = try S2.posToIJ(orientation: Int(orientation), position: pos)
+			let ij = S2.posToIJ(orientation: Int(orientation), position: pos)
 
 			for d in 0 ..< 2 {
 				// The dimension 0 index (i/u) is in bit 1 of ij.
 				let m = 1 - ((ij >> (1 - d)) & 1)
 				_uv[d][m] = uvMid.get(index: d)
 				_uv[d][1 - m] = uv[d][1 - m]
 			}
-			let child = try S2Cell(cellId: id, face: face, level: level + 1, orientation: orientation ^ UInt8(S2.posToOrientation(position: pos)), uv: _uv)
+			let child = S2Cell(cellId: id, face: face, level: level + 1, orientation: orientation ^ UInt8(S2.posToOrientation(position: pos)), uv: _uv)
 			children.append(child)
 
 			id = id.next()
@@ -197,6 +197,57 @@ public struct S2Cell: S2Region, Equatable {
 		return uvPoint.x >= uv[0][0] && uvPoint.x <= uv[0][1] && uvPoint.y >= uv[1][0] && uvPoint.y <= uv[1][1]
 	}
 
+	/**
+	* Return the average area for cells at the given level.
+	*/
+	public static func averageArea(level: Int) -> Double {
+		return S2Projections.avgArea.getValue(level: level)
+	}
+	
+	/**
+		Return the average area of cells at this level. This is accurate to within
+		a factor of 1.7 (for S2_QUADRATIC_PROJECTION) and is extremely cheap to compute.
+	*/
+	public var averageArea: Double {
+		return S2Cell.averageArea(level: Int(level))
+	}
+	
+	/**
+		Return the approximate area of this cell. This method is accurate to within
+		3% percent for all cell sizes and accurate to within 0.1% for cells at
+		level 5 or higher (i.e. 300km square or smaller). It is moderately cheap to compute.
+	*/
+	public var approxArea: Double {
+		// All cells at the first two levels have the same area.
+		if level < 2 { return averageArea }
+		
+		// First, compute the approximate area of the cell when projected
+		// perpendicular to its normal. The cross product of its diagonals gives
+		// the normal, and the length of the normal is twice the projected area.
+		let flatArea = 0.5 * (getVertex(2) - getVertex(0)).crossProd((getVertex(3) - getVertex(1))).norm
+		
+		// Now, compensate for the curvature of the cell surface by pretending
+		// that the cell is shaped like a spherical cap. The ratio of the
+		// area of a spherical cap to the area of its projected disc turns out
+		// to be 2 / (1 + sqrt(1 - r*r)) where "r" is the radius of the disc.
+		// For example, when r=0 the ratio is 1, and when r=1 the ratio is 2.
+		// Here we set Pi*r*r == flat_area to find the equivalent disc.
+		return flatArea * 2 / (1 + sqrt(1 - min(M_1_PI * flatArea, 1.0)))
+	}
+	
+	/**
+		Return the area of this cell as accurately as possible. This method is more
+		expensive but it is accurate to 6 digits of precision even for leaf cells
+		(whose area is approximately 1e-18).
+	*/
+	public var exactArea: Double {
+		let v0 = getVertex(0)
+		let v1 = getVertex(1)
+		let v2 = getVertex(2)
+		let v3 = getVertex(3)
+		return S2.area(a: v0, b: v1, c: v2) + S2.area(a: v0, b: v2, c: v3)
+	}
+	
 	////////////////////////////////////////////////////////////////////////
 	// MARK: S2Region
 	////////////////////////////////////////////////////////////////////////
@@ -236,7 +287,7 @@ public struct S2Cell: S2Region, Equatable {
 			let i = S2Projections.getUAxis(face: Int(face)).z == 0 ? (u < 0 ? 1 : 0) : (u > 0 ? 1 : 0)
 			let j = S2Projections.getVAxis(face: Int(face)).z == 0 ? (v < 0 ? 1 : 0) : (v > 0 ? 1 : 0)
 
-			var lat = R1Interval.fromPointPair(p1: getLatitude(i: i, j: j), p2: getLatitude(i: 1 - i, j: 1 - j))
+			var lat = R1Interval(p1: getLatitude(i: i, j: j), p2: getLatitude(i: 1 - i, j: 1 - j))
 			lat = lat.expanded(radius: S2Cell.maxError).intersection(with: S2LatLngRect.fullLat)
 			if (lat.lo == -M_PI_2 || lat.hi == M_PI_2) {
 				return S2LatLngRect(lat: lat, lng: S1Interval.full)
@@ -264,11 +315,11 @@ public struct S2Cell: S2Region, Equatable {
 	}
 
 	public func contains(cell: S2Cell) -> Bool {
-		return false
+		return cellId.contains(other: cell.cellId)
 	}
 
 	public func mayIntersect(cell: S2Cell) -> Bool {
-		return false
+		return cellId.intersects(other: cell.cellId)
 	}
 
 	// Return the latitude or longitude of the cell vertex given by (i,j),