From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Wed, 4 Nov 2009 00:13:04 +0000 (+0000)
Subject: Remove unit test failure cases in 32-bit architectures. Now have to test correctness
X-Git-Tag: 1.5.0b1~301
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=e0381459b7ebd05794657f371fe45f9cfe507de7;p=postgis
Remove unit test failure cases in 32-bit architectures. Now have to test correctness
of algorithms on test data in 64-bit environment.
git-svn-id: http://svn.osgeo.org/postgis/trunk@4732 b70326c6-7e19-0410-871a-916f4a2858ee
---
diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index bbe2be6b7..b61e4a6ba 100644
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -283,7 +283,7 @@ void test_edge_intersection(void)
rv = edge_intersection(e1, e2, &g);
CU_ASSERT_EQUAL(rv, LW_FALSE);
- /* Second Medford case, very short segment vs very long one
+ /* Second Medford case, very short segment vs very long one
e1.start.lat = 0.73826546728290887156;
e1.start.lon = -2.14426380171833042;
e1.end.lat = 0.73826545883786642843;
@@ -368,34 +368,48 @@ void test_edge_intersection(void)
CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
CU_ASSERT_EQUAL(rv, LW_TRUE);
- /* End touches arc at north pole */
- edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
- edge_set(90.0, 80.0, -90.0, 90.0, &e2);
+ /* Equal edges return true */
+ edge_set(45.0, 10.0, 50.0, 20.0, &e1);
+ edge_set(45.0, 10.0, 50.0, 20.0, &e2);
rv = edge_intersection(e1, e2, &g);
point_rad2deg(&g);
- CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
CU_ASSERT_EQUAL(rv, LW_TRUE);
- /* End touches end at north pole */
- edge_set(-180.0, 80.0, 0.0, 90.0, &e1);
- edge_set(90.0, 80.0, -90.0, 90.0, &e2);
+ /* Parallel edges (same great circle, different end points) return true */
+ edge_set(40.0, 0.0, 70.0, 0.0, &e1);
+ edge_set(60.0, 0.0, 50.0, 0.0, &e2);
rv = edge_intersection(e1, e2, &g);
point_rad2deg(&g);
- CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
CU_ASSERT_EQUAL(rv, LW_TRUE);
- /* Equal edges return true */
- edge_set(45.0, 10.0, 50.0, 20.0, &e1);
- edge_set(45.0, 10.0, 50.0, 20.0, &e2);
+ /* End touches arc at north pole */
+ edge_set(-180.0, 80.0, 0.0, 80.0, &e1);
+ edge_set(90.0, 80.0, -90.0, 90.0, &e2);
rv = edge_intersection(e1, e2, &g);
point_rad2deg(&g);
+#if 0
+ printf("\n");
+ printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
+ printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
+ printf("g = (%.15g %.15g)\n", g.lon, g.lat);
+ printf("rv = %d\n", rv);
+#endif
+ CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
CU_ASSERT_EQUAL(rv, LW_TRUE);
- /* Parallel edges (same great circle, different end points) return true */
- edge_set(40.0, 0.0, 70.0, 0.0, &e1);
- edge_set(60.0, 0.0, 50.0, 0.0, &e2);
+ /* End touches end at north pole */
+ edge_set(-180.0, 80.0, 0.0, 90.0, &e1);
+ edge_set(90.0, 80.0, -90.0, 90.0, &e2);
rv = edge_intersection(e1, e2, &g);
point_rad2deg(&g);
+#if 0
+ printf("\n");
+ printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e1.start.lon, e1.start.lat, e1.end.lon, e1.end.lat);
+ printf("LINESTRING(%.15g %.15g, %.15g %.15g)\n", e2.start.lon, e2.start.lat, e2.end.lon, e2.end.lat);
+ printf("g = (%.15g %.15g)\n", g.lon, g.lat);
+ printf("rv = %d\n", rv);
+#endif
+ CU_ASSERT_DOUBLE_EQUAL(g.lat, 90.0, 0.00001);
CU_ASSERT_EQUAL(rv, LW_TRUE);
}
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 75584edc5..eee0d8c1d 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -438,7 +438,19 @@ int edge_point_in_cone(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
vp_dot_vcp = dot_product(vp, vcp);
LWDEBUGF(4,"vp_dot_vcp %.19g",vp_dot_vcp);
/* If p is more similar than start then p is inside the cone */
- if ( vp_dot_vcp >= vs_dot_vcp )
+ LWDEBUGF(4,"fabs(vp_dot_vcp - vs_dot_vcp) %.39g",fabs(vp_dot_vcp - vs_dot_vcp));
+
+ /*
+ ** We want to test that vp_dot_vcp is >= vs_dot_vcp but there are
+ ** numerical stability issues for values that are very very nearly
+ ** equal. Unfortunately there are also values of vp_dot_vcp that are legitimately
+ ** very close to but still less than vs_dot_vcp which we also need to catch.
+ ** The tolerance of 10-17 seems to do the trick on 32-bit and 64-bit architectures,
+ ** for the test cases here.
+ ** However, tuning the tolerance value feels like a dangerous hack.
+ ** Fundamentally, the problem is that this test is so sensitive.
+ */
+ if ( vp_dot_vcp > vs_dot_vcp || fabs(vp_dot_vcp - vs_dot_vcp) < 1e-17 )
{
LWDEBUG(4, "point is in cone");
return LW_TRUE;
@@ -774,10 +786,10 @@ int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *
}
}
unit_normal(ea, eb, &v);
- LWDEBUGF(4, "v == POINT(%.8g %.8g %.8g)", v.x, v.y, v.z);
+ LWDEBUGF(4, "v == POINT(%.12g %.12g %.12g)", v.x, v.y, v.z);
g->lat = atan2(v.z, sqrt(v.x * v.x + v.y * v.y));
g->lon = atan2(v.y, v.x);
- LWDEBUGF(4, "g == GPOINT(%.6g %.6g)", g->lat, g->lon);
+ LWDEBUGF(4, "g == GPOINT(%.12g %.12g)", g->lat, g->lon);
LWDEBUGF(4, "g == POINT(%.12g %.12g)", rad2deg(g->lon), rad2deg(g->lat));
if ( edge_contains_point(e1, *g) && edge_contains_point(e2, *g) )
{
@@ -785,6 +797,7 @@ int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *
}
else
{
+ LWDEBUG(4, "flipping point to other side of sphere");
g->lat = -1.0 * g->lat;
g->lon = g->lon + M_PI;
if ( g->lon > M_PI )