GRASS Programmer's Manual  6.4.2(2012)
Vlib/snap.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines