break_lines.c

Go to the documentation of this file.
00001 /***************************************************************
00002  *
00003  * MODULE:       vector library
00004  * 
00005  * AUTHOR(S):    Radim Blazek
00006  *               
00007  * PURPOSE:      Clean lines
00008  *               
00009  * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00010  *
00011  *               This program is free software under the 
00012  *               GNU General Public License (>=v2). 
00013  *               Read the file COPYING that comes with GRASS
00014  *               for details.
00015  *
00016  **************************************************************/
00017 #include <stdlib.h> 
00018 #include <grass/gis.h>
00019 #include <grass/Vect.h>
00020 
00036 void 
00037 Vect_break_lines ( struct Map_info *Map, int type, struct Map_info *Err, FILE *msgout )
00038 {
00039     struct line_pnts *APoints, *BPoints, *Points;
00040     struct line_pnts **AXLines, **BXLines;
00041     struct line_cats *ACats, *BCats, *Cats;
00042     int    j, k, l, ret, atype, btype, aline, bline, found;
00043     int    nlines, naxlines, nbxlines, nx;
00044     double *xx=NULL, *yx=NULL, *zx=NULL;
00045     BOUND_BOX  ABox, BBox; 
00046     struct ilist *List; 
00047     int nbreaks;
00048     int touch1_n=0, touch1_s=0, touch1_e=0, touch1_w=0; /* other vertices except node1 touching box */ 
00049     int touch2_n=0, touch2_s=0, touch2_e=0, touch2_w=0; /* other vertices except node2 touching box */ 
00050     int is3d;
00051     int node, anode1, anode2, bnode1, bnode2;
00052     double nodex, nodey;
00053     int printed;
00054 
00055     APoints = Vect_new_line_struct ();
00056     BPoints = Vect_new_line_struct ();
00057     Points = Vect_new_line_struct ();
00058     ACats = Vect_new_cats_struct ();
00059     BCats = Vect_new_cats_struct ();
00060     Cats = Vect_new_cats_struct ();
00061     List = Vect_new_list ();
00062     
00063     nlines = Vect_get_num_lines (Map);
00064     is3d = Vect_is_3d ( Map );
00065 
00066     G_debug (3, "nlines =  %d", nlines );
00067     /* To find intersection of two lines (Vect_line_intersection) is quite slow.
00068      * Fortunately usual lines/boundaries in GIS often forms a network where lines
00069      * are connected by end points, and touch by MBR. This function checks and occasionaly
00070      * skips such cases. This is currently done for 2D only */
00071     
00072     /* Go through all lines in vector, for each select lines which overlap MBR of
00073     *  this line exclude those connected by one endpoint (see above)
00074     *  and try to intersect, if lines intersect write new lines at the end of 
00075     *  the file, and process next line (remainining lines overlaping box are skipped) */
00076     nbreaks = 0;
00077     printed = 0;
00078     if (msgout) fprintf (msgout, "Intersections: %5d", nbreaks ); 
00079     for ( aline = 1; aline <= nlines; aline++ ){ 
00080         G_debug (3, "aline =  %d", aline);
00081         if ( !Vect_line_alive ( Map, aline ) ) continue;
00082 
00083         atype = Vect_read_line (Map, APoints, ACats, aline);
00084         if ( !(atype & type) ) continue;
00085 
00086         Vect_get_line_box ( Map, aline, &ABox );
00087 
00088         /* Find which sides of the box are touched by intermediate (non-end) points of line */
00089         if ( !is3d ) {
00090             touch1_n = touch1_s = touch1_e = touch1_w = 0;
00091             for ( j = 1; j < APoints->n_points; j++ ){
00092                 if ( APoints->y[j] == ABox.N ) touch1_n = 1;
00093                 if ( APoints->y[j] == ABox.S ) touch1_s = 1;
00094                 if ( APoints->x[j] == ABox.E ) touch1_e = 1;  
00095                 if ( APoints->x[j] == ABox.W ) touch1_w = 1;  
00096             }
00097             G_debug(3, "touch1: n = %d s = %d e = %d w = %d", touch1_n, touch1_s, touch1_e, touch1_w );
00098             touch2_n = touch2_s = touch2_e = touch2_w = 0;
00099             for ( j = 0; j < APoints->n_points-1; j++ ){
00100                 if ( APoints->y[j] == ABox.N ) touch2_n = 1;
00101                 if ( APoints->y[j] == ABox.S ) touch2_s = 1;
00102                 if ( APoints->x[j] == ABox.E ) touch2_e = 1;  
00103                 if ( APoints->x[j] == ABox.W ) touch2_w = 1;  
00104             }
00105             G_debug(3, "touch2: n = %d s = %d e = %d w = %d", touch2_n, touch2_s, touch2_e, touch2_w );
00106         }
00107         
00108         Vect_select_lines_by_box ( Map, &ABox, type, List);
00109         G_debug (3, "  %d lines selected by box", List->n_values);
00110         
00111         for ( j = 0; j <  List->n_values; j++ ){ 
00112             bline = List->value[j];
00113             G_debug (3, "  j = %d bline = %d", j, bline);
00114     
00115             /* Check if thouch by end node only */
00116             if ( !is3d ) {
00117                 Vect_get_line_nodes ( Map, aline, &anode1, &anode2 );
00118                 Vect_get_line_nodes ( Map, bline, &bnode1, &bnode2 );
00119                 Vect_get_line_box ( Map, bline, &BBox );
00120 
00121                 if ( anode1 == bnode1 || anode1 == bnode2 ) 
00122                     node = anode1; 
00123                 else if ( anode2 == bnode1 || anode2 == bnode2 )
00124                     node = anode2;
00125                 else
00126                     node = 0;
00127                 
00128                 if ( node ) {
00129                     Vect_get_node_coor ( Map, node, &nodex, &nodey, NULL);
00130                     if ( ( node == anode1 && nodey == ABox.N && !touch1_n && nodey == BBox.S  ) ||
00131                          ( node == anode2 && nodey == ABox.N && !touch2_n && nodey == BBox.S  ) ||
00132                          ( node == anode1 && nodey == ABox.S && !touch1_s && nodey == BBox.N ) ||
00133                          ( node == anode2 && nodey == ABox.S && !touch2_s && nodey == BBox.N ) ||
00134                          ( node == anode1 && nodex == ABox.E && !touch1_e && nodex == BBox.W ) ||
00135                          ( node == anode2 && nodex == ABox.E && !touch2_e && nodex == BBox.W ) ||
00136                          ( node == anode1 && nodex == ABox.W && !touch1_w && nodex == BBox.E ) ||
00137                          ( node == anode2 && nodex == ABox.W && !touch2_w && nodex == BBox.E ) )
00138                     {
00139                         G_debug(3, "lines %d and %d touching by end nodes only -> no intersection", aline, bline);
00140                         continue;
00141                     }
00142                 }
00143             }
00144                 
00145             btype = Vect_read_line (Map, BPoints, BCats, bline);
00146     
00147             Vect_line_intersection(APoints, BPoints, &AXLines, &BXLines, &naxlines, &nbxlines, 0);
00148             G_debug(3, "  naxlines = %d nbxlines = %d",  naxlines, nbxlines);
00149 
00150             /* This part handles a special case when aline == bline, no other intersection was found
00151              * and the line is forming collapsed loop, for example  0,0;1,0;0,0 should be broken at 1,0.
00152              * ---> */
00153             if ( aline == bline && naxlines == 0 && nbxlines == 0 && APoints->n_points >= 3  ) {
00154                 int centre;
00155                 int i;
00156                 
00157                 G_debug(3, "  Check collapsed loop" );
00158                 if ( APoints->n_points % 2 ) { /* odd number of vertices */
00159                     centre = APoints->n_points / 2; /* index of centre */
00160                     if ( APoints->x[centre-1] == APoints->x[centre+1] &&
00161                          APoints->y[centre-1] == APoints->y[centre+1] &&
00162                          APoints->z[centre-1] == APoints->z[centre+1] ) /* -> break */
00163                     {
00164                         AXLines = (struct line_pnts **) G_malloc ( 2 * sizeof(struct line_pnts *) );
00165                         AXLines[0] = Vect_new_line_struct ();
00166                         AXLines[1] = Vect_new_line_struct ();
00167                         
00168                         for ( i = 0; i <= centre; i++ ) 
00169                             Vect_append_point ( AXLines[0], APoints->x[i], APoints->y[i], APoints->z[i]);
00170 
00171                         for ( i = centre; i < APoints->n_points; i++ ) 
00172                             Vect_append_point ( AXLines[1], APoints->x[i], APoints->y[i], APoints->z[i]);
00173 
00174                         naxlines = 2;
00175                     }
00176                 }
00177             }
00178             /* <--- */
00179             
00180             if ( Err ) { /* array for intersections (more than needed */
00181                 xx = (double *) G_malloc ( (naxlines + nbxlines) * sizeof ( double ) ); 
00182                 yx = (double *) G_malloc ( (naxlines + nbxlines) * sizeof ( double ) ); 
00183                 zx = (double *) G_malloc ( (naxlines + nbxlines) * sizeof ( double ) ); 
00184             }
00185             nx = 0; /* number of intersections to be written to Err */
00186             if ( naxlines > 0 ) { /* intersection -> write out */
00187                 Vect_delete_line (Map, aline); 
00188                 for ( k = 0; k < naxlines; k++ ){ 
00189                     /* Write new line segments */
00190                     /* line may collapse, don't write zero length lines */
00191                     Vect_line_prune ( AXLines[k] );
00192                     if ( (atype & GV_POINTS) || AXLines[k]->n_points > 1 ) {
00193                         ret = Vect_write_line ( Map, atype, AXLines[k], ACats );  
00194                         G_debug (3, "Line %d written, npoints = %d", ret, AXLines[k]->n_points);
00195                     }
00196                     
00197                     /* Write intersection points */
00198                     if ( Err ) {
00199                         if ( k > 0 ) {
00200                             xx[nx] = AXLines[k]->x[0];
00201                             yx[nx] = AXLines[k]->y[0];
00202                             zx[nx] = AXLines[k]->z[0];
00203                             nx++;
00204                         }
00205                     }   
00206                     Vect_destroy_line_struct (  AXLines[k] );
00207                 }
00208                 G_free ( AXLines );
00209                 nbreaks += naxlines - 1;
00210             }
00211                 
00212             if ( nbxlines > 0 ) { 
00213                 if ( aline != bline ) { /* Self intersection, do not write twice, TODO: is it OK? */
00214                     Vect_delete_line (Map, bline); 
00215                     for ( k = 0; k < nbxlines; k++ ){ 
00216                         /* Write new line segments */
00217                         /* line may collapse, don't write zero length lines */
00218                         Vect_line_prune ( BXLines[k] );
00219                         if ( (btype & GV_POINTS) || BXLines[k]->n_points > 1 ) {
00220                             ret = Vect_write_line ( Map, btype, BXLines[k], BCats );  
00221                             G_debug (5, "Line %d written", ret);
00222                         }
00223                         
00224                         /* Write intersection points */
00225                         if ( Err ) {
00226                             if ( k > 0 ) {
00227                                 found = 0;
00228                                 for ( l = 0; l < nx; l++ ) {
00229                                     if ( xx[l] == BXLines[k]->x[0] && yx[l] == BXLines[k]->y[0] &&
00230                                          zx[l] == BXLines[k]->z[0] )
00231                                     {
00232                                         found = 1;
00233                                         break;
00234                                     }
00235                                 }
00236                                 if ( !found ) {
00237                                     xx[nx] = BXLines[k]->x[0];
00238                                     yx[nx] = BXLines[k]->y[0];
00239                                     zx[nx] = BXLines[k]->z[0];
00240                                     nx++;
00241                                 }
00242                             }
00243                         }       
00244                         Vect_destroy_line_struct (  BXLines[k] );
00245                     }
00246                     nbreaks += nbxlines - 1;
00247                 } else {
00248                     for ( k = 0; k < nbxlines; k++ ) 
00249                         Vect_destroy_line_struct (  BXLines[k] );
00250                 }
00251                 G_free ( BXLines );
00252             }   
00253             if ( Err ) {
00254                 for ( l = 0; l < nx; l++ ) { /* Write out errors */
00255                     Vect_reset_line ( Points );
00256                     Vect_append_point ( Points, xx[l], yx[l], zx[l] );                            
00257                     ret = Vect_write_line ( Err, GV_POINT, Points, Cats );
00258                 }
00259                 
00260                 G_free ( xx );
00261                 G_free ( yx );
00262                 G_free ( zx );
00263             }
00264             
00265             if (msgout && printed > 1000) {
00266                 fprintf (msgout, "\rIntersections: %5d (line %d)", nbreaks, aline ); 
00267                 fflush ( msgout );
00268                 printed = 0;
00269             }
00270             printed++;
00271             if ( naxlines > 0 ) break; /* first line was broken and deleted -> take the next one */
00272         }
00273         nlines = Vect_get_num_lines (Map);
00274         G_debug (3, "nlines =  %d\n", nlines );
00275     }
00276     if (msgout) fprintf (msgout, "\rIntersections: %5d                         \n", nbreaks ); 
00277     Vect_destroy_list ( List );
00278 }
00279 

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5