Added new point-on-edge routine still no joy on bad test case.

git-svn-id: http://svn.osgeo.org/postgis/trunk@4589 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index f6ac7937af6387aa491bfaa6c19b7d477afce143..03ba57c3eeb1bde4312f0ab4df750307b3e7865f 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -282,7 +282,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;
@@ -294,7 +294,7 @@ void test_edge_intersection(void)
        rv = edge_intersection(e1, e2, &g);
        CU_ASSERT_EQUAL(rv, LW_FALSE);
        */
-
+       
        /* Intersection at (0 0) */
        edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
        edge_set(0.0, -1.0, 0.0, 1.0, &e2);
@@ -305,11 +305,13 @@ void test_edge_intersection(void)
        CU_ASSERT_EQUAL(rv, LW_TRUE);
 
        /*  No intersection at (0 0) */
+       edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
        edge_set(0.0, -1.0, 0.0, -2.0, &e2);
        rv = edge_intersection(e1, e2, &g);
        CU_ASSERT_EQUAL(rv, LW_FALSE);
 
        /*  End touches middle of segment at (0 0) */
+       edge_set(-1.0, 0.0, 1.0, 0.0, &e1);
        edge_set(0.0, -1.0, 0.0, 0.0, &e2);
        rv = edge_intersection(e1, e2, &g);
        point_rad2deg(&g);
@@ -318,6 +320,7 @@ void test_edge_intersection(void)
 
        /*  End touches end of segment at (0 0) */
        edge_set(0.0, 0.0, 1.0, 0.0, &e1);
+       edge_set(0.0, -1.0, 0.0, 0.0, &e2);
        rv = edge_intersection(e1, e2, &g);
        point_rad2deg(&g);
 #if 0

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 89869ede77de938f0018cd19c61e3fb12bf6ba26..56d16c49db5e9559bee2485aeffdf8c27d13af75 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -352,6 +352,7 @@ void y_to_z(POINT3D *p)
        p->y = tmp;
 }
 
+
 /**
 * Returns true if the point p is on the great circle plane.
 * Forms the scalar triple product of A,B,p and if the volume of the
@@ -371,7 +372,7 @@ int edge_point_on_plane(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
        LWDEBUGF(4,"dot product %.9g",w);
        if ( FP_IS_ZERO(w) )
        {
-               LWDEBUG(4, "point is on plane");
+               LWDEBUG(4, "point is on plane (dot product is zero)");
                return LW_TRUE;
        }
        LWDEBUG(4, "point is not on plane");
@@ -414,22 +415,112 @@ int edge_point_in_cone(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
 /**
 * True if the longitude of p is within the range of the longitude of the ends of e
 */
