GRASS Programmer's Manual 6.4.1(2011)
break_lines.c
Go to the documentation of this file.
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 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines