GRASS Programmer's Manual  6.4.2(2012)
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             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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines