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 */
GS_geodistance(double * from,double * to,const char * units)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 */
GS_distance(float * from,float * to)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 */
GS_P2distance(float * from,float * to)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 */
GS_v3eq(float * v1,float * v2)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 */
GS_v3add(float * v1,float * v2)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 */
GS_v3sub(float * v1,float * v2)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 */
GS_v3mult(float * v1,float k)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 */
GS_v3norm(float * v1)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 */
GS_v2norm(float * v1)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 */
GS_dv3norm(double * dv1)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 */
GS_v3normalize(float * v1,float * v2)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 */
GS_v3dir(float * v1,float * v2,float * v3)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 */
GS_v2dir(float * v1,float * v2,float * v3)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 */
GS_v3cross(float * v1,float * v2,float * v3)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 */
GS_v3mag(float * v1,float * mag)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 */
GS_coordpair_repeats(float * p1,float * p2,int nhist)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 }
487