GRASS Programmer's Manual  6.5.svn(2014)-r66266
gis/distance.c
Go to the documentation of this file.
1
21 #include <math.h>
22 #include <grass/gis.h>
23 #include <grass/glocale.h>
24
25 static double min4(double, double, double, double);
26 static double min2(double, double);
27
28
29 static int projection = 0;
30 static double factor = 1.0;
31
32
45 {
46  double a, e2;
47
48  factor = 1.0;
49  switch (projection = G_projection()) {
50  case PROJECTION_LL:
53  return 2;
54  default:
56  if (factor <= 0.0) {
57  factor = 1.0; /* assume meter grid */
58  return 0;
59  }
60  return 1;
61  }
62 }
63
64
78 double G_distance(double e1, double n1, double e2, double n2)
79 {
80  if (projection == PROJECTION_LL)
81  return G_geodesic_distance(e1, n1, e2, n2);
82  else
83  return factor * hypot(e1 - e2, n1 - n2);
84 }
85
86
95 double
96 G_distance_between_line_segments(double ax1, double ay1,
97  double ax2, double ay2,
98  double bx1, double by1,
99  double bx2, double by2)
100 {
101  double ra, rb;
102  double x, y;
103
104  /* if the segments intersect, then the distance is zero */
105  if (G_intersect_line_segments(ax1, ay1, ax2, ay2,
106  bx1, by1, bx2, by2, &ra, &rb, &x, &y) > 0)
107  return 0.0;
108  return
109  min4(G_distance_point_to_line_segment(ax1, ay1, bx1, by1, bx2, by2),
110  G_distance_point_to_line_segment(ax2, ay2, bx1, by1, bx2, by2),
111  G_distance_point_to_line_segment(bx1, by1, ax1, ay1, ax2, ay2),
112  G_distance_point_to_line_segment(bx2, by2, ax1, ay1, ax2, ay2)
113  );
114 }
115
116
126 double G_distance_point_to_line_segment(double xp, double yp, /* the point */
127  double x1, double y1, double x2,
128  double y2)
129 { /* the line segment */
130  double dx, dy;
131  double x, y;
132  double xq, yq, ra, rb;
133  int t;
134
135  /* define the perpendicular to the segment thru the point */
136  dx = x1 - x2;
137  dy = y1 - y2;
138
139  if (dx == 0.0 && dy == 0.0)
140  return G_distance(x1, y1, xp, yp);
141
142  if (fabs(dy) > fabs(dx)) {
143  xq = xp + dy;
144  yq = (dx / dy) * (xp - xq) + yp;
145  }
146  else {
147  yq = yp + dx;
148  xq = (dy / dx) * (yp - yq) + xp;
149  }
150
151  /* find the intersection of the perpendicular with the segment */
152  switch (t =
153  G_intersect_line_segments(xp, yp, xq, yq, x1, y1, x2, y2, &ra,
154  &rb, &x, &y)) {
155  case 0:
156  case 1:
157  break;
158  default:
159  /* parallel/colinear cases shouldn't occur with perpendicular lines */
160  G_warning(_("G_distance_point_to_line_segment: shouldn't happen: "
161  "code=%d P=(%f,%f) S=(%f,%f)(%f,%f)"),
162  t, xp, yp, x1, y1, x2, y2);
163  return -1.0;
164  }
165
166  /* if x,y lies on the segment, then the distance is from to x,y */
167  if (rb >= 0 && rb <= 1.0)
168  return G_distance(x, y, xp, yp);
169
170  /* otherwise the distance is the short of the distances to the endpoints
171  * of the segment
172  */
173  return min2(G_distance(x1, y1, xp, yp), G_distance(x2, y2, xp, yp));
174 }
175
176 static double min4(double a, double b, double c, double d)
177 {
178  return min2(min2(a, b), min2(c, d));
179 }
180
181 static double min2(double a, double b)
182 {
183  return a < b ? a : b;
184 }
float b
Definition: named_colr.c:8
int G_intersect_line_segments(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2, double *ra, double *rb, double *x, double *y)
Definition: gis/intersect.c:80
int G_begin_distance_calculations(void)
Begin distance calculations.
Definition: gis/distance.c:44
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition: get_ellipse.c:66
int y
Definition: plot.c:34
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
double G_distance(double e1, double n1, double e2, double n2)
Returns distance in meters.
Definition: gis/distance.c:78
int G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition: geodist.c:52
double G_distance_between_line_segments(double ax1, double ay1, double ax2, double ay2, double bx1, double by1, double bx2, double by2)
Returns distance between two line segments in meters.
Definition: gis/distance.c:96
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition: geodist.c:212
double G_distance_point_to_line_segment(double xp, double yp, double x1, double y1, double x2, double y2)
Returns distance between a point and line segment in meters.
Definition: gis/distance.c:126
double G_database_units_to_meters_factor(void)
conversion to meters
Definition: proj3.c:80
int G_projection(void)
query cartographic projection
Definition: proj1.c:33