From: Paul Ramsey Date: Tue, 29 Sep 2009 20:41:04 +0000 (+0000) Subject: Flip the clairaut calculations to return both top and bottom in one go. X-Git-Tag: 1.5.0b1~459 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=a71f3b2d8341d2e6009c46025b41f6b9c8d18ef4;p=postgis 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 356c0a1ad..f4a6c002c 100644 --- 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 5c980adff..a4867b75d 100644 --- 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 c0ec1c001..5b004cc15 100644 --- 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 cffcca198..47a8b6252 100644 --- 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);