Flip the clairaut calculations to return both top and bottom in one go.

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

diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c
index 356c0a1ad5e463ee2e9d0f6d467feda43925a957..f4a6c002c3275e7a1b062039b2207b9a88373f53 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic.c
+++ b/liblwgeom/cunit/cu_geodetic.c
@@ -68,8 +68,7 @@ void test_signum(void)
 void test_gbox_from_spherical_coordinates(void)
 {
     const double gtolerance = 0.000001;
-    const int loops = 0;
-//    const int loops = 0;
+    const int loops = 20;
     int i;
        double ll[64];
        GBOX *gbox;
@@ -121,6 +120,7 @@ void test_gbox_from_spherical_coordinates(void)
             ( fabs( gbox->zmax - gbox_slow->zmax ) > gtolerance ) )
         {
             printf("\n-------\n");
+            printf("If you are seeing this, cut and paste it, it is a randomly generated test case!\n");
             printf("LOOP: %d\n", i);
             printf("SEGMENT (Lon Lat): (%.9g %.9g) (%.9g %.9g)\n", ll[0], ll[1], ll[2], ll[3]);
             printf("CALC: %s\n", gbox_to_string(gbox));
@@ -143,7 +143,7 @@ void test_clairaut(void)
   
     GEOGRAPHIC_POINT gs, ge;
     POINT3D vs, ve;
-    GEOGRAPHIC_POINT g_out, v_out;  
+    GEOGRAPHIC_POINT g_out_top, g_out_bottom, v_out_top, v_out_bottom;  
 
     gs.lat = deg2rad(60.0);
     gs.lon = deg2rad(-45.0);
@@ -153,17 +153,13 @@ void test_clairaut(void)
     geog2cart(gs, &vs);
     geog2cart(ge, &ve);
     
-    clairaut_cartesian(vs, ve, LW_TRUE, &v_out);
-    clairaut_geographic(gs, ge, LW_TRUE, &g_out);    
+    clairaut_cartesian(vs, ve, &v_out_top, &v_out_bottom);
+    clairaut_geographic(gs, ge, &g_out_top, &g_out_bottom);    
 
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lat, g_out.lat, 0.000001);
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lon, g_out.lon, 0.000001);
-
-    clairaut_cartesian(vs, ve, LW_FALSE, &v_out);
-    clairaut_geographic(gs, ge, LW_FALSE, &g_out);    
-
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lat, g_out.lat, 0.000001);
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lon, g_out.lon, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lat, g_out_bottom.lat, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lon, g_out_bottom.lon, 0.000001);
 
     gs.lat = 1.3021240033804449;
     ge.lat = 1.3021240033804449;
@@ -173,17 +169,13 @@ void test_clairaut(void)
     geog2cart(gs, &vs);
     geog2cart(ge, &ve);
 
-    clairaut_cartesian(vs, ve, LW_TRUE, &v_out);
-    clairaut_geographic(gs, ge, LW_TRUE, &g_out);    
-
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lat, g_out.lat, 0.000001);
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lon, g_out.lon, 0.000001);
+    clairaut_cartesian(vs, ve, &v_out_top, &v_out_bottom);
+    clairaut_geographic(gs, ge, &g_out_top, &g_out_bottom);    
 
-    clairaut_cartesian(vs, ve, LW_FALSE, &v_out);
-    clairaut_geographic(gs, ge, LW_FALSE, &g_out);    
-
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lat, g_out.lat, 0.000001);
-       CU_ASSERT_DOUBLE_EQUAL(v_out.lon, g_out.lon, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_top.lat, g_out_top.lat, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_top.lon, g_out_top.lon, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lat, g_out_bottom.lat, 0.000001);
+       CU_ASSERT_DOUBLE_EQUAL(v_out_bottom.lon, g_out_bottom.lon, 0.000001);
 }
     
 
