00001
00020 #include <math.h>
00021 #include <grass/gis.h>
00022 #include <grass/Vect.h>
00023 #include <grass/glocale.h>
00024
00025
00026 static int sort_new(const void *pa, const void *pb);
00027
00028
00029 typedef struct
00030 {
00031 double x, y;
00032 int anchor;
00033
00034
00035 } XPNT;
00036
00037 typedef struct
00038 {
00039 int anchor;
00040 double along;
00041 } NEW;
00042
00043
00044 int add_item(int id, struct ilist *list)
00045 {
00046 dig_list_add(list, id);
00047 return 1;
00048 }
00049
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079 void
00080 Vect_snap_lines_list(struct Map_info *Map, struct ilist *List_lines,
00081 double thresh, struct Map_info *Err)
00082 {
00083 struct line_pnts *Points, *NPoints;
00084 struct line_cats *Cats;
00085 int line, ltype, line_idx;
00086 double thresh2;
00087
00088 struct Node *RTree;
00089 int point;
00090 int nanchors, ntosnap;
00091 int nsnapped, ncreated;
00092 int apoints, npoints, nvertices;
00093 XPNT *XPnts;
00094 NEW *New = NULL;
00095 int anew = 0, nnew;
00096 struct Rect rect;
00097 struct ilist *List;
00098 int *Index = NULL;
00099 int aindex = 0;
00100
00101 if (List_lines->n_values < 1)
00102 return;
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
00112
00113 apoints = 0;
00114 point = 1;
00115 nvertices = 0;
00116 XPnts = NULL;
00117
00118 for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
00119 int v;
00120
00121 line = List_lines->value[line_idx];
00122
00123 G_debug(3, "line = %d", line);
00124 if (!Vect_line_alive(Map, line))
00125 continue;
00126
00127 ltype = Vect_read_line(Map, Points, Cats, line);
00128
00129 for (v = 0; v < Points->n_points; v++) {
00130 G_debug(3, " vertex v = %d", v);
00131 nvertices++;
00132
00133
00134 rect.boundary[0] = Points->x[v];
00135 rect.boundary[3] = Points->x[v];
00136 rect.boundary[1] = Points->y[v];
00137 rect.boundary[4] = Points->y[v];
00138 rect.boundary[2] = 0;
00139 rect.boundary[5] = 0;
00140
00141
00142 Vect_reset_list(List);
00143 RTreeSearch(RTree, &rect, (void *)add_item, List);
00144 G_debug(3, "List : nvalues = %d", List->n_values);
00145
00146 if (List->n_values == 0) {
00147
00148 RTreeInsertRect(&rect, point, &RTree, 0);
00149 if ((point - 1) == apoints) {
00150 apoints += 10000;
00151 XPnts =
00152 (XPNT *) G_realloc(XPnts,
00153 (apoints + 1) * sizeof(XPNT));
00154 }
00155 XPnts[point].x = Points->x[v];
00156 XPnts[point].y = Points->y[v];
00157 XPnts[point].anchor = -1;
00158 point++;
00159 }
00160 }
00161 }
00162 npoints = point - 1;
00163
00164
00165
00166 nanchors = ntosnap = 0;
00167 for (point = 1; point <= npoints; point++) {
00168 int i;
00169
00170 G_debug(3, " point = %d", point);
00171
00172 if (XPnts[point].anchor >= 0)
00173 continue;
00174
00175 XPnts[point].anchor = 0;
00176 nanchors++;
00177
00178
00179 rect.boundary[0] = XPnts[point].x - thresh;
00180 rect.boundary[3] = XPnts[point].x + thresh;
00181 rect.boundary[1] = XPnts[point].y - thresh;
00182 rect.boundary[4] = XPnts[point].y + thresh;
00183 rect.boundary[2] = 0;
00184 rect.boundary[5] = 0;
00185
00186 Vect_reset_list(List);
00187 RTreeSearch(RTree, &rect, (void *)add_item, List);
00188 G_debug(4, " %d points in threshold box", List->n_values);
00189
00190 for (i = 0; i < List->n_values; i++) {
00191 int pointb;
00192 double dx, dy, dist2;
00193
00194 pointb = List->value[i];
00195 if (pointb == point)
00196 continue;
00197
00198 dx = XPnts[pointb].x - XPnts[point].x;
00199 dy = XPnts[pointb].y - XPnts[point].y;
00200 dist2 = dx * dx + dy * dy;
00201
00202 if (dist2 <= thresh2) {
00203 XPnts[pointb].anchor = point;
00204 ntosnap++;
00205 }
00206 }
00207 }
00208
00209
00210
00211
00212
00213 nsnapped = ncreated = 0;
00214
00215 for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) {
00216 int v, spoint, anchor;
00217 int changed = 0;
00218
00219 line = List_lines->value[line_idx];
00220
00221 G_debug(3, "line = %d", line);
00222 if (!Vect_line_alive(Map, line))
00223 continue;
00224
00225 ltype = Vect_read_line(Map, Points, Cats, line);
00226
00227 if (Points->n_points >= aindex) {
00228 aindex = Points->n_points;
00229 Index = (int *)G_realloc(Index, aindex * sizeof(int));
00230 }
00231
00232
00233 for (v = 0; v < Points->n_points; v++) {
00234
00235 rect.boundary[0] = Points->x[v];
00236 rect.boundary[3] = Points->x[v];
00237 rect.boundary[1] = Points->y[v];
00238 rect.boundary[4] = Points->y[v];
00239 rect.boundary[2] = 0;
00240 rect.boundary[5] = 0;
00241
00242
00243 Vect_reset_list(List);
00244
00245 RTreeSearch(RTree, &rect, (void *)add_item, List);
00246
00247 spoint = List->value[0];
00248 anchor = XPnts[spoint].anchor;
00249
00250 if (anchor > 0) {
00251 Points->x[v] = XPnts[anchor].x;
00252 Points->y[v] = XPnts[anchor].y;
00253 nsnapped++;
00254 changed = 1;
00255 Index[v] = anchor;
00256 }
00257 else {
00258 Index[v] = spoint;
00259 }
00260 }
00261
00262
00263 Vect_reset_line(NPoints);
00264
00265
00266 for (v = 0; v < Points->n_points - 1; v++) {
00267 int i;
00268 double x1, x2, y1, y2, xmin, xmax, ymin, ymax;
00269
00270 G_debug(3, " segment = %d end anchors : %d %d", v, Index[v],
00271 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],
00279 Points->z[v]);
00280
00281
00282 if (x1 <= x2) {
00283 xmin = x1;
00284 xmax = x2;
00285 }
00286 else {
00287 xmin = x2;
00288 xmax = x1;
00289 }
00290 if (y1 <= y2) {
00291 ymin = y1;
00292 ymax = y2;
00293 }
00294 else {
00295 ymin = y2;
00296 ymax = y1;
00297 }
00298
00299 rect.boundary[0] = xmin - thresh;
00300 rect.boundary[3] = xmax + thresh;
00301 rect.boundary[1] = ymin - thresh;
00302 rect.boundary[4] = ymax + thresh;
00303 rect.boundary[2] = 0;
00304 rect.boundary[5] = 0;
00305
00306
00307 Vect_reset_list(List);
00308 RTreeSearch(RTree, &rect, (void *)add_item, List);
00309
00310 G_debug(3, " %d points in box", List->n_values);
00311
00312
00313 nnew = 0;
00314 for (i = 0; i < List->n_values; i++) {
00315 double dist2, along;
00316
00317 spoint = List->value[i];
00318 G_debug(4, " spoint = %d anchor = %d", spoint,
00319 XPnts[spoint].anchor);
00320
00321 if (spoint == Index[v] || spoint == Index[v + 1])
00322 continue;
00323 if (XPnts[spoint].anchor > 0)
00324 continue;
00325
00326
00327 dist2 =
00328 dig_distance2_point_to_line(XPnts[spoint].x,
00329 XPnts[spoint].y, 0, x1, y1, 0,
00330 x2, y2, 0, 0, NULL, NULL,
00331 NULL, &along, NULL);
00332
00333 G_debug(4, " distance = %lf", sqrt(dist2));
00334
00335 if (dist2 <= thresh2) {
00336 G_debug(4, " anchor in thresh, along = %lf", along);
00337
00338 if (nnew == anew) {
00339 anew += 100;
00340 New = (NEW *) G_realloc(New, anew * sizeof(NEW));
00341 }
00342 New[nnew].anchor = spoint;
00343 New[nnew].along = along;
00344 nnew++;
00345 }
00346 }
00347 G_debug(3, " nnew = %d", nnew);
00348
00349 if (nnew > 0) {
00350
00351 qsort(New, sizeof(char) * nnew, sizeof(NEW), sort_new);
00352
00353 for (i = 0; i < nnew; i++) {
00354 anchor = New[i].anchor;
00355
00356 Vect_append_point(NPoints, XPnts[anchor].x,
00357 XPnts[anchor].y, 0);
00358 ncreated++;
00359 }
00360 changed = 1;
00361 }
00362 }
00363
00364
00365 v = Points->n_points - 1;
00366 Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
00367
00368 if (changed) {
00369 Vect_line_prune(NPoints);
00370 if (NPoints->n_points > 1 || ltype & GV_LINES) {
00371 Vect_rewrite_line(Map, line, ltype, NPoints, Cats);
00372 }
00373 else {
00374 Vect_delete_line(Map, line);
00375 }
00376 if (Err) {
00377 Vect_write_line(Err, ltype, Points, Cats);
00378 }
00379 }
00380 }
00381
00382 Vect_destroy_line_struct(Points);
00383 Vect_destroy_line_struct(NPoints);
00384 Vect_destroy_cats_struct(Cats);
00385 G_free(XPnts);
00386 G_free(Index);
00387 G_free(New);
00388 RTreeDestroyNode(RTree);
00389 }
00390
00391
00392 static int sort_new(const void *pa, const void *pb)
00393 {
00394 NEW *p1 = (NEW *) pa;
00395 NEW *p2 = (NEW *) pb;
00396
00397 if (p1->along < p2->along)
00398 return -1;
00399 if (p1->along > p2->along)
00400 return 1;
00401 return 1;
00402 }
00403
00416 void
00417 Vect_snap_lines(struct Map_info *Map, int type, double thresh,
00418 struct Map_info *Err)
00419 {
00420 int line, nlines, ltype;
00421
00422 struct ilist *List;
00423
00424 List = Vect_new_list();
00425
00426 nlines = Vect_get_num_lines(Map);
00427
00428 for (line = 1; line <= nlines; line++) {
00429 G_debug(3, "line = %d", line);
00430
00431 if (!Vect_line_alive(Map, line))
00432 continue;
00433
00434 ltype = Vect_read_line(Map, NULL, NULL, line);
00435
00436 if (!(ltype & type))
00437 continue;
00438
00439 Vect_list_append(List, line);
00440 }
00441
00442 Vect_snap_lines_list(Map, List, thresh, Err);
00443
00444 Vect_destroy_list(List);
00445
00446 return;
00447 }