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;
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
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,
00127 double x1, double y1, double x2,
00128 double y2)
00129 {
00130 double dx, dy;
00131 double x, y;
00132 double xq, yq, ra, rb;
00133 int t;
00134
00135
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
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
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
00167 if (rb >= 0 && rb <= 1.0)
00168 return G_distance(x, y, xp, yp);
00169
00170
00171
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 }