00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdlib.h>
00018 #include <math.h>
00019 #include "gis.h"
00020 #include "Vect.h"
00021
00022
00023 typedef struct {
00024 double x, y;
00025 int anchor;
00026
00027
00028 } XPNT;
00029
00030 typedef struct {
00031 int anchor;
00032 double along;
00033 } NEW;
00034
00035
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
00065
00066
00067
00068
00069
00070
00071
00072
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;
00088 int nanchors, ntosnap;
00089 int nsnapped, ncreated;
00090 int apoints, npoints, nvertices;
00091 XPNT *XPnts;
00092 NEW *New = NULL;
00093 int anew = 0, nnew;
00094 struct Rect rect;
00095 struct ilist *List;
00096 int *Index = NULL;
00097 int aindex = 0;
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
00111 apoints = 0;
00112 point = 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
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
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 ) {
00143
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
00170
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;
00179 nanchors++;
00180
00181
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
00215
00216
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
00238 for ( v = 0; v < Points->n_points; v++ ){
00239
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
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 ) {
00253 Points->x[v] = XPnts[anchor].x;
00254 Points->y[v] = XPnts[anchor].y;
00255 nsnapped++;
00256 changed = 1;
00257 Index[v] = anchor;
00258 } else {
00259 Index[v] = spoint;
00260 }
00261 }
00262
00263
00264 Vect_reset_line (NPoints);
00265
00266
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
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
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
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;
00305 if ( XPnts[spoint].anchor > 0 ) continue;
00306
00307
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
00327 if ( nnew > 0 ) {
00328
00329 qsort ( New, nnew, sizeof ( NEW), sort_new);
00330
00331 for ( i = 0; i < nnew; i++ ) {
00332 anchor = New[i].anchor;
00333
00334 Vect_append_point ( NPoints, XPnts[anchor].x, XPnts[anchor].y, 0 );
00335 ncreated++;
00336 }
00337 changed = 1;
00338 }
00339 }
00340
00341 v = Points->n_points-1;
00342 Vect_append_point ( NPoints, Points->x[v], Points->y[v], Points->z[v] );
00343
00344 if ( changed ) {
00345 Vect_line_prune ( Points );
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
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