GRASS Programmer's Manual
6.4.1(2011)
|
00001 00020 #include <math.h> 00021 #include <grass/gis.h> 00022 #include <grass/Vect.h> 00023 #include <grass/glocale.h> 00024 00025 /* function prototypes */ 00026 static int sort_new(const void *pa, const void *pb); 00027 00028 /* Vertex */ 00029 typedef struct 00030 { 00031 double x, y; 00032 int anchor; /* 0 - anchor, do not snap this point, that means snap others to this */ 00033 /* >0 - index of anchor to which snap this point */ 00034 /* -1 - init value */ 00035 } XPNT; 00036 00037 typedef struct 00038 { 00039 int anchor; 00040 double along; 00041 } NEW; 00042 00043 /* This function is called by RTreeSearch() to add selected node/line/area/isle to thelist */ 00044 int add_item(int id, struct ilist *list) 00045 { 00046 dig_list_add(list, id); 00047 return 1; 00048 } 00049 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 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; /* index in points array */ 00090 int nanchors, ntosnap; /* number of anchors and number of points to be snapped */ 00091 int nsnapped, ncreated; /* number of snapped verices, number of new vertices (on segments) */ 00092 int apoints, npoints, nvertices; /* number of allocated points, registered points, vertices */ 00093 XPNT *XPnts; /* Array of points */ 00094 NEW *New = NULL; /* Array of new points */ 00095 int anew = 0, nnew; /* allocated new points , number of new points */ 00096 struct Rect rect; 00097 struct ilist *List; 00098 int *Index = NULL; /* indexes of anchors for vertices */ 00099 int aindex = 0; /* allocated Index */ 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 /* Go through all lines in vector, and add each point to structure of points */ 00113 apoints = 0; 00114 point = 1; /* index starts from 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 G_percent(line_idx, List_lines->n_values, 2); 00122 00123 line = List_lines->value[line_idx]; 00124 00125 G_debug(3, "line = %d", line); 00126 if (!Vect_line_alive(Map, line)) 00127 continue; 00128 00129 ltype = Vect_read_line(Map, Points, Cats, line); 00130 00131 for (v = 0; v < Points->n_points; v++) { 00132 G_debug(3, " vertex v = %d", v); 00133 nvertices++; 00134 00135 /* Box */ 00136 rect.boundary[0] = Points->x[v]; 00137 rect.boundary[3] = Points->x[v]; 00138 rect.boundary[1] = Points->y[v]; 00139 rect.boundary[4] = Points->y[v]; 00140 rect.boundary[2] = 0; 00141 rect.boundary[5] = 0; 00142 00143 /* Already registered ? */ 00144 Vect_reset_list(List); 00145 RTreeSearch(RTree, &rect, (void *)add_item, List); 00146 G_debug(3, "List : nvalues = %d", List->n_values); 00147 00148 if (List->n_values == 0) { /* Not found */ 00149 /* Add to tree and to structure */ 00150 RTreeInsertRect(&rect, point, &RTree, 0); 00151 if ((point - 1) == apoints) { 00152 apoints += 10000; 00153 XPnts = 00154 (XPNT *) G_realloc(XPnts, 00155 (apoints + 1) * sizeof(XPNT)); 00156 } 00157 XPnts[point].x = Points->x[v]; 00158 XPnts[point].y = Points->y[v]; 00159 XPnts[point].anchor = -1; 00160 point++; 00161 } 00162 } 00163 } 00164 G_percent(line_idx, List_lines->n_values, 2); /* finish it */ 00165 00166 npoints = point - 1; 00167 00168 /* Go through all registered points and if not yet marked mark it as anchor and assign this anchor 00169 * to all not yet marked points in threshold */ 00170 nanchors = ntosnap = 0; 00171 for (point = 1; point <= npoints; point++) { 00172 int i; 00173 00174 G_percent(point, npoints, 2); 00175 00176 G_debug(3, " point = %d", point); 00177 00178 if (XPnts[point].anchor >= 0) 00179 continue; 00180 00181 XPnts[point].anchor = 0; /* make it anchor */ 00182 nanchors++; 00183 00184 /* Find points in threshold */ 00185 rect.boundary[0] = XPnts[point].x - thresh; 00186 rect.boundary[3] = XPnts[point].x + thresh; 00187 rect.boundary[1] = XPnts[point].y - thresh; 00188 rect.boundary[4] = XPnts[point].y + thresh; 00189 rect.boundary[2] = 0; 00190 rect.boundary[5] = 0; 00191 00192 Vect_reset_list(List); 00193 RTreeSearch(RTree, &rect, (void *)add_item, List); 00194 G_debug(4, " %d points in threshold box", List->n_values); 00195 00196 for (i = 0; i < List->n_values; i++) { 00197 int pointb; 00198 double dx, dy, dist2; 00199 00200 pointb = List->value[i]; 00201 if (pointb == point) 00202 continue; 00203 00204 dx = XPnts[pointb].x - XPnts[point].x; 00205 dy = XPnts[pointb].y - XPnts[point].y; 00206 dist2 = dx * dx + dy * dy; 00207 00208 if (dist2 > thresh2) /* outside threshold */ 00209 continue; 00210 00211 /* doesn't have an anchor yet */ 00212 if (XPnts[pointb].anchor == -1) { 00213 XPnts[pointb].anchor = point; 00214 ntosnap++; 00215 } 00216 else if (XPnts[pointb].anchor > 0) { /* check distance to previously assigned anchor */ 00217 double dist2_a; 00218 00219 dx = XPnts[XPnts[pointb].anchor].x - XPnts[pointb].x; 00220 dy = XPnts[XPnts[pointb].anchor].y - XPnts[pointb].y; 00221 dist2_a = dx * dx + dy * dy; 00222 00223 /* replace old anchor */ 00224 if (dist2 < dist2_a) { 00225 XPnts[pointb].anchor = point; 00226 } 00227 } 00228 } 00229 } 00230 00231 /* Go through all lines and: 00232 * 1) for all vertices: if not anchor snap it to its anchor 00233 * 2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course) */ 00234 00235 nsnapped = ncreated = 0; 00236 00237 for (line_idx = 0; line_idx < List_lines->n_values; line_idx++) { 00238 int v, spoint, anchor; 00239 int changed = 0; 00240 00241 G_percent(line_idx, List_lines->n_values, 2); 00242 00243 line = List_lines->value[line_idx]; 00244 00245 G_debug(3, "line = %d", line); 00246 if (!Vect_line_alive(Map, line)) 00247 continue; 00248 00249 ltype = Vect_read_line(Map, Points, Cats, line); 00250 00251 if (Points->n_points >= aindex) { 00252 aindex = Points->n_points; 00253 Index = (int *)G_realloc(Index, aindex * sizeof(int)); 00254 } 00255 00256 /* Snap all vertices */ 00257 for (v = 0; v < Points->n_points; v++) { 00258 /* Box */ 00259 rect.boundary[0] = Points->x[v]; 00260 rect.boundary[3] = Points->x[v]; 00261 rect.boundary[1] = Points->y[v]; 00262 rect.boundary[4] = Points->y[v]; 00263 rect.boundary[2] = 0; 00264 rect.boundary[5] = 0; 00265 00266 /* Find point ( should always find one point ) */ 00267 Vect_reset_list(List); 00268 00269 RTreeSearch(RTree, &rect, (void *)add_item, List); 00270 00271 spoint = List->value[0]; 00272 anchor = XPnts[spoint].anchor; 00273 00274 if (anchor > 0) { /* to be snapped */ 00275 Points->x[v] = XPnts[anchor].x; 00276 Points->y[v] = XPnts[anchor].y; 00277 nsnapped++; 00278 changed = 1; 00279 Index[v] = anchor; /* point on new location */ 00280 } 00281 else { 00282 Index[v] = spoint; /* old point */ 00283 } 00284 } 00285 00286 /* New points */ 00287 Vect_reset_line(NPoints); 00288 00289 /* Snap all segments to anchors in threshold */ 00290 for (v = 0; v < Points->n_points - 1; v++) { 00291 int i; 00292 double x1, x2, y1, y2, xmin, xmax, ymin, ymax; 00293 00294 G_debug(3, " segment = %d end anchors : %d %d", v, Index[v], 00295 Index[v + 1]); 00296 00297 x1 = Points->x[v]; 00298 x2 = Points->x[v + 1]; 00299 y1 = Points->y[v]; 00300 y2 = Points->y[v + 1]; 00301 00302 Vect_append_point(NPoints, Points->x[v], Points->y[v], 00303 Points->z[v]); 00304 00305 /* Box */ 00306 if (x1 <= x2) { 00307 xmin = x1; 00308 xmax = x2; 00309 } 00310 else { 00311 xmin = x2; 00312 xmax = x1; 00313 } 00314 if (y1 <= y2) { 00315 ymin = y1; 00316 ymax = y2; 00317 } 00318 else { 00319 ymin = y2; 00320 ymax = y1; 00321 } 00322 00323 rect.boundary[0] = xmin - thresh; 00324 rect.boundary[3] = xmax + thresh; 00325 rect.boundary[1] = ymin - thresh; 00326 rect.boundary[4] = ymax + thresh; 00327 rect.boundary[2] = 0; 00328 rect.boundary[5] = 0; 00329 00330 /* Find points */ 00331 Vect_reset_list(List); 00332 RTreeSearch(RTree, &rect, (void *)add_item, List); 00333 00334 G_debug(3, " %d points in box", List->n_values); 00335 00336 /* Snap to anchor in threshold different from end points */ 00337 nnew = 0; 00338 for (i = 0; i < List->n_values; i++) { 00339 double dist2, along; 00340 00341 spoint = List->value[i]; 00342 G_debug(4, " spoint = %d anchor = %d", spoint, 00343 XPnts[spoint].anchor); 00344 00345 if (spoint == Index[v] || spoint == Index[v + 1]) 00346 continue; /* end point */ 00347 if (XPnts[spoint].anchor > 0) 00348 continue; /* point is not anchor */ 00349 00350 /* Check the distance */ 00351 dist2 = 00352 dig_distance2_point_to_line(XPnts[spoint].x, 00353 XPnts[spoint].y, 0, x1, y1, 0, 00354 x2, y2, 0, 0, NULL, NULL, 00355 NULL, &along, NULL); 00356 00357 G_debug(4, " distance = %lf", sqrt(dist2)); 00358 00359 if (dist2 <= thresh2) { 00360 G_debug(4, " anchor in thresh, along = %lf", along); 00361 00362 if (nnew == anew) { 00363 anew += 100; 00364 New = (NEW *) G_realloc(New, anew * sizeof(NEW)); 00365 } 00366 New[nnew].anchor = spoint; 00367 New[nnew].along = along; 00368 nnew++; 00369 } 00370 } 00371 G_debug(3, " nnew = %d", nnew); 00372 /* insert new vertices */ 00373 if (nnew > 0) { 00374 /* sort by distance along the segment */ 00375 qsort(New, sizeof(char) * nnew, sizeof(NEW), sort_new); 00376 00377 for (i = 0; i < nnew; i++) { 00378 anchor = New[i].anchor; 00379 /* Vect_line_insert_point ( Points, ++v, XPnts[anchor].x, XPnts[anchor].y, 0); */ 00380 Vect_append_point(NPoints, XPnts[anchor].x, 00381 XPnts[anchor].y, 0); 00382 ncreated++; 00383 } 00384 changed = 1; 00385 } 00386 } 00387 00388 /* append end point */ 00389 v = Points->n_points - 1; 00390 Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]); 00391 00392 if (changed) { /* rewrite the line */ 00393 Vect_line_prune(NPoints); /* remove duplicates */ 00394 if (NPoints->n_points > 1 || ltype & GV_LINES) { 00395 Vect_rewrite_line(Map, line, ltype, NPoints, Cats); 00396 } 00397 else { 00398 Vect_delete_line(Map, line); 00399 } 00400 if (Err) { 00401 Vect_write_line(Err, ltype, Points, Cats); 00402 } 00403 } 00404 } /* for each line */ 00405 G_percent(line_idx, List_lines->n_values, 2); /* finish it */ 00406 00407 Vect_destroy_line_struct(Points); 00408 Vect_destroy_line_struct(NPoints); 00409 Vect_destroy_cats_struct(Cats); 00410 G_free(XPnts); 00411 G_free(Index); 00412 G_free(New); 00413 RTreeDestroyNode(RTree); 00414 } 00415 00416 /* for qsort */ 00417 static int sort_new(const void *pa, const void *pb) 00418 { 00419 NEW *p1 = (NEW *) pa; 00420 NEW *p2 = (NEW *) pb; 00421 00422 if (p1->along < p2->along) 00423 return -1; 00424 if (p1->along > p2->along) 00425 return 1; 00426 return 1; 00427 } 00428 00441 void 00442 Vect_snap_lines(struct Map_info *Map, int type, double thresh, 00443 struct Map_info *Err) 00444 { 00445 int line, nlines, ltype; 00446 00447 struct ilist *List; 00448 00449 List = Vect_new_list(); 00450 00451 nlines = Vect_get_num_lines(Map); 00452 00453 for (line = 1; line <= nlines; line++) { 00454 G_debug(3, "line = %d", line); 00455 00456 if (!Vect_line_alive(Map, line)) 00457 continue; 00458 00459 ltype = Vect_read_line(Map, NULL, NULL, line); 00460 00461 if (!(ltype & type)) 00462 continue; 00463 00464 /* Vect_list_append(List, line); */ 00465 dig_list_add(List, line); 00466 } 00467 00468 Vect_snap_lines_list(Map, List, thresh, Err); 00469 00470 Vect_destroy_list(List); 00471 00472 return; 00473 }