-int edge_contains_longitude(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
+int edge_contains_coplanar_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
 {
+       GEOGRAPHIC_EDGE g;
+       GEOGRAPHIC_POINT q;
+       double slon = fabs(e.start.lon) + fabs(e.end.lon);
+       double dlon = fabs(fabs(e.start.lon) - fabs(e.end.lon));
+       double slat = e.start.lat + e.end.lat;
+
        LWDEBUGF(4, "e.start == GPOINT(%.6g %.6g) ", e.start.lat, e.start.lon);
        LWDEBUGF(4, "e.end == GPOINT(%.6g %.6g) ", e.end.lat, e.end.lon);
        LWDEBUGF(4, "p == GPOINT(%.6g %.6g) ", p.lat, p.lon);
-       if ( e.start.lon < p.lon && p.lon < e.end.lon )
+       
+       /* Copy values into working registers */
+       g = e;
+       q = p;
+
+       /* Vertical plane, we need to do this calculation in latitude */
+       if( FP_EQUALS( g.start.lon, g.end.lon ) )
        {
-               LWDEBUG(4, "returning TRUE");
-               return LW_TRUE;
+               LWDEBUG(4, "vertical plane, we need to do this calculation in latitude");
+               /* Supposed to be co-planar... */
+               if ( ! FP_EQUALS( q.lon, g.start.lon ) )
+                       return LW_FALSE;
+
+               if ( ( g.start.lat <= q.lat && q.lat <= g.end.lat ) || 
+                    ( g.end.lat <= q.lat && q.lat <= g.start.lat ) )
+               {
+                       return LW_TRUE;
+               }
+               else
+               {
+                       return LW_FALSE;
+               }
        }
-       if ( e.end.lon < p.lon && p.lon < e.start.lon )
+
+       /* Over the pole, we need normalize latitude and do this calculation in latitude */
+       if ( FP_EQUALS( slon, M_PI ) && ( signum(g.start.lon) != signum(g.end.lon) || FP_EQUALS(dlon, M_PI) ) )
        {
-               LWDEBUG(4, "returning TRUE");
+               LWDEBUG(4, "over the pole...");
+               /* Antipodal, everything (or nothing?) is inside */
+               if ( FP_EQUALS( slat, 0.0 ) )
+                       return LW_TRUE;
+
+               /* Point *is* the north pole */
+               if ( slat > 0.0 && FP_EQUALS(q.lat, M_PI_2 ) )
+                       return LW_TRUE;
+               
+               /* Point *is* the south pole */
+               if ( slat < 0.0 && FP_EQUALS(q.lat, -1.0 * M_PI_2) )
+                       return LW_TRUE;
+
+               LWDEBUG(4, "coplanar?...");
+
+               /* Supposed to be co-planar... */
+               if ( ! FP_EQUALS( q.lon, g.start.lon ) )
+                       return LW_FALSE;
+
+               LWDEBUG(4, "north or south?...");
+                               
+               /* Over north pole, test based on south pole */
+               if ( slat > 0.0 )
+               {
+                       LWDEBUG(4, "over the north pole...");
+                       if( q.lat > FP_MIN(g.start.lat, g.end.lat) )
+                               return LW_TRUE;
+                       else 
+                               return LW_FALSE;
+               }
+               else
+               /* Over south pole, test based on north pole */
+               {
+                       LWDEBUG(4, "over the south pole...");
+                       if( q.lat < FP_MAX(g.start.lat, g.end.lat) )
+                               return LW_TRUE;
+                       else 
+                               return LW_FALSE;
+               }                       
+       }
+
+       /* Dateline crossing, flip everything to the opposite hemisphere */
+       else if( slon > M_PI && ( signum(g.start.lon) != signum(g.end.lon) ) )
+       {
+               LWDEBUG(4, "crosses dateline, flip longitudes...");
+               if ( g.start.lon > 0.0 )
+                       g.start.lon -= M_PI;
+               else
+                       g.start.lon += M_PI;
+               if ( g.end.lon > 0.0 )
+                       g.end.lon -= M_PI;
+               else
+                       g.end.lon += M_PI;
+
+               if ( q.lon > 0.0 )
+                       q.lon -= M_PI;
+               else
+                       q.lon += M_PI;                  
+       }
+
+       if ( ( g.start.lon <= q.lon && q.lon <= g.end.lon ) || 
+            ( g.end.lon <= q.lon && q.lon <= g.start.lon ) )
+       {
+               LWDEBUG(4, "true, this edge contains point");
                return LW_TRUE;
        }
-       LWDEBUG(4, "returning FALSE");
+
+       LWDEBUG(4, "false, this edge does not contain point");
        return LW_FALSE;
 }
 
@@ -483,6 +574,7 @@ double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
 int edge_contains_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
 {
        if ( edge_point_in_cone(e, p) && edge_point_on_plane(e, p) )
+/*     if ( edge_contains_coplanar_point(e, p) && edge_point_on_plane(e, p) ) */
        {
                LWDEBUG(4, "point is on edge");
                return LW_TRUE;
@@ -773,9 +865,7 @@ int edge_calculate_gbox_slow(GEOGRAPHIC_EDGE e, GBOX *gbox)
        }
 
        /* Walk along the chord between start and end incrementally, 
-          normalizing at each step. Only corner case occurs if path
-          passes through 0,0,0 (as the normalization will return zero)
-          but the box should always have extrema beyond that. */
+          normalizing at each step. */
        geog2cart(e.start, &start);
        geog2cart(e.end, &end);
        dx = (end.x - start.x)/steps;

diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index 46f5aef0340dc7782cbd54e43ea2b577e4343b77..31d7e4704571197a3ccdbf86d14c878e6f6ddc30 100644 (file)
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -53,7 +53,7 @@ void x_to_z(POINT3D *p);
 void y_to_z(POINT3D *p);
 int edge_point_on_plane(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
 int edge_point_in_cone(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
-int edge_contains_longitude(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
+int edge_contains_coplanar_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
 int edge_contains_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
 double z_to_latitude(double z, int top);
 int clairaut_cartesian(POINT3D start, POINT3D end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);