Area: new implementation avoiding isCCW, and taking radius.

library/src/com/google/maps/android/SphericalUtil.java

@@ -189,100 +189,62 @@ public static double computeLength(List<LatLng> path) {
     }
 
     /**
-     * Returns the area of a closed path. The computed area uses the same units as
-     * the radius.
+     * Returns the area of a closed path on Earth.
      * @param path A closed path.
-     * @return The loop's area in square meters.
+     * @return The path's area in square meters.
      */
     public static double computeArea(List<LatLng> path) {
         return abs(computeSignedArea(path));
     }
 
     /**
-     * Returns the signed area of a closed path. The signed area may be used to
-     * determine the orientation of the path. The computed area uses the same
-     * units as the radius.
+     * Returns the signed area of a closed path on Earth. The sign of the area may be used to
+     * determine the orientation of the path.
+     * "inside" is the surface that does not contain the South Pole.
      * @param loop A closed path.
      * @return The loop's area in square meters.
      */
-    public static double computeSignedArea(List<LatLng> loop) {
-        // For each edge, accumulate the signed area of the triangle formed with
-        // the first point. We can skip the first and last edge as they form
-        // triangles of zero area.
-        LatLng origin = loop.get(0);
-        double total = 0;
-        for (int i = 1, I = loop.size() - 1; i < I; ++i) {
-            total += computeSignedTriangleArea(origin, loop.get(i), loop.get(i + 1));
-        }
-        return total * EARTH_RADIUS * EARTH_RADIUS;
-    }
-
-    /**
-     * Compute the signed area of the triangle [a, b, c] on the unit sphere.
-     */
-    static double computeSignedTriangleArea(LatLng a, LatLng b, LatLng c) {
-        return computeTriangleArea(a, b, c) * isCCW(a, b, c);
+    public static double computeSignedArea(List<LatLng> path) {
+        return computeSignedArea(path, EARTH_RADIUS);
     }
 
     /**
-     * Compute the area of the triangle [a, b, c] on the unit sphere.
-     * We use l'Huilier's theorem, which is the Spherical analogue of Heron's
-     * theorem for the area of a triangle in R2.
-     * @return The area.
+     * Returns the signed area of a closed path on a sphere of given radius.
+     * The computed area uses the same units as the radius squared.
+     * Used by SphericalUtilTest.
      */
-    static double computeTriangleArea(LatLng a, LatLng b, LatLng c) {
-        LatLng[] points = new LatLng[]{a, b, c, a};  // Simplify cyclic indexing
-
-        // Compute the length of each edge, and s which is half the perimeter.
-        double[] angles = new double[3];
-        double s = 0;
-        for (int i = 0; i < 3; ++i) {
-            angles[i] = computeAngleBetween(points[i], points[i + 1]);
-            s += angles[i];
-        }
-        s /= 2;
-
-        // Apply l'Huilier's theorem
-        double product = tan(s / 2);
-        for (int i = 0; i < 3; ++i) {
-            product *= tan((s - angles[i]) / 2);
+    static double computeSignedArea(List<LatLng> path, double radius) {
+        int size = path.size();
+        if (size < 3) { return 0; }
+        double total = 0;
+        LatLng prev = path.get(size - 1);
+        double prevTanLat = tan((PI / 2 - toRadians(prev.latitude)) / 2);
+        double prevLng = toRadians(prev.longitude);
+        // For each edge, accumulate the signed area of the triangle formed by the North Pole
+        // and that edge ("polar triangle").
+        for (LatLng point : path) {
+            double tanLat = tan((PI / 2 - toRadians(point.latitude)) / 2);
+            double lng = toRadians(point.longitude);
+            total += polarTriangleArea(tanLat, lng, prevTanLat, prevLng);
+            prevTanLat = tanLat;
+            prevLng = lng;
         }
-        return 4 * atan(sqrt(abs(product)));
+        return total * (radius * radius);
     }
 
     /**
-     * Compute the ordering of 3 points in a triangle:
-     * Counter ClockWise (CCW) vs ClockWise (CW).
-     * Results are indeterminate for triangles of area 0.
-     * @return +1 if CCW, -1 if CW.
+     * Returns the signed area of a triangle which has North Pole as a vertex.
+     * Formula derived from "Area of a spherical triangle given two edges and the included angle"
+     * as per "Spherical Trigonometry" by Todhunter, page 71, section 103, point 2.
+     * See http://books.google.com/books?id=3uBHAAAAIAAJ&pg=PA71
+     * The arguments named "tan" are tan((pi/2 - latitude)/2).
      */
-    static int isCCW(LatLng a, LatLng b, LatLng c) {
-        // Convert the 3 points to 3 unit vectors on the sphere.
-        LatLng[] points = new LatLng[]{a, b, c};
-        double[][] pointsR3 = new double[3][];
-        for (int i = 0; i < 3; ++i) {
-            LatLng latLng = points[i];
-            double lat = toRadians(latLng.latitude);
-            double lng = toRadians(latLng.longitude);
-            double[] r3 = new double[3];
-            r3[0] = cos(lat) * cos(lng);
-            r3[1] = cos(lat) * sin(lng);
-            r3[2] = sin(lat);
-            pointsR3[i] = r3;
-        }
-
-        // Compute the determinant of the matrix formed by the 3 unit vectors.
-        double det = pointsR3[0][0] * pointsR3[1][1] * pointsR3[2][2] +
-                pointsR3[1][0] * pointsR3[2][1] * pointsR3[0][2] +
-                pointsR3[2][0] * pointsR3[0][1] * pointsR3[1][2] -
-                pointsR3[0][0] * pointsR3[2][1] * pointsR3[1][2] -
-                pointsR3[1][0] * pointsR3[0][1] * pointsR3[2][2] -
-                pointsR3[2][0] * pointsR3[1][1] * pointsR3[0][2];
-
-        // Threshold to sign
-        return det > 0 ? 1 : -1;
+    private static double polarTriangleArea(double tan1, double lng1, double tan2, double lng2) {
+        double deltaLng = lng1 - lng2;
+        double t = tan1 * tan2;
+        return 2 * atan2(t * sin(deltaLng), 1 + t * cos(deltaLng));
     }
-
+    
     /**
      * Wraps the given value into the inclusive-exclusive interval between min and max.
      * @param n   The value to wrap.

library/tests/src/com/google/maps/android/SphericalUtilTest.java

@@ -33,7 +33,7 @@ private static void expectLatLngApproxEquals(LatLng actual, LatLng expected) {
     }
 
     private static void expectNearNumber(double actual, double expected, double epsilon) {
-        Assert.assertTrue(String.format("Expected %f to be near %f", actual, expected),
+        Assert.assertTrue(String.format("Expected %g to be near %g", actual, expected),
                 Math.abs(expected - actual) <= epsilon);
     }
 
@@ -269,23 +269,35 @@ public void testComputeLength() {
         expectNearNumber(SphericalUtil.computeLength(latLngs), Math.PI * EARTH_RADIUS, 1e-6);
     }
 
+    private static double computeSignedTriangleArea(LatLng a, LatLng b, LatLng c) {
+        List<LatLng> path = Arrays.asList(a, b, c);
+        return SphericalUtil.computeSignedArea(path, 1);
+    }
+
+    private static double computeTriangleArea(LatLng a, LatLng b, LatLng c) {
+        return Math.abs(computeSignedTriangleArea(a, b, c));
+    }
+
+    private static int isCCW(LatLng a, LatLng b, LatLng c) {
+        return computeSignedTriangleArea(a, b, c) > 0 ? 1 : -1;
+    }
+
     public void testIsCCW() {
         // One face of the octahedron
-        expectEq(1, SphericalUtil.isCCW(right, up, front));
-        expectEq(1, SphericalUtil.isCCW(up, front, right));
-        expectEq(1, SphericalUtil.isCCW(front, right, up));
-        expectEq(-1, SphericalUtil.isCCW(front, up, right));
-        expectEq(-1, SphericalUtil.isCCW(up, right, front));
-        expectEq(-1, SphericalUtil.isCCW(right, front, up));
+        expectEq(1, isCCW(right, up, front));
+        expectEq(1, isCCW(up, front, right));
+        expectEq(1, isCCW(front, right, up));
+        expectEq(-1, isCCW(front, up, right));
+        expectEq(-1, isCCW(up, right, front));
+        expectEq(-1, isCCW(right, front, up));
     }
 
     public void testComputeTriangleArea() {
-        expectNearNumber(SphericalUtil.computeTriangleArea(right, up, front), Math.PI / 2, 1e-6);
-
-        expectNearNumber(SphericalUtil.computeTriangleArea(front, up, right), Math.PI / 2, 1e-6);
+        expectNearNumber(computeTriangleArea(right, up, front), Math.PI / 2, 1e-6);
+        expectNearNumber(computeTriangleArea(front, up, right), Math.PI / 2, 1e-6);
 
         // computeArea returns area of zero on small polys
-        double area = SphericalUtil.computeTriangleArea(
+        double area = computeTriangleArea(
                 new LatLng(0, 0),
                 new LatLng(0, Math.toDegrees(1E-6)),
                 new LatLng(Math.toDegrees(1E-6), 0));
@@ -296,14 +308,14 @@ public void testComputeTriangleArea() {
 
     public void testComputeSignedTriangleArea() {
         expectNearNumber(
-                SphericalUtil.computeSignedTriangleArea(
+                computeSignedTriangleArea(
                         new LatLng(0, 0), new LatLng(0, 0.1), new LatLng(0.1, 0.1)),
                 Math.toRadians(0.1) * Math.toRadians(0.1) / 2, 1e-6);
 
-        expectNearNumber(SphericalUtil.computeSignedTriangleArea(right, up, front),
+        expectNearNumber(computeSignedTriangleArea(right, up, front),
                 Math.PI / 2, 1e-6);
 
-        expectNearNumber(SphericalUtil.computeSignedTriangleArea(front, up, right),
+        expectNearNumber(computeSignedTriangleArea(front, up, right),
                 -Math.PI / 2, 1e-6);
     }