break_polygons.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 <math.h>
00019 #include <grass/gis.h>
00020 #include <grass/Vect.h>
00021 
00022 typedef struct {
00023     double x, y;
00024     double a1, a2; /* angles */
00025     char cross; /* 0 - do not break, 1 - break */
00026     char used; /* 0 - was not used to break line, 1 - was used to break line
00027                *   this is stored because points are automaticaly marked as cross, even if not used 
00028                *   later to break lines */       
00029 } XPNT;
00030 
00031 static int fpoint;
00032 /* Function called from RTreeSearch for point found */
00033 void srch (int id, int *arg)
00034 {
00035     fpoint = id;
00036 }
00037 
00056 void 
00057 Vect_break_polygons ( struct Map_info *Map, int type, struct Map_info *Err, FILE *msgout )
00058 {
00059     struct line_pnts *BPoints, *Points;
00060     struct line_cats *Cats;
00061     int    i, j, k, ret, ltype, broken, last, nlines;
00062     int nbreaks;
00063     struct Node *RTree;
00064     int    apoints, npoints, nallpoints, nmarks;
00065     XPNT   *XPnts;
00066     struct Rect rect;
00067     double dx, dy, a1=0, a2=0;
00068     int closed, last_point, cross;
00069     int printed;
00070 
00071     RTree = RTreeNewIndex();
00072 
00073     BPoints = Vect_new_line_struct ();
00074     Points = Vect_new_line_struct ();
00075     Cats = Vect_new_cats_struct ();
00076     
00077     nlines = Vect_get_num_lines (Map);
00078 
00079     G_debug (3, "nlines =  %d", nlines );
00080     /* Go through all lines in vector, and add each point to structure of points,
00081      * if such point already exists check angles of segments and if differ mark for break */
00082 
00083     apoints = 0;
00084     nmarks = 0;
00085     npoints = 1; /* index starts from 1 !*/
00086     nallpoints = 0;
00087     XPnts = NULL;
00088     
00089     printed = 0;
00090     for ( i = 1; i <= nlines; i++ ){ 
00091         G_debug (3, "i =  %d", i);
00092         if ( !Vect_line_alive ( Map, i ) ) continue;
00093 
00094         ltype = Vect_read_line (Map, Points, Cats, i);
00095         if ( !(ltype & type) ) continue;
00096 
00097         /* This would be confused by duplicate coordinates (angle cannot be calculated) ->
00098          * prune line first */
00099         Vect_line_prune ( Points );
00100         
00101         /* If first and last point are identical it is close polygon, we dont need to register last point
00102          * and we can calculate angle for first.
00103          * If first and last point are not identical we have to mark for break both */
00104         last_point = Points->n_points-1;
00105         if ( Points->x[0] == Points->x[last_point] && Points->y[0] == Points->y[last_point] )
00106             closed = 1;
00107         else
00108             closed = 0;
00109         
00110         for ( j = 0; j <  Points->n_points; j++ ){ 
00111             G_debug (3, "j =  %d", j);
00112             nallpoints++;
00113 
00114             if ( j == last_point && closed ) continue; /* do not register last of close polygon */
00115 
00116             /* Box */
00117             rect.boundary[0] = Points->x[j];  rect.boundary[3] = Points->x[j];      
00118             rect.boundary[1] = Points->y[j];  rect.boundary[4] = Points->y[j];      
00119             rect.boundary[2] = 0;  rect.boundary[5] = 0;            
00120 
00121             /* Already in DB? */
00122             fpoint = -1;
00123             RTreeSearch(RTree, &rect, (void *)srch, 0);
00124             G_debug (3, "fpoint =  %d", fpoint);
00125 
00126             if ( Points->n_points <= 2 || (!closed && (j == 0 || j == last_point) ) ) {
00127                 cross = 1; /* mark for cross in any case */
00128             } else { /* calculate angles */
00129                 cross = 0;
00130                 if ( j == 0 && closed ) { /* closed polygon */
00131                     dx = Points->x[last_point] - Points->x[0]; dy = Points->y[last_point] - Points->y[0]; 
00132                     a1 = atan2 ( dy, dx );
00133                     dx = Points->x[1] - Points->x[0]; dy = Points->y[1] - Points->y[0]; 
00134                     a2 = atan2 ( dy, dx );
00135                 } else {
00136                     dx = Points->x[j-1] - Points->x[j]; dy = Points->y[j-1] - Points->y[j]; 
00137                     a1 = atan2 ( dy, dx );
00138                     dx = Points->x[j+1] - Points->x[j]; dy = Points->y[j+1] - Points->y[j]; 
00139                     a2 = atan2 ( dy, dx );
00140                 }
00141             }
00142 
00143             if ( fpoint > 0 ) { /* Found */
00144                 if ( XPnts[fpoint].cross == 1 ) continue; /* already marked */
00145                 
00146                 /* Check angles */
00147                 if ( cross ) {
00148                     XPnts[fpoint].cross = 1;
00149                     nmarks++;
00150                 } else { 
00151                     G_debug (3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1, XPnts[fpoint].a1,
00152                                         a2, XPnts[fpoint].a2);
00153                     if ( ( a1 == XPnts[fpoint].a1 && a2 == XPnts[fpoint].a2 ) ||
00154                          ( a1 == XPnts[fpoint].a2 && a2 == XPnts[fpoint].a1 ) ) 
00155                     { /* identical */
00156 
00157                     } else {
00158                         XPnts[fpoint].cross = 1;
00159                         nmarks++;
00160                     }
00161                 }
00162             } else { 
00163                 /* Add to tree and to structure */
00164                 RTreeInsertRect( &rect, npoints, &RTree, 0);
00165                 if ( npoints >= apoints ) {
00166                     apoints += 10000;
00167                     XPnts = (XPNT *) G_realloc ( XPnts, (apoints + 1) * sizeof (XPNT) );
00168                 }
00169                 XPnts[npoints].x = Points->x[j];
00170                 XPnts[npoints].y = Points->y[j];
00171                 XPnts[npoints].used = 0;
00172                 if ( j == 0 || j == (Points->n_points - 1) || Points->n_points < 3 ) {
00173                     XPnts[npoints].a1 = 0; 
00174                     XPnts[npoints].a2 = 0; 
00175                     XPnts[npoints].cross = 1;
00176                     nmarks++;
00177                 } else { 
00178                     XPnts[npoints].a1 = a1;
00179                     XPnts[npoints].a2 = a2;
00180                     XPnts[npoints].cross = 0;
00181                 }
00182                     
00183                 npoints++;                  
00184             }
00185         }
00186         if ( msgout && printed > 5000 ) {
00187             fprintf (msgout, "\rRegistering points ... %d", npoints - 1); 
00188             fflush ( msgout );
00189             printed = 0;
00190         }
00191         printed++;
00192     }
00193     if ( msgout ) {
00194         fprintf (msgout, "\rRegistering points ... %d", npoints - 1 ); 
00195         fprintf ( msgout, "\n" ); 
00196         fprintf ( msgout, "All points (vertices): %5d\n", nallpoints ); 
00197         fprintf ( msgout, "Registered points (unique coordinates): %5d\n", npoints - 1 ); 
00198         fprintf ( msgout, "Points marked for break: %5d\n", nmarks ); 
00199     }
00200 
00201     /* sleep (10); */
00202     
00203     nbreaks = 0;
00204     if ( msgout ) fprintf (msgout, "Breaks: %5d", nbreaks ); 
00205     printed = 0;
00206 
00207     /* Second loop through lines (existing when loop is started, no need to process lines written again)
00208      * and break at points marked for break */
00209     for ( i = 1; i <= nlines; i++ ){ 
00210         int n_orig_points;
00211         
00212         G_debug (3, "i =  %d", i);
00213         if ( !Vect_line_alive ( Map, i ) ) continue;
00214 
00215         ltype = Vect_read_line (Map, Points, Cats, i);
00216         if ( !(ltype & type) ) continue;
00217         if ( !(ltype & GV_LINES ) ) continue; /* Nonsense to break points */
00218 
00219         /* Duplicates would result in zero length lines -> prune line first */
00220         n_orig_points = Points->n_points;
00221         Vect_line_prune ( Points );
00222         
00223         broken = 0; last = 0;
00224         G_debug (3, "n_points =  %d", Points->n_points);
00225         for ( j = 1; j <  Points->n_points; j++ ){ 
00226             G_debug (3, "j =  %d", j);
00227             nallpoints++;
00228 
00229             /* Box */
00230             rect.boundary[0] = Points->x[j];  rect.boundary[3] = Points->x[j];      
00231             rect.boundary[1] = Points->y[j];  rect.boundary[4] = Points->y[j];      
00232             rect.boundary[2] = 0;  rect.boundary[5] = 0;            
00233 
00234             if ( Points->n_points <= 1 || (j == (Points->n_points-1) && !broken) ) break; 
00235               /* One point only or 
00236                * last point and line is not broken, do nothing */
00237 
00238             RTreeSearch(RTree, &rect, (void *)srch, 0);
00239             G_debug (3, "fpoint =  %d", fpoint);
00240 
00241             if ( XPnts[fpoint].cross ) { /* realy use to break line */
00242                 XPnts[fpoint].used = 1;
00243             }
00244                 
00245             /* break or write last segment of broken line */
00246             if ( (j == (Points->n_points-1) && broken ) || XPnts[fpoint].cross  ) { 
00247                 Vect_reset_line ( BPoints );
00248                 for ( k = last; k <= j; k++ ){ 
00249                     Vect_append_point ( BPoints, Points->x[k], Points->y[k], Points->z[k] ); 
00250                 }
00251 
00252                 /* Result may collapse to one point */
00253                 Vect_line_prune ( BPoints );
00254                 if ( BPoints->n_points > 1 ) {
00255                     ret = Vect_write_line ( Map, ltype, BPoints, Cats );  
00256                     G_debug (3, "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d", 
00257                                          ret, j, Points->n_points, BPoints->n_points);
00258                 }
00259                 
00260                 if ( !broken ) Vect_delete_line (Map, i); /* not yet deleted */
00261 
00262                 last = j;
00263                 broken = 1;
00264                 nbreaks++;
00265             }
00266                 
00267             if ( msgout && printed > 1000 ) {
00268                 fprintf ( msgout, "\rBreaks: %5d", nbreaks ); 
00269                 fflush ( msgout );
00270                 printed = 0;
00271             }
00272             printed++;
00273         }
00274         if ( !broken && n_orig_points > Points->n_points ) { /* was pruned before -> rewrite */
00275             if ( Points->n_points > 1 ) {
00276                 Vect_rewrite_line ( Map, i, ltype, Points, Cats );  
00277                 G_debug (3, "Line %d pruned, npoints = %d", i, Points->n_points );
00278             } else { 
00279                 Vect_delete_line (Map, i);
00280                 G_debug (3, "Line %d was deleted", i);
00281             }
00282         } else {
00283             G_debug (3, "Line %d was not changed", i);
00284         }
00285     }
00286 
00287     if ( msgout )
00288         fprintf ( msgout, "\rBreaks: %5d", nbreaks ); 
00289 
00290     /* Write points on breaks */
00291     if ( Err ) {
00292         Vect_reset_cats ( Cats );
00293         for ( i = 1; i < npoints; i++ ){ 
00294             if ( XPnts[i].used ) {
00295                 Vect_reset_line ( Points );
00296                 Vect_append_point ( Points, XPnts[i].x, XPnts[i].y, 0 ); 
00297                 Vect_write_line ( Err, GV_POINT, Points, Cats );
00298             }
00299         }
00300     }
00301 
00302     if ( msgout ) fprintf ( msgout, "\n" ); 
00303 
00304     G_free ( XPnts );
00305     RTreeDestroyNode ( RTree);
00306 }
00307 

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