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         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             /* Box */
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             /* Already registered ? */
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) {  /* Not found */
00147                 /* Add to tree and to structure */
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     /* Go through all registered points and if not yet marked mark it as anchor and assign this anchor
00165      * to all not yet marked points in threshold */
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;        /* make it anchor */
00176         nanchors++;
00177 
00178         /* Find points in threshold */
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     /* Go through all lines and: 
00210      *   1) for all vertices: if not anchor snap it to its anchor
00211      *   2) for all segments: snap it to all anchors in threshold (except anchors of vertices of course) */
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         /* Snap all vertices */
00233         for (v = 0; v < Points->n_points; v++) {
00234             /* Box */
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             /* Find point ( should always find one point ) */
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) {   /* to be snapped */
00251                 Points->x[v] = XPnts[anchor].x;
00252                 Points->y[v] = XPnts[anchor].y;
00253                 nsnapped++;
00254                 changed = 1;
00255                 Index[v] = anchor;      /* point on new location */
00256             }
00257             else {
00258                 Index[v] = spoint;      /* old point */
00259             }
00260         }
00261 
00262         /* New points */
00263         Vect_reset_line(NPoints);
00264 
00265         /* Snap all segments to anchors in threshold */
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             /* Box */
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             /* Find points */
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             /* Snap to anchor in threshold different from end points */
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;   /* end point */
00323                 if (XPnts[spoint].anchor > 0)
00324                     continue;   /* point is not anchor */
00325 
00326                 /* Check the distance */
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             /* insert new vertices */
00349             if (nnew > 0) {
00350                 /* sort by distance along the segment */
00351                 qsort(New, sizeof(char) * nnew, sizeof(NEW), sort_new);
00352 
00353                 for (i = 0; i < nnew; i++) {
00354                     anchor = New[i].anchor;
00355                     /* Vect_line_insert_point ( Points, ++v, XPnts[anchor].x, XPnts[anchor].y, 0); */
00356                     Vect_append_point(NPoints, XPnts[anchor].x,
00357                                       XPnts[anchor].y, 0);
00358                     ncreated++;
00359                 }
00360                 changed = 1;
00361             }
00362         }
00363 
00364         /* append end point */
00365         v = Points->n_points - 1;
00366         Vect_append_point(NPoints, Points->x[v], Points->y[v], Points->z[v]);
00367 
00368         if (changed) {          /* rewrite the line */
00369             Vect_line_prune(NPoints);   /* remove duplicates */
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     }                           /* for each line */
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 /* for qsort */
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 }

Generated on Thu Jul 16 13:21:17 2009 for GRASS Programmer's Manual by  doxygen 1.5.6