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

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