diff --git a/library/src/com/google/maps/android/SphericalUtil.java b/library/src/com/google/maps/android/SphericalUtil.java index 429749f82..6bdf6ff2e 100644 --- a/library/src/com/google/maps/android/SphericalUtil.java +++ b/library/src/com/google/maps/android/SphericalUtil.java @@ -189,100 +189,62 @@ public static double computeLength(List 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 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 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 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 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. diff --git a/library/tests/src/com/google/maps/android/SphericalUtilTest.java b/library/tests/src/com/google/maps/android/SphericalUtilTest.java index 515d20962..345ab594b 100644 --- a/library/tests/src/com/google/maps/android/SphericalUtilTest.java +++ b/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 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); }