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     G_verbose_message(_("Intersections: %5d"), nbreaks);
00116 
00117     for (iline = 0; iline < nlines; iline++) {
00118         if (List_break) {
00119             aline = List_break->value[iline];
00120         }
00121         else {
00122             aline = iline + 1;
00123         }
00124 
00125         if (List_ref && !Vect_val_in_list(List_ref, aline))
00126             continue;
00127 
00128         G_debug(3, "aline =  %d", aline);
00129         if (!Vect_line_alive(Map, aline))
00130             continue;
00131 
00132         atype = Vect_read_line(Map, APoints, ACats, aline);
00133         if (!(atype & type))
00134             continue;
00135 
00136         Vect_get_line_box(Map, aline, &ABox);
00137 
00138         /* Find which sides of the box are touched by intermediate (non-end) points of line */
00139         if (!is3d) {
00140             touch1_n = touch1_s = touch1_e = touch1_w = 0;
00141             for (j = 1; j < APoints->n_points; j++) {
00142                 if (APoints->y[j] == ABox.N)
00143                     touch1_n = 1;
00144                 if (APoints->y[j] == ABox.S)
00145                     touch1_s = 1;
00146                 if (APoints->x[j] == ABox.E)
00147                     touch1_e = 1;
00148                 if (APoints->x[j] == ABox.W)
00149                     touch1_w = 1;
00150             }
00151             G_debug(3, "touch1: n = %d s = %d e = %d w = %d", touch1_n,
00152                     touch1_s, touch1_e, touch1_w);
00153             touch2_n = touch2_s = touch2_e = touch2_w = 0;
00154             for (j = 0; j < APoints->n_points - 1; j++) {
00155                 if (APoints->y[j] == ABox.N)
00156                     touch2_n = 1;
00157                 if (APoints->y[j] == ABox.S)
00158                     touch2_s = 1;
00159                 if (APoints->x[j] == ABox.E)
00160                     touch2_e = 1;
00161                 if (APoints->x[j] == ABox.W)
00162                     touch2_w = 1;
00163             }
00164             G_debug(3, "touch2: n = %d s = %d e = %d w = %d", touch2_n,
00165                     touch2_s, touch2_e, touch2_w);
00166         }
00167 
00168         Vect_select_lines_by_box(Map, &ABox, type, List);
00169         G_debug(3, "  %d lines selected by box", List->n_values);
00170 
00171         for (j = 0; j < List->n_values; j++) {
00172             bline = List->value[j];
00173             if (List_break && !Vect_val_in_list(List_break, bline)) {
00174                 continue;
00175             }
00176             G_debug(3, "  j = %d bline = %d", j, bline);
00177 
00178             /* Check if thouch by end node only */
00179             if (!is3d) {
00180                 Vect_get_line_nodes(Map, aline, &anode1, &anode2);
00181                 Vect_get_line_nodes(Map, bline, &bnode1, &bnode2);
00182                 Vect_get_line_box(Map, bline, &BBox);
00183 
00184                 if (anode1 == bnode1 || anode1 == bnode2)
00185                     node = anode1;
00186                 else if (anode2 == bnode1 || anode2 == bnode2)
00187                     node = anode2;
00188                 else
00189                     node = 0;
00190 
00191                 if (node) {
00192                     Vect_get_node_coor(Map, node, &nodex, &nodey, NULL);
00193                     if ((node == anode1 && nodey == ABox.N && !touch1_n &&
00194                          nodey == BBox.S) || (node == anode2 &&
00195                                               nodey == ABox.N && !touch2_n &&
00196                                               nodey == BBox.S) ||
00197                         (node == anode1 && nodey == ABox.S && !touch1_s &&
00198                          nodey == BBox.N) || (node == anode2 &&
00199                                               nodey == ABox.S && !touch2_s &&
00200                                               nodey == BBox.N) ||
00201                         (node == anode1 && nodex == ABox.E && !touch1_e &&
00202                          nodex == BBox.W) || (node == anode2 &&
00203                                               nodex == ABox.E && !touch2_e &&
00204                                               nodex == BBox.W) ||
00205                         (node == anode1 && nodex == ABox.W && !touch1_w &&
00206                          nodex == BBox.E) || (node == anode2 &&
00207                                               nodex == ABox.W && !touch2_w &&
00208                                               nodex == BBox.E)) {
00209                         G_debug(3,
00210                                 "lines %d and %d touching by end nodes only -> no intersection",
00211                                 aline, bline);
00212                         continue;
00213                     }
00214                 }
00215             }
00216 
00217             btype = Vect_read_line(Map, BPoints, BCats, bline);
00218 
00219             Vect_line_intersection(APoints, BPoints, &AXLines, &BXLines,
00220                                    &naxlines, &nbxlines, 0);
00221             G_debug(3, "  naxlines = %d nbxlines = %d", naxlines, nbxlines);
00222 
00223             /* This part handles a special case when aline == bline, no other intersection was found
00224              * and the line is forming collapsed loop, for example  0,0;1,0;0,0 should be broken at 1,0.
00225              * ---> */
00226             if (aline == bline && naxlines == 0 && nbxlines == 0 &&
00227                 APoints->n_points >= 3) {
00228                 int centre;
00229                 int i;
00230 
00231                 G_debug(3, "  Check collapsed loop");
00232                 if (APoints->n_points % 2) {    /* odd number of vertices */
00233                     centre = APoints->n_points / 2;     /* index of centre */
00234                     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 */
00235                         AXLines =
00236                             (struct line_pnts **)G_malloc(2 *
00237                                                           sizeof(struct
00238                                                                  line_pnts
00239                                                                  *));
00240                         AXLines[0] = Vect_new_line_struct();
00241                         AXLines[1] = Vect_new_line_struct();
00242 
00243                         for (i = 0; i <= centre; i++)
00244                             Vect_append_point(AXLines[0], APoints->x[i],
00245                                               APoints->y[i], APoints->z[i]);
00246 
00247                         for (i = centre; i < APoints->n_points; i++)
00248                             Vect_append_point(AXLines[1], APoints->x[i],
00249                                               APoints->y[i], APoints->z[i]);
00250 
00251                         naxlines = 2;
00252                     }
00253                 }
00254             }
00255             /* <--- */
00256 
00257             if (Err) {          /* array for intersections (more than needed */
00258                 xx = (double *)G_malloc((naxlines + nbxlines) *
00259                                         sizeof(double));
00260                 yx = (double *)G_malloc((naxlines + nbxlines) *
00261                                         sizeof(double));
00262                 zx = (double *)G_malloc((naxlines + nbxlines) *
00263                                         sizeof(double));
00264             }
00265             nx = 0;             /* number of intersections to be written to Err */
00266             if (naxlines > 0) { /* intersection -> write out */
00267                 Vect_delete_line(Map, aline);
00268                 for (k = 0; k < naxlines; k++) {
00269                     /* Write new line segments */
00270                     /* line may collapse, don't write zero length lines */
00271                     Vect_line_prune(AXLines[k]);
00272                     if ((atype & GV_POINTS) || AXLines[k]->n_points > 1) {
00273                         ret = Vect_write_line(Map, atype, AXLines[k], ACats);
00274                         if (List_ref) {
00275                             Vect_list_append(List_ref, ret);
00276                         }
00277                         G_debug(3, "Line %d written, npoints = %d", ret,
00278                                 AXLines[k]->n_points);
00279                         if (List_break) {
00280                             Vect_list_append(List_break, ret);
00281                         }
00282                     }
00283 
00284                     /* Write intersection points */
00285                     if (Err) {
00286                         if (k > 0) {
00287                             xx[nx] = AXLines[k]->x[0];
00288                             yx[nx] = AXLines[k]->y[0];
00289                             zx[nx] = AXLines[k]->z[0];
00290                             nx++;
00291                         }
00292                     }
00293                     Vect_destroy_line_struct(AXLines[k]);
00294                 }
00295                 G_free(AXLines);
00296                 nbreaks += naxlines - 1;
00297             }
00298 
00299             if (nbxlines > 0) {
00300                 if (aline != bline) {   /* Self intersection, do not write twice, TODO: is it OK? */
00301                     Vect_delete_line(Map, bline);
00302                     for (k = 0; k < nbxlines; k++) {
00303                         /* Write new line segments */
00304                         /* line may collapse, don't write zero length lines */
00305                         Vect_line_prune(BXLines[k]);
00306                         if ((btype & GV_POINTS) || BXLines[k]->n_points > 1) {
00307                             ret =
00308                                 Vect_write_line(Map, btype, BXLines[k],
00309                                                 BCats);
00310                             G_debug(5, "Line %d written", ret);
00311                             if (List_break) {
00312                                 Vect_list_append(List_break, ret);
00313                             }
00314                         }
00315 
00316                         /* Write intersection points */
00317                         if (Err) {
00318                             if (k > 0) {
00319                                 found = 0;
00320                                 for (l = 0; l < nx; l++) {
00321                                     if (xx[l] == BXLines[k]->x[0] &&
00322                                         yx[l] == BXLines[k]->y[0] &&
00323                                         zx[l] == BXLines[k]->z[0]) {
00324                                         found = 1;
00325                                         break;
00326                                     }
00327                                 }
00328                                 if (!found) {
00329                                     xx[nx] = BXLines[k]->x[0];
00330                                     yx[nx] = BXLines[k]->y[0];
00331                                     zx[nx] = BXLines[k]->z[0];
00332                                     nx++;
00333                                 }
00334                             }
00335                         }
00336                         Vect_destroy_line_struct(BXLines[k]);
00337                     }
00338                     nbreaks += nbxlines - 1;
00339                 }
00340                 else {
00341                     for (k = 0; k < nbxlines; k++)
00342                         Vect_destroy_line_struct(BXLines[k]);
00343                 }
00344                 G_free(BXLines);
00345             }
00346             if (Err) {
00347                 for (l = 0; l < nx; l++) {      /* Write out errors */
00348                     Vect_reset_line(Points);
00349                     Vect_append_point(Points, xx[l], yx[l], zx[l]);
00350                     ret = Vect_write_line(Err, GV_POINT, Points, Cats);
00351                 }
00352 
00353                 G_free(xx);
00354                 G_free(yx);
00355                 G_free(zx);
00356             }
00357             if (naxlines > 0)
00358                 break;          /* first line was broken and deleted -> take the next one */
00359         }
00360 
00361         if (List_break) {
00362             nlines = List_break->n_values;
00363         }
00364         else {
00365             nlines = Vect_get_num_lines(Map);
00366         }
00367         G_debug(3, "nlines =  %d", nlines);
00368     }                           /* for each line */
00369 
00370     G_verbose_message(_("Intersections: %5d"), nbreaks);
00371 
00372     Vect_destroy_list(List);
00373 
00374     return nbreaks;
00375 }