GRASS Programmer's Manual
6.4.2(2012)
|
00001 00018 #include <stdlib.h> 00019 #include <grass/gis.h> 00020 #include <grass/Vect.h> 00021 #include <grass/glocale.h> 00022 00035 void 00036 Vect_break_lines(struct Map_info *Map, int type, struct Map_info *Err) 00037 { 00038 Vect_break_lines_list(Map, NULL, NULL, type, Err); 00039 00040 return; 00041 } 00042 00065 int 00066 Vect_break_lines_list(struct Map_info *Map, struct ilist *List_break, 00067 struct ilist *List_ref, int type, struct Map_info *Err) 00068 { 00069 struct line_pnts *APoints, *BPoints, *Points; 00070 struct line_pnts **AXLines, **BXLines; 00071 struct line_cats *ACats, *BCats, *Cats; 00072 int j, k, l, ret, atype, btype, aline, bline, found, iline, nlines; 00073 int naxlines, nbxlines, nx; 00074 double *xx = NULL, *yx = NULL, *zx = NULL; 00075 BOUND_BOX ABox, BBox; 00076 struct ilist *List; 00077 int nbreaks; 00078 int touch1_n = 0, touch1_s = 0, touch1_e = 0, touch1_w = 0; /* other vertices except node1 touching box */ 00079 int touch2_n = 0, touch2_s = 0, touch2_e = 0, touch2_w = 0; /* other vertices except node2 touching box */ 00080 int is3d; 00081 int node, anode1, anode2, bnode1, bnode2; 00082 double nodex, nodey; 00083 00084 APoints = Vect_new_line_struct(); 00085 BPoints = Vect_new_line_struct(); 00086 Points = Vect_new_line_struct(); 00087 ACats = Vect_new_cats_struct(); 00088 BCats = Vect_new_cats_struct(); 00089 Cats = Vect_new_cats_struct(); 00090 List = Vect_new_list(); 00091 00092 is3d = Vect_is_3d(Map); 00093 00094 if (List_break) { 00095 nlines = List_break->n_values; 00096 } 00097 else { 00098 nlines = Vect_get_num_lines(Map); 00099 } 00100 G_debug(3, "nlines = %d", nlines); 00101 00102 /* To find intersection of two lines (Vect_line_intersection) is quite slow. 00103 * Fortunately usual lines/boundaries in GIS often forms a network where lines 00104 * are connected by end points, and touch by MBR. This function checks and occasionaly 00105 * skips such cases. This is currently done for 2D only 00106 */ 00107 00108 /* Go through all lines in vector, for each select lines which overlap MBR of 00109 * this line exclude those connected by one endpoint (see above) 00110 * and try to intersect, if lines intersect write new lines at the end of 00111 * the file, and process next line (remaining lines overlapping box are skipped) 00112 */ 00113 nbreaks = 0; 00114 00115 for (iline = 0; iline < nlines; iline++) { 00116 G_percent(iline, nlines, 1); 00117 if (List_break) { 00118 aline = List_break->value[iline]; 00119 } 00120 else { 00121 aline = iline + 1; 00122 } 00123 00124 if (List_ref && !Vect_val_in_list(List_ref, aline)) 00125 continue; 00126 00127 G_debug(3, "aline = %d", aline); 00128 if (!Vect_line_alive(Map, aline)) 00129 continue; 00130 00131 atype = Vect_read_line(Map, APoints, ACats, aline); 00132 if (!(atype & type)) 00133 continue; 00134 00135 Vect_get_line_box(Map, aline, &ABox); 00136 00137 /* Find which sides of the box are touched by intermediate (non-end) points of line */ 00138 if (!is3d) { 00139 touch1_n = touch1_s = touch1_e = touch1_w = 0; 00140 for (j = 1; j < APoints->n_points; j++) { 00141 if (APoints->y[j] == ABox.N) 00142 touch1_n = 1; 00143 if (APoints->y[j] == ABox.S) 00144 touch1_s = 1; 00145 if (APoints->x[j] == ABox.E) 00146 touch1_e = 1; 00147 if (APoints->x[j] == ABox.W) 00148 touch1_w = 1; 00149 } 00150 G_debug(3, "touch1: n = %d s = %d e = %d w = %d", touch1_n, 00151 touch1_s, touch1_e, touch1_w); 00152 touch2_n = touch2_s = touch2_e = touch2_w = 0; 00153 for (j = 0; j < APoints->n_points - 1; j++) { 00154 if (APoints->y[j] == ABox.N) 00155 touch2_n = 1; 00156 if (APoints->y[j] == ABox.S) 00157 touch2_s = 1; 00158 if (APoints->x[j] == ABox.E) 00159 touch2_e = 1; 00160 if (APoints->x[j] == ABox.W) 00161 touch2_w = 1; 00162 } 00163 G_debug(3, "touch2: n = %d s = %d e = %d w = %d", touch2_n, 00164 touch2_s, touch2_e, touch2_w); 00165 } 00166 00167 Vect_select_lines_by_box(Map, &ABox, type, List); 00168 G_debug(3, " %d lines selected by box", List->n_values); 00169 00170 for (j = 0; j < List->n_values; j++) { 00171 bline = List->value[j]; 00172 if (List_break && !Vect_val_in_list(List_break, bline)) { 00173 continue; 00174 } 00175 G_debug(3, " j = %d bline = %d", j, bline); 00176 00177 /* Check if thouch by end node only */ 00178 if (!is3d) { 00179 Vect_get_line_nodes(Map, aline, &anode1, &anode2); 00180 Vect_get_line_nodes(Map, bline, &bnode1, &bnode2); 00181 Vect_get_line_box(Map, bline, &BBox); 00182 00183 if (anode1 == bnode1 || anode1 == bnode2) 00184 node = anode1; 00185 else if (anode2 == bnode1 || anode2 == bnode2) 00186 node = anode2; 00187 else 00188 node = 0; 00189 00190 if (node) { 00191 Vect_get_node_coor(Map, node, &nodex, &nodey, NULL); 00192 if ((node == anode1 && nodey == ABox.N && !touch1_n && 00193 nodey == BBox.S) || (node == anode2 && 00194 nodey == ABox.N && !touch2_n && 00195 nodey == BBox.S) || 00196 (node == anode1 && nodey == ABox.S && !touch1_s && 00197 nodey == BBox.N) || (node == anode2 && 00198 nodey == ABox.S && !touch2_s && 00199 nodey == BBox.N) || 00200 (node == anode1 && nodex == ABox.E && !touch1_e && 00201 nodex == BBox.W) || (node == anode2 && 00202 nodex == ABox.E && !touch2_e && 00203 nodex == BBox.W) || 00204 (node == anode1 && nodex == ABox.W && !touch1_w && 00205 nodex == BBox.E) || (node == anode2 && 00206 nodex == ABox.W && !touch2_w && 00207 nodex == BBox.E)) { 00208 G_debug(3, 00209 "lines %d and %d touching by end nodes only -> no intersection", 00210 aline, bline); 00211 continue; 00212 } 00213 } 00214 } 00215 00216 btype = Vect_read_line(Map, BPoints, BCats, bline); 00217 00218 AXLines = NULL; 00219 BXLines = NULL; 00220 Vect_line_intersection(APoints, BPoints, &AXLines, &BXLines, 00221 &naxlines, &nbxlines, 0); 00222 G_debug(3, " naxlines = %d nbxlines = %d", naxlines, nbxlines); 00223 00224 /* This part handles a special case when aline == bline, no other intersection was found 00225 * and the line is forming collapsed loop, for example 0,0;1,0;0,0 should be broken at 1,0. 00226 * ---> */ 00227 if (aline == bline && naxlines == 0 && nbxlines == 0 && 00228 APoints->n_points >= 3) { 00229 int centre; 00230 int i; 00231 00232 G_debug(3, " Check collapsed loop"); 00233 if (APoints->n_points % 2) { /* odd number of vertices */ 00234 centre = APoints->n_points / 2; /* index of centre */ 00235 if (APoints->x[centre - 1] == APoints->x[centre + 1] && APoints->y[centre - 1] == APoints->y[centre + 1] && APoints->z[centre - 1] == APoints->z[centre + 1]) { /* -> break */ 00236 AXLines = 00237 (struct line_pnts **)G_malloc(2 * 00238 sizeof(struct 00239 line_pnts 00240 *)); 00241 AXLines[0] = Vect_new_line_struct(); 00242 AXLines[1] = Vect_new_line_struct(); 00243 00244 for (i = 0; i <= centre; i++) 00245 Vect_append_point(AXLines[0], APoints->x[i], 00246 APoints->y[i], APoints->z[i]); 00247 00248 for (i = centre; i < APoints->n_points; i++) 00249 Vect_append_point(AXLines[1], APoints->x[i], 00250 APoints->y[i], APoints->z[i]); 00251 00252 naxlines = 2; 00253 } 00254 } 00255 } 00256 /* <--- */ 00257 00258 if (Err) { /* array for intersections (more than needed */ 00259 xx = (double *)G_malloc((naxlines + nbxlines) * 00260 sizeof(double)); 00261 yx = (double *)G_malloc((naxlines + nbxlines) * 00262 sizeof(double)); 00263 zx = (double *)G_malloc((naxlines + nbxlines) * 00264 sizeof(double)); 00265 } 00266 nx = 0; /* number of intersections to be written to Err */ 00267 if (naxlines > 0) { /* intersection -> write out */ 00268 Vect_delete_line(Map, aline); 00269 for (k = 0; k < naxlines; k++) { 00270 /* Write new line segments */ 00271 /* line may collapse, don't write zero length lines */ 00272 Vect_line_prune(AXLines[k]); 00273 if ((atype & GV_POINTS) || AXLines[k]->n_points > 1) { 00274 ret = Vect_write_line(Map, atype, AXLines[k], ACats); 00275 if (List_ref) { 00276 Vect_list_append(List_ref, ret); 00277 } 00278 G_debug(3, "Line %d written, npoints = %d", ret, 00279 AXLines[k]->n_points); 00280 if (List_break) { 00281 Vect_list_append(List_break, ret); 00282 } 00283 } 00284 00285 /* Write intersection points */ 00286 if (Err) { 00287 if (k > 0) { 00288 xx[nx] = AXLines[k]->x[0]; 00289 yx[nx] = AXLines[k]->y[0]; 00290 zx[nx] = AXLines[k]->z[0]; 00291 nx++; 00292 } 00293 } 00294 Vect_destroy_line_struct(AXLines[k]); 00295 } 00296 nbreaks += naxlines - 1; 00297 } 00298 if (AXLines) 00299 G_free(AXLines); 00300 00301 if (nbxlines > 0) { 00302 if (aline != bline) { /* Self intersection, do not write twice, TODO: is it OK? */ 00303 Vect_delete_line(Map, bline); 00304 for (k = 0; k < nbxlines; k++) { 00305 /* Write new line segments */ 00306 /* line may collapse, don't write zero length lines */ 00307 Vect_line_prune(BXLines[k]); 00308 if ((btype & GV_POINTS) || BXLines[k]->n_points > 1) { 00309 ret = 00310 Vect_write_line(Map, btype, BXLines[k], 00311 BCats); 00312 G_debug(5, "Line %d written", ret); 00313 if (List_break) { 00314 Vect_list_append(List_break, ret); 00315 } 00316 } 00317 00318 /* Write intersection points */ 00319 if (Err) { 00320 if (k > 0) { 00321 found = 0; 00322 for (l = 0; l < nx; l++) { 00323 if (xx[l] == BXLines[k]->x[0] && 00324 yx[l] == BXLines[k]->y[0] && 00325 zx[l] == BXLines[k]->z[0]) { 00326 found = 1; 00327 break; 00328 } 00329 } 00330 if (!found) { 00331 xx[nx] = BXLines[k]->x[0]; 00332 yx[nx] = BXLines[k]->y[0]; 00333 zx[nx] = BXLines[k]->z[0]; 00334 nx++; 00335 } 00336 } 00337 } 00338 Vect_destroy_line_struct(BXLines[k]); 00339 } 00340 nbreaks += nbxlines - 1; 00341 } 00342 else { 00343 for (k = 0; k < nbxlines; k++) 00344 Vect_destroy_line_struct(BXLines[k]); 00345 } 00346 } 00347 if (BXLines) 00348 G_free(BXLines); 00349 if (Err) { 00350 for (l = 0; l < nx; l++) { /* Write out errors */ 00351 Vect_reset_line(Points); 00352 Vect_append_point(Points, xx[l], yx[l], zx[l]); 00353 ret = Vect_write_line(Err, GV_POINT, Points, Cats); 00354 } 00355 00356 G_free(xx); 00357 G_free(yx); 00358 G_free(zx); 00359 } 00360 if (naxlines > 0) 00361 break; /* first line was broken and deleted -> take the next one */ 00362 } 00363 00364 if (List_break) { 00365 nlines = List_break->n_values; 00366 } 00367 else { 00368 nlines = Vect_get_num_lines(Map); 00369 } 00370 G_debug(3, "nlines = %d", nlines); 00371 } /* for each line */ 00372 G_percent(nlines, nlines, 1); /* finish it */ 00373 00374 G_verbose_message(_("Intersections: %5d"), nbreaks); 00375 00376 Vect_destroy_list(List); 00377 00378 return nbreaks; 00379 }