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 8f03c5ba62916771c5cee5964f2046c66eb35a28..9966731f8e3271c5196983e15809a1539139f17b 100644 (file)
--- 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 295fa75617ef9456dce3b254ed784fa637882ec3..ebc77b9eb24a5ac40c7f8c287322d532f5e02b81 100644 (file)
--- 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 627f673e14392a51cd6aaab97094486727b11042..acd10251fb5da1aba0bd4e39331621ee00716d08 100644 (file)
--- 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.