1 /*!
2 \file lib/vector/Vlib/line.c
3
4 \brief Vector library - vector feature geometry
5
6 (C) 2001-2009 by the GRASS Development Team
7
8 This program is free software under the GNU General Public License
9 (>=v2). Read the file COPYING that comes with GRASS for details.
10
11 \author Original author CERL, probably Dave Gerdes or Mike Higgins.
12 \author Update to GRASS 5.7 Radim Blazek and David D. Gray.
13 \author Some updates for GRASS 7 by Martin Landa <landa.martin gmail.com>
14 */
15
16 #include <stdlib.h>
17 #include <math.h>
18 #include <grass/vector.h>
19 #include <grass/glocale.h>
20
21 /*!
22 \brief Creates and initializes a struct line_pnts (internal use only)
23
24 Use Vect_new_line_struct() instead.
25
26 \return pointer to allocated line_pnts structure
27 \return NULL on error
28 */
29 struct line_pnts *Vect__new_line_struct(void);
30
31 /*!
32 \brief Creates and initializes a line_pnts structure
33
34 This structure is used for reading and writing vector lines and
35 polygons. The library routines handle all memory allocation. If 3
36 lines in memory are needed at the same time, then simply 3 line_pnts
37 structures have to be used.
38
39 To free allocated memory call Vect_destroy_line_struct().
40
41 Calls G_fatal_error() on error.
42
43 \return pointer to line_pnts
44 */
Vect_new_line_struct()45 struct line_pnts *Vect_new_line_struct()
46 {
47 struct line_pnts *p;
48
49 if (NULL == (p = Vect__new_line_struct()))
50 G_fatal_error("Vect_new_line_struct(): %s", _("Out of memory"));
51
52 return p;
53 }
54
Vect__new_line_struct()55 struct line_pnts *Vect__new_line_struct()
56 {
57 struct line_pnts *p;
58
59 p = (struct line_pnts *)malloc(sizeof(struct line_pnts));
60
61 /* alloc_points MUST be initialized to zero */
62 if (p)
63 p->alloc_points = p->n_points = 0;
64
65 if (p)
66 p->x = p->y = p->z = NULL;
67
68 return p;
69 }
70
71 /*!
72 \brief Frees all memory associated with a line_pnts structure,
73 including the structure itself
74
75 \param p pointer to line_pnts structure
76 */
Vect_destroy_line_struct(struct line_pnts * p)77 void Vect_destroy_line_struct(struct line_pnts *p)
78 {
79 if (p) { /* probably a moot test */
80 if (p->alloc_points) {
81 G_free((char *)p->x);
82 G_free((char *)p->y);
83 G_free((char *)p->z);
84 }
85 G_free((char *)p);
86 }
87 }
88
89 /*!
90 \brief Copy points from array to line_pnts structure
91
92 \param Points pointer to line_ptns structure
93 \param x,y,z array of coordinates
94 \param n number of points to be copied
95
96 \return 0 on success
97 \return -1 on out of memory
98 */
Vect_copy_xyz_to_pnts(struct line_pnts * Points,const double * x,const double * y,const double * z,int n)99 int Vect_copy_xyz_to_pnts(struct line_pnts *Points, const double *x, const double *y,
100 const double *z, int n)
101 {
102 int i;
103
104 if (0 > dig_alloc_points(Points, n))
105 return -1;
106
107 for (i = 0; i < n; i++) {
108 Points->x[i] = x[i];
109 Points->y[i] = y[i];
110 if (z != NULL)
111 Points->z[i] = z[i];
112 else
113 Points->z[i] = 0;
114 Points->n_points = n;
115 }
116
117 return 0;
118 }
119
120
121 /*!
122 \brief Reset line
123
124 Make sure line structure is clean to be re-used, i.e. it has no
125 points associated with it Points must have previously been created
126 with Vect_new_line_struct().
127
128 \param Points pointer to line_pnts structure to be reset
129 */
Vect_reset_line(struct line_pnts * Points)130 void Vect_reset_line(struct line_pnts *Points)
131 {
132 Points->n_points = 0;
133 }
134
135 /*!
136 \brief Appends one point to the end of a line.
137
138 If you are re-using a line struct, be sure to clear out old data
139 first by calling Vect_reset_line().
140
141 Calls G_fatal_error() when out of memory.
142
143 \param Points pointer to line_pnts structure
144 \param x,y,z point coordinates to be added
145
146 \return number of points
147 \return -1 on error (out of memory)
148 */
Vect_append_point(struct line_pnts * Points,double x,double y,double z)149 int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
150 {
151 register int n;
152
153 if (0 > dig_alloc_points(Points, Points->n_points + 1)) {
154 G_fatal_error(_("Out of memory"));
155 return -1;
156 }
157
158 n = Points->n_points;
159 Points->x[n] = x;
160 Points->y[n] = y;
161 Points->z[n] = z;
162
163 return ++(Points->n_points);
164 }
165
166 /*!
167 \brief Insert new point at index position and move all old points
168 at that position and above up
169
170 \param Points pointer to line_pnts structure
171 \param index (from 0 to Points->n_points-1)
172 \param x,y,z point coordinates
173
174 \return number of points
175 \return -1 on error (alocation)
176 */
Vect_line_insert_point(struct line_pnts * Points,int index,double x,double y,double z)177 int Vect_line_insert_point(struct line_pnts *Points, int index, double x,
178 double y, double z)
179 {
180 int n;
181
182 if (index < 0 || index > Points->n_points - 1)
183 G_fatal_error("Vect_line_insert_point(): %s",
184 _("Index out of range in"));
185
186 if (0 > dig_alloc_points(Points, Points->n_points + 1))
187 return -1;
188
189 /* move up */
190 for (n = Points->n_points; n > index; n--) {
191 Points->x[n] = Points->x[n - 1];
192 Points->y[n] = Points->y[n - 1];
193 Points->z[n] = Points->z[n - 1];
194 }
195
196 Points->x[index] = x;
197 Points->y[index] = y;
198 Points->z[index] = z;
199
200 return ++(Points->n_points);
201 }
202
203 /*!
204 \brief Delete point at given index and move all points above down
205
206 \param Points pointer to line_pnts structure
207 \param index point (from 0 to Points->n_points-1)
208
209 \return number of points
210 */
Vect_line_delete_point(struct line_pnts * Points,int index)211 int Vect_line_delete_point(struct line_pnts *Points, int index)
212 {
213 int n;
214
215 if (index < 0 || index > Points->n_points - 1)
216 G_fatal_error("Vect_line_insert_point(): %s",
217 _("Index out of range in"));
218
219 if (Points->n_points == 0)
220 return 0;
221
222 /* move down */
223 for (n = index; n < Points->n_points - 1; n++) {
224 Points->x[n] = Points->x[n + 1];
225 Points->y[n] = Points->y[n + 1];
226 Points->z[n] = Points->z[n + 1];
227 }
228
229 return --(Points->n_points);
230 }
231
232 /*!
233 \brief Get line point of given index
234
235 Calls G_fatal_error() when index is not range in.
236
237 \param Points pointer to line_pnts structure
238 \param index point index (from 0 to Points->n_points-1)
239 \param x pointer to x coordinate or NULL
240 \param y pointer to y coordinate or NULL
241 \param z pointer to z coordinate or NULL
242
243 \return number of points
244 */
Vect_line_get_point(const struct line_pnts * Points,int index,double * x,double * y,double * z)245 int Vect_line_get_point(const struct line_pnts *Points, int index,
246 double *x, double *y, double *z)
247 {
248 if (index < 0 || index > Points->n_points - 1)
249 G_fatal_error("Vect_line_get_point(): %s",
250 _("Index out of range in"));
251
252 if (x)
253 *x = Points->x[index];
254 if (y)
255 *y = Points->y[index];
256 if (z)
257 *z = Points->z[index];
258
259 return Points->n_points;
260 }
261
262 /*!
263 \brief Get number of line points
264
265 \param Points pointer to line_pnts structure
266
267 \return number of points
268 */
Vect_get_num_line_points(const struct line_pnts * Points)269 int Vect_get_num_line_points(const struct line_pnts *Points)
270 {
271 return Points->n_points;
272 }
273
274 /*!
275 \brief Remove duplicate points, i.e. zero length segments
276
277 \param Points pointer to line_pnts structure
278
279 \return number of points
280 */
Vect_line_prune(struct line_pnts * Points)281 int Vect_line_prune(struct line_pnts *Points)
282 {
283 int i, j;
284
285 if (Points->n_points > 0) {
286 j = 1;
287 for (i = 1; i < Points->n_points; i++) {
288 if (Points->x[i] != Points->x[j - 1] ||
289 Points->y[i] != Points->y[j - 1]
290 || Points->z[i] != Points->z[j - 1]) {
291 Points->x[j] = Points->x[i];
292 Points->y[j] = Points->y[i];
293 Points->z[j] = Points->z[i];
294 j++;
295 }
296 }
297 Points->n_points = j;
298 }
299
300 return (Points->n_points);
301 }
302
303 /*!
304 \brief Remove points in threshold
305
306 \param Points pointer to line_pnts structure
307 \param threshold threshold value
308
309 \return number of points in result
310 */
Vect_line_prune_thresh(struct line_pnts * Points,double threshold)311 int Vect_line_prune_thresh(struct line_pnts *Points, double threshold)
312 {
313 int ret;
314
315 ret = dig_prune(Points, threshold);
316
317 if (ret < Points->n_points)
318 Points->n_points = ret;
319
320 return (Points->n_points);
321 }
322
323 /*!
324 \brief Appends points to the end of a line.
325
326 Note, this will append to whatever is in line_pnts structure. If you
327 are re-using a line struct, be sure to clear out old data first by
328 calling Vect_reset_line().
329
330 \param Points pointer to line_pnts structure
331 \param APoints points to be included
332 \param direction direction (GV_FORWARD, GV_BACKWARD)
333
334 \return new number of points
335 \return -1 on out of memory
336 */
Vect_append_points(struct line_pnts * Points,const struct line_pnts * APoints,int direction)337 int Vect_append_points(struct line_pnts *Points, const struct line_pnts *APoints,
338 int direction)
339 {
340 int i, n, on, an;
341
342 on = Points->n_points;
343 an = APoints->n_points;
344 n = on + an;
345
346 /* Should be OK, dig_alloc_points calls realloc */
347 if (0 > dig_alloc_points(Points, n))
348 return (-1);
349
350 if (direction == GV_FORWARD) {
351 for (i = 0; i < an; i++) {
352 Points->x[on + i] = APoints->x[i];
353 Points->y[on + i] = APoints->y[i];
354 Points->z[on + i] = APoints->z[i];
355 }
356 }
357 else {
358 for (i = 0; i < an; i++) {
359 Points->x[on + i] = APoints->x[an - i - 1];
360 Points->y[on + i] = APoints->y[an - i - 1];
361 Points->z[on + i] = APoints->z[an - i - 1];
362 }
363 }
364
365 Points->n_points = n;
366 return n;
367 }
368
369
370 /*!
371 \brief Copy points from line structure to array
372
373 x/y/z arrays MUST be at least as large as Points->n_points
374
375 Also note that n is a pointer to int.
376
377 \param Points pointer to line_pnts structure
378 \param x,y,z coordinates arrays
379 \param n number of points
380
381 \return number of points copied
382 */
Vect_copy_pnts_to_xyz(const struct line_pnts * Points,double * x,double * y,double * z,int * n)383 int Vect_copy_pnts_to_xyz(const struct line_pnts *Points, double *x, double *y,
384 double *z, int *n)
385 {
386 int i;
387
388 for (i = 0; i < *n; i++) {
389 x[i] = Points->x[i];
390 y[i] = Points->y[i];
391 if (z != NULL)
392 z[i] = Points->z[i];
393 *n = Points->n_points;
394 }
395
396 return (Points->n_points);
397 }
398
399 /*!
400 \brief Find point on line in the specified distance.
401
402 From the beginning, measured along line.
403
404 If the distance is greater than line length or negative, error is returned.
405
406 \param Points pointer to line_pnts structure
407 \param distance distance value
408 \param x,y,z pointers to point coordinates or NULL
409 \param angle pointer to angle of line in that point (radians, counter clockwise from x axis) or NULL
410 \param slope pointer to slope angle in radians (positive up)
411
412 \return number of segment the point is on (first is 1),
413 \return 0 error when point is outside the line
414 */
Vect_point_on_line(const struct line_pnts * Points,double distance,double * x,double * y,double * z,double * angle,double * slope)415 int Vect_point_on_line(const struct line_pnts *Points, double distance,
416 double *x, double *y, double *z, double *angle,
417 double *slope)
418 {
419 int j, np, seg = 0;
420 double dist = 0, length;
421 double xp = 0, yp = 0, zp = 0, dx = 0, dy = 0, dz = 0, dxy =
422 0, dxyz, k, rest;
423
424 G_debug(3, "Vect_point_on_line(): distance = %f", distance);
425 if ((distance < 0) || (Points->n_points < 2))
426 return 0;
427
428 /* Check if first or last */
429 length = Vect_line_length(Points);
430 G_debug(3, " length = %f", length);
431 if (distance < 0 || distance > length) {
432 G_debug(3, " -> outside line");
433 return 0;
434 }
435
436 np = Points->n_points;
437 if (distance == 0) {
438 G_debug(3, " -> first point");
439 xp = Points->x[0];
440 yp = Points->y[0];
441 zp = Points->z[0];
442 dx = Points->x[1] - Points->x[0];
443 dy = Points->y[1] - Points->y[0];
444 dz = Points->z[1] - Points->z[0];
445 dxy = hypot(dx, dy);
446 seg = 1;
447 }
448 else if (distance == length) {
449 G_debug(3, " -> last point");
450 xp = Points->x[np - 1];
451 yp = Points->y[np - 1];
452 zp = Points->z[np - 1];
453 dx = Points->x[np - 1] - Points->x[np - 2];
454 dy = Points->y[np - 1] - Points->y[np - 2];
455 dz = Points->z[np - 1] - Points->z[np - 2];
456 dxy = hypot(dx, dy);
457 seg = np - 1;
458 }
459 else {
460 for (j = 0; j < Points->n_points - 1; j++) {
461 /* dxyz = G_distance(Points->x[j], Points->y[j],
462 Points->x[j+1], Points->y[j+1]); */
463 dx = Points->x[j + 1] - Points->x[j];
464 dy = Points->y[j + 1] - Points->y[j];
465 dz = Points->z[j + 1] - Points->z[j];
466 dxy = hypot(dx, dy);
467 dxyz = hypot(dxy, dz);
468
469 dist += dxyz;
470 if (dist >= distance) { /* point is on the current line part */
471 rest = distance - dist + dxyz; /* from first point of segment to point */
472 k = rest / dxyz;
473
474 xp = Points->x[j] + k * dx;
475 yp = Points->y[j] + k * dy;
476 zp = Points->z[j] + k * dz;
477 seg = j + 1;
478 break;
479 }
480 }
481 }
482
483 if (x != NULL)
484 *x = xp;
485 if (y != NULL)
486 *y = yp;
487 if (z != NULL)
488 *z = zp;
489
490 /* calculate angle */
491 if (angle != NULL)
492 *angle = atan2(dy, dx);
493
494 /* calculate slope */
495 if (slope != NULL)
496 *slope = atan2(dz, dxy);
497
498 return seg;
499 }
500
501 /*!
502 \brief Create line segment.
503
504 Creates segment of InPoints from start to end measured along the
505 line and write it to OutPoints.
506
507 If the distance is greater than line length or negative, error is
508 returned.
509
510 \param InPoints input line
511 \param start segment number
512 \param end segment number
513 \param OutPoints output line
514
515 \return 1 success
516 \return 0 error when start > length or end < 0 or start < 0 or end > length
517 */
Vect_line_segment(const struct line_pnts * InPoints,double start,double end,struct line_pnts * OutPoints)518 int Vect_line_segment(const struct line_pnts *InPoints, double start, double end,
519 struct line_pnts *OutPoints)
520 {
521 int i, seg1, seg2;
522 double length, tmp;
523 double x1, y1, z1, x2, y2, z2;
524
525 G_debug(3, "Vect_line_segment(): start = %f, end = %f, n_points = %d",
526 start, end, InPoints->n_points);
527
528 Vect_reset_line(OutPoints);
529
530 if (start > end) {
531 tmp = start;
532 start = end;
533 end = tmp;
534 }
535
536 /* Check start/end */
537 if (end < 0)
538 return 0;
539 length = Vect_line_length(InPoints);
540 if (start > length)
541 return 0;
542
543 /* Find coordinates and segments of start/end */
544 seg1 = Vect_point_on_line(InPoints, start, &x1, &y1, &z1, NULL, NULL);
545 seg2 = Vect_point_on_line(InPoints, end, &x2, &y2, &z2, NULL, NULL);
546
547 G_debug(3, " -> seg1 = %d seg2 = %d", seg1, seg2);
548
549 if (seg1 == 0 || seg2 == 0) {
550 G_warning(_("Segment outside line, no segment created"));
551 return 0;
552 }
553
554 Vect_append_point(OutPoints, x1, y1, z1);
555
556 for (i = seg1; i < seg2; i++) {
557 Vect_append_point(OutPoints, InPoints->x[i], InPoints->y[i],
558 InPoints->z[i]);
559 };
560
561 Vect_append_point(OutPoints, x2, y2, z2);
562 Vect_line_prune(OutPoints);
563
564 return 1;
565 }
566
567 /*!
568 \brief Calculate line length, 3D-length in case of 3D vector line
569
570 For Lat-Long use Vect_line_geodesic_length() instead.
571
572 \param Points pointer to line_pnts structure geometry
573
574 \return line length
575 */
Vect_line_length(const struct line_pnts * Points)576 double Vect_line_length(const struct line_pnts *Points)
577 {
578 int j;
579 double dx, dy, dz, len = 0;
580
581 if (Points->n_points < 2)
582 return 0;
583
584 for (j = 0; j < Points->n_points - 1; j++) {
585 dx = Points->x[j + 1] - Points->x[j];
586 dy = Points->y[j + 1] - Points->y[j];
587 dz = Points->z[j + 1] - Points->z[j];
588 len += hypot(hypot(dx, dy), dz);
589 }
590
591 return len;
592 }
593
594
595 /*!
596 \brief Calculate line length.
597
598 If projection is LL, the length is measured along the geodesic.
599
600 \param Points pointer to line_pnts structure geometry
601
602 \return line length
603 */
Vect_line_geodesic_length(const struct line_pnts * Points)604 double Vect_line_geodesic_length(const struct line_pnts *Points)
605 {
606 int j, dc;
607 double dx, dy, dz, dxy, len = 0;
608
609 dc = G_begin_distance_calculations();
610
611 if (Points->n_points < 2)
612 return 0;
613
614 for (j = 0; j < Points->n_points - 1; j++) {
615 if (dc == 2)
616 dxy =
617 G_geodesic_distance(Points->x[j], Points->y[j],
618 Points->x[j + 1], Points->y[j + 1]);
619 else {
620 dx = Points->x[j + 1] - Points->x[j];
621 dy = Points->y[j + 1] - Points->y[j];
622 dxy = hypot(dx, dy);
623 }
624
625 dz = Points->z[j + 1] - Points->z[j];
626 len += hypot(dxy, dz);
627 }
628
629 return len;
630 }
631
632 /*!
633 \brief Calculate distance of point to line.
634
635 Sets (if not null):
636 - px, py - point on line,
637 - dist - distance to line,
638 - spdist - distance to point on line from segment beginning,
639 - lpdist - distance to point on line from line beginning along line
640
641 \param points pointer to line_pnts structure
642 \param ux,uy,uz point coordinates
643 \param with_z flag if to use z coordinate (3D calculation)
644 \param[out] px,py,pz point on line
645 \param[out] dist distance to line
646 \param[out] spdist distance to point on line from segment beginning
647 \param[out] lpdist distance to point on line from line beginning along line
648
649 \return nearest segment (first is 1)
650 */
Vect_line_distance(const struct line_pnts * points,double ux,double uy,double uz,int with_z,double * px,double * py,double * pz,double * dist,double * spdist,double * lpdist)651 int Vect_line_distance(const struct line_pnts *points,
652 double ux, double uy, double uz,
653 int with_z,
654 double *px, double *py, double *pz,
655 double *dist, double *spdist, double *lpdist)
656 {
657 int i;
658 double distance;
659 double new_dist;
660 double tpx, tpy, tpz, tdist, tspdist, tlpdist = 0;
661 double dx, dy, dz;
662 int n_points;
663 int segment;
664
665 n_points = points->n_points;
666
667 if (n_points == 1) {
668 distance =
669 dig_distance2_point_to_line(ux, uy, uz, points->x[0],
670 points->y[0], points->z[0],
671 points->x[0], points->y[0],
672 points->z[0], with_z, NULL, NULL,
673 NULL, NULL, NULL);
674 tpx = points->x[0];
675 tpy = points->y[0];
676 tpz = points->z[0];
677 tdist = sqrt(distance);
678 tspdist = 0;
679 tlpdist = 0;
680 segment = 0;
681
682 }
683 else {
684
685 distance =
686 dig_distance2_point_to_line(ux, uy, uz, points->x[0],
687 points->y[0], points->z[0],
688 points->x[1], points->y[1],
689 points->z[1], with_z, NULL, NULL,
690 NULL, NULL, NULL);
691 segment = 1;
692
693 for (i = 1; i < n_points - 1; i++) {
694 new_dist = dig_distance2_point_to_line(ux, uy, uz,
695 points->x[i], points->y[i],
696 points->z[i],
697 points->x[i + 1],
698 points->y[i + 1],
699 points->z[i + 1], with_z,
700 NULL, NULL, NULL, NULL,
701 NULL);
702 if (new_dist < distance) {
703 distance = new_dist;
704 segment = i + 1;
705 }
706 }
707
708 /* we have segment and now we can recalculate other values (speed) */
709 new_dist = dig_distance2_point_to_line(ux, uy, uz,
710 points->x[segment - 1],
711 points->y[segment - 1],
712 points->z[segment - 1],
713 points->x[segment],
714 points->y[segment],
715 points->z[segment], with_z,
716 &tpx, &tpy, &tpz, &tspdist,
717 NULL);
718
719 /* calculate distance from beginning of line */
720 if (lpdist) {
721 tlpdist = 0;
722 for (i = 0; i < segment - 1; i++) {
723 dx = points->x[i + 1] - points->x[i];
724 dy = points->y[i + 1] - points->y[i];
725 if (with_z)
726 dz = points->z[i + 1] - points->z[i];
727 else
728 dz = 0;
729
730 tlpdist += hypot(hypot(dx, dy), dz);
731 }
732 tlpdist += tspdist;
733 }
734 tdist = sqrt(distance);
735 }
736
737 if (px)
738 *px = tpx;
739 if (py)
740 *py = tpy;
741 if (pz && with_z)
742 *pz = tpz;
743 if (dist)
744 *dist = tdist;
745 if (spdist)
746 *spdist = tspdist;
747 if (lpdist)
748 *lpdist = tlpdist;
749
750 return (segment);
751 }
752
753 /*!
754 \brief Calculate geodesic distance of point to line in meters.
755
756 Sets (if not null):
757 - px, py - point on line,
758 - dist - distance to line,
759 - spdist - distance to point on line from segment beginning,
760 - lpdist - distance to point on line from line beginning along line
761
762 \param points pointer to line_pnts structure
763 \param ux,uy,uz point coordinates
764 \param with_z flag if to use z coordinate (3D calculation)
765 \param[out] px,py,pz point on line
766 \param[out] dist distance to line
767 \param[out] spdist distance to point on line from segment beginning
768 \param[out] lpdist distance to point on line from line beginning along line
769
770 \return nearest segment (first is 1)
771 */
Vect_line_geodesic_distance(const struct line_pnts * points,double ux,double uy,double uz,int with_z,double * px,double * py,double * pz,double * dist,double * spdist,double * lpdist)772 int Vect_line_geodesic_distance(const struct line_pnts *points,
773 double ux, double uy, double uz,
774 int with_z,
775 double *px, double *py, double *pz,
776 double *dist, double *spdist,
777 double *lpdist)
778 {
779 int i;
780 double distance;
781 double new_dist;
782 double tpx, tpy, tpz, ttpx, ttpy, ttpz;
783 double tdist, tspdist, tlpdist = 0, tlpdistseg;
784 double dz;
785 int n_points;
786 int segment;
787
788 G_begin_distance_calculations();
789
790 n_points = points->n_points;
791
792 if (n_points == 1) {
793 distance = G_distance(ux, uy, points->x[0], points->y[0]);
794 if (with_z)
795 distance = hypot(distance, uz - points->z[0]);
796
797 tpx = points->x[0];
798 tpy = points->y[0];
799 tpz = points->z[0];
800 tdist = distance;
801 tspdist = 0;
802 tlpdist = 0;
803 segment = 0;
804 }
805 else {
806 distance =
807 dig_distance2_point_to_line(ux, uy, uz, points->x[0],
808 points->y[0], points->z[0],
809 points->x[1], points->y[1],
810 points->z[1], with_z,
811 &tpx, &tpy, &tpz,
812 NULL, NULL);
813
814 distance = G_distance(ux, uy, tpx, tpy);
815 if (with_z)
816 distance = hypot(distance, uz - tpz);
817
818 segment = 1;
819
820 for (i = 1; i < n_points - 1; i++) {
821 new_dist = dig_distance2_point_to_line(ux, uy, uz,
822 points->x[i], points->y[i],
823 points->z[i],
824 points->x[i + 1],
825 points->y[i + 1],
826 points->z[i + 1], with_z,
827 &ttpx, &ttpy, &ttpz,
828 NULL, NULL);
829
830 new_dist = G_distance(ux, uy, ttpx, ttpy);
831 if (with_z)
832 new_dist = hypot(new_dist, uz - ttpz);
833
834 if (new_dist < distance) {
835 distance = new_dist;
836 segment = i + 1;
837 tpx = ttpx;
838 tpy = ttpy;
839 tpz = ttpz;
840 }
841 }
842
843 /* calculate distance from beginning of segment */
844 tspdist = G_distance(points->x[segment - 1], points->y[segment - 1],
845 tpx, tpy);
846 if (with_z) {
847 dz = points->z[segment - 1] - tpz;
848 tspdist += hypot(tspdist, dz);
849 }
850
851 /* calculate distance from beginning of line */
852 if (lpdist) {
853 tlpdist = 0;
854 for (i = 0; i < segment - 1; i++) {
855 tlpdistseg = G_distance(points->x[i], points->y[i],
856 points->x[i + 1], points->y[i + 1]);
857
858 if (with_z) {
859 dz = points->z[i + 1] - points->z[i];
860 tlpdistseg += hypot(tlpdistseg, dz);
861 }
862
863 tlpdist += tlpdistseg;
864 }
865 tlpdist += tspdist;
866 }
867 tdist = distance;
868 }
869
870 if (px)
871 *px = tpx;
872 if (py)
873 *py = tpy;
874 if (pz && with_z)
875 *pz = tpz;
876 if (dist)
877 *dist = tdist;
878 if (spdist)
879 *spdist = tspdist;
880 if (lpdist)
881 *lpdist = tlpdist;
882
883 return (segment);
884 }
885
886
887 /*!
888 \brief Calculate distance of 2 points
889
890 Simply uses Pythagoras.
891
892 \param x1,y1,z1 first point
893 \param x2,y2,z2 second point
894 \param with_z use z coordinate
895
896 \return distance
897 */
Vect_points_distance(double x1,double y1,double z1,double x2,double y2,double z2,int with_z)898 double Vect_points_distance(double x1, double y1, double z1, /* point 1 */
899 double x2, double y2, double z2, /* point 2 */
900 int with_z)
901 {
902 double dx, dy, dz;
903
904
905 dx = x2 - x1;
906 dy = y2 - y1;
907 dz = z2 - z1;
908
909 if (with_z)
910 return hypot(hypot(dx, dy), dz);
911 else
912 return hypot(dx, dy);
913
914 }
915
916 /*!
917 \brief Get bounding box of line
918
919 \param Points pointer to line_pnts structure
920 \param[out] Box bounding box
921 */
Vect_line_box(const struct line_pnts * Points,struct bound_box * Box)922 void Vect_line_box(const struct line_pnts *Points, struct bound_box *Box)
923 {
924 dig_line_box(Points, Box);
925 }
926
927 /*!
928 \brief Reverse the order of vertices
929
930 \param Points pointer to line_pnts structure to be changed
931 */
Vect_line_reverse(struct line_pnts * Points)932 void Vect_line_reverse(struct line_pnts *Points)
933 {
934 int i, j, np;
935 double x, y, z;
936
937 np = (int)Points->n_points / 2;
938
939 for (i = 0; i < np; i++) {
940 j = Points->n_points - i - 1;
941 x = Points->x[i];
942 y = Points->y[i];
943 z = Points->z[i];
944 Points->x[i] = Points->x[j];
945 Points->y[i] = Points->y[j];
946 Points->z[i] = Points->z[j];
947 Points->x[j] = x;
948 Points->y[j] = y;
949 Points->z[j] = z;
950 }
951 }
952
953 /*!
954 \brief Fetches FIRST category number for given vector line and field
955
956 \param Map pointer to Map_info structure
957 \param line line id
958 \param field layer number
959
960 \return -1 no category
961 \return category number (>=0)
962 */
Vect_get_line_cat(const struct Map_info * Map,int line,int field)963 int Vect_get_line_cat(const struct Map_info *Map, int line, int field)
964 {
965
966 static struct line_cats *cats = NULL;
967 int cat, ltype;
968
969 if (cats == NULL)
970 cats = Vect_new_cats_struct();
971
972 ltype = Vect_read_line(Map, NULL, cats, line);
973 Vect_cat_get(cats, field, &cat);
974 G_debug(3, "Vect_get_line_cat: display line %d, ltype %d, cat %d", line,
975 ltype, cat);
976
977 return cat;
978 }
979