GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-2653d7b383
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  \brief Change v2 so that v1v2 is a unit vector
314
315  \param v1 first vector
316  \param v2[in,out] second vector
317
318  \return 0 if magnitude of dx is zero
319  \return 1 if magnitude of dx > 0
320  */
321 int GS_v3normalize(float *v1, float *v2)
322 {
323  float n, dx, dy, dz;
324
325  dx = v2[X] - v1[X];
326  dy = v2[Y] - v1[Y];
327  dz = v2[Z] - v1[Z];
328  n = sqrt(dx * dx + dy * dy + dz * dz);
329
330  if (n == 0.0) {
331  return (0);
332  }
333
334  v2[X] = v1[X] + dx / n;
335  v2[Y] = v1[Y] + dy / n;
336  v2[Z] = v1[Z] + dz / n;
337
338  return (1);
339 }
340
341 /*!
342  \brief Get a normalized direction from v1 to v2, store in v3
343
344  \param v1 first vector
345  \param v2 second vector
346  \param[out] v3 output vector
347
348  \return 0 if magnitude of dx is zero
349  \return 1 if magnitude of dx > 0
350  */
351 int GS_v3dir(float *v1, float *v2, float *v3)
352 {
353  float n, dx, dy, dz;
354
355  dx = v2[X] - v1[X];
356  dy = v2[Y] - v1[Y];
357  dz = v2[Z] - v1[Z];
358  n = sqrt(dx * dx + dy * dy + dz * dz);
359
360  if (n == 0.0) {
361  v3[X] = v3[Y] = v3[Z] = 0.0;
362  return (0);
363  }
364
365  v3[X] = dx / n;
366  v3[Y] = dy / n;
367  v3[Z] = dz / n;
368
369  return (1);
370 }
371
372 /*!
373  \brief Get a normalized direction from v1 to v2, store in v3 (2D)
374
375  \param v1 first vector
376  \param v2 second vector
377  \param[out] v3 output vector
378
379  \return 0 if magnitude of dx is zero
380  \return 1 if magnitude of dx > 0
381  */
382 void GS_v2dir(float *v1, float *v2, float *v3)
383 {
384  float n, dx, dy;
385
386  dx = v2[X] - v1[X];
387  dy = v2[Y] - v1[Y];
388  n = sqrt(dx * dx + dy * dy);
389
390  v3[X] = dx / n;
391  v3[Y] = dy / n;
392
393  return;
394 }
395
396 /*!
397  \brief Get the cross product v3 = v1 cross v2
398
399  \param v1 first vector
400  \param v2 second vector
401  \param[out] v3 output vector
402  */
403 void GS_v3cross(float *v1, float *v2, float *v3)
404 {
405  v3[X] = (v1[Y] * v2[Z]) - (v1[Z] * v2[Y]);
406  v3[Y] = (v1[Z] * v2[X]) - (v1[X] * v2[Z]);
407  v3[Z] = (v1[X] * v2[Y]) - (v1[Y] * v2[X]);
408
409  return;
410 }
411
412 /*!
413  \brief Magnitude of vector
414
415  \param v1 vector
416  \param[out] mag magnitude value
417  */
418 void GS_v3mag(float *v1, float *mag)
419 {
420  *mag = sqrt(v1[X] * v1[X] + v1[Y] * v1[Y] + v1[Z] * v1[Z]);
421
422  return;
423 }
424
425 /*!
427
428  Initialize by calling with a number nhist to represent number of
429  previous entries to check, then call with zero as nhist
430
431  \param p1 first point
432  \param p2 second point
433  \param nhist ?
434
435  \return -1 on error
436  \return -2
437  \return 1
438  \return 9
439  */
440 int GS_coordpair_repeats(float *p1, float *p2, int nhist)
441 {
442  static float *entrys = NULL;
443  static int next = 0;
444  static int len = 0;
445  int i;
446
447  if (nhist) {
448  if (entrys) {
449  G_free(entrys);
450  }
451
452  entrys = (float *)G_malloc(4 * nhist * sizeof(float));
453
454  if (!entrys)
455  return (-1);
456
457  len = nhist;
458  next = 0;
459  }
460
461  if (!len) {
462  return (-2);
463  }
464
465  for (i = 0; i < next; i += 4) {
466  if (entrys[i] == p1[0] && entrys[i + 1] == p1[1] &&
467  entrys[i + 2] == p2[0] && entrys[i + 3] == p2[1]) {
468  return (1);
469  }
470  }
471
472  if (len == next / 4) {
473  next = 0;
474  }
475
476  entrys[next] = p1[0];
477  entrys[next + 1] = p1[1];
478  entrys[next + 2] = p2[0];
479  entrys[next + 3] = p2[1];
480  next += 4;
481
482  return (0);
483 }
#define NULL
Definition: ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_malloc(n)
Definition: defs/gis.h:94
double Gs_distance(double *, double *)
Calculates distance in METERS between two points in current projection (2D)
Definition: gs3.c:79
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:382
void GS_v3sub(float *v1, float *v2)
Subtract vectors.
Definition: gs_util.c:212
void GS_v3mult(float *v1, float k)
Multiple vectors.
Definition: gs_util.c:229
void GS_v3mag(float *v1, float *mag)
Magnitude of vector.
Definition: gs_util.c:418
int GS_dv3norm(double *dv1)
Changes v1 so that it is a unit vector.
Definition: gs_util.c:295
int GS_coordpair_repeats(float *p1, float *p2, int nhist)
Definition: gs_util.c:440
float GS_distance(float *from, float *to)
Calculate distance.
Definition: gs_util.c:141
int GS_v2norm(float *v1)
Change v1 so that it is a unit vector (3D)
Definition: gs_util.c:271
void GS_v3add(float *v1, float *v2)
Sum vectors.
Definition: gs_util.c:195
int GS_v3norm(float *v1)
Change v1 so that it is a unit vector (2D)
Definition: gs_util.c:246
void GS_v3eq(float *v1, float *v2)
Copy vector values.
Definition: gs_util.c:178
double GS_geodistance(double *from, double *to, const char *units)
Calculate distance between 2 coordinates.
Definition: gs_util.c:55
float GS_P2distance(float *from, float *to)
Calculate distance in plane.
Definition: gs_util.c:160
int GS_v3dir(float *v1, float *v2, float *v3)
Get a normalized direction from v1 to v2, store in v3.
Definition: gs_util.c:351
int GS_v3normalize(float *v1, float *v2)
Change v2 so that v1v2 is a unit vector.
Definition: gs_util.c:321
void GS_v3cross(float *v1, float *v2, float *v3)
Get the cross product v3 = v1 cross v2.
Definition: gs_util.c:403
#define X
Definition: ogsf.h:140
#define Z
Definition: ogsf.h:142
#define Y
Definition: ogsf.h:141
#define x