GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72729
gis/distance.c
Go to the documentation of this file.
1 /*!
2  \file lib/gis/distance.c
3
4  \brief GIS Library - Distance calculation functions.
5
6  WARNING: this code is preliminary and may be changed,
7  including calling sequences to any of the functions
8  defined here.
9
10  (C) 2001-2009, 2011 by the GRASS Development Team
11
12  This program is free software under the GNU General Public License
13  (>=v2). Read the file COPYING that comes with GRASS for details.
14
15  \author Original author CERL
16 */
17
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
32 /*!
33  \brief Begin distance calculations.
34
35  Initializes the distance calculations. It is used both for the
36  planimetric and latitude-longitude projections.
37
38  \return 0 if projection has no metrix (ie. imagery)
39  \return 1 if projection is planimetric
40  \return 2 if projection is latitude-longitude
41 */
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
62 /*!
63  \brief Returns distance in meters.
64
65  This routine computes the distance, in meters, from
66  <i>x1</i>,<i>y1</i> to <i>x2</i>,<i>y2</i>. If the projection is
67  latitude-longitude, this distance is measured along the
68  geodesic. Two routines perform geodesic distance calculations.
69
70  \param e1,n1 east-north coordinates of first point
71  \param e2,n2 east-north coordinates of second point
72
73  \return distance
74 */
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
83 /*!
84  \brief Returns distance between two line segments in meters.
85
86  \param ax1,ay1,ax2,ay2 first segment
87  \param bx1,by1,bx2,by2 second segment
88
89  \return distance value
90 */
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
111 /*!
112  \brief Returns distance between a point and line segment in meters.
113
114  \param xp,yp point coordinates
115  \param x1,y1 segment point coordinates
116  \param x2,y2 segment point coordinates
117
118  \return distance
119 */
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 through 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 }
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
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
#define x
struct state * st
Definition: parser.c:103
double t
Definition: r_raster.c:39
double b
Definition: r_raster.c:39
#define PROJECTION_LL
Projection code - Latitude-Longitude.
Definition: gis.h:99
double G_distance(double e1, double n1, double e2, double n2)
Returns distance in meters.
Definition: gis/distance.c:75
#define _(str)
Definition: glocale.h:13
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:135
struct state state
Definition: parser.c:102
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:204