snap.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 "gis.h"
00020 #include "Vect.h"
00021 
00022 /* Vertex */
00023 typedef struct {
00024     double x, y;
00025     int    anchor; /* 0 - anchor, do not snap this point, that means snap others to this */
00026                    /* >0  - index of anchor to which snap this point */
00027                    /* -1  - init value */
00028 } XPNT;
00029 
00030 typedef struct {
00031     int    anchor; 
00032     double along;
00033 } NEW;
00034 
00035 /* This function is called by  RTreeSearch() to add selected node/line/area/isle to thelist */
00036 int add_item(int id, struct ilist *list)
00037 {
00038         dig_list_add ( list, id );
00039             return 1;
00040 }
00041 
00042 static int sort_new( const void *, const void * );
00043 
00064 /* As mentioned above, lines are not necessarily snapped to nearest vertex! For example:
00065      |                    
00066      | 1         line 3 is snapped to line 1,
00067      |           then line 2 is not snapped to common node at lines 1 and 3,
00068                  because it is already outside of threshold
00069 ----------- 3   
00070 
00071      |
00072      | 2
00073      |  
00074 */
00075 
00076 
00077 void 
00078 Vect_snap_lines ( struct Map_info *Map, int type, double thresh, struct Map_info *Err, FILE *msgout )
00079 {
00080     struct line_pnts *Points, *NPoints;
00081     struct line_cats *Cats;
00082     int    nlines, line, ltype;
00083     double thresh2;
00084     int    printed;
00085 
00086     struct Node *RTree;
00087     int    point;   /* index in points array */
00088     int    nanchors, ntosnap; /* number of anchors and number of points to be snapped */
00089     int    nsnapped, ncreated; /* number of snapped verices, number of new vertices (on segments) */
00090     int    apoints, npoints, nvertices; /* number of allocated points, registered points, vertices */
00091     XPNT   *XPnts;  /* Array of points */
00092     NEW    *New = NULL;    /* Array of new points */
00093     int    anew = 0, nnew;   /* allocated new points , number of new points */
00094     struct Rect rect;
00095     struct ilist *List;
00096     int *Index = NULL;  /* indexes of anchors for vertices */
00097     int aindex = 0; /* allocated Index */
00098 
00099     Points = Vect_new_line_struct ();
00100     NPoints = Vect_new_line_struct ();
00101     Cats = Vect_new_cats_struct ();
00102     List = Vect_new_list();
00103     RTree = RTreeNewIndex();
00104     
00105     thresh2 = thresh * thresh;
00106     nlines = Vect_get_num_lines (Map);
00107 
00108     G_debug (3, "nlines =  %d", nlines );
00109     
00110     /* Go through all lines in vector, and add each point to structure of points */
00111     apoints = 0;
00112     point = 1; /* index starts from 1 !*/
00113     nvertices = 0;
00114     XPnts = NULL;
00115     printed = 0;
00116 
00117     if ( msgout ) fprintf (msgout, "Registering points ..."); 
00118     
00119     for ( line = 1; line <= nlines; line++ ){ 
00120         int v;
00121         
00122         G_debug (3, "line =  %d", line);
00123         if ( !Vect_line_alive ( Map, line ) ) continue;
00124 
00125         ltype = Vect_read_line (Map, Points, Cats, line);
00126         if ( !(ltype & type) ) continue;
00127 
00128         for ( v = 0; v <  Points->n_points; v++ ){ 
00129             G_debug (3, "  vertex v = %d", v);
00130             nvertices++;
00131 
00132             /* Box */
00133             rect.boundary[0] = Points->x[v];  rect.boundary[3] = Points->x[v];      
00134             rect.boundary[1] = Points->y[v];  rect.boundary[4] = Points->y[v];      
00135             rect.boundary[2] = 0;  rect.boundary[5] = 0;            
00136 
00137             /* Already registered ? */
00138             Vect_reset_list ( List );
00139             RTreeSearch(RTree, &rect, (void *)add_item, List);
00140             G_debug (3, "List : nvalues =  %d", List->n_values);
00141 
00142             if ( List->n_values == 0 ) { /* Not found */
00143                 /* Add to tree and to structure */
00144                 RTreeInsertRect( &rect, point, &RTree, 0);
00145                 if ( (point - 1) == apoints ) {
00146                     apoints += 10000;
00147                     XPnts = (XPNT *) G_realloc ( XPnts, (apoints + 1) * sizeof (XPNT) );
00148                 }
00149                 XPnts[point].x = Points->x[v];
00150                 XPnts[point].y = Points->y[v];
00151                 XPnts[point].anchor = -1;
00152                 point++;                    
00153             }
00154         }
00155         if ( msgout && printed > 1000 ) {
00156             fprintf (msgout, "\rRegistering points ... %d", point - 1); 
00157             fflush ( msgout );
00158             printed = 0;
00159         }
00160         printed++;
00161     }
00162     npoints = point - 1;
00163     if ( msgout ) {
00164         fprintf (msgout, "\r                                               \r" ); 
00165         fprintf ( msgout, "All vertices: %5d\n", nvertices ); 
00166         fprintf ( msgout, "Registered points (unique coordinates): %5d\n", npoints ); 
00167     }
00168 
00169     /* Go through all registered points and if not yet marked mark it as anchor and assign this anchor
00170      * to all not yet marked points in threshold */
00171     nanchors = ntosnap = 0;
00172     for ( point = 1; point <= npoints; point++ ) { 
00173         int i;
00174         G_debug (3, "  point = %d", point);
00175         
00176         if ( XPnts[point].anchor >= 0 ) continue;
00177 
00178         XPnts[point].anchor = 0; /* make it anchor */
00179         nanchors++;
00180 
00181         /* Find points in threshold */
00182         rect.boundary[0] = XPnts[point].x - thresh;  
00183         rect.boundary[3] = XPnts[point].x + thresh;         
00184         rect.boundary[1] = XPnts[point].y - thresh;  
00185         rect.boundary[4] = XPnts[point].y + thresh;
00186         rect.boundary[2] = 0;  rect.boundary[5] = 0;        
00187 
00188         Vect_reset_list ( List );
00189         RTreeSearch(RTree, &rect, (void *)add_item, List);
00190         G_debug (4, "  %d points in threshold box", List->n_values);
00191         
00192         for ( i = 0; i < List->n_values; i++ ) {
00193             int    pointb;
00194             double dx, dy, dist2;
00195 
00196             pointb = List->value[i];
00197             if ( pointb == point ) continue;
00198 
00199             dx = XPnts[pointb].x - XPnts[point].x;
00200             dy = XPnts[pointb].y - XPnts[point].y;
00201             dist2 = dx * dx + dy * dy;
00202 
00203             if ( dist2 <= thresh2 ) {
00204                 XPnts[pointb].anchor = point;
00205                 ntosnap++;
00206             }
00207         }
00208     }
00209     if ( msgout ) {
00210         fprintf ( msgout, "Nodes marked as anchor     : %5d\n", nanchors ); 
00211         fprintf ( msgout, "Nodes marked to be snapped : %5d\n", ntosnap ); 
00212     }
00213 
00214     /* Go through all lines and: 
00215      *   1) for all vertices: if not anchor snap it to its anchor
00216      *   2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course) */
00217     
00218     printed = 0;
00219     nsnapped = ncreated = 0;
00220     if ( msgout ) fprintf (msgout, "Snaps: %5d", nsnapped + ncreated ); 
00221     
00222     for ( line = 1; line <= nlines; line++ ){ 
00223         int v, spoint, anchor;
00224         int changed = 0;
00225         
00226         G_debug (3, "line =  %d", line);
00227         if ( !Vect_line_alive ( Map, line ) ) continue;
00228 
00229         ltype = Vect_read_line (Map, Points, Cats, line);
00230         if ( !(ltype & type) ) continue;
00231 
00232         if ( Points->n_points >= aindex ) {
00233             aindex = Points->n_points;
00234             Index = (int *) G_realloc ( Index, aindex * sizeof(int) );
00235         }
00236 
00237         /* Snap all vertices */
00238         for ( v = 0; v <  Points->n_points; v++ ){ 
00239             /* Box */
00240             rect.boundary[0] = Points->x[v];  rect.boundary[3] = Points->x[v];      
00241             rect.boundary[1] = Points->y[v];  rect.boundary[4] = Points->y[v];      
00242             rect.boundary[2] = 0;  rect.boundary[5] = 0;            
00243 
00244             /* Find point ( should always find one point )*/
00245             Vect_reset_list ( List );
00246             
00247             RTreeSearch(RTree, &rect, (void *)add_item, List);
00248 
00249             spoint = List->value[0];
00250             anchor = XPnts[spoint].anchor;
00251 
00252             if ( anchor > 0 ) { /* to be snapped */
00253                 Points->x[v] = XPnts[anchor].x;
00254                 Points->y[v] = XPnts[anchor].y;
00255                 nsnapped++;                 
00256                 changed = 1;
00257                 Index[v] = anchor; /* point on new location */
00258             } else {
00259                 Index[v] = spoint; /* old point */
00260             }
00261         }
00262         
00263         /* New points */
00264         Vect_reset_line (NPoints);
00265 
00266         /* Snap all segments to anchors in threshold */
00267         for ( v = 0; v < Points->n_points - 1; v++ ){ 
00268             int    i;
00269             double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
00270             
00271             G_debug (3, "  segment = %d end anchors : %d  %d", v, Index[v], Index[v+1]);
00272             
00273             x1 = Points->x[v];
00274             x2 = Points->x[v+1];
00275             y1 = Points->y[v];
00276             y2 = Points->y[v+1];
00277 
00278             Vect_append_point ( NPoints, Points->x[v], Points->y[v], Points->z[v] );
00279 
00280             /* Box */
00281             if ( x1 <= x2 ) { xmin = x1; xmax = x2; } else { xmin = x2; xmax = x1; }
00282             if ( y1 <= y2 ) { ymin = y1; ymax = y2; } else { ymin = y2; ymax = y1; }
00283             
00284             rect.boundary[0] = xmin - thresh;  
00285             rect.boundary[3] = xmax + thresh;       
00286             rect.boundary[1] = ymin - thresh;  
00287             rect.boundary[4] = ymax + thresh;       
00288             rect.boundary[2] = 0;  rect.boundary[5] = 0;            
00289 
00290             /* Find points */
00291             Vect_reset_list ( List );
00292             RTreeSearch(RTree, &rect, (void *)add_item, List);
00293         
00294             G_debug (3, "  %d points in box", List->n_values);
00295 
00296             /* Snap to anchor in threshold different from end points */
00297             nnew = 0;
00298             for ( i = 0; i < List->n_values; i++ ) {
00299                 double dist2, along;
00300                 
00301                 spoint = List->value[i];
00302                 G_debug (4, "    spoint = %d anchor = %d", spoint, XPnts[spoint].anchor);
00303 
00304                 if ( spoint == Index[v] || spoint == Index[v+1] ) continue; /* end point */
00305                 if ( XPnts[spoint].anchor > 0 ) continue; /* point is not anchor */
00306 
00307                 /* Check the distance */
00308                 dist2 = dig_distance2_point_to_line ( XPnts[spoint].x, XPnts[spoint].y, 0,
00309                         x1, y1, 0, x2, y2, 0, 0, NULL, NULL, NULL, &along, NULL );
00310                     
00311                 G_debug (4, "      distance = %lf", sqrt(dist2));
00312 
00313                 if ( dist2 <= thresh2 ) {
00314                     G_debug (4, "      anchor in thresh, along = %lf", along);
00315 
00316                     if ( nnew == anew ) {
00317                         anew += 100;
00318                         New = (NEW *) G_realloc ( New, anew * sizeof (NEW) );
00319                     }
00320                     New[nnew].anchor = spoint;
00321                     New[nnew].along = along;
00322                     nnew++;                 
00323                 }
00324             }
00325             G_debug (3, "  nnew = %d", nnew);
00326             /* insert new vertices */
00327             if ( nnew > 0 ) {
00328                 /* sort by distance along the segment */
00329                 qsort ( New, nnew, sizeof ( NEW), sort_new);
00330                 
00331                 for ( i = 0; i < nnew; i++ ) {
00332                     anchor = New[i].anchor;
00333                     /* Vect_line_insert_point ( Points, ++v, XPnts[anchor].x, XPnts[anchor].y, 0); */
00334                     Vect_append_point ( NPoints, XPnts[anchor].x, XPnts[anchor].y, 0 );
00335                     ncreated++;
00336                 }
00337                 changed = 1;
00338             }
00339         }
00340         /* append end point */
00341         v = Points->n_points-1; 
00342         Vect_append_point ( NPoints, Points->x[v], Points->y[v], Points->z[v] );
00343 
00344         if ( changed ) { /* rewrite the line */
00345             Vect_line_prune ( Points );  /* remove duplicates */
00346             if ( Points->n_points > 1 || ltype & GV_LINES )
00347                 Vect_rewrite_line ( Map, line, ltype, NPoints, Cats );  
00348             else
00349                 Vect_delete_line ( Map, line);
00350         }
00351         if ( msgout && printed > 1000 ) {
00352             fprintf (msgout, "\rSnaps: %5d  (line = %d)", nsnapped + ncreated, line ); 
00353             fflush ( msgout );
00354             printed = 0;
00355         }
00356         printed++;
00357 
00358     }
00359     if ( msgout ) {
00360         fprintf ( msgout, "\rSnapped vertices : %5d                             \n", nsnapped ); 
00361         fprintf ( msgout, "New vertices     : %5d\n", ncreated ); 
00362     }
00363     
00364     Vect_destroy_line_struct ( Points );
00365     Vect_destroy_line_struct ( NPoints );
00366     Vect_destroy_cats_struct ( Cats );
00367     G_free ( XPnts );
00368     G_free (Index);
00369     G_free ( New );
00370     RTreeDestroyNode ( RTree);
00371 }
00372 
00373 /* for qsort */
00374 static int sort_new (  const void *pa, const void *pb )
00375 {
00376     NEW *p1 = (NEW *) pa;
00377     NEW *p2 = (NEW *) pb;
00378 
00379     if ( p1->along < p2->along ) return -1;
00380     if ( p1->along > p2->along ) return 1;
00381     return 1;
00382 }
00383 

Generated on Sat Jul 22 22:05:57 2006 for GRASS by  doxygen 1.4.7