GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
line.c
Go to the documentation of this file.
1 
19 #include <stdlib.h>
20 #include <math.h>
21 #include <grass/gis.h>
22 #include <grass/Vect.h>
23 #include <grass/glocale.h>
24 
40 struct line_pnts *Vect__new_line_struct(void);
41 
55 struct line_pnts *Vect_new_line_struct()
56 {
57  struct line_pnts *p;
58 
59  if (NULL == (p = Vect__new_line_struct()))
60  G_fatal_error("Vect_new_line_struct(): %s", _("Out of memory"));
61 
62  return p;
63 }
64 
65 struct line_pnts *Vect__new_line_struct()
66 {
67  struct line_pnts *p;
68 
69  p = (struct line_pnts *)malloc(sizeof(struct line_pnts));
70 
71  /* alloc_points MUST be initialized to zero */
72  if (p)
73  p->alloc_points = p->n_points = 0;
74 
75  if (p)
76  p->x = p->y = p->z = NULL;
77 
78  return p;
79 }
80 
88 int Vect_destroy_line_struct(struct line_pnts *p)
89 {
90  if (p) { /* probably a moot test */
91  if (p->alloc_points) {
92  G_free((char *)p->x);
93  G_free((char *)p->y);
94  G_free((char *)p->z);
95  }
96  G_free((char *)p);
97  }
98 
99  return 0;
100 }
101 
112 int
113 Vect_copy_xyz_to_pnts(struct line_pnts *Points, double *x, double *y,
114  double *z, int n)
115 {
116  register int i;
117 
118  if (0 > dig_alloc_points(Points, n))
119  return (-1);
120 
121  for (i = 0; i < n; i++) {
122  Points->x[i] = x[i];
123  Points->y[i] = y[i];
124  if (z != NULL)
125  Points->z[i] = z[i];
126  else
127  Points->z[i] = 0;
128  Points->n_points = n;
129  }
130 
131  return (0);
132 }
133 
134 
146 int Vect_reset_line(struct line_pnts *Points)
147 {
148  Points->n_points = 0;
149 
150  return 0;
151 }
152 
166 int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
167 {
168  register int n;
169 
170  if (0 > dig_alloc_points(Points, Points->n_points + 1))
171  return (-1);
172 
173  n = Points->n_points;
174  Points->x[n] = x;
175  Points->y[n] = y;
176  Points->z[n] = z;
177 
178  return ++(Points->n_points);
179 }
180 
191 int
192 Vect_line_insert_point(struct line_pnts *Points, int index, double x,
193  double y, double z)
194 {
195  register int n;
196 
197  if (index < 0 || index > Points->n_points - 1)
198  G_fatal_error("%s Vect_line_insert_point()",
199  _("Index out of range in"));
200 
201  if (0 > dig_alloc_points(Points, Points->n_points + 1))
202  return (-1);
203 
204  /* move up */
205  for (n = Points->n_points; n > index; n--) {
206  Points->x[n] = Points->x[n - 1];
207  Points->y[n] = Points->y[n - 1];
208  Points->z[n] = Points->z[n - 1];
209  }
210 
211  Points->x[index] = x;
212  Points->y[index] = y;
213  Points->z[index] = z;
214  return ++(Points->n_points);
215 }
216 
225 int Vect_line_delete_point(struct line_pnts *Points, int index)
226 {
227  register int n;
228 
229  if (index < 0 || index > Points->n_points - 1)
230  G_fatal_error("%s Vect_line_insert_point()",
231  _("Index out of range in"));
232 
233  if (Points->n_points == 0)
234  return 0;
235 
236  /* move down */
237  for (n = index; n < Points->n_points - 1; n++) {
238  Points->x[n] = Points->x[n + 1];
239  Points->y[n] = Points->y[n + 1];
240  Points->z[n] = Points->z[n + 1];
241  }
242 
243  return --(Points->n_points);
244 }
245 
253 int Vect_line_prune(struct line_pnts *Points)
254 {
255  int i, j;
256 
257  if (Points->n_points > 0) {
258  j = 1;
259  for (i = 1; i < Points->n_points; i++) {
260  if (Points->x[i] != Points->x[j - 1] ||
261  Points->y[i] != Points->y[j - 1]
262  || Points->z[i] != Points->z[j - 1]) {
263  Points->x[j] = Points->x[i];
264  Points->y[j] = Points->y[i];
265  Points->z[j] = Points->z[i];
266  j++;
267  }
268  }
269  Points->n_points = j;
270  }
271 
272  return (Points->n_points);
273 }
274 
283 int Vect_line_prune_thresh(struct line_pnts *Points, double threshold)
284 {
285  int ret;
286 
287  ret = dig_prune(Points, threshold);
288 
289  if (ret < Points->n_points)
290  Points->n_points = ret;
291 
292  return (Points->n_points);
293 }
294 
309 int
310 Vect_append_points(struct line_pnts *Points, struct line_pnts *APoints,
311  int direction)
312 {
313  int i, n, on, an;
314 
315  on = Points->n_points;
316  an = APoints->n_points;
317  n = on + an;
318 
319  /* Should be OK, dig_alloc_points calls realloc */
320  if (0 > dig_alloc_points(Points, n))
321  return (-1);
322 
323  if (direction == GV_FORWARD) {
324  for (i = 0; i < an; i++) {
325  Points->x[on + i] = APoints->x[i];
326  Points->y[on + i] = APoints->y[i];
327  Points->z[on + i] = APoints->z[i];
328  }
329  }
330  else {
331  for (i = 0; i < an; i++) {
332  Points->x[on + i] = APoints->x[an - i - 1];
333  Points->y[on + i] = APoints->y[an - i - 1];
334  Points->z[on + i] = APoints->z[an - i - 1];
335  }
336  }
337 
338  Points->n_points = n;
339  return n;
340 }
341 
342 
356 int
357 Vect_copy_pnts_to_xyz(struct line_pnts *Points, double *x, double *y,
358  double *z, int *n)
359 {
360  register int i;
361 
362  for (i = 0; i < *n; i++) {
363  x[i] = Points->x[i];
364  y[i] = Points->y[i];
365  if (z != NULL)
366  z[i] = Points->z[i];
367  *n = Points->n_points;
368  }
369 
370  return (Points->n_points);
371 }
372 
390 int
391 Vect_point_on_line(struct line_pnts *Points, double distance,
392  double *x, double *y, double *z, double *angle,
393  double *slope)
394 {
395  int j, np, seg = 0;
396  double dist = 0, length;
397  double xp = 0, yp = 0, zp = 0, dx = 0, dy = 0, dz = 0, dxy =
398  0, dxyz, k, rest;
399 
400  G_debug(3, "Vect_point_on_line(): distance = %f", distance);
401  if ((distance < 0) || (Points->n_points < 2))
402  return 0;
403 
404  /* Check if first or last */
405  length = Vect_line_length(Points);
406  G_debug(3, " length = %f", length);
407  if (distance < 0 || distance > length) {
408  G_debug(3, " -> outside line");
409  return 0;
410  }
411 
412  np = Points->n_points;
413  if (distance == 0) {
414  G_debug(3, " -> first point");
415  xp = Points->x[0];
416  yp = Points->y[0];
417  zp = Points->z[0];
418  dx = Points->x[1] - Points->x[0];
419  dy = Points->y[1] - Points->y[0];
420  dz = Points->z[1] - Points->z[0];
421  dxy = hypot(dx, dy);
422  seg = 1;
423  }
424  else if (distance == length) {
425  G_debug(3, " -> last point");
426  xp = Points->x[np - 1];
427  yp = Points->y[np - 1];
428  zp = Points->z[np - 1];
429  dx = Points->x[np - 1] - Points->x[np - 2];
430  dy = Points->y[np - 1] - Points->y[np - 2];
431  dz = Points->z[np - 1] - Points->z[np - 2];
432  dxy = hypot(dx, dy);
433  seg = np - 1;
434  }
435  else {
436  for (j = 0; j < Points->n_points - 1; j++) {
437  /* dxyz = G_distance(Points->x[j], Points->y[j],
438  Points->x[j+1], Points->y[j+1]); */
439  dx = Points->x[j + 1] - Points->x[j];
440  dy = Points->y[j + 1] - Points->y[j];
441  dz = Points->z[j + 1] - Points->z[j];
442  dxy = hypot(dx, dy);
443  dxyz = hypot(dxy, dz);
444 
445  dist += dxyz;
446  if (dist >= distance) { /* point is on the current line part */
447  rest = distance - dist + dxyz; /* from first point of segment to point */
448  k = rest / dxyz;
449 
450  xp = Points->x[j] + k * dx;
451  yp = Points->y[j] + k * dy;
452  zp = Points->z[j] + k * dz;
453  seg = j + 1;
454  break;
455  }
456  }
457  }
458 
459  if (x != NULL)
460  *x = xp;
461  if (y != NULL)
462  *y = yp;
463  if (z != NULL)
464  *z = zp;
465 
466  /* calculate angle */
467  if (angle != NULL)
468  *angle = atan2(dy, dx);
469 
470  /* calculate slope */
471  if (slope != NULL)
472  *slope = atan2(dz, dxy);
473 
474  return seg;
475 }
476 
494 int
495 Vect_line_segment(struct line_pnts *InPoints, double start, double end,
496  struct line_pnts *OutPoints)
497 {
498  int i, seg1, seg2;
499  double length, tmp;
500  double x1, y1, z1, x2, y2, z2;
501 
502  G_debug(3, "Vect_line_segment(): start = %f, end = %f, n_points = %d",
503  start, end, InPoints->n_points);
504 
505  Vect_reset_line(OutPoints);
506 
507  if (start > end) {
508  tmp = start;
509  start = end;
510  end = tmp;
511  }
512 
513  /* Check start/end */
514  if (end < 0)
515  return 0;
516  length = Vect_line_length(InPoints);
517  if (start > length)
518  return 0;
519 
520  /* Find coordinates and segments of start/end */
521  seg1 = Vect_point_on_line(InPoints, start, &x1, &y1, &z1, NULL, NULL);
522  seg2 = Vect_point_on_line(InPoints, end, &x2, &y2, &z2, NULL, NULL);
523 
524  G_debug(3, " -> seg1 = %d seg2 = %d", seg1, seg2);
525 
526  if (seg1 == 0 || seg2 == 0) {
527  G_warning(_("Segment outside line, no segment created"));
528  return 0;
529  }
530 
531  Vect_append_point(OutPoints, x1, y1, z1);
532 
533  for (i = seg1; i < seg2; i++) {
534  Vect_append_point(OutPoints, InPoints->x[i], InPoints->y[i],
535  InPoints->z[i]);
536  };
537 
538  Vect_append_point(OutPoints, x2, y2, z2);
539 
540  return 1;
541 }
542 
552 double Vect_line_length(struct line_pnts *Points)
553 {
554  int j;
555  double dx, dy, dz, len = 0;
556 
557  if (Points->n_points < 2)
558  return 0;
559 
560  for (j = 0; j < Points->n_points - 1; j++) {
561  dx = Points->x[j + 1] - Points->x[j];
562  dy = Points->y[j + 1] - Points->y[j];
563  dz = Points->z[j + 1] - Points->z[j];
564  len += hypot(hypot(dx, dy), dz);
565  }
566 
567  return len;
568 }
569 
570 
580 double Vect_line_geodesic_length(struct line_pnts *Points)
581 {
582  int j, dc;
583  double dx, dy, dz, dxy, len = 0;
584 
586 
587  if (Points->n_points < 2)
588  return 0;
589 
590  for (j = 0; j < Points->n_points - 1; j++) {
591  if (dc == 2)
592  dxy =
593  G_geodesic_distance(Points->x[j], Points->y[j],
594  Points->x[j + 1], Points->y[j + 1]);
595  else {
596  dx = Points->x[j + 1] - Points->x[j];
597  dy = Points->y[j + 1] - Points->y[j];
598  dxy = hypot(dx, dy);
599  }
600 
601  dz = Points->z[j + 1] - Points->z[j];
602  len += hypot(dxy, dz);
603  }
604 
605  return len;
606 }
607 
627 int
628 Vect_line_distance(struct line_pnts *points,
629  double ux, double uy, double uz,
630  int with_z,
631  double *px, double *py, double *pz,
632  double *dist, double *spdist, double *lpdist)
633 {
634  register int i;
635  register double distance;
636  register double new_dist;
637  double tpx, tpy, tpz, tdist, tspdist, tlpdist = 0;
638  double dx, dy, dz;
639  register int n_points;
640  int segment;
641 
642  n_points = points->n_points;
643 
644  if (n_points == 1) {
645  distance =
646  dig_distance2_point_to_line(ux, uy, uz, points->x[0],
647  points->y[0], points->z[0],
648  points->x[0], points->y[0],
649  points->z[0], with_z, NULL, NULL,
650  NULL, NULL, NULL);
651  tpx = points->x[0];
652  tpy = points->y[0];
653  tpz = points->z[0];
654  tdist = sqrt(distance);
655  tspdist = 0;
656  tlpdist = 0;
657  segment = 0;
658 
659  }
660  else {
661 
662  distance =
663  dig_distance2_point_to_line(ux, uy, uz, points->x[0],
664  points->y[0], points->z[0],
665  points->x[1], points->y[1],
666  points->z[1], with_z, NULL, NULL,
667  NULL, NULL, NULL);
668  segment = 1;
669 
670  for (i = 1; i < n_points - 1; i++) {
671  new_dist = dig_distance2_point_to_line(ux, uy, uz,
672  points->x[i], points->y[i],
673  points->z[i],
674  points->x[i + 1],
675  points->y[i + 1],
676  points->z[i + 1], with_z,
677  NULL, NULL, NULL, NULL,
678  NULL);
679  if (new_dist < distance) {
680  distance = new_dist;
681  segment = i + 1;
682  }
683  }
684 
685  /* we have segment and now we can recalculate other values (speed) */
686  new_dist = dig_distance2_point_to_line(ux, uy, uz,
687  points->x[segment - 1],
688  points->y[segment - 1],
689  points->z[segment - 1],
690  points->x[segment],
691  points->y[segment],
692  points->z[segment], with_z,
693  &tpx, &tpy, &tpz, &tspdist,
694  NULL);
695 
696  /* calculate distance from beginning of line */
697  if (lpdist) {
698  tlpdist = 0;
699  for (i = 0; i < segment - 1; i++) {
700  dx = points->x[i + 1] - points->x[i];
701  dy = points->y[i + 1] - points->y[i];
702  if (with_z)
703  dz = points->z[i + 1] - points->z[i];
704  else
705  dz = 0;
706 
707  tlpdist += hypot(hypot(dx, dy), dz);
708  }
709  tlpdist += tspdist;
710  }
711  tdist = sqrt(distance);
712  }
713 
714  if (px)
715  *px = tpx;
716  if (py)
717  *py = tpy;
718  if (pz && with_z)
719  *pz = tpz;
720  if (dist)
721  *dist = tdist;
722  if (spdist)
723  *spdist = tspdist;
724  if (lpdist)
725  *lpdist = tlpdist;
726 
727  return (segment);
728 }
729 
730 
742 double Vect_points_distance(double x1, double y1, double z1, /* point 1 */
743  double x2, double y2, double z2, /* point 2 */
744  int with_z)
745 {
746  double dx, dy, dz;
747 
748 
749  dx = x2 - x1;
750  dy = y2 - y1;
751  dz = z2 - z1;
752 
753  if (with_z)
754  return hypot(hypot(dx, dy), dz);
755  else
756  return hypot(dx, dy);
757 
758 }
759 
768 int Vect_line_box(struct line_pnts *Points, BOUND_BOX * Box)
769 {
770  dig_line_box(Points, Box);
771  return 0;
772 }
773 
781 void Vect_line_reverse(struct line_pnts *Points)
782 {
783  int i, j, np;
784  double x, y, z;
785 
786  np = (int)Points->n_points / 2;
787 
788  for (i = 0; i < np; i++) {
789  j = Points->n_points - i - 1;
790  x = Points->x[i];
791  y = Points->y[i];
792  z = Points->z[i];
793  Points->x[i] = Points->x[j];
794  Points->y[i] = Points->y[j];
795  Points->z[i] = Points->z[j];
796  Points->x[j] = x;
797  Points->y[j] = y;
798  Points->z[j] = z;
799  }
800 }
801 
812 int Vect_get_line_cat(struct Map_info *Map, int line, int field)
813 {
814 
815  static struct line_cats *cats = NULL;
816  int cat, ltype;
817 
818  if (cats == NULL)
819  cats = Vect_new_cats_struct();
820 
821  ltype = Vect_read_line(Map, NULL, cats, line);
822  Vect_cat_get(cats, field, &cat);
823  G_debug(3, "Vect_get_line_cat: display line %d, ltype %d, cat %d", line,
824  ltype, cat);
825 
826  return cat;
827 }