GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-3a95b7a4e6
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 =
124  Vect_find_line_list(Map, x, y, z, GV_LINES, thresh,
125  WITHOUT_Z, 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,
130  found);
131  if (extend_lines(Map, !j, line, found, parallel, thresh,
132  List)) {
133  G_debug(3,
134  "Vedit_extend_lines(): lines=%d,%d -> extended",
135  line, found);
136  nlines_modified += 2;
137  extended = 1;
138  }
139  }
140 
141  Vect_list_append(List_exclude, found);
142  } while (List_found->n_values > 0 && !extended);
143  }
144  }
145 
146  Vect_destroy_list(List_exclude);
147  Vect_destroy_list(List_found);
148 
149  return nlines_modified;
150 }
151 
152 int extend_lines(struct Map_info *Map, int first, int line_from, int line_to,
153  int parallel, double thresh, struct ilist *List UNUSED)
154 {
155  /* TODO: If line_from extends to the i'th segment of line_to but the
156  * line_from node is closest to the j'th segment of line_to, this function
157  * wouldn't work because it only checks intersection of the start/end
158  * segment of line_from and the closest segment of line_to (i'th segment).
159  */
160  int line_new;
161  int type_from, type_to;
162 
163  struct line_pnts *Points_from, *Points_to, *Points_final;
164  struct line_cats *Cats_from, *Cats_to;
165 
166  Points_from = Vect_new_line_struct();
167  Points_to = Vect_new_line_struct();
168  Points_final = Vect_new_line_struct();
169  Cats_from = Vect_new_cats_struct();
170  Cats_to = Vect_new_cats_struct();
171 
172  type_from = Vect_read_line(Map, Points_from, Cats_from, line_from);
173  type_to = Vect_read_line(Map, Points_to, Cats_to, line_to);
174 
175  line_new = 0;
176  if (!(type_from & GV_LINES) || !(type_to & GV_LINES))
177  line_new = -1;
178 
179  /* avoid too much indentation */
180  do {
181  int n_points, seg, is, line_to_extended;
182  double x, y, px, py, x1, y1;
183  double dist, spdist, lpdist, length;
184  double angle_t, angle_f;
185 
186  if (line_new == -1)
187  break;
188 
189  n_points = Points_from->n_points - 1;
190 
191  if (first) {
192  x = Points_from->x[0];
193  y = Points_from->y[0];
194  }
195  else {
196  x = Points_from->x[n_points];
197  y = Points_from->y[n_points];
198  }
199  seg = Vect_line_distance(Points_to, x, y, 0.0, WITHOUT_Z, &px, &py,
200  NULL, &dist, &spdist, &lpdist);
201 
202  if (!(seg > 0 && dist > 0.0 && (thresh < 0. || dist <= thresh)))
203  break;
204 
205  /* lines in threshold */
206  length = first ? 0 : Vect_line_length(Points_from);
207 
208  /* find angles */
209  if (!Vect_point_on_line(Points_from, length, NULL, NULL, NULL, &angle_f,
210  NULL) ||
211  !Vect_point_on_line(Points_to, lpdist, NULL, NULL, NULL, &angle_t,
212  NULL))
213  break;
214 
215  line_to_extended = 0;
216 
217  /* extend both lines and find intersection */
218  if (!find_extended_intersection(x, y, angle_f, px, py, angle_t, &x1,
219  &y1)) {
220  /* parallel lines */
221  if (!parallel)
222  break;
223 
224  x1 = px;
225  y1 = py;
226  if (first)
227  Vect_line_insert_point(Points_from, 0, x1, y1, 0.0);
228  else
229  Vect_append_point(Points_from, x1, y1, 0.0);
230  }
231  else {
232  /* skip if extended into the wrong direction */
233  if (!check_extended_direction(x, y, angle_f, first, x1, y1))
234  break;
235 
236  /* skip if extended too far from line_from */
237  if (!Vect_line_distance(Points_from, x1, y1, 0.0, WITHOUT_Z, NULL,
238  NULL, NULL, &dist, NULL, NULL) ||
239  dist > thresh)
240  break;
241 
242  Vect_line_distance(Points_to, x1, y1, 0.0, WITHOUT_Z, NULL, NULL,
243  NULL, &dist, NULL, NULL);
244  /* if intersection point is not on line_to */
245  if (dist > TOL) {
246  double x2, y2;
247 
248  /* skip if not extended from a line_to node */
249  if (seg > 1 && seg < Points_to->n_points - 1)
250  break;
251 
252  if (seg == 1) {
253  line_to_extended = 1;
254 
255  x2 = Points_to->x[0];
256  y2 = Points_to->y[0];
257  }
258  else {
259  line_to_extended = 2;
260 
261  x2 = Points_to->x[Points_to->n_points - 1];
262  y2 = Points_to->y[Points_to->n_points - 1];
263  }
264 
265  /* skip if extended into the wrong direction */
266  if (!check_extended_direction(x2, y2, angle_t, seg == 1, x1,
267  y1))
268  break;
269  }
270  /* otherwise, split line_to later */
271 
272  /* lines extended -> extend/split line_to */
273  /* update line_from */
274  if (first) {
275  Points_from->x[0] = x1;
276  Points_from->y[0] = y1;
277  }
278  else {
279  Points_from->x[n_points] = x1;
280  Points_from->y[n_points] = y1;
281  }
282  }
283 
284  line_new = Vect_rewrite_line(Map, line_from, type_from, Points_from,
285  Cats_from);
286  /* Vect_list_append(List, line_new); */
287 
288  Vect_reset_line(Points_final);
289  if (line_to_extended == 1) {
290  /* extend line_to start node */
291  Vect_append_point(Points_final, x1, y1, 0.0);
292  for (is = 0; is < Points_to->n_points; is++)
293  Vect_append_point(Points_final, Points_to->x[is],
294  Points_to->y[is], Points_to->z[is]);
295  line_new =
296  Vect_rewrite_line(Map, line_to, type_to, Points_final, Cats_to);
297  }
298  else if (line_to_extended == 2) {
299  /* extend line_to end node */
300  for (is = 0; is < Points_to->n_points; is++)
301  Vect_append_point(Points_final, Points_to->x[is],
302  Points_to->y[is], Points_to->z[is]);
303  Vect_append_point(Points_final, x1, y1, 0.0);
304  line_new =
305  Vect_rewrite_line(Map, line_to, type_to, Points_final, Cats_to);
306  }
307  else {
308  int n_parts = 0;
309 
310  /* break line_to */
311  /* update line_to -- first part */
312  for (is = 0; is < seg; is++)
313  Vect_append_point(Points_final, Points_to->x[is],
314  Points_to->y[is], Points_to->z[is]);
315  Vect_append_point(Points_final, x1, y1, 0.0);
316 
317  if (Vect_line_length(Points_final) > 0) {
318  n_parts++;
319  line_new = Vect_rewrite_line(Map, line_to, type_to,
320  Points_final, Cats_to);
321  /* Vect_list_append(List, line_new); */
322  }
323 
324  /* write second part */
325  Vect_reset_line(Points_final);
326  Vect_append_point(Points_final, x1, y1, 0.0);
327  for (is = seg; is < Points_to->n_points; is++)
328  Vect_append_point(Points_final, Points_to->x[is],
329  Points_to->y[is], Points_to->z[is]);
330 
331  if (Vect_line_length(Points_final) > 0) {
332  if (n_parts > 0)
333  line_new =
334  Vect_write_line(Map, type_to, Points_final, Cats_to);
335  else
336  line_new = Vect_rewrite_line(Map, line_to, type_to,
337  Points_final, Cats_to);
338  /* Vect_list_append(List, line_new); */
339  }
340  }
341  } while (0);
342 
343  Vect_destroy_line_struct(Points_from);
344  Vect_destroy_line_struct(Points_to);
345  Vect_destroy_line_struct(Points_final);
346  Vect_destroy_cats_struct(Cats_from);
347  Vect_destroy_cats_struct(Cats_to);
348 
349  return line_new > 0 ? 1 : 0;
350 }
351 
352 static int find_extended_intersection(double x1, double y1, double angle1,
353  double x2, double y2, double angle2,
354  double *x, double *y)
355 {
356  double c1, s1, c2, s2, d, a;
357 
358  if (fabs(sin(angle1 - angle2)) <= TOL) {
359  /* two lines are parallel */
360  double angle;
361 
362  angle = atan2(y2 - y1, x2 - x1);
363  if (fabs(sin(angle - angle1)) <= TOL) {
364  /* they are on the same line */
365  *x = x2;
366  *y = y2;
367 
368  return 1;
369  }
370 
371  /* two lines don't intersect */
372  return 0;
373  }
374 
375  c1 = cos(angle1);
376  s1 = sin(angle1);
377  c2 = cos(angle2);
378  s2 = sin(angle2);
379  d = -c1 * s2 + c2 * s1;
380  if (d == 0.0)
381  /* shouldn't happen again */
382  return 0;
383 
384  a = (-s2 * (x2 - x1) + c2 * (y2 - y1)) / d;
385  *x = x1 + a * c1;
386  *y = y1 + a * s1;
387 
388  return 1;
389 }
390 
391 static int check_extended_direction(double x, double y, double angle,
392  int start_node, double extx, double exty)
393 {
394  double tmp;
395  int xdir, ydir, xext, yext;
396 
397  if (start_node)
398  angle += M_PI;
399 
400  /* expected directions */
401  tmp = cos(angle);
402  xdir = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
403  tmp = sin(angle);
404  ydir = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
405 
406  /* extended directions */
407  tmp = extx - x;
408  xext = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
409  tmp = exty - y;
410  yext = (fabs(tmp) <= TOL ? 0 : (tmp > 0 ? 1 : -1));
411 
412  if (xext != 0 && yext != 0) {
413  /* no if extended into the wrong direction */
414  if (xdir / xext <= 0 || ydir / yext <= 0)
415  return 0;
416  /* otherwise, ok */
417  }
418  else if (xext == 0 && yext == 0) {
419  /* snapped to the node, ok */
420  }
421  else if (xext == 0) {
422  /* vertical extension */
423  /* no if not expected or extended into the wrong direction */
424  if (xdir != 0 || ydir / yext <= 0)
425  return 0;
426  /* otherwise, ok */
427  }
428  else {
429  /* horizontal extension */
430  /* no if not expected or extended into the wrong direction */
431  if (ydir != 0 || xdir / xext <= 0)
432  return 0;
433  /* otherwise, ok */
434  }
435 
436  return 1;
437 }
#define NULL
Definition: ccmath.h:32
int G_debug(int, const char *,...) __attribute__((format(printf
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
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 Vect_get_line_nodes(struct Map_info *, int, int *, int *)
Get line nodes.
Definition: level_two.c:303
int Vect_get_node_coor(struct Map_info *, int, double *, double *, double *)
Get node coordinates.
Definition: level_two.c:274
double Vect_line_length(const struct line_pnts *)
Calculate line length, 3D-length in case of 3D vector line.
Definition: line.c:575
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:413
int Vect_get_line_type(struct Map_info *, int)
Get line type.
Definition: level_two.c:254
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_list(struct ilist *)
Frees all memory associated with a struct ilist, including the struct itself.
void Vect_destroy_cats_struct(struct line_cats *)
Frees all memory associated with line_cats structure, including the struct itself.
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
int Vect_read_line(struct Map_info *, struct line_pnts *, struct line_cats *, int)
Read vector feature (topological level required)
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:648
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
int Vect_line_alive(struct Map_info *, int)
Check if feature is alive or dead (topological level required)
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:176
off_t Vect_write_line(struct Map_info *, int, const struct line_pnts *, const struct line_cats *)
Writes a new feature.
struct line_pnts * Vect_new_line_struct(void)
Creates and initializes a line_pnts structure.
Definition: line.c:45
void Vect_reset_line(struct line_pnts *)
Reset line.
Definition: line.c:129
struct line_cats * Vect_new_cats_struct(void)
Creates and initializes line_cats structure.
int Vect_reset_list(struct ilist *)
Reset ilist structure.
int Vect_append_point(struct line_pnts *, double, double, double)
Appends one point to the end of a line.
Definition: line.c:148
#define GV_LINES
Definition: dig_defines.h:193
#define WITHOUT_Z
2D/3D vector data
Definition: dig_defines.h:171
#define GV_POINTS
Definition: dig_defines.h:192
#define TOL
Definition: extend.c:17
int extend_lines(struct Map_info *Map, int first, int line_from, int line_to, int parallel, double thresh, struct ilist *List UNUSED)
Definition: extend.c:152
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 UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition: gis.h:47
#define M_PI
Definition: gis.h:158
Vector map info.
Definition: dig_structs.h:1243
List of integers.
Definition: gis.h:708
int n_values
Number of values in the list.
Definition: gis.h:716
int * value
Array of values.
Definition: gis.h:712
Feature category info.
Definition: dig_structs.h:1677
Feature geometry info - coordinates.
Definition: dig_structs.h:1651
double * y
Array of Y coordinates.
Definition: dig_structs.h:1659
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