GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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, 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  */
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  */
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  */
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  */
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 */
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 */
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  */
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  */
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 */
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 */
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  */
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 */
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  */
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 */
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 
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  */
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  */
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 
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 */
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 */
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 */
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  */
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 }
Bounding box.
Definition: dig_structs.h:65
#define GV_FORWARD
Line direction indicator forward/backward.
Definition: dig_defines.h:178
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int alloc_points
Allocated space for points.
Definition: dig_structs.h:1696
void Vect_line_box(const struct line_pnts *Points, struct bound_box *Box)
Get bounding box of line.
Definition: line.c:922
int Vect_get_num_line_points(const struct line_pnts *Points)
Get number of line points.
Definition: line.c:269
float Box[8][3]
Vertices for box.
Definition: gsd_objs.c:1421
struct line_pnts * Vect_new_line_struct()
Creates and initializes a line_pnts structure.
Definition: line.c:45
int dig_line_box(const struct line_pnts *, struct bound_box *)
int n_points
Number of points.
Definition: dig_structs.h:1692
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
int Vect_line_prune_thresh(struct line_pnts *Points, double threshold)
Remove points in threshold.
Definition: line.c:311
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
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
int Vect_line_segment(const struct line_pnts *InPoints, double start, double end, struct line_pnts *OutPoints)
Create line segment.
Definition: line.c:518
#define NULL
Definition: ccmath.h:32
#define x
double Vect_line_length(const struct line_pnts *Points)
Calculate line length, 3D-length in case of 3D vector line.
Definition: line.c:576
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:149
Feature category info.
Definition: dig_structs.h:1702
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:337
void * malloc(YYSIZE_T)
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
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:772
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
void Vect_line_reverse(struct line_pnts *Points)
Reverse the order of vertices.
Definition: line.c:932
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:177
double G_geodesic_distance(double, double, double, double)
Calculates geodesic distance.
Definition: geodist.c:196
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:651
int * cat
Array of categories.
Definition: dig_structs.h:1711
double dig_distance2_point_to_line(double, double, double, double, double, double, double, double, double, int, double *, double *, double *, double *, int *)
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:245
int dig_alloc_points(struct line_pnts *, int)
allocate room for &#39;num&#39; X and Y arrays in struct line_pnts
Definition: struct_alloc.c:336
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:383
struct line_pnts * Vect__new_line_struct(void)
Creates and initializes a struct line_pnts (internal use only)
Definition: line.c:55
Vector map info.
Definition: dig_structs.h:1259
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
int int G_begin_distance_calculations(void)
Begin distance calculations.
Definition: gis/distance.c:42
int Vect_line_prune(struct line_pnts *Points)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:281
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
void G_warning(const char *,...) __attribute__((format(printf
int Vect_get_line_cat(const struct Map_info *Map, int line, int field)
Fetches FIRST category number for given vector line and field.
Definition: line.c:963
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
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:898
#define _(str)
Definition: glocale.h:10
double G_distance(double, double, double, double)
Returns distance in meters.
Definition: gis/distance.c:75
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:211
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:415
void Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:130
double Vect_line_geodesic_length(const struct line_pnts *Points)
Calculate line length.
Definition: line.c:604
int Vect_read_line(const struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int dig_prune(struct line_pnts *, double)
Definition: prune.c:74
int G_debug(int, const char *,...) __attribute__((format(printf
int Vect_cat_get(const struct line_cats *, int, int *)
Get first found category of given field.