GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72846
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gs_util.c
Go to the documentation of this file.
1 /*!
2  \file lib/ogsf/gs_util.c
3 
4  \brief OGSF library - loading and manipulating surfaces
5 
6  GRASS OpenGL gsurf OGSF Library
7 
8  (C) 1999-2008 by the GRASS Development Team
9 
10  This program is free software under the
11  GNU General Public License (>=v2).
12  Read the file COPYING that comes with GRASS
13  for details.
14 
15  \author Bill Brown USACERL, GMSL/University of Illinois
16  \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
17  */
18 
19 #include <stdlib.h>
20 #include <math.h>
21 #include <string.h>
22 
23 #include <grass/gis.h>
24 #include <grass/ogsf.h>
25 
26 /*!
27  \brief Calculate distance between 2 coordinates
28 
29  Units is one of:
30  - "meters",
31  - "miles",
32  - "kilometers",
33  - "feet",
34  - "yards",
35  - "nmiles" (nautical miles),
36  - "rods",
37  - "inches",
38  - "centimeters",
39  - "millimeters",
40  - "micron",
41  - "nanometers",
42  - "cubits",
43  - "hands",
44  - "furlongs",
45  - "chains"
46 
47  Default is meters.
48 
49  \param from starting point
50  \param to ending point
51  \param units map units
52 
53  \return distance between two geographic coordinates in current projection
54  */
55 double GS_geodistance(double *from, double *to, const char *units)
56 {
57  double meters;
58 
59  meters = Gs_distance(from, to);
60 
61  if (!units) {
62  return (meters);
63  }
64 
65  if (strcmp(units, "meters") == 0) {
66  return (meters);
67  }
68 
69  if (strcmp(units, "miles") == 0) {
70  return (meters * .0006213712);
71  }
72 
73  if (strcmp(units, "kilometers") == 0) {
74  return (meters * .001);
75  }
76 
77  if (strcmp(units, "feet") == 0) {
78  return (meters * 3.280840);
79  }
80 
81  if (strcmp(units, "yards") == 0) {
82  return (meters * 1.093613);
83  }
84 
85  if (strcmp(units, "rods") == 0) {
86  return (meters * .1988388);
87  }
88 
89  if (strcmp(units, "inches") == 0) {
90  return (meters * 39.37008);
91  }
92 
93  if (strcmp(units, "centimeters") == 0) {
94  return (meters * 100.0);
95  }
96 
97  if (strcmp(units, "millimeters") == 0) {
98  return (meters * 1000.0);
99  }
100 
101  if (strcmp(units, "micron") == 0) {
102  return (meters * 1000000.0);
103  }
104 
105  if (strcmp(units, "nanometers") == 0) {
106  return (meters * 1000000000.0);
107  }
108 
109  if (strcmp(units, "cubits") == 0) {
110  return (meters * 2.187227);
111  }
112 
113  if (strcmp(units, "hands") == 0) {
114  return (meters * 9.842520);
115  }
116 
117  if (strcmp(units, "furlongs") == 0) {
118  return (meters * .004970970);
119  }
120 
121  if (strcmp(units, "nmiles") == 0) {
122  /* nautical miles */
123  return (meters * .0005399568);
124  }
125 
126  if (strcmp(units, "chains") == 0) {
127  return (meters * .0497097);
128  }
129 
130  return (meters);
131 }
132 
133 /*!
134  \brief Calculate distance
135 
136  \param from 'from' point (X,Y,Z)
137  \param to 'to' point (X,Y,Z)
138 
139  \return distance
140  */
141 float GS_distance(float *from, float *to)
142 {
143  float x, y, z;
144 
145  x = from[X] - to[X];
146  y = from[Y] - to[Y];
147  z = from[Z] - to[Z];
148 
149  return (float)sqrt(x * x + y * y + z * z);
150 }
151 
152 /*!
153  \brief Calculate distance in plane
154 
155  \param from 'from' point (X,Y)
156  \param to 'to' point (X,Y)
157 
158  \return distance
159  */
160 float GS_P2distance(float *from, float *to)
161 {
162  float x, y;
163 
164  x = from[X] - to[X];
165  y = from[Y] - to[Y];
166 
167  return (float)sqrt(x * x + y * y);
168 }
169 
170 /*!
171  \brief Copy vector values
172 
173  v1 = v2
174 
175  \param[out] v1 first vector
176  \param v2 second vector
177  */
178 void GS_v3eq(float *v1, float *v2)
179 {
180  v1[X] = v2[X];
181  v1[Y] = v2[Y];
182  v1[Z] = v2[Z];
183 
184  return;
185 }
186 
187 /*!
188  \brief Sum vectors
189 
190  v1 += v2
191 
192  \param[in,out] v1 first vector
193  \param v2 second vector
194  */
195 void GS_v3add(float *v1, float *v2)
196 {
197  v1[X] += v2[X];
198  v1[Y] += v2[Y];
199  v1[Z] += v2[Z];
200 
201  return;
202 }
203 
204 /*!
205  \brief Subtract vectors
206 
207  v1 -= v2
208 
209  \param[in,out] v1 first vector
210  \param v2 second vector
211  */
212 void GS_v3sub(float *v1, float *v2)
213 {
214  v1[X] -= v2[X];
215  v1[Y] -= v2[Y];
216  v1[Z] -= v2[Z];
217 
218  return;
219 }
220 
221 /*!
222  \brief Multiple vectors
223 
224  v1 *= k
225 
226  \param[in,out] v1 vector
227  \param k multiplicator
228  */
229 void GS_v3mult(float *v1, float k)
230 {
231  v1[X] *= k;
232  v1[Y] *= k;
233  v1[Z] *= k;
234 
235  return;
236 }
237 
238 /*!
239  \brief Change v1 so that it is a unit vector (2D)
240 
241  \param[in,out] v1 vector
242 
243  \return 0 if magnitude of v1 is zero
244  \return 1 if magnitude of v1 > 0
245  */
246 int GS_v3norm(float *v1)
247 {
248  float n;
249 
250  n = sqrt(v1[X] * v1[X] + v1[Y] * v1[Y] + v1[Z] * v1[Z]);
251 
252  if (n == 0.0) {
253  return (0);
254  }
255 
256  v1[X] /= n;
257  v1[Y] /= n;
258  v1[Z] /= n;
259 
260  return (1);
261 }
262 
263 /*!
264  \brief Change v1 so that it is a unit vector (3D)
265 
266  \param[in,out] v1 vector
267 
268  \return 0 if magnitude of v1 is zero
269  \return 1 if magnitude of v1 > 0
270  */
271 int GS_v2norm(float *v1)
272 {
273  float n;
274 
275  n = sqrt(v1[X] * v1[X] + v1[Y] * v1[Y]);
276 
277  if (n == 0.0) {
278  return (0);
279  }
280 
281  v1[X] /= n;
282  v1[Y] /= n;
283 
284  return (1);
285 }
286 
287 /*!
288  \brief Changes v1 so that it is a unit vector
289 
290  \param dv1 vector
291 
292  \return 0 if magnitude of dv1 is zero
293  \return 1 if magnitude of dv1 > 0
294  */
295 int GS_dv3norm(double *dv1)
296 {
297  double n;
298 
299  n = sqrt(dv1[X] * dv1[X] + dv1[Y] * dv1[Y] + dv1[Z] * dv1[Z]);
300 
301  if (n == 0.0) {
302  return (0);
303  }
304 
305  dv1[X] /= n;
306  dv1[Y] /= n;
307  dv1[Z] /= n;
308 
309  return (1);
310 }
311 
312 
313 /*!
314  \brief Change v2 so that v1v2 is a unit vector
315 
316  \param v1 first vector
317  \param v2[in,out] second vector
318 
319  \return 0 if magnitude of dx is zero
320  \return 1 if magnitude of dx > 0
321  */
322 int GS_v3normalize(float *v1, float *v2)
323 {
324  float n, dx, dy, dz;
325 
326  dx = v2[X] - v1[X];
327  dy = v2[Y] - v1[Y];
328  dz = v2[Z] - v1[Z];
329  n = sqrt(dx * dx + dy * dy + dz * dz);
330 
331  if (n == 0.0) {
332  return (0);
333  }
334 
335  v2[X] = v1[X] + dx / n;
336  v2[Y] = v1[Y] + dy / n;
337  v2[Z] = v1[Z] + dz / n;
338 
339  return (1);
340 }
341 
342 
343 /*!
344  \brief Get a normalized direction from v1 to v2, store in v3
345 
346  \param v1 first vector
347  \param v2 second vector
348  \param[out] v3 output vector
349 
350  \return 0 if magnitude of dx is zero
351  \return 1 if magnitude of dx > 0
352  */
353 int GS_v3dir(float *v1, float *v2, float *v3)
354 {
355  float n, dx, dy, dz;
356 
357  dx = v2[X] - v1[X];
358  dy = v2[Y] - v1[Y];
359  dz = v2[Z] - v1[Z];
360  n = sqrt(dx * dx + dy * dy + dz * dz);
361 
362  if (n == 0.0) {
363  v3[X] = v3[Y] = v3[Z] = 0.0;
364  return (0);
365  }
366 
367  v3[X] = dx / n;
368  v3[Y] = dy / n;
369  v3[Z] = dz / n;
370 
371  return (1);
372 }
373 
374 
375 /*!
376  \brief Get a normalized direction from v1 to v2, store in v3 (2D)
377 
378  \param v1 first vector
379  \param v2 second vector
380  \param[out] v3 output vector
381 
382  \return 0 if magnitude of dx is zero
383  \return 1 if magnitude of dx > 0
384  */
385 void GS_v2dir(float *v1, float *v2, float *v3)
386 {
387  float n, dx, dy;
388 
389  dx = v2[X] - v1[X];
390  dy = v2[Y] - v1[Y];
391  n = sqrt(dx * dx + dy * dy);
392 
393  v3[X] = dx / n;
394  v3[Y] = dy / n;
395 
396  return;
397 }
398 
399 /*!
400  \brief Get the cross product v3 = v1 cross v2
401 
402  \param v1 first vector
403  \param v2 second vector
404  \param[out] v3 output vector
405  */
406 void GS_v3cross(float *v1, float *v2, float *v3)
407 {
408  v3[X] = (v1[Y] * v2[Z]) - (v1[Z] * v2[Y]);
409  v3[Y] = (v1[Z] * v2[X]) - (v1[X] * v2[Z]);
410  v3[Z] = (v1[X] * v2[Y]) - (v1[Y] * v2[X]);
411 
412  return;
413 }
414 
415 /*!
416  \brief Magnitude of vector
417 
418  \param v1 vector
419  \param[out] mag magnitude value
420  */
421 void GS_v3mag(float *v1, float *mag)
422 {
423  *mag = sqrt(v1[X] * v1[X] + v1[Y] * v1[Y] + v1[Z] * v1[Z]);
424 
425  return;
426 }
427 
428 /*!
429  \brief ADD
430 
431  Initialize by calling with a number nhist to represent number of
432  previous entrys to check, then call with zero as nhist
433 
434  \param p1 first point
435  \param p2 second point
436  \param nhist ?
437 
438  \return -1 on error
439  \return -2
440  \return 1
441  \return 9
442  */
443 int GS_coordpair_repeats(float *p1, float *p2, int nhist)
444 {
445  static float *entrys = NULL;
446  static int next = 0;
447  static int len = 0;
448  int i;
449 
450  if (nhist) {
451  if (entrys) {
452  G_free(entrys);
453  }
454 
455  entrys = (float *)G_malloc(4 * nhist * sizeof(float));
456 
457  if (!entrys)
458  return (-1);
459 
460  len = nhist;
461  next = 0;
462  }
463 
464  if (!len) {
465  return (-2);
466  }
467 
468  for (i = 0; i < next; i += 4) {
469  if (entrys[i] == p1[0] && entrys[i + 1] == p1[1]
470  && entrys[i + 2] == p2[0] && entrys[i + 3] == p2[1]) {
471  return (1);
472  }
473  }
474 
475  if (len == next / 4) {
476  next = 0;
477  }
478 
479  entrys[next] = p1[0];
480  entrys[next + 1] = p1[1];
481  entrys[next + 2] = p2[0];
482  entrys[next + 3] = p2[1];
483  next += 4;
484 
485  return (0);
486 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
int GS_v3normalize(float *v1, float *v2)
Change v2 so that v1v2 is a unit vector.
Definition: gs_util.c:322
void GS_v3sub(float *v1, float *v2)
Subtract vectors.
Definition: gs_util.c:212
int GS_v3dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3.
Definition: gs_util.c:353
int GS_v3norm(float *v1)
Change v1 so that it is a unit vector (2D)
Definition: gs_util.c:246
void GS_v3add(float *v1, float *v2)
Sum vectors.
Definition: gs_util.c:195
#define NULL
Definition: ccmath.h:32
#define x
void GS_v3cross(float *v1, float *v2, float *v3)
Get the cross product v3 = v1 cross v2.
Definition: gs_util.c:406
int GS_v2norm(float *v1)
Change v1 so that it is a unit vector (3D)
Definition: gs_util.c:271
void GS_v2dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3 (2D)
Definition: gs_util.c:385
double GS_geodistance(double *from, double *to, const char *units)
Calculate distance between 2 coordinates.
Definition: gs_util.c:55
double Gs_distance(double *from, double *to)
Calculates distance in METERS between two points in current projection (2D)
Definition: gs3.c:84
void GS_v3mult(float *v1, float k)
Multiple vectors.
Definition: gs_util.c:229
#define Z
Definition: ogsf.h:139
#define Y
Definition: ogsf.h:138
int GS_dv3norm(double *dv1)
Changes v1 so that it is a unit vector.
Definition: gs_util.c:295
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition: gs_util.c:178
#define X
Definition: ogsf.h:137
float GS_distance(float *from, float *to)
Calculate distance.
Definition: gs_util.c:141
int GS_coordpair_repeats(float *p1, float *p2, int nhist)
ADD.
Definition: gs_util.c:443
void GS_v3mag(float *v1, float *mag)
Magnitude of vector.
Definition: gs_util.c:421
float GS_P2distance(float *from, float *to)
Calculate distance in plane.
Definition: gs_util.c:160