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 <grass/gis.h>
00020 #include <grass/Vect.h>
00021 #include <grass/glocale.h>
00022
00023
00024 static int sort_new(const void *pa, const void *pb);
00025
00026
00027
00028 typedef struct {
00029 double x, y;
00030 int anchor;
00031
00032
00033 } XPNT;
00034
00035 typedef struct {
00036 int anchor;
00037 double along;
00038 } NEW;
00039
00040
00041 int add_item(int id, struct ilist *list)
00042 {
00043 dig_list_add ( list, id );
00044 return 1;
00045 }
00046
00047
00068
00069
00070
00071
00072
00073
00074
00075
00076
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;
00092 int nanchors, ntosnap;
00093 int nsnapped, ncreated;
00094 int apoints, npoints, nvertices;
00095 XPNT *XPnts;
00096 NEW *New = NULL;
00097 int anew = 0, nnew;
00098 struct Rect rect;
00099 struct ilist *List;
00100 int *Index = NULL;
00101 int aindex = 0;
00102 int width = 26;
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
00116 apoints = 0;
00117 point = 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
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
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 ) {
00148
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
00176
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;
00185 nanchors++;
00186
00187
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
00221
00222
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
00244 for ( v = 0; v < Points->n_points; v++ ){
00245
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
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 ) {
00259 Points->x[v] = XPnts[anchor].x;
00260 Points->y[v] = XPnts[anchor].y;
00261 nsnapped++;
00262 changed = 1;
00263 Index[v] = anchor;
00264 } else {
00265 Index[v] = spoint;
00266 }
00267 }
00268
00269
00270 Vect_reset_line (NPoints);
00271
00272
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
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
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
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;
00311 if ( XPnts[spoint].anchor > 0 ) continue;
00312
00313
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
00333 if ( nnew > 0 ) {
00334
00335 qsort ( New, sizeof(char)*nnew, sizeof ( NEW), sort_new);
00336
00337 for ( i = 0; i < nnew; i++ ) {
00338 anchor = New[i].anchor;
00339
00340 Vect_append_point ( NPoints, XPnts[anchor].x, XPnts[anchor].y, 0 );
00341 ncreated++;
00342 }
00343 changed = 1;
00344 }
00345 }
00346
00347
00348 v = Points->n_points-1;
00349 Vect_append_point ( NPoints, Points->x[v], Points->y[v], Points->z[v] );
00350
00351 if ( changed ) {
00352 Vect_line_prune ( NPoints );
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 }
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
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