From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Fri, 9 Oct 2009 04:16:44 +0000 (+0000)
Subject: Add ST_PointOutside() function for testing purposes.
X-Git-Tag: 1.5.0b1~385
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=3657335216222395cf4ee68c45b26c63066fb3b8;p=postgis

Add ST_PointOutside() function for testing purposes.


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

diff --git a/liblwgeom/libgeom.h b/liblwgeom/libgeom.h
index f234e589b..2a0bc6321 100644
--- a/liblwgeom/libgeom.h
+++ b/liblwgeom/libgeom.h
@@ -393,6 +393,11 @@ extern int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point);
 */
 extern int ptarray_calculate_gbox_geodetic(POINTARRAY *pa, GBOX *gbox);
 
+/**
+* Calculate a spherical point that falls outside the geocentric gbox
+*/
+void gbox_pt_outside(GBOX gbox, POINT2D *pt_outside);
+
 /**
 * Create a new gbox with the dimensionality indicated by the flags. Caller
 * is responsible for freeing.
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index a0f60575b..3ac393b25 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -1198,16 +1198,18 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
 }
 
 /**
-* Given a gbox, return a cartesian unit vector to a point that is 
+* Given a unit geocentric gbox, return a lon/lat (degrees) coordinate point point that is 
 * guaranteed to be outside the box (and therefore anything it contains).
 */
-static void gbox_pt_outside(GBOX gbox, POINT3D *pt)
+void gbox_pt_outside(GBOX gbox, POINT2D *pt_outside)
 {
 	static double grow = M_PI / 180.0 / 60.0; /* one arc-minute */
 	double d;
 	int i;
 	GBOX ge;
 	POINT3D corners[8];
+	POINT3D pt;
+	GEOGRAPHIC_POINT g;
 	
 	/* Assign our box and expand it slightly. */
 	ge = gbox;
@@ -1256,40 +1258,49 @@ static void gbox_pt_outside(GBOX gbox, POINT3D *pt)
 		normalize(&(corners[i]));
 		if( ! gbox_contains_point3d(gbox, corners[i]) )
 		{
-			*pt = corners[i];
+			pt = corners[i];
+			normalize(&pt);
+			cart2geog(pt, &g);
+			pt_outside->x = rad2deg(g.lon);
+			pt_outside->y = rad2deg(g.lat);
 			return;
 		}	
 	}
 
-	pt->x = 1.0;
-	pt->y = 0.0;
-	pt->z = 0.0;
+	pt.x = 1.0;
+	pt.y = 0.0;
+	pt.z = 0.0;
 	
 	if((1.0 - gbox.xmax) > 0.1) 
 	{
-		pt->x = gbox.xmax + (1.0 - gbox.xmax) * 0.01;
-		d = sqrt((1.0 - pow(pt->x, 2.0))/2.0);
-		pt->y = d;
-		pt->z = d;
+		pt.x = gbox.xmax + (1.0 - gbox.xmax) * 0.01;
+		d = sqrt((1.0 - pow(pt.x, 2.0))/2.0);
+		pt.y = d;
+		pt.z = d;
 	} 
 	else if((1.0 - gbox.ymax) > 0.1) 
 	{
-		pt->y = gbox.ymax + (1.0 - gbox.ymax) * 0.01;
-		d = sqrt((1.0 - pow(pt->y, 2.0))/2.0);
-		pt->x = d;
-		pt->z = d;
+		pt.y = gbox.ymax + (1.0 - gbox.ymax) * 0.01;
+		d = sqrt((1.0 - pow(pt.y, 2.0))/2.0);
+		pt.x = d;
+		pt.z = d;
 	} 
 	else if((1.0 - gbox.zmax) > 0.1) 
 	{
-		pt->z = gbox.zmax + (1.0 - gbox.zmax) * 0.01;
-		d = sqrt((1.0 - pow(pt->z, 2.0))/2.0);
-		pt->x = d;
-		pt->y = d;
+		pt.z = gbox.zmax + (1.0 - gbox.zmax) * 0.01;
+		d = sqrt((1.0 - pow(pt.z, 2.0))/2.0);
+		pt.x = d;
+		pt.y = d;
 	}
-	normalize(pt);
+	normalize(&pt);
+	cart2geog(pt, &g);
+	pt_outside->x = rad2deg(g.lon);
+	pt_outside->y = rad2deg(g.lat);
 	return;
 }
 
+
+
 /**
 * Returns the area of the ring (ring must be closed) in square radians (surface of
 * the sphere is 4*PI).
@@ -1514,8 +1525,6 @@ static double ptarray_distance_sphere(POINTARRAY *pa1, POINTARRAY *pa2, double t
 double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox)
 {
 	int type;
-	POINT3D p;
-	GEOGRAPHIC_POINT g;
 	POINT2D pt_outside;
 	
 	assert(lwgeom);
@@ -1531,10 +1540,9 @@ double lwgeom_area_sphere(LWGEOM *lwgeom, GBOX gbox)
 	if( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
 		return 0.0;
 
-	gbox_pt_outside(gbox, &p);
-	cart2geog(p, &g);
-	pt_outside.x = rad2deg(g.lon);
-	pt_outside.y = rad2deg(g.lat);
+	gbox_pt_outside(gbox, &pt_outside);
+	
+	LWDEBUGF(2, "pt_outside = POINT(%.20g %.20g)", pt_outside.x, pt_outside.y);
 	
 	/* Actually calculate area */
 	if( type == POLYGONTYPE )
