GRASS Programmer's Manual 6.4.1(2011)
|
00001 00021 #include <math.h> 00022 #include <grass/gis.h> 00023 #include <grass/glocale.h> 00024 00025 static double min4(double, double, double, double); 00026 static double min2(double, double); 00027 00028 00029 static int projection = 0; 00030 static double factor = 1.0; 00031 00032 00044 int G_begin_distance_calculations(void) 00045 { 00046 double a, e2; 00047 00048 factor = 1.0; 00049 switch (projection = G_projection()) { 00050 case PROJECTION_LL: 00051 G_get_ellipsoid_parameters(&a, &e2); 00052 G_begin_geodesic_distance(a, e2); 00053 return 2; 00054 default: 00055 factor = G_database_units_to_meters_factor(); 00056 if (factor <= 0.0) { 00057 factor = 1.0; /* assume meter grid */ 00058 return 0; 00059 } 00060 return 1; 00061 } 00062 } 00063 00064 00078 double G_distance(double e1, double n1, double e2, double n2) 00079 { 00080 if (projection == PROJECTION_LL) 00081 return G_geodesic_distance(e1, n1, e2, n2); 00082 else 00083 return factor * hypot(e1 - e2, n1 - n2); 00084 } 00085 00086 00095 double 00096 G_distance_between_line_segments(double ax1, double ay1, 00097 double ax2, double ay2, 00098 double bx1, double by1, 00099 double bx2, double by2) 00100 { 00101 double ra, rb; 00102 double x, y; 00103 00104 /* if the segments intersect, then the distance is zero */ 00105 if (G_intersect_line_segments(ax1, ay1, ax2, ay2, 00106 bx1, by1, bx2, by2, &ra, &rb, &x, &y) > 0) 00107 return 0.0; 00108 return 00109 min4(G_distance_point_to_line_segment(ax1, ay1, bx1, by1, bx2, by2), 00110 G_distance_point_to_line_segment(ax2, ay2, bx1, by1, bx2, by2), 00111 G_distance_point_to_line_segment(bx1, by1, ax1, ay1, ax2, ay2), 00112 G_distance_point_to_line_segment(bx2, by2, ax1, ay1, ax2, ay2) 00113 ); 00114 } 00115 00116 00126 double G_distance_point_to_line_segment(double xp, double yp, /* the point */ 00127 double x1, double y1, double x2, 00128 double y2) 00129 { /* the line segment */ 00130 double dx, dy; 00131 double x, y; 00132 double xq, yq, ra, rb; 00133 int t; 00134 00135 /* define the perpendicular to the segment thru the point */ 00136 dx = x1 - x2; 00137 dy = y1 - y2; 00138 00139 if (dx == 0.0 && dy == 0.0) 00140 return G_distance(x1, y1, xp, yp); 00141 00142 if (fabs(dy) > fabs(dx)) { 00143 xq = xp + dy; 00144 yq = (dx / dy) * (xp - xq) + yp; 00145 } 00146 else { 00147 yq = yp + dx; 00148 xq = (dy / dx) * (yp - yq) + xp; 00149 } 00150 00151 /* find the intersection of the perpendicular with the segment */ 00152 switch (t = 00153 G_intersect_line_segments(xp, yp, xq, yq, x1, y1, x2, y2, &ra, 00154 &rb, &x, &y)) { 00155 case 0: 00156 case 1: 00157 break; 00158 default: 00159 /* parallel/colinear cases shouldn't occur with perpendicular lines */ 00160 G_warning(_("G_distance_point_to_line_segment: shouldn't happen: " 00161 "code=%d P=(%f,%f) S=(%f,%f)(%f,%f)"), 00162 t, xp, yp, x1, y1, x2, y2); 00163 return -1.0; 00164 } 00165 00166 /* if x,y lies on the segment, then the distance is from to x,y */ 00167 if (rb >= 0 && rb <= 1.0) 00168 return G_distance(x, y, xp, yp); 00169 00170 /* otherwise the distance is the short of the distances to the endpoints 00171 * of the segment 00172 */ 00173 return min2(G_distance(x1, y1, xp, yp), G_distance(x2, y2, xp, yp)); 00174 } 00175 00176 static double min4(double a, double b, double c, double d) 00177 { 00178 return min2(min2(a, b), min2(c, d)); 00179 } 00180 00181 static double min2(double a, double b) 00182 { 00183 return a < b ? a : b; 00184 }