From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Thu, 25 Oct 2012 22:17:39 +0000 (+0000)
Subject: Expose some geodetic functions a little higher
X-Git-Tag: 2.1.0beta2~452
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=4252d6b53f89ab0a114c583f949bb05b7b3059d0;p=postgis

Expose some geodetic functions a little higher


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

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 8f03c5ba6..9966731f8 100644
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -149,6 +149,58 @@ static void test_sphere_project(void)
 	CU_ASSERT_DOUBLE_EQUAL(dir2, -2.35612, 0.0001);	
 }
 
+#if 0
+/**
+* Tests the relative numerical stability of the "robust" and
+* naive cross product calculation methods.
+*/
+static void cross_product_stability(void)
+{
+	POINT2D p1, p2;
+	int i;
+	GEOGRAPHIC_POINT g1, g2;
+	POINT3D A1, A2;
+	POINT3D Nr, Nc;
+	POINT3D Or, Oc;
+
+	p1.x = 10.0;
+	p1.y = 45.0;
+	p2.x = 10.0;
+	p2.y = 50.0;
+	
+	geographic_point_init(p1.x, p1.y, &g1);
+	ll2cart(&p1, &A1);
+
+	for ( i = 0; i < 40; i++ )
+	{
+		geographic_point_init(p2.x, p2.y, &g2);
+		ll2cart(&p2, &A2);
+		
+		/* Skea */
+		robust_cross_product(&g1, &g2, &Nr);
+		normalize(&Nr);
+		
+		/* Ramsey */
+		unit_normal(&A1, &A2, &Nc);
+
+		if ( i > 0 ) 
+		{
+			printf("\n- %d -------------------- %.24g ------------------------\n", i, p2.y);
+			printf("Skea:         %.24g,%.24g,%.24g\n", Nr.x, Nr.y, Nr.z);
+			printf("Skea Diff:    %.24g,%.24g,%.24g\n", Or.x-Nr.x, Or.y-Nr.y, Or.z-Nr.z);
+			printf("Ramsey:       %.24g,%.24g,%.24g\n", Nc.x, Nc.y, Nc.z);
+			printf("Ramsey Diff:  %.24g,%.24g,%.24g\n", Oc.x-Nc.x, Oc.y-Nc.y, Oc.z-Nc.z);
+			printf("Diff:         %.24g,%.24g,%.24g\n", Nr.x-Nc.x, Nr.y-Nc.y, Nr.z-Nc.z);
+		}
+		
+		Or = Nr;
+		Oc = Nc;
+		
+		p2.y += (p1.y - p2.y)/2.0;
+	}
+}
+#endif
+
 static void test_gbox_from_spherical_coordinates(void)
 {
 #if RANDOM_TEST
@@ -247,12 +299,12 @@ static void test_gserialized_get_gbox_geocentric(void)
 		printf("line %d: diff %.9g\n", i, fabs(gbox.xmin - gbox_slow.xmin)+fabs(gbox.ymin - gbox_slow.ymin)+fabs(gbox.zmin - gbox_slow.zmin));
 		printf("------------\n");
 #endif
-		CU_ASSERT_DOUBLE_EQUAL(gbox.xmin, gbox_slow.xmin, 0.000001);
-		CU_ASSERT_DOUBLE_EQUAL(gbox.ymin, gbox_slow.ymin, 0.000001);
-		CU_ASSERT_DOUBLE_EQUAL(gbox.zmin, gbox_slow.zmin, 0.000001);
-		CU_ASSERT_DOUBLE_EQUAL(gbox.xmax, gbox_slow.xmax, 0.000001);
-		CU_ASSERT_DOUBLE_EQUAL(gbox.ymax, gbox_slow.ymax, 0.000001);
-		CU_ASSERT_DOUBLE_EQUAL(gbox.zmax, gbox_slow.zmax, 0.000001);
+		CU_ASSERT_DOUBLE_EQUAL(gbox.xmin, gbox_slow.xmin, 0.00000001);
+		CU_ASSERT_DOUBLE_EQUAL(gbox.ymin, gbox_slow.ymin, 0.00000001);
+		CU_ASSERT_DOUBLE_EQUAL(gbox.zmin, gbox_slow.zmin, 0.00000001);
+		CU_ASSERT_DOUBLE_EQUAL(gbox.xmax, gbox_slow.xmax, 0.00000001);
+		CU_ASSERT_DOUBLE_EQUAL(gbox.ymax, gbox_slow.ymax, 0.00000001);
+		CU_ASSERT_DOUBLE_EQUAL(gbox.zmax, gbox_slow.zmax, 0.00000001);
 	}
 
 }
diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index 295fa7561..ebc77b9eb 100644
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -367,7 +367,7 @@ void cart2geog(const POINT3D *p, GEOGRAPHIC_POINT *g)
 /**
 * Convert lon/lat coordinates to cartesion coordinates on unit sphere
 */
-static void ll2cart(const POINT2D *g, POINT3D *p)
+void ll2cart(const POINT2D *g, POINT3D *p)
 {
 	double x_rad = M_PI * g->x / 180.0;
 	double y_rad = M_PI * g->y / 180.0;
@@ -379,12 +379,12 @@ static void ll2cart(const POINT2D *g, POINT3D *p)
 
 /**
 * Convert cartesion coordinates on unit sphere to lon/lat coordinates 
-*/
 static void cart2ll(const POINT3D *p, POINT2D *g)
 {
 	g->x = longitude_degrees_normalize(180.0 * atan2(p->y, p->x) / M_PI);
 	g->y = latitude_degrees_normalize(180.0 * asin(p->z) / M_PI);
 }
+*/
 
 /**
 * Calculate the dot product of two unit vectors
@@ -483,7 +483,7 @@ static void normalize2d(POINT2D *p)
 * Calculates the unit normal to two vectors, trying to avoid
 * problems with over-narrow or over-wide cases.
 */
-static void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
+void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal)
 {
 	double p_dot = dot_product(P1, P2);
 	POINT3D P3;
diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index 627f673e1..acd10251f 100644
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -119,7 +119,9 @@ void vector_sum(const POINT3D *a, const POINT3D *b, POINT3D *n);
 double vector_angle(const POINT3D* v1, const POINT3D* v2);
 void vector_rotate(const POINT3D* v1, const POINT3D* v2, double angle, POINT3D* n);
 void normalize(POINT3D *p);
+void unit_normal(const POINT3D *P1, const POINT3D *P2, POINT3D *normal);
 double sphere_direction(const GEOGRAPHIC_POINT *s, const GEOGRAPHIC_POINT *e, double d);
+void ll2cart(const POINT2D *g, POINT3D *p);
 
 /*
 ** Prototypes for spheroid functions.