GRASS Programmer's Manual  6.4.2(2012)
gis/distance.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines