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 bbe2be6b7f9bfc7e0d19623fafe05cb5df770dad..b61e4a6baf7434bf11eed300d1dbb0e0c3ffe31f 100644 (file)
--- 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 75584edc5ca4a6842a379a5facd2a85a352fd99f..eee0d8c1dc2ec61977bc1bbf58274838a38f5ffc 100644 (file)
--- 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 )