GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-381f661b3f
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:
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, double ax2,
92  double ay2, double bx1, double by1,
93  double bx2, double by2)
94 {
95  double ra, rb;
96  double x, y;
97 
98  /* if the segments intersect, then the distance is zero */
99  if (G_intersect_line_segments(ax1, ay1, ax2, ay2, bx1, by1, bx2, by2, &ra,
100  &rb, &x, &y) > 0)
101  return 0.0;
102  return min4(G_distance_point_to_line_segment(ax1, ay1, bx1, by1, bx2, by2),
103  G_distance_point_to_line_segment(ax2, ay2, bx1, by1, bx2, by2),
104  G_distance_point_to_line_segment(bx1, by1, ax1, ay1, ax2, ay2),
105  G_distance_point_to_line_segment(bx2, by2, ax1, ay1, ax2, ay2));
106 }
107 
108 /*!
109  \brief Returns distance between a point and line segment in meters.
110 
111  \param xp,yp point coordinates
112  \param x1,y1 segment point coordinates
113  \param x2,y2 segment point coordinates
114 
115  \return distance
116  */
117 double G_distance_point_to_line_segment(double xp, double yp, double x1,
118  double y1, double x2, double y2)
119 {
120  double dx, dy;
121  double x, y;
122  double xq, yq, ra, rb;
123  int t;
124 
125  /* define the perpendicular to the segment through the point */
126  dx = x1 - x2;
127  dy = y1 - y2;
128 
129  if (dx == 0.0 && dy == 0.0)
130  return G_distance(x1, y1, xp, yp);
131 
132  if (fabs(dy) > fabs(dx)) {
133  xq = xp + dy;
134  yq = (dx / dy) * (xp - xq) + yp;
135  }
136  else {
137  yq = yp + dx;
138  xq = (dy / dx) * (yp - yq) + xp;
139  }
140 
141  /* find the intersection of the perpendicular with the segment */
142  t = G_intersect_line_segments(xp, yp, xq, yq, x1, y1, x2, y2, &ra, &rb, &x,
143  &y);
144  switch (t) {
145  case 0:
146  case 1:
147  break;
148  default:
149  /* parallel/colinear cases shouldn't occur with perpendicular lines */
150  G_warning(_("%s: shouldn't happen: "
151  "code=%d P=(%f,%f) S=(%f,%f)(%f,%f)"),
152  "G_distance_point_to_line_segment", t, xp, yp, x1, y1, x2,
153  y2);
154  return -1.0;
155  }
156 
157  /* if x,y lies on the segment, then the distance is from to x,y */
158  if (rb >= 0 && rb <= 1.0)
159  return G_distance(x, y, xp, yp);
160 
161  /* otherwise the distance is the short of the distances to the endpoints
162  * of the segment
163  */
164  return min2(G_distance(x1, y1, xp, yp), G_distance(x2, y2, xp, yp));
165 }
166 
167 static double min4(double a, double b, double c, double d)
168 {
169  return min2(min2(a, b), min2(c, d));
170 }
171 
172 static double min2(double a, double b)
173 {
174  return a < b ? a : b;
175 }
void G_begin_geodesic_distance(double, double)
Begin geodesic distance.
Definition: geodist.c:50
void G_warning(const char *,...) __attribute__((format(printf
int G_get_ellipsoid_parameters(double *, double *)
get ellipsoid parameters
Definition: get_ellipse.c:66
double G_database_units_to_meters_factor(void)
Conversion to meters.
Definition: proj3.c:146
double G_geodesic_distance(double, double, double, double)
Calculates geodesic distance.
Definition: geodist.c:195
int G_intersect_line_segments(double, double, double, double, double, double, double, double, double *, double *, double *, double *)
Definition: gis/intersect.c:84
int G_projection(void)
Query cartographic projection.
Definition: proj1.c:32
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:117
int G_begin_distance_calculations(void)
Begin distance calculations.
Definition: gis/distance.c:42
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_distance(double e1, double n1, double e2, double n2)
Returns distance in meters.
Definition: gis/distance.c:75
#define PROJECTION_LL
Projection code - Latitude-Longitude.
Definition: gis.h:130
#define _(str)
Definition: glocale.h:10
struct state state
Definition: parser.c:103
struct state * st
Definition: parser.c:104
double b
Definition: r_raster.c:39
double t
Definition: r_raster.c:39
#define x