@@ -198,11 +190,11 @@ void test_gserialized_get_gbox_geocentric(void)
        
        for ( i = 0; i < gbox_data_length; i++ )
        {
-//             if( i != 25 && i != 35 ) continue; /* skip our one bad case */
-//             printf("\n------------\n");
-//             printf("%s\n", gbox_data[2*i]);
-//             printf("%s\n", gbox_data[2*i+1]);
-               lwg = lwgeom_from_ewkt(gbox_data[i], PARSER_CHECK_NONE);
+               if( i != 0 ) continue; /* skip our bad case */
+/* 
+               printf("\n\n------------\n");
+               printf("%s\n", gbox_data[i]);
+*/             lwg = lwgeom_from_ewkt(gbox_data[i], PARSER_CHECK_NONE);
                g = gserialized_from_lwgeom(lwg, 1, 0);
                g->flags = FLAGS_SET_GEODETIC(g->flags, 1);
                lwgeom_free(lwg);
@@ -211,11 +203,12 @@ void test_gserialized_get_gbox_geocentric(void)
                gbox_geocentric_slow = LW_TRUE;
                gbox_slow = gserialized_calculate_gbox_geocentric(g);
                gbox_geocentric_slow = LW_FALSE;
-//             printf("\nCALC: %s\n", gbox_to_string(gbox));
-//             printf("GOOD: %s\n", gbox_to_string(gbox_good));
-//             printf("GOOD CALC: %s\n", gbox_to_string(gbox_good_calc));
-//             printf("line %d: diff %.9g\n", i, fabs(gbox->xmin - gbox_good->xmin)+fabs(gbox->ymin - gbox_good->ymin)+fabs(gbox->zmin - gbox_good->zmin));
-               CU_ASSERT_DOUBLE_EQUAL(gbox->xmin, gbox_slow->xmin, 0.000001);
+/* 
+               printf("\nCALC: %s\n", gbox_to_string(gbox));
+               printf("GOOD: %s\n", gbox_to_string(gbox_slow));
+               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");
+*/             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);

