GRASS Programmer's Manual  6.5.svn(2014)-r66266
 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 
57 struct line_pnts *Vect_new_line_struct()
58 {
59  struct line_pnts *p;
60 
61  if (NULL == (p = Vect__new_line_struct()))
62  G_fatal_error("Vect_new_line_struct(): %s", _("Out of memory"));
63 
64  return p;
65 }
66 
67 struct line_pnts *Vect__new_line_struct()
68 {
69  struct line_pnts *p;
70 
71  p = (struct line_pnts *)malloc(sizeof(struct line_pnts));
72 
73  /* alloc_points MUST be initialized to zero */
74  if (p)
75  p->alloc_points = p->n_points = 0;
76 
77  if (p)
78  p->x = p->y = p->z = NULL;
79 
80  return p;
81 }
82 
90 int Vect_destroy_line_struct(struct line_pnts *p)
91 {
92  if (p) { /* probably a moot test */
93  if (p->alloc_points) {
94  G_free((char *)p->x);
95  G_free((char *)p->y);
96  G_free((char *)p->z);
97  }
98  G_free((char *)p);
99  }
100 
101  return 0;
102 }
103 
114 int
115 Vect_copy_xyz_to_pnts(struct line_pnts *Points, double *x, double *y,
116  double *z, int n)
117 {
118  register int i;
119 
120  if (0 > dig_alloc_points(Points, n))
121  return (-1);
122 
123  for (i = 0; i < n; i++) {
124  Points->x[i] = x[i];
125  Points->y[i] = y[i];
126  if (z != NULL)
127  Points->z[i] = z[i];
128  else
129  Points->z[i] = 0;
130  Points->n_points = n;
131  }
132 
133  return (0);
134 }
135 
136 
148 int Vect_reset_line(struct line_pnts *Points)
149 {
150  Points->n_points = 0;
151 
152  return 0;
153 }
154 
168 int Vect_append_point(struct line_pnts *Points, double x, double y, double z)
169 {
170  register int n;
171 
172  if (0 > dig_alloc_points(Points, Points->n_points + 1))
173  return (-1);
174 
175  n = Points->n_points;
176  Points->x[n] = x;
177  Points->y[n] = y;
178  Points->z[n] = z;
179 
180  return ++(Points->n_points);
181 }
182 
193 int
194 Vect_line_insert_point(struct line_pnts *Points, int index, double x,
195  double y, double z)
196 {
197  register int n;
198 
199  if (index < 0 || index > Points->n_points - 1)
200  G_fatal_error("%s Vect_line_insert_point()",
201  _("Index out of range in"));
202 
203  if (0 > dig_alloc_points(Points, Points->n_points + 1))
204  return (-1);
205 
206  /* move up */
207  for (n = Points->n_points; n > index; n--) {
208  Points->x[n] = Points->x[n - 1];
209  Points->y[n] = Points->y[n - 1];
210  Points->z[n] = Points->z[n - 1];
211  }
212 
213  Points->x[index] = x;
214  Points->y[index] = y;
215  Points->z[index] = z;
216  return ++(Points->n_points);
217 }
218 
227 int Vect_line_delete_point(struct line_pnts *Points, int index)
228 {
229  register int n;
230 
231  if (index < 0 || index > Points->n_points - 1)
232  G_fatal_error("%s Vect_line_insert_point()",
233  _("Index out of range in"));
234 
235  if (Points->n_points == 0)
236  return 0;
237 
238  /* move down */
239  for (n = index; n < Points->n_points - 1; n++) {
240  Points->x[n] = Points->x[n + 1];
241  Points->y[n] = Points->y[n + 1];
242  Points->z[n] = Points->z[n + 1];
243  }
244 
245  return --(Points->n_points);
246 }
247 
255 int Vect_line_prune(struct line_pnts *Points)
256 {
257  int i, j;
258 
259  if (Points->n_points > 0) {
260  j = 1;
261  for (i = 1; i < Points->n_points; i++) {
262  if (Points->x[i] != Points->x[j - 1] ||
263  Points->y[i] != Points->y[j - 1]
264  || Points->z[i] != Points->z[j - 1]) {
265  Points->x[j] = Points->x[i];
266  Points->y[j] = Points->y[i];
267  Points->z[j] = Points->z[i];
268  j++;
269  }
270  }
271  Points->n_points = j;
272  }
273 
274  return (Points->n_points);
275 }
276 
285 int Vect_line_prune_thresh(struct line_pnts *Points, double threshold)
286 {
287  int ret;
288 
289  ret = dig_prune(Points, threshold);
290 
291  if (ret < Points->n_points)
292  Points->n_points = ret;
293 
294  return (Points->n_points);
295 }
296 
311 int
312 Vect_append_points(struct line_pnts *Points, struct line_pnts *APoints,
313  int direction)
314 {
315  int i, n, on, an;
316 
317  on = Points->n_points;
318  an = APoints->n_points;
319  n = on + an;
320 
321  /* Should be OK, dig_alloc_points calls realloc */
322  if (0 > dig_alloc_points(Points, n))
323  return (-1);
324 
325  if (direction == GV_FORWARD) {
326  for (i = 0; i < an; i++) {
327  Points->x[on + i] = APoints->x[i];
328  Points->y[on + i] = APoints->y[i];
329  Points->z[on + i] = APoints->z[i];
330  }
331  }
332  else {
333  for (i = 0; i < an; i++) {
334  Points->x[on + i] = APoints->x[an - i - 1];
335  Points->y[on + i] = APoints->y[an - i - 1];
336  Points->z[on + i] = APoints->z[an - i - 1];
337  }
338  }
339 
340  Points->n_points = n;
341  return n;
342 }
343 
344 
358 int
359 Vect_copy_pnts_to_xyz(struct line_pnts *Points, double *x, double *y,
360  double *z, int *n)
361 {
362  register int i;
363 
364  for (i = 0; i < *n; i++) {
365  x[i] = Points->x[i];
366  y[i] = Points->y[i];
367  if (z != NULL)
368  z[i] = Points->z[i];
369  *n = Points->n_points;
370  }
371 
372  return (Points->n_points);
373 }
374 
392 int
393 Vect_point_on_line(struct line_pnts *Points, double distance,
394  double *x, double *y, double *z, double *angle,
395  double *slope)
396 {
397  int j, np, seg = 0;
398  double dist = 0, length;
399  double xp = 0, yp = 0, zp = 0, dx = 0, dy = 0, dz = 0, dxy =
400  0, dxyz, k, rest;
401 
402  G_debug(3, "Vect_point_on_line(): distance = %f", distance);
403  if ((distance < 0) || (Points->n_points < 2))
404  return 0;
405 
406  /* Check if first or last */
407  length = Vect_line_length(Points);
408  G_debug(3, " length = %f", length);
409  if (distance < 0 || distance > length) {
410  G_debug(3, " -> outside line");
411  return 0;
412  }
413 
414  np = Points->n_points;
415  if (distance == 0) {
416  G_debug(3, " -> first point");
417  xp = Points->x[0];
418  yp = Points->y[0];
419  zp = Points->z[0];
420  dx = Points->x[1] - Points->x[0];
421  dy = Points->y[1] - Points->y[0];
422  dz = Points->z[1] - Points->z[0];
423  dxy = hypot(dx, dy);
424  seg = 1;
425  }
426  else if (distance == length) {
427  G_debug(3, " -> last point");
428  xp = Points->x[np - 1];
429  yp = Points->y[np - 1];
430  zp = Points->z[np - 1];
431  dx = Points->x[np - 1] - Points->x[np - 2];
432  dy = Points->y[np - 1] - Points->y[np - 2];
433  dz = Points->z[np - 1] - Points->z[np - 2];
434  dxy = hypot(dx, dy);
435  seg = np - 1;
436  }
437  else {
438  for (j = 0; j < Points->n_points - 1; j++) {
439  /* dxyz = G_distance(Points->x[j], Points->y[j],
440  Points->x[j+1], Points->y[j+1]); */
441  dx = Points->x[j + 1] - Points->x[j];
442  dy = Points->y[j + 1] - Points->y[j];
443  dz = Points->z[j + 1] - Points->z[j];
444  dxy = hypot(dx, dy);
445  dxyz = hypot(dxy, dz);
446 
447  dist += dxyz;
448  if (dist >= distance) { /* point is on the current line part */
449  rest = distance - dist + dxyz; /* from first point of segment to point */
450  k = rest / dxyz;
451 
452  xp = Points->x[j] + k * dx;
453  yp = Points->y[j] + k * dy;
454  zp = Points->z[j] + k * dz;
455  seg = j + 1;
456  break;
457  }
458  }
459  }
460 
461  if (x != NULL)
462  *x = xp;
463  if (y != NULL)
464  *y = yp;
465  if (z != NULL)
466  *z = zp;
467 
468  /* calculate angle */
469  if (angle != NULL)
470  *angle = atan2(dy, dx);
471 
472  /* calculate slope */
473  if (slope != NULL)
474  *slope = atan2(dz, dxy);
475 
476  return seg;
477 }
478 
496 int
497 Vect_line_segment(struct line_pnts *InPoints, double start, double end,
498  struct line_pnts *OutPoints)
499 {
500  int i, seg1, seg2;
501  double length, tmp;
502  double x1, y1, z1, x2, y2, z2;
503 
504  G_debug(3, "Vect_line_segment(): start = %f, end = %f, n_points = %d",
505  start, end, InPoints->n_points);
506 
507  Vect_reset_line(OutPoints);
508 
509  if (start > end) {
510  tmp = start;
511  start = end;
512  end = tmp;
513  }
514 
515  /* Check start/end */
516  if (end < 0)
517  return 0;
518  length = Vect_line_length(InPoints);
519  if (start > length)
520  return 0;
521 
522  /* Find coordinates and segments of start/end */
523  seg1 = Vect_point_on_line(InPoints, start, &x1, &y1, &z1, NULL, NULL);
524  seg2 = Vect_point_on_line(InPoints, end, &x2, &y2, &z2, NULL, NULL);
525 
526  G_debug(3, " -> seg1 = %d seg2 = %d", seg1, seg2);
527 
528  if (seg1 == 0 || seg2 == 0) {
529  G_warning(_("Segment outside line, no segment created"));
530  return 0;
531  }
532 
533  Vect_append_point(OutPoints, x1, y1, z1);
534 
535  for (i = seg1; i < seg2; i++) {
536  Vect_append_point(OutPoints, InPoints->x[i], InPoints->y[i],
537  InPoints->z[i]);
538  };
539 
540  Vect_append_point(OutPoints, x2, y2, z2);
541  Vect_line_prune(OutPoints);
542 
543  return 1;
544 }
545 
555 double Vect_line_length(struct line_pnts *Points)
556 {
557  int j;
558  double dx, dy, dz, len = 0;
559 
560  if (Points->n_points < 2)
561  return 0;
562 
563  for (j = 0; j < Points->n_points - 1; j++) {
564  dx = Points->x[j + 1] - Points->x[j];
565  dy = Points->y[j + 1] - Points->y[j];
566  dz = Points->z[j + 1] - Points->z[j];
567  len += hypot(hypot(dx, dy), dz);
568  }
569 
570  return len;
571 }
572 
573 
583 double Vect_line_geodesic_length(struct line_pnts *Points)
584 {
585  int j, dc;
586  double dx, dy, dz, dxy, len = 0;
587 
589 
590  if (Points->n_points < 2)
591  return 0;
592 
593  for (j = 0; j < Points->n_points - 1; j++) {
594  if (dc == 2)
595  dxy =
596  G_geodesic_distance(Points->x[j], Points->y[j],
597  Points->x[j + 1], Points->y[j + 1]);
598  else {
599  dx = Points->x[j + 1] - Points->x[j];
600  dy = Points->y[j + 1] - Points->y[j];
601  dxy = hypot(dx, dy);
602  }
603 
604  dz = Points->z[j + 1] - Points->z[j];
605  len += hypot(dxy, dz);
606  }
607 
608  return len;
609 }
610 
630 int
631 Vect_line_distance(struct line_pnts *points,
632  double ux, double uy, double uz,
633  int with_z,
634  double *px, double *py, double *pz,
635  double *dist, double *spdist, double *lpdist)
636 {
637  register int i;
638  register double distance;
639  register double new_dist;
640  double tpx, tpy, tpz, tdist, tspdist, tlpdist = 0;
641  double dx, dy, dz;
642  register int n_points;
643  int segment;
644 
645  n_points = points->n_points;
646 
647  if (n_points == 1) {
648  distance =
649  dig_distance2_point_to_line(ux, uy, uz, points->x[0],
650  points->y[0], points->z[0],
651  points->x[0], points->y[0],
652  points->z[0], with_z, NULL, NULL,
653  NULL, NULL, NULL);
654  tpx = points->x[0];
655  tpy = points->y[0];
656  tpz = points->z[0];
657  tdist = sqrt(distance);
658  tspdist = 0;
659  tlpdist = 0;
660  segment = 0;
661 
662  }
663  else {
664 
665  distance =
666  dig_distance2_point_to_line(ux, uy, uz, points->x[0],
667  points->y[0], points->z[0],
668  points->x[1], points->y[1],
669  points->z[1], with_z, NULL, NULL,
670  NULL, NULL, NULL);
671  segment = 1;
672 
673  for (i = 1; i < n_points - 1; i++) {
674  new_dist = dig_distance2_point_to_line(ux, uy, uz,
675  points->x[i], points->y[i],
676  points->z[i],
677  points->x[i + 1],
678  points->y[i + 1],
679  points->z[i + 1], with_z,
680  NULL, NULL, NULL, NULL,
681  NULL);
682  if (new_dist < distance) {
683  distance = new_dist;
684  segment = i + 1;
685  }
686  }
687 
688  /* we have segment and now we can recalculate other values (speed) */
689  new_dist = dig_distance2_point_to_line(ux, uy, uz,
690  points->x[segment - 1],
691  points->y[segment - 1],
692  points->z[segment - 1],
693  points->x[segment],
694  points->y[segment],
695  points->z[segment], with_z,
696  &tpx, &tpy, &tpz, &tspdist,
697  NULL);
698 
699  /* calculate distance from beginning of line */
700  if (lpdist) {
701  tlpdist = 0;
702  for (i = 0; i < segment - 1; i++) {
703  dx = points->x[i + 1] - points->x[i];
704  dy = points->y[i + 1] - points->y[i];
705  if (with_z)
706  dz = points->z[i + 1] - points->z[i];
707  else
708  dz = 0;
709 
710  tlpdist += hypot(hypot(dx, dy), dz);
711  }
712  tlpdist += tspdist;
713  }
714  tdist = sqrt(distance);
715  }
716 
717  if (px)
718  *px = tpx;
719  if (py)
720  *py = tpy;
721  if (pz && with_z)
722  *pz = tpz;
723  if (dist)
724  *dist = tdist;
725  if (spdist)
726  *spdist = tspdist;
727  if (lpdist)
728  *lpdist = tlpdist;
729 
730  return (segment);
731 }
732 
733 
745 double Vect_points_distance(double x1, double y1, double z1, /* point 1 */
746  double x2, double y2, double z2, /* point 2 */
747  int with_z)
748 {
749  double dx, dy, dz;
750 
751 
752  dx = x2 - x1;
753  dy = y2 - y1;
754  dz = z2 - z1;
755 
756  if (with_z)
757  return hypot(hypot(dx, dy), dz);
758  else
759  return hypot(dx, dy);
760 
761 }
762 
771 int Vect_line_box(struct line_pnts *Points, BOUND_BOX * Box)
772 {
773  dig_line_box(Points, Box);
774  return 0;
775 }
776 
784 void Vect_line_reverse(struct line_pnts *Points)
785 {
786  int i, j, np;
787  double x, y, z;
788 
789  np = (int)Points->n_points / 2;
790 
791  for (i = 0; i < np; i++) {
792  j = Points->n_points - i - 1;
793  x = Points->x[i];
794  y = Points->y[i];
795  z = Points->z[i];
796  Points->x[i] = Points->x[j];
797  Points->y[i] = Points->y[j];
798  Points->z[i] = Points->z[j];
799  Points->x[j] = x;
800  Points->y[j] = y;
801  Points->z[j] = z;
802  }
803 }
804 
815 int Vect_get_line_cat(struct Map_info *Map, int line, int field)
816 {
817 
818  static struct line_cats *cats = NULL;
819  int cat, ltype;
820 
821  if (cats == NULL)
822  cats = Vect_new_cats_struct();
823 
824  ltype = Vect_read_line(Map, NULL, cats, line);
825  Vect_cat_get(cats, field, &cat);
826  G_debug(3, "Vect_get_line_cat: display line %d, ltype %d, cat %d", line,
827  ltype, cat);
828 
829  return cat;
830 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
int Vect_line_box(struct line_pnts *Points, BOUND_BOX *Box)
Get bounding box of line.
Definition: line.c:771
int G_begin_distance_calculations(void)
Begin distance calculations.
Definition: gis/distance.c:44
int dig_line_box(struct line_pnts *Points, BOUND_BOX *Box)
Definition: diglib/box.c:24
float Box[8][3]
Vertices for box.
Definition: gsd_objs.c:1420
struct line_pnts * Vect_new_line_struct()
Creates and initializes a struct line_pnts.
Definition: line.c:57
int Vect_line_prune_thresh(struct line_pnts *Points, double threshold)
Remove points in threshold.
Definition: line.c:285
int Vect_append_points(struct line_pnts *Points, struct line_pnts *APoints, int direction)
Appends points to the end of a line.
Definition: line.c:312
int Vect_reset_line(struct line_pnts *Points)
Reset line.
Definition: line.c:148
double dig_distance2_point_to_line(double x, double y, double z, double x1, double y1, double z1, double x2, double y2, double z2, int with_z, double *px, double *py, double *pz, double *pdist, int *status)
int y
Definition: plot.c:34
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:168
int Vect_copy_xyz_to_pnts(struct line_pnts *Points, double *x, double *y, double *z, int n)
Copy points from array to line structure.
Definition: line.c:115
int Vect_copy_pnts_to_xyz(struct line_pnts *Points, double *x, double *y, double *z, int *n)
Copy points from line structure to array.
Definition: line.c:359
void Vect_line_reverse(struct line_pnts *Points)
Reverse the order of vertices.
Definition: line.c:784
void * malloc(YYSIZE_T)
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:194
int Vect_point_on_line(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:393
int dig_alloc_points(struct line_pnts *points, int num)
Definition: struct_alloc.c:232
int
Definition: g3dcolor.c:48
struct line_pnts * Vect__new_line_struct(void)
Creates and initializes a struct line_pnts.
Definition: line.c:67
struct line_cats * Vect_new_cats_struct()
Creates and initializes line_cats structure.
return NULL
Definition: dbfopen.c:1394
double Vect_line_length(struct line_pnts *Points)
Calculate line length, 3D-length in case of 3D vector line.
Definition: line.c:555
int Vect_line_prune(struct line_pnts *Points)
Remove duplicate points, i.e. zero length segments.
Definition: line.c:255
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
int Vect_line_segment(struct line_pnts *InPoints, double start, double end, struct line_pnts *OutPoints)
Create line segment.
Definition: line.c:497
tuple Map
Definition: render.py:1310
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
double Vect_line_geodesic_length(struct line_pnts *Points)
Calculate line length.
Definition: line.c:583
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:745
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:227
for(cat=0;;cat++)
Definition: g3dcats.c:140
double G_geodesic_distance(double lon1, double lat1, double lon2, double lat2)
Calculates geodesic distance.
Definition: geodist.c:212
CELL cat
Definition: g3dcats.c:90
int G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
int Vect_line_distance(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 line distance.
Definition: line.c:631
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:815
int dig_prune(struct line_pnts *points, double thresh)
Definition: prune.c:74
int n
Definition: dataquad.c:291
int Vect_destroy_line_struct(struct line_pnts *p)
Frees all memory associated with a struct line_pnts, including the struct itself. ...
Definition: line.c:90
int Vect_read_line(struct Map_info *Map, struct line_pnts *line_p, struct line_cats *line_c, int line)
Read vector feature.
int Vect_cat_get(struct line_cats *Cats, int field, int *cat)
Get first found category of given field.