GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
extend.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/vedit/extend.c
3 
4  \brief Vedit library - extend lines (adopted from break.c)
5 
6  (C) 2017 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 Huidae Cho <grass4u gmail.com>
12 */
13 
14 #include <math.h>
15 #include <grass/vedit.h>
16 
17 #define TOL 1e-9
18 
19 static int extend_lines(struct Map_info *, int, int, int, int, double,
20  struct ilist *);
21 static int find_extended_intersection(double, double, double, double, double,
22  double, double *, double *);
23 static int check_extended_direction(double, double, double, int, double,
24  double);
25 
26 /*!
27  \brief Extend lines in given threshold
28 
29  \code
30  1. Extend first line only
31  \ \
32  id1 \ -> \
33  \
34  id2 ---------- -----+----
35 
36 
37  2. Extend both lines
38  \ \
39  id1 \ -> \
40  \
41  id2 --- +----
42 
43 
44  3. Extend first line when both are on the same line
45  id1 --- --- id2 -> -----+----
46 
47 
48  4. Connect two parallel lines (parallel=1)
49  id1 ------ -------
50  -> /
51  id2 ------ +-----
52 
53 
54  5. Don't connect two parallel lines (parallel=0)
55  id1 ------ ------
56  ->
57  id2 ------ ------
58  \endcode
59 
60  \param Map pointer to Map_info
61  \param List list of selected lines
62  \param nodes 1 for start node, 2 for end node, other for both
63  \param parallel connect parallel lines
64  \param thresh threshold value
65 
66  \return number of modified lines
67  */
68 int Vedit_extend_lines(struct Map_info *Map, struct ilist *List, int nodes,
69  int parallel, double thresh)
70 {
71  int nlines_modified;
72  int i, j, first_node, n_nodes;
73 
74  struct ilist *List_exclude, *List_found;
75 
76  nlines_modified = 0;
77 
78  List_exclude = Vect_new_list();
79  List_found = Vect_new_list();
80 
81  first_node = 0;
82  n_nodes = 2;
83 
84  switch (nodes) {
85  case 1:
86  n_nodes = 1;
87  break;
88  case 2:
89  first_node = 1;
90  break;
91  }
92 
93  /* collect lines to be modified */
94  for (i = 0; i < List->n_values; i++) {
95  int line, extended, node[2];
96 
97  line = List->value[i];
98 
99  if (!Vect_line_alive(Map, line))
100  continue;
101 
102  if (Vect_get_line_type(Map, line) & GV_POINTS)
103  continue;
104 
105  node[0] = node[1] = -1;
106  Vect_get_line_nodes(Map, line, &(node[0]), &(node[1]));
107  if (node[0] < 0 || node[1] < 0)
108  continue;
109 
110  extended = 0;
111  Vect_reset_list(List_exclude);
112  Vect_list_append(List_exclude, line);
113  for (j = first_node; j < n_nodes && !extended; j++) {
114  double x, y, z;
115 
116  /* for each line node find lines in threshold */
117  Vect_get_node_coor(Map, node[j], &x, &y, &z);
118 
119  do {
120  int found;
121 
122  /* find first nearest line */
123  found = Vect_find_line_list(Map, x, y, z,
124  GV_LINES, thresh, WITHOUT_Z,
125  List_exclude, List_found);
126 
127  if (found > 0 && Vect_line_alive(Map, found)) {
128  /* try to extend lines (given node) */
129  G_debug(3, "Vedit_extend_lines(): lines=%d,%d", line, found);
130  if (extend_lines(Map, !j, line, found, parallel, thresh,
131  List)) {
132  G_debug(3, "Vedit_extend_lines(): lines=%d,%d -> extended",
133  line, found);
134  nlines_modified += 2;
135  extended = 1;
136  }
137  }
138 
139  Vect_list_append(List_exclude, found);
140  } while(List_found->n_values > 0 && !extended);
141  }
142  }
143 
144  Vect_destroy_list(List_exclude);
145  Vect_destroy_list(List_found);
146 
147  return nlines_modified;
148 }
149 
150 int extend_lines(struct Map_info *Map, int first, int line_from, int line_to,
151  int parallel, double thresh, struct ilist *List)
152 {
153  /* TODO: If line_from extends to the i'th segment of line_to but the
154  * line_from node is closest to the j'th segment of line_to, this function
155  * wouldn't work because it only checks intersection of the start/end
156  * segment of line_from and the closest segment of line_to (i'th segment).
157  */
158  int line_new;
159  int type_from, type_to;
160 
161  struct line_pnts *Points_from, *Points_to, *Points_final;
162  struct line_cats *Cats_from, *Cats_to;
163 
164  Points_from = Vect_new_line_struct();
165  Points_to = Vect_new_line_struct();
166  Points_final = Vect_new_line_struct();
167  Cats_from = Vect_new_cats_struct();
168  Cats_to = Vect_new_cats_struct();
169 
170  type_from = Vect_read_line(Map, Points_from, Cats_from, line_from);
171  type_to = Vect_read_line(Map, Points_to, Cats_to, line_to);
172 
173  line_new = 0;
174  if (!(type_from & GV_LINES) || !(type_to & GV_LINES))
175  line_new = -1;
176 
177  /* avoid too much indentation */
178  do {
179  int n_points, seg, is, line_to_extended;
180  double x, y, px, py, x1, y1;
181  double dist, spdist, lpdist, length;
182  double angle_t, angle_f;
183 
184  if (line_new == -1)
185  break;
186 
187  n_points = Points_from->n_points - 1;
188 
189  if (first) {
190  x = Points_from->x[0];
191  y = Points_from->y[0];
192  }
193  else {
194  x = Points_from->x[n_points];
195  y = Points_from->y[n_points];
196  }
197  seg = Vect_line_distance(Points_to, x, y, 0.0, WITHOUT_Z,
198  &px, &py, NULL, &dist, &spdist, &lpdist);
199 
200  if (!(seg > 0 && dist > 0.0 && (thresh < 0. || dist <= thresh)))
201  break;
202 
203  /* lines in threshold */
204  length = first ? 0 : Vect_line_length(Points_from);
205 
206  /* find angles */
207  if (!Vect_point_on_line(Points_from, length, NULL, NULL, NULL, &angle_f,
208  NULL) ||
209  !Vect_point_on_line(Points_to, lpdist, NULL, NULL, NULL, &angle_t,
210  NULL))
211  break;
212 
213  line_to_extended = 0;
214 
215  /* extend both lines and find intersection */
216  if (!find_extended_intersection(x, y, angle_f, px, py, angle_t,
217  &x1, &y1)) {
218  /* parallel lines */
219  if (!parallel)
220  break;
221 
222  x1 = px;
223  y1 = py;
224  if (first)
225  Vect_line_insert_point(Points_from, 0, x1, y1, 0.0);
226  else
227  Vect_append_point(Points_from, x1, y1, 0.0);
228  } else {
229  /* skip if extended into the wrong direction */
230  if (!check_extended_direction(x, y, angle_f, first, x1, y1))
231  break;
232 
233  /* skip if extended too far from line_from */
234  if (!Vect_line_distance(Points_from, x1, y1, 0.0, WITHOUT_Z, NULL,
235  NULL, NULL, &dist, NULL, NULL) ||
236  dist > thresh)
237  break;
238 
239  Vect_line_distance(Points_to, x1, y1, 0.0, WITHOUT_Z, NULL, NULL,
240  NULL, &dist, NULL, NULL);
241  /* if intersection point is not on line_to */
242  if (dist > TOL) {
243  double x2, y2;
244 
245  /* skip if not extended from a line_to node */
246  if (seg > 1 && seg < Points_to->n_points - 1)
247  break;
248 
249  if (seg == 1) {
250  line_to_extended = 1;
251 
252  x2 = Points_to->x[0];
253  y2 = Points_to->y[0];
254  } else {
255  line_to_extended = 2;
256 
257  x2 = Points_to->x[Points_to->n_points - 1];
258  y2 = Points_to->y[Points_to->n_points - 1];
259  }
260 
261  /* skip if extended into the wrong direction */
262  if (!check_extended_direction(x2, y2, angle_t, seg == 1, x1,
263  y1))
264  break;
265  }
266  /* otherwise, split line_to later */
267 
268  /* lines extended -> extend/split line_to */
269  /* update line_from */
270  if (first) {
271  Points_from->x[0] = x1;
272  Points_from->y[0] = y1;
273  } else {
274  Points_from->x[n_points] = x1;
275  Points_from->y[n_points] = y1;
276  }
277  }
278 
279  line_new = Vect_rewrite_line(Map, line_from, type_from, Points_from,
280  Cats_from);
281  /* Vect_list_append(List, line_new); */
282 
283  Vect_reset_line(Points_final);
284  if (line_to_extended == 1) {
285  /* extend line_to start node */
286  Vect_append_point(Points_final, x1, y1, 0.0);
287  for (is = 0; is < Points_to->n_points; is++)
288  Vect_append_point(Points_final, Points_to->x[is],
289  Points_to->y[is], Points_to->z[is]);
290  line_new = Vect_rewrite_line(Map, line_to, type_to, Points_final,
291  Cats_to);
292  } else if (line_to_extended == 2) {
293  /* extend line_to end node */
294  for (is = 0; is < Points_to->n_points; is++)
295  Vect_append_point(Points_final, Points_to->x[is],
296  Points_to->y[is], Points_to->z[is]);
297  Vect_append_point(Points_final, x1, y1, 0.0);
298  line_new = Vect_rewrite_line(Map, line_to, type_to, Points_final,
299  Cats_to);
300  } else {
301  int n_parts = 0;
302 
303  /* break line_to */
304  /* update line_to -- first part */
305  for (is = 0; is < seg; is++)
306  Vect_append_point(Points_final, Points_to->x[is],
307  Points_to->y[is], Points_to->z[is]);
308  Vect_append_point(Points_final, x1, y1, 0.0);
309 
310  if (Vect_line_length(Points_final) > 0) {
311  n_parts++;
312  line_new = Vect_rewrite_line(Map, line_to, type_to,
313  Points_final, Cats_to);
314  /* Vect_list_append(List, line_new); */
315  }
316 
317  /* write second part */
318  Vect_reset_line(Points_final);
319  Vect_append_point(Points_final, x1, y1, 0.0);
320  for (is = seg; is < Points_to->n_points; is++)
321  Vect_append_point(Points_final, Points_to->x[is],
322  Points_to->y[is], Points_to->z[is]);
323 
324  if (Vect_line_length(Points_final) > 0) {
325  if (n_parts > 0)
326  line_new = Vect_write_line(Map, type_to, Points_final,
327  Cats_to);
328  else
329  line_new = Vect_rewrite_line(Map, line_to, type_to,
330  Points_final, Cats_to);
331  /* Vect_list_append(List, line_new); */
332  }
333  }
334  } while(0);
335 
336  Vect_destroy_line_struct(Points_from);
337  Vect_destroy_line_struct(Points_to);
338  Vect_destroy_line_struct(Points_final);
339  Vect_destroy_cats_struct(Cats_from);
340  Vect_destroy_cats_struct(Cats_to);
341 
342  return line_new > 0 ? 1 : 0;
343 }
344 
345 static int find_extended_intersection(double x1, double y1, double angle1,
346  double x2, double y2, double angle2,
347  double *x, double *y)
348 {
349  double c1, s1, c2, s2, d, a;
350 
351  if (fabs(sin(angle1 - angle2)) <= TOL) {
352  /* two lines are parallel */
353  double angle;
354 
355  angle = atan2(y2 - y1, x2 - x1);
356  if (fabs(sin(angle - angle1)) <= TOL) {
357  /* they are on the same line */
358  *x = x2;
359  *y = y2;
360 
361  return 1;
362  }
363 
364  /* two lines don't intersect */
365  return 0;
366  }
367 
368  c1 = cos(angle1);
369  s1 = sin(angle1);
370  c2 = cos(angle2);
371  s2 = sin(angle2);
372  d = -c1 * s2 + c2 * s1;
373  if (d == 0.0)
374  /* shouldn't happen again */
375  return 0;
376 
377  a = (-s2 * (x2 - x1) + c2 * (y2 - y1)) / d;
378  *x = x1 + a * c1;
379  *y = y1 + a * s1;
380 
381  return 1;
382 }
383 
384 static int check_extended_direction(double x, double y, double angle,
385  int start_node, double extx, double exty)
386 {
387  double tmp;
388  int xdir, ydir, xext, yext;
389 
390  if (start_node)
391  angle += M_PI;
392 
393  /* expected directions */
394  tmp = cos(angle);
395  xdir = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
396  tmp = sin(angle);
397  ydir = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
398 
399  /* extended directions */
400  tmp = extx - x;
401  xext = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
402  tmp = exty - y;
403  yext = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
404 
405  if (xext != 0 && yext != 0) {
406  /* no if extended into the wrong direction */
407  if (xdir / xext <= 0 || ydir / yext <= 0)
408  return 0;
409  /* otherwise, ok */
410  } else if (xext == 0 && yext == 0) {
411  /* snapped to the node, ok */
412  } else if (xext == 0) {
413  /* vertical extension */
414  /* no if not expected or extended into the wrong direction */
415  if (xdir != 0 || ydir / yext <= 0)
416  return 0;
417  /* otherwise, ok */
418  } else {
419  /* horizontal extension */
420  /* no if not expected or extended into the wrong direction */
421  if (ydir != 0 || xdir / xext <= 0)
422  return 0;
423  /* otherwise, ok */
424  }
425 
426  return 1;
427 }
int Vect_reset_list(struct ilist *)
Reset ilist structure.
#define TOL
Definition: extend.c:17
off_t Vect_rewrite_line(struct Map_info *, off_t, int, const struct line_pnts *, const struct line_cats *)
Rewrites existing feature (topological level required)
int n_points
Number of points.
Definition: dig_structs.h:1692
int n_values
Number of values in the list.
Definition: gis.h:698
#define GV_POINTS
Definition: dig_defines.h:191
int Vect_line_distance(const struct line_pnts *, double, double, double, int, double *, double *, double *, double *, double *, double *)
Calculate distance of point to line.
Definition: line.c:651
#define M_PI
Definition: gis.h:134
double Vect_line_length(const struct line_pnts *)
Calculate line length, 3D-length in case of 3D vector line.
Definition: line.c:576
#define NULL
Definition: ccmath.h:32
#define x
Feature category info.
Definition: dig_structs.h:1702
int Vect_get_node_coor(const struct Map_info *, int, double *, double *, double *)
Get node coordinates.
Definition: level_two.c:278
double * x
Array of X coordinates.
Definition: dig_structs.h:1680
Feature geometry info - coordinates.
Definition: dig_structs.h:1675
void Vect_destroy_list(struct ilist *)
Frees all memory associated with a struct ilist, including the struct itself.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
void Vect_destroy_cats_struct(struct line_cats *)
Frees all memory associated with line_cats structure, including the struct itself.
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:149
int Vedit_extend_lines(struct Map_info *Map, struct ilist *List, int nodes, int parallel, double thresh)
Extend lines in given threshold.
Definition: extend.c:68
#define WITHOUT_Z
2D/3D vector data
Definition: dig_defines.h:170
int Vect_get_line_nodes(const struct Map_info *, int, int *, int *)
Get line nodes.
Definition: level_two.c:307
Vector map info.
Definition: dig_structs.h:1259
double * y
Array of Y coordinates.
Definition: dig_structs.h:1684
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
int Vect_line_insert_point(struct line_pnts *, int, double, double, double)
Insert new point at index position and move all old points at that position and above up...
Definition: line.c:177
int Vect_line_alive(const struct Map_info *, int)
Check if feature is alive or dead (topological level required)
double * z
Array of Z coordinates.
Definition: dig_structs.h:1688
List of integers.
Definition: gis.h:689
int * value
Array of values.
Definition: gis.h:694
#define GV_LINES
Definition: dig_defines.h:192
int Vect_point_on_line(const struct line_pnts *, double, double *, double *, double *, double *, double *)
Find point on line in the specified distance.
Definition: line.c:415
int Vect_find_line_list(struct Map_info *, double, double, double, int, double, int, const struct ilist *, struct ilist *)
Find the nearest line(s).
void Vect_destroy_line_struct(struct line_pnts *)
Frees all memory associated with a line_pnts structure, including the structure itself.
Definition: line.c:77
int Vect_read_line(const struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
int G_debug(int, const char *,...) __attribute__((format(printf
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:130
int Vect_get_line_type(const struct Map_info *, int)
Get line type.
Definition: level_two.c:258