From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Tue, 6 Oct 2009 04:19:26 +0000 (+0000)
Subject: New point-outside routine, and allow distances against empty geometries.
X-Git-Tag: 1.5.0b1~412
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=fab491e7f54123f7e222f606c1b29895da763722;p=postgis

New point-outside routine, and allow distances against empty geometries.


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

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 5fbf4bd1a..825f4bfc7 100644
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -555,8 +555,4 @@ void test_lwgeom_distance_sphere(void)
 	lwgeom_free(lwg1);
 	lwgeom_free(lwg2);
 
-}
-
-void test_ptarray_point_in_ring_winding(void)
-{
-}
+}
\ No newline at end of file
diff --git a/liblwgeom/g_box.c b/liblwgeom/g_box.c
index 3d671ac4f..37bc380e1 100644
--- a/liblwgeom/g_box.c
+++ b/liblwgeom/g_box.c
@@ -31,6 +31,15 @@ int gbox_merge_point3d(POINT3D p, GBOX *gbox)
 	return G_SUCCESS;
 }
 
+int gbox_contains_point3d(GBOX gbox, POINT3D pt)
+{
+	if( gbox.xmin > pt.x || gbox.ymin > pt.y || gbox.zmin > pt.z ||
+	    gbox.xmax < pt.x || gbox.ymax < pt.y || gbox.zmax < pt.z )
+	{
+		return LW_FALSE;
+	}
+	return LW_TRUE;
+}
 
 int gbox_merge(GBOX new_box, GBOX *merge_box)
 {
diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h
index a80783581..ff1602f81 100644
--- a/liblwgeom/libgeom.h
+++ b/liblwgeom/libgeom.h
@@ -405,6 +405,11 @@ extern int gbox_merge(GBOX new_box, GBOX *merged_box);
 */
 extern int gbox_merge_point3d(POINT3D p, GBOX *gbox);
 
+/**
+* Return true if the point is inside the gbox
+*/
+extern int gbox_contains_point3d(GBOX gbox, POINT3D pt);
+
 /**
 * Allocate a string representation of the #GBOX, based on dimensionality of flags.
 */
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 02a30c857..599f9c4ee 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1182,7 +1182,64 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
 */
 static void gbox_pt_outside(GBOX gbox, POINT3D *pt)
 {
+	static double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
 	double d;
+	int i;
+	GBOX ge;
+	POINT3D corners[8];
+	
+	/* Assign our box and expand it slightly. */
+	ge = gbox;
+	ge.xmin -= grow;
+	ge.ymin -= grow;
+	ge.zmin -= grow;
+	ge.xmax += grow;
+	ge.ymax += grow;
+	ge.zmax += grow;
+	
+	/* Build our eight corner points */
+	corners[0].x = ge.xmin;
+	corners[0].y = ge.ymin;
+	corners[0].z = ge.zmin;
+
+	corners[1].x = ge.xmin;
+	corners[1].y = ge.ymax;
+	corners[1].z = ge.zmin;
+
+	corners[2].x = ge.xmin;
+	corners[2].y = ge.ymin;
+	corners[2].z = ge.zmax;
+
+	corners[3].x = ge.xmax;
+	corners[3].y = ge.ymin;
+	corners[3].z = ge.zmin;
+
+	corners[4].x = ge.xmax;
+	corners[4].y = ge.ymax;
+	corners[4].z = ge.zmin;
+
+	corners[5].x = ge.xmax;
+	corners[5].y = ge.ymin;
+	corners[5].z = ge.zmax;
+
+	corners[6].x = ge.xmin;
+	corners[6].y = ge.ymax;
+	corners[6].z = ge.zmax;
+
+	corners[7].x = ge.xmax;
+	corners[7].y = ge.ymax;
+	corners[7].z = ge.zmax;
+
+	for( i = 0; i < 8; i++ )
+	{
+		normalize(&(corners[i]));
+		if( ! gbox_contains_point3d(gbox, corners[i]) )
+		{
+			*pt = corners[i];
+			return;
+		}	
+	}
+
 	pt->x = 1.0;
 	pt->y = 0.0;
 	pt->z = 0.0;
@@ -1272,57 +1329,6 @@ int ptarray_point_in_ring(POINTARRAY *pa, POINT2D pt_outside, POINT2D pt_to_test
 	return LW_FALSE;
 }
 
-/**
-* This routine returns LW_TRUE if the point is inside the ring or on the boundary, LW_FALSE otherwise.
-* The pt_outside must be guaranteed to be outside the ring (use the geography_pt_outside() function
-* to derive one in postgis, or the gbox_pt_outside() function if you don't mind burning CPU cycles
-* building a gbox first).
-*/
-int ptarray_point_in_ring_winding(POINTARRAY *pa, POINT2D pt_to_test) 
-{
-	GEOGRAPHIC_POINT gpt;
-	GEOGRAPHIC_EDGE edge1, edge2;
-	POINT3D norm1, norm2;
-	POINT2D p;
-	double wind_number = 0;
-	int i;
-	
-	/* Null input, not enough points for a ring? You ain't closed! */
-	if( ! pa || pa->npoints < 4 ) 
-		return LW_FALSE;
-
-	/* Set up our test point */
-	geographic_point_init(pt_to_test.x, pt_to_test.y, &gpt);
-	edge1.start = gpt;
-	edge2.start = gpt;
-
-	/* Walk every edge and see if the stab line hits it */
-	for( i = 1; i < pa->npoints; i++ )
-	{
-		getPoint2d_p(pa, i-1, &p);
-		geographic_point_init(p.x, p.y, &(edge1.end));
-		getPoint2d_p(pa, i, &p);
-		geographic_point_init(p.x, p.y, &(edge2.end));
-
-		/* Calculate normals to the great circle planes */
-		robust_cross_product(edge1.start, edge1.end, &norm1);
-		robust_cross_product(edge2.start, edge2.end, &norm2);
-		normalize(&norm1);
-		normalize(&norm2);
-		
-		/* Add the wind */
-		wind_number += acos(dot_product(norm1, norm2));
-		
-	}
-
-	if( wind_number > M_PI )
-	{
-		return LW_TRUE;
-	}
-	
-	return LW_FALSE;
-}
-
 
 static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double tolerance)
 {
@@ -1447,7 +1453,7 @@ double lwgeom_distance_sphere(LWGEOM *lwgeom1, LWGEOM *lwgeom2, GBOX *gbox1, GBO
 	/* What's the distance to an empty geometry? We don't know. */
 	if( lwgeom_is_empty(lwgeom1) || lwgeom_is_empty(lwgeom2) )
 	{
-		return -1.0;
+		return 0.0;
 	}
 		
 	type1 = TYPE_GETTYPE(lwgeom1->type);