From 6d5e19110046ea2aefc645f40dfce4b3dadb8109 Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Thu, 22 Feb 2018 20:59:57 -0500 Subject: Fix issue #812 Implement Polygon AddEdge fix in C library (will be version 1.49.2 of C library for proj 5.0.0). Still to do: add tests to expand code coverage. This will only affect geodtest.c which is not part of the compiled library. --- src/geodesic.c | 14 +++++++++----- src/geodesic.h | 2 +- src/geodtest.c | 18 ++++++++++++++++++ 3 files changed, 28 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/geodesic.c b/src/geodesic.c index 233dc34c..6a383b72 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -89,10 +89,14 @@ static void Init() { tolb = tol0 * tol2; xthresh = 1000 * tol2; degree = pi/180; + #if defined(NAN) + NaN = NAN; + #else { real minus1 = -1; NaN = sqrt(minus1); } + #endif init = 1; } } @@ -1809,13 +1813,13 @@ int transitdirect(real lon1, real lon2) { #if HAVE_C99_MATH lon1 = remainder(lon1, (real)(720)); lon2 = remainder(lon2, (real)(720)); - return ( (lon2 >= 0 && lon2 < 360 ? 0 : 1) - - (lon1 >= 0 && lon1 < 360 ? 0 : 1) ); + return ( (lon2 <= 0 && lon2 > -360 ? 1 : 0) - + (lon1 <= 0 && lon1 > -360 ? 1 : 0) ); #else lon1 = fmod(lon1, (real)(720)); lon2 = fmod(lon2, (real)(720)); - return ( ((lon2 >= 0 && lon2 < 360) || lon2 < -360 ? 0 : 1) - - ((lon1 >= 0 && lon1 < 360) || lon1 < -360 ? 0 : 1) ); + return ( ((lon2 <= 0 && lon2 > -360) || lon2 > 360 ? 1 : 0) - + ((lon1 <= 0 && lon1 > -360) || lon1 > 360 ? 1 : 0) ); #endif } @@ -1888,7 +1892,7 @@ void geod_polygon_addpoint(const struct geod_geodesic* g, void geod_polygon_addedge(const struct geod_geodesic* g, struct geod_polygon* p, real azi, real s) { - if (p->num) { /* Do nothing is num is zero */ + if (p->num) { /* Do nothing is num is zero */ real lat, lon, S12 = 0; /* Initialize S12 to stop Visual Studio warning */ geod_gendirect(g, p->lat, p->lon, azi, GEOD_LONG_UNROLL, s, &lat, &lon, 0, diff --git a/src/geodesic.h b/src/geodesic.h index 43fd0d1f..536ba7c4 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -132,7 +132,7 @@ * The patch level of the geodesic library. (This tracks the version of * GeographicLib.) **********************************************************************/ -#define GEODESIC_VERSION_PATCH 1 +#define GEODESIC_VERSION_PATCH 2 /** * Pack the version components into a single integer. Users should not rely on diff --git a/src/geodtest.c b/src/geodtest.c index 0f2c0ac2..517e5bed 100644 --- a/src/geodtest.c +++ b/src/geodtest.c @@ -787,6 +787,23 @@ static int Planimeter13() { return result; } +static int AddEdge1() { + /* Check fix to transitdirect vs transit zero handling inconsistency */ + struct geod_geodesic g; + struct geod_polygon p; + double area; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_polygon_init(&p, 0); + geod_polygon_addpoint(&g, &p, 0, 0); + geod_polygon_addedge(&g, &p, 90, 1000); + geod_polygon_addedge(&g, &p, 0, 1000); + geod_polygon_addedge(&g, &p, -90, 1000); + geod_polygon_compute(&g, &p, 0, 1, &area, 0); + result += assertEquals(area, 1000000.0, 0.01); + return result; +} + int main() { int n = 0, i; if ((i = testinverse())) {++n; printf("testinverse fail: %d\n", i);} @@ -823,6 +840,7 @@ int main() { if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);} if ((i = Planimeter12())) {++n; printf("Planimeter12 fail: %d\n", i);} if ((i = Planimeter13())) {++n; printf("Planimeter13 fail: %d\n", i);} + if ((i = AddEdge1())) {++n; printf("AddEdge1 fail: %d\n", i);} return n; } -- cgit v1.2.3 From d83971aa739e421c8d73e08189f274dc0eb29f3d Mon Sep 17 00:00:00 2001 From: Charles Karney Date: Fri, 23 Feb 2018 11:46:47 -0500 Subject: Add tests to improve the code coverage. Now all of geodesic.c is covered except for 3 lines (and 2 of these are intentionally dead code). This corresponds to tag v1.49.2-c in the GeographicLib code base. While testing polygons which encircle the globe multiple times, I uncovered a problem where the range of the area was not reduced to the allowed range (either [0, area0) or (-area0/2, area0/2]) correctly. Since the documentation explicity restricted the calculation of polygon areas to simple polygons, we'll defer fixing this for now. (However the intention was always to handle the area "algebraically" so that, for example, a "bowtie" has zero area. So I will plan on fixing this for 1.50.) Update copyright dates + NEWS. --- src/geodesic.c | 2 +- src/geodesic.h | 2 +- src/geodtest.c | 214 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 215 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/geodesic.c b/src/geodesic.c index 6a383b72..23557b5c 100644 --- a/src/geodesic.c +++ b/src/geodesic.c @@ -18,7 +18,7 @@ * * See the comments in geodesic.h for documentation. * - * Copyright (c) Charles Karney (2012-2017) and licensed + * Copyright (c) Charles Karney (2012-2018) and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ */ diff --git a/src/geodesic.h b/src/geodesic.h index 536ba7c4..cd7c2b70 100644 --- a/src/geodesic.h +++ b/src/geodesic.h @@ -107,7 +107,7 @@ * twice about restructuring the internals of the C code since this may make * porting fixes from the C++ code more difficult. * - * Copyright (c) Charles Karney (2012-2017) and licensed + * Copyright (c) Charles Karney (2012-2018) and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ * diff --git a/src/geodtest.c b/src/geodtest.c index 517e5bed..76c0d14b 100644 --- a/src/geodtest.c +++ b/src/geodtest.c @@ -4,7 +4,7 @@ * * Run these tests by configuring with cmake and running "make test". * - * Copyright (c) Charles Karney (2015-2017) and licensed + * Copyright (c) Charles Karney (2015-2018) and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ @@ -688,6 +688,48 @@ static int GeodSolve78() { return result; } +static int GeodSolve80() { + /* Some tests to add code coverage: computing scale in special cases + zero + length geodesic (includes GeodSolve80 - GeodSolve83) + using an incapable + line. */ + double a12, s12, azi1, azi2, m12, M12, M21, S12; + struct geod_geodesic g; + struct geod_geodesicline l; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_geninverse(&g, 0, 0, 0, 90, 0, 0, 0, 0, &M12, &M21, 0); + result += assertEquals(M12, -0.0052842753, 0.5e-10); + result += assertEquals(M21, -0.0052842753, 0.5e-10); + geod_geninverse(&g, 0, 0, 1e-6, 1e-6, 0, 0, 0, 0, &M12, &M21, 0); + result += assertEquals(M12, 1, 0.5e-10); + result += assertEquals(M21, 1, 0.5e-10); + a12 = geod_geninverse(&g, 20.001, 0, 20.001, 0, + &s12, &azi1, &azi2, &m12, &M12, &M21, &S12); + result += a12 == 0 ? 0 : 1; + result += s12 == 0 ? 0 : 1; + result += azi1 == 180 ? 0 : 1; + result += azi2 == 180 ? 0 : 1; + result += m12 == 0 ? 0 : 1; + result += assertEquals(M12, 1, 1e-15); + result += assertEquals(M21, 1, 1e-15); + result += S12 == 0 ? 0 : 1; + a12 = geod_geninverse(&g, 90, 0, 90, 180, + &s12, &azi1, &azi2, &m12, &M12, &M21, &S12); + result += a12 == 0 ? 0 : 1; + result += s12 == 0 ? 0 : 1; + result += azi1 == 0 ? 0 : 1; + result += azi2 == 180 ? 0 : 1; + result += m12 == 0 ? 0 : 1; + result += assertEquals(M12, 1, 1e-15); + result += assertEquals(M21, 1, 1e-15); + result += assertEquals(S12, 127516405431022, 0.5); + /* An incapable line which can't take distance as input */ + geod_lineinit(&l, &g, 1, 2, 90, GEOD_LATITUDE); + a12 = geod_genposition(&l, 0, 1000, 0, 0, 0, 0, 0, 0, 0, 0); + result += a12 != a12 ? 0 : 1; + return result; +} + static int Planimeter0() { /* Check fix for pole-encircling bug found 2011-03-16 */ double pa[4][2] = {{89, 0}, {89, 90}, {89, 180}, {89, 270}}; @@ -787,6 +829,140 @@ static int Planimeter13() { return result; } +static int Planimeter15() { + /* Coverage tests, includes Planimeter15 - Planimeter18 (combinations of + reverse and sign) + calls to testpoint, testedge, geod_polygonarea. */ + struct geod_geodesic g; + struct geod_polygon p; + double lat[] = {2, 1, 3}, lon[] = {1, 2, 3}; + double area, s12, azi1; + double r = 18454562325.45119, + a0 = 510065621724088.5093; /* ellipsoid area */ + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_polygon_init(&p, 0); + geod_polygon_addpoint(&g, &p, lat[0], lon[0]); + geod_polygon_addpoint(&g, &p, lat[1], lon[1]); + geod_polygon_testpoint(&g, &p, lat[2], lon[2], 0, 1, &area, 0); + result += assertEquals(area, r, 0.5); + geod_polygon_testpoint(&g, &p, lat[2], lon[2], 0, 0, &area, 0); + result += assertEquals(area, r, 0.5); + geod_polygon_testpoint(&g, &p, lat[2], lon[2], 1, 1, &area, 0); + result += assertEquals(area, -r, 0.5); + geod_polygon_testpoint(&g, &p, lat[2], lon[2], 1, 0, &area, 0); + result += assertEquals(area, a0-r, 0.5); + geod_inverse(&g, lat[1], lon[1], lat[2], lon[2], &s12, &azi1, 0); + geod_polygon_testedge(&g, &p, azi1, s12, 0, 1, &area, 0); + result += assertEquals(area, r, 0.5); + geod_polygon_testedge(&g, &p, azi1, s12, 0, 0, &area, 0); + result += assertEquals(area, r, 0.5); + geod_polygon_testedge(&g, &p, azi1, s12, 1, 1, &area, 0); + result += assertEquals(area, -r, 0.5); + geod_polygon_testedge(&g, &p, azi1, s12, 1, 0, &area, 0); + result += assertEquals(area, a0-r, 0.5); + geod_polygon_addpoint(&g, &p, lat[2], lon[2]); + geod_polygon_compute(&g, &p, 0, 1, &area, 0); + result += assertEquals(area, r, 0.5); + geod_polygon_compute(&g, &p, 0, 0, &area, 0); + result += assertEquals(area, r, 0.5); + geod_polygon_compute(&g, &p, 1, 1, &area, 0); + result += assertEquals(area, -r, 0.5); + geod_polygon_compute(&g, &p, 1, 0, &area, 0); + result += assertEquals(area, a0-r, 0.5); + geod_polygonarea(&g, lat, lon, 3, &area, 0); + result += assertEquals(area, r, 0.5); + return result; +} + +static int Planimeter19() { + /* Coverage tests, includes Planimeter19 - Planimeter20 (degenerate + polygons) + extra cases. */ + struct geod_geodesic g; + struct geod_polygon p; + double area, perim; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_polygon_init(&p, 0); + geod_polygon_compute(&g, &p, 0, 1, &area, &perim); + result += area == 0 ? 0 : 1; + result += perim == 0 ? 0 : 1; + geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, &area, &perim); + result += area == 0 ? 0 : 1; + result += perim == 0 ? 0 : 1; + geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, &area, &perim); + result += area != area ? 0 : 1; + result += perim != perim ? 0 : 1; + geod_polygon_addpoint(&g, &p, 1, 1); + geod_polygon_compute(&g, &p, 0, 1, &area, &perim); + result += area == 0 ? 0 : 1; + result += perim == 0 ? 0 : 1; + geod_polygon_init(&p, 1); + geod_polygon_compute(&g, &p, 0, 1, 0, &perim); + result += perim == 0 ? 0 : 1; + geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, 0, &perim); + result += perim == 0 ? 0 : 1; + geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, 0, &perim); + result += perim != perim ? 0 : 1; + geod_polygon_addpoint(&g, &p, 1, 1); + geod_polygon_compute(&g, &p, 0, 1, 0, &perim); + result += perim == 0 ? 0 : 1; + return result; +} + +static int Planimeter21() { + /* Some test to add code coverage: multiple circlings of pole (includes + Planimeter21 - Planimeter28) + invocations via testpoint and testedge. + Some of the results for i = 4 in the loop are wrong because we don't + reduce the area to the allowed range correctly. However these cases are + not "simple" polygons, so we'll defer fixing the problem for now. + */ + struct geod_geodesic g; + struct geod_polygon p; + double area, lat = 45, + a = 39.2144607176828184218, s = 8420705.40957178156285, + r = 39433884866571.4277, /* Area for one circuit */ + a0 = 510065621724088.5093; /* Ellipsoid area */ + int result = 0, i; + geod_init(&g, wgs84_a, wgs84_f); + geod_polygon_init(&p, 0); + geod_polygon_addpoint(&g, &p, lat, 60); + geod_polygon_addpoint(&g, &p, lat, 180); + geod_polygon_addpoint(&g, &p, lat, -60); + geod_polygon_addpoint(&g, &p, lat, 60); + geod_polygon_addpoint(&g, &p, lat, 180); + geod_polygon_addpoint(&g, &p, lat, -60); + for (i = 3; i <= 4; ++i) { + geod_polygon_addpoint(&g, &p, lat, 60); + geod_polygon_addpoint(&g, &p, lat, 180); + geod_polygon_testpoint(&g, &p, lat, -60, 0, 1, &area, 0); + if (i != 4) result += assertEquals(area, i*r, 0.5); + geod_polygon_testpoint(&g, &p, lat, -60, 0, 0, &area, 0); + if (i != 4) result += assertEquals(area, i*r, 0.5); + geod_polygon_testpoint(&g, &p, lat, -60, 1, 1, &area, 0); + if (i != 4) result += assertEquals(area, -i*r, 0.5); + geod_polygon_testpoint(&g, &p, lat, -60, 1, 0, &area, 0); + result += assertEquals(area, -i*r + a0, 0.5); + geod_polygon_testedge(&g, &p, a, s, 0, 1, &area, 0); + if (i != 4) result += assertEquals(area, i*r, 0.5); + geod_polygon_testedge(&g, &p, a, s, 0, 0, &area, 0); + if (i != 4) result += assertEquals(area, i*r, 0.5); + geod_polygon_testedge(&g, &p, a, s, 1, 1, &area, 0); + if (i != 4) result += assertEquals(area, -i*r, 0.5); + geod_polygon_testedge(&g, &p, a, s, 1, 0, &area, 0); + result += assertEquals(area, -i*r + a0, 0.5); + geod_polygon_addpoint(&g, &p, lat, -60); + geod_polygon_compute(&g, &p, 0, 1, &area, 0); + if (i != 4) result += assertEquals(area, i*r, 0.5); + geod_polygon_compute(&g, &p, 0, 0, &area, 0); + if (i != 4) result += assertEquals(area, i*r, 0.5); + geod_polygon_compute(&g, &p, 1, 1, &area, 0); + if (i != 4) result += assertEquals(area, -i*r, 0.5); + geod_polygon_compute(&g, &p, 1, 0, &area, 0); + result += assertEquals(area, -i*r + a0, 0.5); + } + return result; +} + static int AddEdge1() { /* Check fix to transitdirect vs transit zero handling inconsistency */ struct geod_geodesic g; @@ -804,6 +980,37 @@ static int AddEdge1() { return result; } +static int EmptyPoly() { + struct geod_geodesic g; + struct geod_polygon p; + double perim, area; + int result = 0; + geod_init(&g, wgs84_a, wgs84_f); + geod_polygon_init(&p, 0); + geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, &area, &perim); + result += area == 0 ? 0 : 1; + result += perim == 0 ? 0 : 1; + geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, &area, &perim); + result += area != area ? 0 : 1; + result += perim != perim ? 0 : 1; + geod_polygon_compute(&g, &p, 0, 1, &area, &perim); + result += area == 0 ? 0 : 1; + result += perim == 0 ? 0 : 1; + geod_polygon_init(&p, 1); + geod_polygon_testpoint(&g, &p, 1, 1, 0, 1, 0, &perim); + result += perim == 0 ? 0 : 1; + geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, 0, &perim); + result += perim != perim ? 0 : 1; + geod_polygon_compute(&g, &p, 0, 1, 0, &perim); + result += perim == 0 ? 0 : 1; + geod_polygon_addpoint(&g, &p, 1, 1); + geod_polygon_testedge(&g, &p, 90, 1000, 0, 1, 0, &perim); + result += assertEquals(perim, 1000, 1e-10); + geod_polygon_testpoint(&g, &p, 2, 2, 0, 1, 0, &perim); + result += assertEquals(perim, 156876.149, 0.5e-3); + return result; +} + int main() { int n = 0, i; if ((i = testinverse())) {++n; printf("testinverse fail: %d\n", i);} @@ -835,12 +1042,17 @@ int main() { if ((i = GeodSolve74())) {++n; printf("GeodSolve74 fail: %d\n", i);} if ((i = GeodSolve76())) {++n; printf("GeodSolve76 fail: %d\n", i);} if ((i = GeodSolve78())) {++n; printf("GeodSolve78 fail: %d\n", i);} + if ((i = GeodSolve80())) {++n; printf("GeodSolve80 fail: %d\n", i);} if ((i = Planimeter0())) {++n; printf("Planimeter0 fail: %d\n", i);} if ((i = Planimeter5())) {++n; printf("Planimeter5 fail: %d\n", i);} if ((i = Planimeter6())) {++n; printf("Planimeter6 fail: %d\n", i);} if ((i = Planimeter12())) {++n; printf("Planimeter12 fail: %d\n", i);} if ((i = Planimeter13())) {++n; printf("Planimeter13 fail: %d\n", i);} + if ((i = Planimeter15())) {++n; printf("Planimeter15 fail: %d\n", i);} + if ((i = Planimeter19())) {++n; printf("Planimeter19 fail: %d\n", i);} + if ((i = Planimeter21())) {++n; printf("Planimeter21 fail: %d\n", i);} if ((i = AddEdge1())) {++n; printf("AddEdge1 fail: %d\n", i);} + if ((i = EmptyPoly())) {++n; printf("EmptyPoly fail: %d\n", i);} return n; } -- cgit v1.2.3