diff --git a/liblwgeom/cunit/cu_geodetic_data.h b/liblwgeom/cunit/cu_geodetic_data.h
index 5c980adff39afe8f98d64d09a68377a35b96aae4..a4867b75d96e2e092f8d887c14df37b34aeead94 100644 (file)
--- a/liblwgeom/cunit/cu_geodetic_data.h
+++ b/liblwgeom/cunit/cu_geodetic_data.h
@@ -1,5 +1,6 @@
-int gbox_data_length = 47;
+int gbox_data_length = 48;
 char gbox_data[][512] = { 
+"LINESTRING(158 -85,-57 86)",
 "LINESTRING(-3.083333333333333333333333333333333 9.83333333333333333333333333333333,15.5 -5.25)",
 "LINESTRING(-35.0 52.5,50.0 60.0)",
 "LINESTRING(-122.5 56.25,-123.5 69.166666)",

diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c
index c0ec1c00184583e5a27a064b2e567eda76b7665a..5b004cc1550d2f169ab849f8eab16990d53e6910 100644 (file)
--- a/liblwgeom/lwgeodetic.c
+++ b/liblwgeom/lwgeodetic.c
@@ -217,6 +217,7 @@ int edge_point_on_plane(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
        robust_cross_product(e.start, e.end, &normal);
        geog2cart(p, &pt);
        w = dot_product(normal, pt);
+    LWDEBUGF(4,"dot product %.9g",w);
        if( FP_IS_ZERO(w) )
        {
                LWDEBUG(4, "point is on plane");
@@ -240,8 +241,11 @@ int edge_point_in_cone(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p)
        vector_sum(vs, ve, &vcp);
        normalize(&vcp);
        vs_dot_vcp = dot_product(vs, vcp);
+    LWDEBUGF(4,"vs_dot_vcp %.9g",vs_dot_vcp);
        ve_dot_vcp = dot_product(ve, vcp);
+    LWDEBUGF(4,"ve_dot_vcp %.9g",ve_dot_vcp);
        vp_dot_vcp = dot_product(vp, vcp);
+    LWDEBUGF(4,"vp_dot_vcp %.9g",vp_dot_vcp);
 
        if( vp_dot_vcp > ve_dot_vcp && vp_dot_vcp > ve_dot_vcp )
        {
@@ -362,7 +366,7 @@ double z_to_latitude(double z, int top)
 * the great circle with the line of maximum/minimum gradiant that lies on
 * the great circle plane.
 */
-int clairaut_cartesian(POINT3D start, POINT3D end, int top, GEOGRAPHIC_POINT *g)
+int clairaut_cartesian(POINT3D start, POINT3D end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
 {
        POINT3D t1, t2;
        GEOGRAPHIC_POINT vN1, vN2;
@@ -373,17 +377,12 @@ int clairaut_cartesian(POINT3D start, POINT3D end, int top, GEOGRAPHIC_POINT *g)
     LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
        cart2geog(t1, &vN1);
        cart2geog(t2, &vN2);
-       if( top )
-       {
-               g->lat = z_to_latitude(t1.z,top);
-               g->lon = vN2.lon;
-       }
-       else
-       {
-               g->lat = z_to_latitude(t2.z,top);
-               g->lon = vN1.lon;
-       }
-    LWDEBUGF(4, "clairaut == GPOINT(%.6g %.6g)", g->lat, g->lon);
+       g_top->lat = z_to_latitude(t1.z,LW_TRUE);
+       g_top->lon = vN2.lon;
+       g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
+       g_bottom->lon = vN1.lon;
+    LWDEBUGF(4, "clairaut top == GPOINT(%.6g %.6g)", g_top->lat, g_top->lon);
+    LWDEBUGF(4, "clairaut bottom == GPOINT(%.6g %.6g)", g_bottom->lat, g_bottom->lon);
        return G_SUCCESS;
 }
 
@@ -392,29 +391,23 @@ int clairaut_cartesian(POINT3D start, POINT3D end, int top, GEOGRAPHIC_POINT *g)
 * the great circle with the line of maximum/minimum gradiant that lies on
 * the great circle plane.
 */
-int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, int top, GEOGRAPHIC_POINT *g)
+int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom)
 {
        POINT3D t1, t2;
        GEOGRAPHIC_POINT vN1, vN2;
     LWDEBUG(4,"entering function");
     robust_cross_product(start, end, &t1);
     robust_cross_product(end, start, &t2);
-    LWDEBUGF(4, "t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
-    LWDEBUGF(4, "t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
+    LWDEBUGF(4, "unit normal t1 == POINT(%.8g %.8g %.8g)", t1.x, t1.y, t1.z);
+    LWDEBUGF(4, "unit normal t2 == POINT(%.8g %.8g %.8g)", t2.x, t2.y, t2.z);
     cart2geog(t1, &vN1);
     cart2geog(t2, &vN2);
-        LWDEBUGF(4, "vN1 == GPOINT(%.6g %.6g) ", vN1.lat, vN1.lon);
-        LWDEBUGF(4, "vN2 == GPOINT(%.6g %.6g) ", vN2.lat, vN2.lon);
-    if( top )
-    {
-        g->lat = z_to_latitude(t1.z,top);
-        g->lon = vN2.lon;
-    }
-    else
-    {
-        g->lat = z_to_latitude(t2.z,top);
-        g->lon = vN1.lon;
-    }
+    g_top->lat = z_to_latitude(t1.z,LW_TRUE);
+    g_top->lon = vN2.lon;
+    g_bottom->lat = z_to_latitude(t2.z,LW_FALSE);
+    g_bottom->lon = vN1.lon;
+    LWDEBUGF(4, "clairaut top == GPOINT(%.6g %.6g)", g_top->lat, g_top->lon);
+    LWDEBUGF(4, "clairaut bottom == GPOINT(%.6g %.6g)", g_bottom->lat, g_bottom->lon);
     return G_SUCCESS;
 }
 
@@ -469,7 +462,7 @@ int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPH
 
 int edge_calculate_gbox_slow(GEOGRAPHIC_EDGE e, GBOX *gbox)
 {
-       int steps = 10000;
+       int steps = 1000000;
        int i;
        double dx, dy, dz;
        double distance = sphere_distance(e.start, e.end);
@@ -543,10 +536,8 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
        /* We're testing, do this the slow way. */
        if(gbox_geocentric_slow) 
        {
-               //printf("\n--- using SLOW calc! ---\n");
                return edge_calculate_gbox_slow(e, gbox);
        }
-       //printf("\n=== using FAST calc! ===\n");
 
        /* Initialize our working copy of the edge */
        g = e;
@@ -555,7 +546,6 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
        LWDEBUGF(4, "edge length: %.8g", distance);
     LWDEBUGF(4, "edge values: (%.6g %.6g, %.6g %.6g)", g.start.lon, g.start.lat, g.end.lon, g.end.lat);
 
-
        /* Edge is zero length, just return the naive box */
        if( FP_IS_ZERO(distance) )
        {
@@ -662,8 +652,7 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
        else
        {
         LWDEBUG(4, "not a pole crossing, calculating clairaut points");
-               clairaut_cartesian(start, end, LW_TRUE, &vT1);
-               clairaut_cartesian(start, end, LW_FALSE, &vT2);
+               clairaut_cartesian(start, end, &vT1, &vT2);
         LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
         LWDEBUGF(4, "vT2 == GPOINT(%.6g %.6g) ", vT2.lat, vT2.lon);
                if( edge_contains_point(g, vT1) )
@@ -690,8 +679,7 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
     endXZ = end;
     x_to_z(&startXZ);
     x_to_z(&endXZ);
-       clairaut_cartesian(startXZ, endXZ, LW_TRUE, &vT1);
-       clairaut_cartesian(startXZ, endXZ, LW_FALSE, &vT2);
+       clairaut_cartesian(startXZ, endXZ, &vT1, &vT2);
     gimbal_lock = LW_FALSE;
     LWDEBUG(4, "vT1/vT2 before flipping back z to x");
     LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
@@ -717,6 +705,12 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
         LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
         LWDEBUGF(4, "vT2 == GPOINT(%.6g %.6g) ", vT2.lat, vT2.lon);
     }
+/* For extra logging if needed
+    geog2cart(vT1, &nT1); 
+    geog2cart(vT2, &nT2);
+    LWDEBUGF(4, "p1 == POINT(%.8g %.8g %.8g)", nT1.x, nT1.y, nT1.z);
+    LWDEBUGF(4, "p2 == POINT(%.8g %.8g %.8g)", nT2.x, nT2.y, nT2.z);
+*/
     if( edge_contains_point(g, vT1) )
     {
                geog2cart(vT1, &p);
@@ -740,8 +734,7 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
     endYZ = end;
     y_to_z(&startYZ);
     y_to_z(&endYZ);
-       clairaut_cartesian(startYZ, endYZ, LW_TRUE, &vT1);
-       clairaut_cartesian(startYZ, endYZ, LW_FALSE, &vT2);
+       clairaut_cartesian(startYZ, endYZ, &vT1, &vT2);
     gimbal_lock = LW_FALSE;
     LWDEBUG(4, "vT1/vT2 before flipping back z to y");
     LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
@@ -767,6 +760,12 @@ int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
         LWDEBUGF(4, "vT1 == GPOINT(%.6g %.6g) ", vT1.lat, vT1.lon);
         LWDEBUGF(4, "vT2 == GPOINT(%.6g %.6g) ", vT2.lat, vT2.lon);
     }
+/* For extra logging if needed
+    geog2cart(vT1, &nT1);
+    geog2cart(vT2, &nT2);
+    LWDEBUGF(4, "p1 == POINT(%.8g %.8g %.8g)", nT1.x, nT1.y, nT1.z);
+    LWDEBUGF(4, "p2 == POINT(%.8g %.8g %.8g)", nT2.x, nT2.y, nT2.z);
+*/
     if( edge_contains_point(g, vT1) )
     {
                geog2cart(vT1, &p);

diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h
index cffcca198bb577eacd9f373e26a0f806fbff24ae..47a8b62522ad4220b391dd56aa7533e26562c0d0 100644 (file)
--- a/liblwgeom/lwgeodetic.h
+++ b/liblwgeom/lwgeodetic.h
@@ -54,8 +54,8 @@ int edge_point_in_cone(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT p);
 int edge_contains_longitude(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, int top, GEOGRAPHIC_POINT *g);
-int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, int top, GEOGRAPHIC_POINT *g);
+int clairaut_cartesian(POINT3D start, POINT3D end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);
+int clairaut_geographic(GEOGRAPHIC_POINT start, GEOGRAPHIC_POINT end, GEOGRAPHIC_POINT *g_top, GEOGRAPHIC_POINT *g_bottom);
 double sphere_distance(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
 double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
 int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n);