From 0c179ca6e8231769fc9c7a2963c01207021a89a7 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Wed, 30 Sep 2009 21:45:58 +0000 Subject: [PATCH] Add return value for point of closest approach on arc/edge distance. git-svn-id: http://svn.osgeo.org/postgis/trunk@4567 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/cunit/cu_geodetic.c | 92 +++++++++++++++++------------------ liblwgeom/lwgeodetic.c | 24 +++++++-- liblwgeom/lwgeodetic.h | 2 +- 3 files changed, 67 insertions(+), 51 deletions(-) diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index 84a34cf1b..f9b160815 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -142,48 +142,6 @@ void test_gbox_from_spherical_coordinates(void) #endif /* RANDOM_TEST */ } - -void test_clairaut(void) -{ - - GEOGRAPHIC_POINT gs, ge; - POINT3D vs, ve; - GEOGRAPHIC_POINT g_out_top, g_out_bottom, v_out_top, v_out_bottom; - - gs.lat = deg2rad(60.0); - gs.lon = deg2rad(-45.0); - ge.lat = deg2rad(60.0); - ge.lon = deg2rad(135.0); - - geog2cart(gs, &vs); - geog2cart(ge, &ve); - - 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_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; - gs.lon = -1.3387392931438733; - ge.lon = 1.80285336044592; - - geog2cart(gs, &vs); - geog2cart(ge, &ve); - - 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_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); -} - - #include "cu_geodetic_data.h" void test_gserialized_get_gbox_geocentric(void) @@ -259,6 +217,43 @@ static void point_set(double lon, double lat, GEOGRAPHIC_POINT *p) point_deg2rad(p); } +void test_clairaut(void) +{ + + GEOGRAPHIC_POINT gs, ge; + POINT3D vs, ve; + GEOGRAPHIC_POINT g_out_top, g_out_bottom, v_out_top, v_out_bottom; + + point_set(-45.0, 60.0, &gs); + point_set(135.0, 60.0, &gs); + + geog2cart(gs, &vs); + geog2cart(ge, &ve); + + 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_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; + gs.lon = -1.3387392931438733; + ge.lon = 1.80285336044592; + + geog2cart(gs, &vs); + geog2cart(ge, &ve); + + 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_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); +} void test_edge_intersection(void) { @@ -367,22 +362,27 @@ void test_edge_distance_to_point(void) { GEOGRAPHIC_EDGE e; GEOGRAPHIC_POINT g; + GEOGRAPHIC_POINT closest; double d; + /* closest point at origin, one degree away */ edge_set(-50.0, 0.0, 50.0, 0.0, &e); point_set(0.0, 1.0, &g); - d = edge_distance_to_point(e, g); - CU_ASSERT_DOUBLE_EQUAL(d, 0.0, M_PI / 180.0); + d = edge_distance_to_point(e, g, 0); + CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 180.0, 0.00001); + /* closest point at origin, one degree away */ edge_set(-50.0, 0.0, 50.0, 0.0, &e); point_set(0.0, 2.0, &g); - d = edge_distance_to_point(e, g); + d = edge_distance_to_point(e, g, &closest); #if 0 printf("LINESTRING(%.8g %.8g, %.8g %.8g)\n", e.start.lon, e.start.lat, e.end.lon, e.end.lat); printf("POINT(%.9g %.9g)\n", g.lon, g.lat); printf("\nDISTANCE == %.8g\n", d); #endif - CU_ASSERT_DOUBLE_EQUAL(d, 0.0, M_PI / 90.0); + CU_ASSERT_DOUBLE_EQUAL(d, M_PI / 90.0, 0.00001); + CU_ASSERT_DOUBLE_EQUAL(closest.lat, 0.0, 0.00001); + CU_ASSERT_DOUBLE_EQUAL(closest.lon, 0.0, 0.00001); } diff --git a/liblwgeom/lwgeodetic.c b/liblwgeom/lwgeodetic.c index 2561022af..218adf879 100644 --- a/liblwgeom/lwgeodetic.c +++ b/liblwgeom/lwgeodetic.c @@ -612,11 +612,11 @@ int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT * return LW_FALSE; } -double edge_distance_to_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT gp) +double edge_distance_to_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT gp, GEOGRAPHIC_POINT *closest) { - double d1 = 1000000000.0, d2, d3; + double d1 = 1000000000.0, d2, d3, d_nearest; POINT3D n, p, k, v1; - GEOGRAPHIC_POINT gk; + GEOGRAPHIC_POINT gk, g_nearest; /* Zero length edge, */ if( point_equal(e.start,e.end) ) @@ -636,7 +636,23 @@ double edge_distance_to_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT gp) d2 = sphere_distance(gp, e.start); d3 = sphere_distance(gp, e.end); - return FP_MIN(FP_MIN(d1,d2),d3); + d_nearest = d1; + g_nearest = gk; + + if( d2 < d_nearest ) + { + d_nearest = d2; + g_nearest = e.start; + } + if( d3 < d_nearest ) + { + d_nearest = d3; + g_nearest = e.end; + } + if(closest) + *closest = g_nearest; + + return d_nearest; } /** diff --git a/liblwgeom/lwgeodetic.h b/liblwgeom/lwgeodetic.h index 02ce86b6d..af5a29403 100644 --- a/liblwgeom/lwgeodetic.h +++ b/liblwgeom/lwgeodetic.h @@ -64,7 +64,7 @@ int sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPH int edge_calculate_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox); int edge_calculate_gbox_slow(GEOGRAPHIC_EDGE e, GBOX *gbox); int edge_intersection(GEOGRAPHIC_EDGE e1, GEOGRAPHIC_EDGE e2, GEOGRAPHIC_POINT *g); -double edge_distance_to_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT gp); +double edge_distance_to_point(GEOGRAPHIC_EDGE e, GEOGRAPHIC_POINT gp, GEOGRAPHIC_POINT *closest); void edge_deg2rad(GEOGRAPHIC_EDGE *e); void edge_rad2deg(GEOGRAPHIC_EDGE *e); void point_deg2rad(GEOGRAPHIC_POINT *p); -- 2.49.0