GRASS GIS 7 Programmer's Manual  7.0.4(2016)-r00000
gis/distance.c
Go to the documentation of this file.
1 
18 #include <math.h>
19 #include <grass/gis.h>
20 #include <grass/glocale.h>
21 
22 static double min4(double, double, double, double);
23 static double min2(double, double);
24 
25 static struct state {
26  int projection;
27  double factor;
28 } state;
29 
30 static struct state *st = &state;
31 
43 {
44  double a, e2;
45 
46  st->factor = 1.0;
47  switch (st->projection = G_projection()) {
48  case PROJECTION_LL:
51  return 2;
52  default:
53  st->factor = G_database_units_to_meters_factor();
54  if (st->factor <= 0.0) {
55  st->factor = 1.0; /* assume meter grid */
56  return 0;
57  }
58  return 1;
59  }
60 }
61 
75 double G_distance(double e1, double n1, double e2, double n2)
76 {
77  if (st->projection == PROJECTION_LL)
78  return G_geodesic_distance(e1, n1, e2, n2);
79  else
80  return st->factor * hypot(e1 - e2, n1 - n2);
81 }
82 
91 double G_distance_between_line_segments(double ax1, double ay1,
92  double ax2, double ay2,
93  double bx1, double by1,
94  double bx2, double by2)
95 {
96  double ra, rb;
97  double x, y;
98 
99  /* if the segments intersect, then the distance is zero */
100  if (G_intersect_line_segments(ax1, ay1, ax2, ay2,
101  bx1, by1, bx2, by2, &ra, &rb, &x, &y) > 0)
102  return 0.0;
103  return
104  min4(G_distance_point_to_line_segment(ax1, ay1, bx1, by1, bx2, by2),
105  G_distance_point_to_line_segment(ax2, ay2, bx1, by1, bx2, by2),
106  G_distance_point_to_line_segment(bx1, by1, ax1, ay1, ax2, ay2),
107  G_distance_point_to_line_segment(bx2, by2, ax1, ay1, ax2, ay2)
108  );
109 }
110 
120 double G_distance_point_to_line_segment(double xp, double yp,
121  double x1, double y1, double x2,
122  double y2)
123 {
124  double dx, dy;
125  double x, y;
126  double xq, yq, ra, rb;
127  int t;
128 
129  /* define the perpendicular to the segment thru the point */
130  dx = x1 - x2;
131  dy = y1 - y2;
132 
133  if (dx == 0.0 && dy == 0.0)
134  return G_distance(x1, y1, xp, yp);
135 
136  if (fabs(dy) > fabs(dx)) {
137  xq = xp + dy;
138  yq = (dx / dy) * (xp - xq) + yp;
139  }
140  else {
141  yq = yp + dx;
142  xq = (dy / dx) * (yp - yq) + xp;
143  }
144 
145  /* find the intersection of the perpendicular with the segment */
146  t = G_intersect_line_segments(xp, yp, xq, yq, x1, y1, x2, y2, &ra,
147  &rb, &x, &y);
148  switch (t) {
149  case 0:
150  case 1:
151  break;
152  default:
153  /* parallel/colinear cases shouldn't occur with perpendicular lines */
154  G_warning(_("%s: shouldn't happen: "
155  "code=%d P=(%f,%f) S=(%f,%f)(%f,%f)"),
156  "G_distance_point_to_line_segment", t, xp, yp, x1, y1, x2, y2);
157  return -1.0;
158  }
159 
160  /* if x,y lies on the segment, then the distance is from to x,y */
161  if (rb >= 0 && rb <= 1.0)
162  return G_distance(x, y, xp, yp);
163 
164  /* otherwise the distance is the short of the distances to the endpoints
165  * of the segment
166  */
167  return min2(G_distance(x1, y1, xp, yp), G_distance(x2, y2, xp, yp));
168 }
169 
170 static double min4(double a, double b, double c, double d)
171 {
172  return min2(min2(a, b), min2(c, d));
173 }
174 
175 static double min2(double a, double b)
176 {
177  return a < b ? a : b;
178 }
void G_begin_geodesic_distance(double a, double e2)
Begin geodesic distance.
Definition: geodist.c:50
int G_begin_distance_calculations(void)
Begin distance calculations.
Definition: gis/distance.c:42
int G_get_ellipsoid_parameters(double *a, double *e2)
get ellipsoid parameters
Definition: get_ellipse.c:67
double b
struct state * st
Definition: parser.c:101
double t
double G_distance(double e1, double n1, double e2, double n2)
Returns distance in meters.
Definition: gis/distance.c:75
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:91
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition: geodist.c:196
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:120
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition: proj3.c:124
struct state state
Definition: parser.c:100
int G_projection(void)
Query cartographic projection.
Definition: proj1.c:32
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
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: intersect.c:80