1 #include <stdlib.h>
2 #include <grass/gis.h>
3 #include <grass/vector.h>
4 #include <grass/glocale.h>
5 #include "proto.h"
6 
7 /**
8  * \brief Consolidate network arcs (edge) based on given point vector map (nodes)
9  *
10  * If there is no connection between network edge and point, new edge
11  * is added, the line broken, and new point added to nfield layer
12  *
13  * \param In,Points input vector maps
14  * \param Out output vector map
15  * \param nfield nodes layer
16  * \param thresh threshold value to find neareast line
17  *
18  * \return number of new arcs
19  */
connect_arcs(struct Map_info * In,struct Map_info * Pnts,struct Map_info * Out,int afield,int nfield,double thresh,int snap)20 int connect_arcs(struct Map_info *In, struct Map_info *Pnts,
21 		 struct Map_info *Out, int afield, int nfield,
22 		 double thresh, int snap)
23 {
24     int narcs;
25     int type, line, seg, i, ltype, broken;
26     double px, py, pz, spdist, dist;
27 
28     struct line_pnts *Points, *Pline, *Pout;
29     struct line_cats *Cats, *Cline, *Cnew;
30     int maxcat, findex, ncats;
31 
32     struct ilist *exclude_list;
33     int pline;
34 
35     exclude_list = Vect_new_list();
36 
37     narcs = 0;
38 
39     Points = Vect_new_line_struct();
40     Pline = Vect_new_line_struct();
41     Pout = Vect_new_line_struct();
42     Cats = Vect_new_cats_struct();
43     Cline = Vect_new_cats_struct();
44     Cnew = Vect_new_cats_struct();
45 
46     /* rewrite all primitives to output file */
47     Vect_copy_map_lines(In, Out);
48     Vect_build_partial(Out, GV_BUILD_BASE);
49 
50     findex = Vect_cidx_get_field_index(In, afield);
51     ncats = Vect_cidx_get_num_cats_by_index(In, findex);
52     Vect_cidx_get_cat_by_index(In, findex, ncats - 1, &maxcat, &type, &line);
53 
54 
55     /* go thorough all points in point map and write a new arcs if missing */
56     pz = 0.0;
57     while ((type = Vect_read_next_line(Pnts, Points, Cats)) >= 0) {
58 	if (type != GV_POINT)
59 	    continue;
60 
61 	/* find the nearest line in given threshold */
62 	line = Vect_find_line_list(Out, Points->x[0], Points->y[0], Points->z[0],
63 				   GV_LINES, thresh, WITHOUT_Z,
64 				   exclude_list, NULL);
65 
66 	if (line < 1 || !Vect_line_alive(Out, line))
67 	    continue;
68 
69 	ltype = Vect_read_line(Out, Pline, Cline, line);
70 
71 	/* find point on the line */
72 	seg = Vect_line_distance(Pline,
73 				 Points->x[0], Points->y[0], Points->z[0],
74 				 WITHOUT_Z, &px, &py, &pz, &dist, &spdist,
75 				 NULL);
76 
77 	if (seg == 0)
78 	    G_fatal_error(_("Failed to find intersection segment"));
79 	/* break the line */
80 	broken = 0;
81 	Vect_reset_line(Pout);
82 	for (i = 0; i < seg; i++) {
83 	    Vect_append_point(Pout, Pline->x[i], Pline->y[i], Pline->z[i]);
84 	}
85 	Vect_append_point(Pout, px, py, pz);
86 	Vect_line_prune(Pout);
87 	if (Pout->n_points > 1) {
88 	    Vect_rewrite_line(Out, line, ltype, Pout, Cline);
89 	    broken++;
90 	}
91 
92 	Vect_reset_line(Pout);
93 	Vect_append_point(Pout, px, py, pz);
94 	for (i = seg; i < Pline->n_points; i++) {
95 	    Vect_append_point(Pout, Pline->x[i], Pline->y[i], Pline->z[i]);
96 	}
97 	Vect_line_prune(Pout);
98 	if (Pout->n_points > 1) {
99 	    if (broken)
100 		Vect_write_line(Out, ltype, Pout, Cline);
101 	    else
102 		Vect_rewrite_line(Out, line, ltype, Pout, Cline);
103 	    broken++;
104 	}
105 	if (broken == 2)
106 	    narcs++;
107 
108 	if (dist > 0.0) {
109 	    if (snap) {
110 		/* snap point */
111 		Points->x[0] = px;
112 		Points->y[0] = py;
113 		Points->z[0] = pz;
114 	    }
115 	    else {
116 		/* write new arc */
117 		Vect_reset_line(Pout);
118 		Vect_append_point(Pout, px, py, pz);
119 		Vect_append_point(Pout, Points->x[0], Points->y[0], Points->z[0]);
120 		maxcat++;
121 		Vect_reset_cats(Cnew);
122 		Vect_cat_set(Cnew, afield, maxcat);
123 		pline = Vect_write_line(Out, ltype, Pout, Cnew);
124 
125 		Vect_list_append(exclude_list, pline);
126 
127 		narcs++;
128 	    }
129 	}
130 
131 	/* add points to 'nfield' layer */
132 	for (i = 0; i < Cats->n_cats; i++) {
133 	    Cats->field[i] = nfield;	/* all points to 'nfield' layer */
134 	}
135 
136 	Vect_write_line(Out, type, Points, Cats);
137     }
138 
139     Vect_destroy_line_struct(Points);
140     Vect_destroy_line_struct(Pline);
141     Vect_destroy_line_struct(Pout);
142     Vect_destroy_cats_struct(Cats);
143     Vect_destroy_cats_struct(Cline);
144     Vect_destroy_cats_struct(Cnew);
145     Vect_destroy_list(exclude_list);
146 
147     return narcs;
148 }
149