GRASS Programmer's Manual 6.4.1(2011)
|
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 Vect_line_intersection(APoints, BPoints, &AXLines, &BXLines, 00219 &naxlines, &nbxlines, 0); 00220 G_debug(3, " naxlines = %d nbxlines = %d", naxlines, nbxlines); 00221 00222 /* This part handles a special case when aline == bline, no other intersection was found 00223 * and the line is forming collapsed loop, for example 0,0;1,0;0,0 should be broken at 1,0. 00224 * ---> */ 00225 if (aline == bline && naxlines == 0 && nbxlines == 0 && 00226 APoints->n_points >= 3) { 00227 int centre; 00228 int i; 00229 00230 G_debug(3, " Check collapsed loop"); 00231 if (APoints->n_points % 2) { /* odd number of vertices */ 00232 centre = APoints->n_points / 2; /* index of centre */ 00233 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 */ 00234 AXLines = 00235 (struct line_pnts **)G_malloc(2 * 00236 sizeof(struct 00237 line_pnts 00238 *)); 00239 AXLines[0] = Vect_new_line_struct(); 00240 AXLines[1] = Vect_new_line_struct(); 00241 00242 for (i = 0; i <= centre; i++) 00243 Vect_append_point(AXLines[0], APoints->x[i], 00244 APoints->y[i], APoints->z[i]); 00245 00246 for (i = centre; i < APoints->n_points; i++) 00247 Vect_append_point(AXLines[1], APoints->x[i], 00248 APoints->y[i], APoints->z[i]); 00249 00250 naxlines = 2; 00251 } 00252 } 00253 } 00254 /* <--- */ 00255 00256 if (Err) { /* array for intersections (more than needed */ 00257 xx = (double *)G_malloc((naxlines + nbxlines) * 00258 sizeof(double)); 00259 yx = (double *)G_malloc((naxlines + nbxlines) * 00260 sizeof(double)); 00261 zx = (double *)G_malloc((naxlines + nbxlines) * 00262 sizeof(double)); 00263 } 00264 nx = 0; /* number of intersections to be written to Err */ 00265 if (naxlines > 0) { /* intersection -> write out */ 00266 Vect_delete_line(Map, aline); 00267 for (k = 0; k < naxlines; k++) { 00268 /* Write new line segments */ 00269 /* line may collapse, don't write zero length lines */ 00270 Vect_line_prune(AXLines[k]); 00271 if ((atype & GV_POINTS) || AXLines[k]->n_points > 1) { 00272 ret = Vect_write_line(Map, atype, AXLines[k], ACats); 00273 if (List_ref) { 00274 Vect_list_append(List_ref, ret); 00275 } 00276 G_debug(3, "Line %d written, npoints = %d", ret, 00277 AXLines[k]->n_points); 00278 if (List_break) { 00279 Vect_list_append(List_break, ret); 00280 } 00281 } 00282 00283 /* Write intersection points */ 00284 if (Err) { 00285 if (k > 0) { 00286 xx[nx] = AXLines[k]->x[0]; 00287 yx[nx] = AXLines[k]->y[0]; 00288 zx[nx] = AXLines[k]->z[0]; 00289 nx++; 00290 } 00291 } 00292 Vect_destroy_line_struct(AXLines[k]); 00293 } 00294 G_free(AXLines); 00295 nbreaks += naxlines - 1; 00296 } 00297 00298 if (nbxlines > 0) { 00299 if (aline != bline) { /* Self intersection, do not write twice, TODO: is it OK? */ 00300 Vect_delete_line(Map, bline); 00301 for (k = 0; k < nbxlines; k++) { 00302 /* Write new line segments */ 00303 /* line may collapse, don't write zero length lines */ 00304 Vect_line_prune(BXLines[k]); 00305 if ((btype & GV_POINTS) || BXLines[k]->n_points > 1) { 00306 ret = 00307 Vect_write_line(Map, btype, BXLines[k], 00308 BCats); 00309 G_debug(5, "Line %d written", ret); 00310 if (List_break) { 00311 Vect_list_append(List_break, ret); 00312 } 00313 } 00314 00315 /* Write intersection points */ 00316 if (Err) { 00317 if (k > 0) { 00318 found = 0; 00319 for (l = 0; l < nx; l++) { 00320 if (xx[l] == BXLines[k]->x[0] && 00321 yx[l] == BXLines[k]->y[0] && 00322 zx[l] == BXLines[k]->z[0]) { 00323 found = 1; 00324 break; 00325 } 00326 } 00327 if (!found) { 00328 xx[nx] = BXLines[k]->x[0]; 00329 yx[nx] = BXLines[k]->y[0]; 00330 zx[nx] = BXLines[k]->z[0]; 00331 nx++; 00332 } 00333 } 00334 } 00335 Vect_destroy_line_struct(BXLines[k]); 00336 } 00337 nbreaks += nbxlines - 1; 00338 } 00339 else { 00340 for (k = 0; k < nbxlines; k++) 00341 Vect_destroy_line_struct(BXLines[k]); 00342 } 00343 G_free(BXLines); 00344 } 00345 if (Err) { 00346 for (l = 0; l < nx; l++) { /* Write out errors */ 00347 Vect_reset_line(Points); 00348 Vect_append_point(Points, xx[l], yx[l], zx[l]); 00349 ret = Vect_write_line(Err, GV_POINT, Points, Cats); 00350 } 00351 00352 G_free(xx); 00353 G_free(yx); 00354 G_free(zx); 00355 } 00356 if (naxlines > 0) 00357 break; /* first line was broken and deleted -> take the next one */ 00358 } 00359 00360 if (List_break) { 00361 nlines = List_break->n_values; 00362 } 00363 else { 00364 nlines = Vect_get_num_lines(Map); 00365 } 00366 G_debug(3, "nlines = %d", nlines); 00367 } /* for each line */ 00368 G_percent(nlines, nlines, 1); /* finish it */ 00369 00370 G_verbose_message(_("Intersections: %5d"), nbreaks); 00371 00372 Vect_destroy_list(List); 00373 00374 return nbreaks; 00375 }