@@ -1800,8 +1808,8 @@ int lwpoly_covers_point2d(const LWPOLY *poly, GBOX gbox, POINT2D pt_to_test)
 {
 	int i;
 	int in_hole_count = 0;
-	POINT3D p, q;
-	GEOGRAPHIC_POINT g, gpt_to_test;
+	POINT3D p;
+	GEOGRAPHIC_POINT gpt_to_test;
 	POINT2D pt_outside;
 	
 	/* Nulls and empties don't contain anything! */
@@ -1813,15 +1821,12 @@ int lwpoly_covers_point2d(const LWPOLY *poly, GBOX gbox, POINT2D pt_to_test)
 
 	/* Point not in box? Done! */
 	geographic_point_init(pt_to_test.x, pt_to_test.y, &gpt_to_test);
-	geog2cart(gpt_to_test, &q);
-	if( ! gbox_contains_point3d(gbox, q) )
+	geog2cart(gpt_to_test, &p);
+	if( ! gbox_contains_point3d(gbox, p) )
 		return LW_FALSE;
 	
 	/* Calculate our outside point from the gbox */
-	gbox_pt_outside(gbox, &p);
-	cart2geog(p, &g);
-	pt_outside.x = rad2deg(g.lon);
-	pt_outside.y = rad2deg(g.lat);
+	gbox_pt_outside(gbox, &pt_outside);
 
 	LWDEBUGF(4, "pt_outside POINT(%.18g %.18g)", pt_outside.x, pt_outside.y);
 	LWDEBUGF(4, "pt_to_test POINT(%.18g %.18g)", pt_to_test.x, pt_to_test.y);
diff --git a/postgis/geography.sql.in.c b/postgis/geography.sql.in.c
index eb2eb1ac9..720790943 100644
--- a/postgis/geography.sql.in.c
+++ b/postgis/geography.sql.in.c
@@ -456,6 +456,12 @@ CREATE OR REPLACE FUNCTION ST_Area(geography)
 	AS 'MODULE_PATHNAME','geography_area_sphere'
 	LANGUAGE 'C' IMMUTABLE STRICT;
 
+-- Availability: 1.5.0
+CREATE OR REPLACE FUNCTION ST_PointOutside(geography)
+	RETURNS geography
+	AS 'MODULE_PATHNAME','geography_point_outside'
+	LANGUAGE 'C' IMMUTABLE STRICT;
+
 -- ---------- ---------- ---------- ---------- ---------- ---------- ----------
 
 COMMIT;
diff --git a/postgis/geography_measurement.c b/postgis/geography_measurement.c
index 6c02ad9f7..45efdd5f1 100644
--- a/postgis/geography_measurement.c
+++ b/postgis/geography_measurement.c
@@ -159,11 +159,11 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS)
 	area = lwgeom_area_sphere(lwgeom, gbox);
 
 	/* Something went wrong... */
-	if( area < 0.0 )
-	{
-		elog(ERROR, "lwgeom_area_sphere returned area < 0.0");
-		PG_RETURN_NULL();
-	}
+//	if( area < 0.0 )
+//	{
+//		elog(ERROR, "lwgeom_area_sphere returned area < 0.0");
+//		PG_RETURN_NULL();
+//	}
 
 	/* Currently normalizing with a fixed WGS84 radius, in future this
 	   should be the average radius of the SRID in play */
@@ -174,4 +174,40 @@ Datum geography_area_sphere(PG_FUNCTION_ARGS)
 	
 	PG_RETURN_FLOAT8(area);
 
+}
+
+/*
+** geography_point_outside(GSERIALIZED *g) 
+** returns point outside the object
+*/
+PG_FUNCTION_INFO_V1(geography_point_outside);
+Datum geography_point_outside(PG_FUNCTION_ARGS)
+{
+	GBOX gbox;
+	GSERIALIZED *g = NULL;
+	GSERIALIZED *g_out = NULL;
+	size_t g_out_size;
+	LWPOINT *lwpoint = NULL;
+	POINT2D pt;
+	
+	/* Get our geometry object loaded into memory. */
+	g = (GSERIALIZED*)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+	
+	/* We need the bounding box to get an outside point for area algorithm */
+	if( ! gbox_from_gserialized(g, &gbox) )
+	{
+		elog(ERROR, "Error in gbox_from_gserialized calculation.");
+		PG_RETURN_NULL();
+	}
+	
+	/* Get an exterior point, based on this gbox */
+	gbox_pt_outside(gbox, &pt);
+	
+	lwpoint = make_lwpoint2d(4326, pt.x, pt.y);
+	
+	g_out = gserialized_from_lwgeom((LWGEOM*)lwpoint, 1, &g_out_size);
+	SET_VARSIZE(g_out, g_out_size);
+	
+	PG_RETURN_POINTER(g_out);
+
 }
\ No newline at end of file