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