GRASS Programmer's Manual  6.4.2(2012)
break_polygons.c
Go to the documentation of this file.
00001 
00020 #include <stdlib.h>
00021 #include <math.h>
00022 #include <grass/gis.h>
00023 #include <grass/Vect.h>
00024 #include <grass/glocale.h>
00025 
00026 typedef struct
00027 {
00028     double x, y;
00029     double a1, a2;              /* angles */
00030     char cross;                 /* 0 - do not break, 1 - break */
00031     char used;                  /* 0 - was not used to break line, 1 - was used to break line
00032                                  *   this is stored because points are automaticaly marked as cross, even if not used 
00033                                  *   later to break lines */
00034 } XPNT;
00035 
00036 static int fpoint;
00037 
00038 /* Function called from RTreeSearch for point found */
00039 void srch(int id, int *arg)
00040 {
00041     fpoint = id;
00042 }
00043 
00061 void
00062 Vect_break_polygons(struct Map_info *Map, int type, struct Map_info *Err)
00063 {
00064     struct line_pnts *BPoints, *Points;
00065     struct line_cats *Cats;
00066     int i, j, k, ret, ltype, broken, last, nlines;
00067     int nbreaks;
00068     struct Node *RTree;
00069     int apoints, npoints, nallpoints, nmarks;
00070     XPNT *XPnts;
00071     struct Rect rect;
00072     double dx, dy, a1 = 0, a2 = 0;
00073     int closed, last_point, cross;
00074 
00075     RTree = RTreeNewIndex();
00076 
00077     BPoints = Vect_new_line_struct();
00078     Points = Vect_new_line_struct();
00079     Cats = Vect_new_cats_struct();
00080 
00081     nlines = Vect_get_num_lines(Map);
00082 
00083     G_debug(3, "nlines =  %d", nlines);
00084     /* Go through all lines in vector, and add each point to structure of points,
00085      * if such point already exists check angles of segments and if differ mark for break */
00086 
00087     apoints = 0;
00088     nmarks = 0;
00089     npoints = 1;                /* index starts from 1 ! */
00090     nallpoints = 0;
00091     XPnts = NULL;
00092 
00093     for (i = 1; i <= nlines; i++) {
00094         G_percent(i, nlines, 1);
00095         G_debug(3, "i =  %d", i);
00096         if (!Vect_line_alive(Map, i))
00097             continue;
00098 
00099         ltype = Vect_read_line(Map, Points, Cats, i);
00100         if (!(ltype & type))
00101             continue;
00102 
00103         /* This would be confused by duplicate coordinates (angle cannot be calculated) ->
00104          * prune line first */
00105         Vect_line_prune(Points);
00106 
00107         /* If first and last point are identical it is close polygon, we dont need to register last point
00108          * and we can calculate angle for first.
00109          * If first and last point are not identical we have to mark for break both */
00110         last_point = Points->n_points - 1;
00111         if (Points->x[0] == Points->x[last_point] &&
00112             Points->y[0] == Points->y[last_point])
00113             closed = 1;
00114         else
00115             closed = 0;
00116 
00117         for (j = 0; j < Points->n_points; j++) {
00118             G_debug(3, "j =  %d", j);
00119             nallpoints++;
00120 
00121             if (j == last_point && closed)
00122                 continue;       /* do not register last of close polygon */
00123 
00124             /* Box */
00125             rect.boundary[0] = Points->x[j];
00126             rect.boundary[3] = Points->x[j];
00127             rect.boundary[1] = Points->y[j];
00128             rect.boundary[4] = Points->y[j];
00129             rect.boundary[2] = 0;
00130             rect.boundary[5] = 0;
00131 
00132             /* Already in DB? */
00133             fpoint = -1;
00134             RTreeSearch(RTree, &rect, (void *)srch, 0);
00135             G_debug(3, "fpoint =  %d", fpoint);
00136 
00137             if (Points->n_points <= 2 ||
00138                 (!closed && (j == 0 || j == last_point))) {
00139                 cross = 1;      /* mark for cross in any case */
00140             }
00141             else {              /* calculate angles */
00142                 cross = 0;
00143                 if (j == 0 && closed) { /* closed polygon */
00144                     dx = Points->x[last_point] - Points->x[0];
00145                     dy = Points->y[last_point] - Points->y[0];
00146                     a1 = atan2(dy, dx);
00147                     dx = Points->x[1] - Points->x[0];
00148                     dy = Points->y[1] - Points->y[0];
00149                     a2 = atan2(dy, dx);
00150                 }
00151                 else {
00152                     dx = Points->x[j - 1] - Points->x[j];
00153                     dy = Points->y[j - 1] - Points->y[j];
00154                     a1 = atan2(dy, dx);
00155                     dx = Points->x[j + 1] - Points->x[j];
00156                     dy = Points->y[j + 1] - Points->y[j];
00157                     a2 = atan2(dy, dx);
00158                 }
00159             }
00160 
00161             if (fpoint > 0) {   /* Found */
00162                 if (XPnts[fpoint].cross == 1)
00163                     continue;   /* already marked */
00164 
00165                 /* Check angles */
00166                 if (cross) {
00167                     XPnts[fpoint].cross = 1;
00168                     nmarks++;
00169                 }
00170                 else {
00171                     G_debug(3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1,
00172                             XPnts[fpoint].a1, a2, XPnts[fpoint].a2);
00173                     if ((a1 == XPnts[fpoint].a1 && a2 == XPnts[fpoint].a2) || (a1 == XPnts[fpoint].a2 && a2 == XPnts[fpoint].a1)) {     /* identical */
00174 
00175                     }
00176                     else {
00177                         XPnts[fpoint].cross = 1;
00178                         nmarks++;
00179                     }
00180                 }
00181             }
00182             else {
00183                 /* Add to tree and to structure */
00184                 RTreeInsertRect(&rect, npoints, &RTree, 0);
00185                 if (npoints >= apoints) {
00186                     apoints += 10000;
00187                     XPnts =
00188                         (XPNT *) G_realloc(XPnts,
00189                                            (apoints + 1) * sizeof(XPNT));
00190                 }
00191                 XPnts[npoints].x = Points->x[j];
00192                 XPnts[npoints].y = Points->y[j];
00193                 XPnts[npoints].used = 0;
00194                 if (j == 0 || j == (Points->n_points - 1) ||
00195                     Points->n_points < 3) {
00196                     XPnts[npoints].a1 = 0;
00197                     XPnts[npoints].a2 = 0;
00198                     XPnts[npoints].cross = 1;
00199                     nmarks++;
00200                 }
00201                 else {
00202                     XPnts[npoints].a1 = a1;
00203                     XPnts[npoints].a2 = a2;
00204                     XPnts[npoints].cross = 0;
00205                 }
00206 
00207                 npoints++;
00208             }
00209         }
00210     }
00211 
00212     /* G_sleep (10); */
00213 
00214     nbreaks = 0;
00215 
00216     /* Second loop through lines (existing when loop is started, no need to process lines written again)
00217      * and break at points marked for break */
00218     for (i = 1; i <= nlines; i++) {
00219         int n_orig_points;
00220 
00221         G_percent(i, nlines, 1);
00222         G_debug(3, "i =  %d", i);
00223         if (!Vect_line_alive(Map, i))
00224             continue;
00225 
00226         ltype = Vect_read_line(Map, Points, Cats, i);
00227         if (!(ltype & type))
00228             continue;
00229         if (!(ltype & GV_LINES))
00230             continue;           /* Nonsense to break points */
00231 
00232         /* Duplicates would result in zero length lines -> prune line first */
00233         n_orig_points = Points->n_points;
00234         Vect_line_prune(Points);
00235 
00236         broken = 0;
00237         last = 0;
00238         G_debug(3, "n_points =  %d", Points->n_points);
00239         for (j = 1; j < Points->n_points; j++) {
00240             G_debug(3, "j =  %d", j);
00241             nallpoints++;
00242 
00243             /* Box */
00244             rect.boundary[0] = Points->x[j];
00245             rect.boundary[3] = Points->x[j];
00246             rect.boundary[1] = Points->y[j];
00247             rect.boundary[4] = Points->y[j];
00248             rect.boundary[2] = 0;
00249             rect.boundary[5] = 0;
00250 
00251             if (Points->n_points <= 1 ||
00252                 (j == (Points->n_points - 1) && !broken))
00253                 break;
00254             /* One point only or 
00255              * last point and line is not broken, do nothing */
00256 
00257             RTreeSearch(RTree, &rect, (void *)srch, 0);
00258             G_debug(3, "fpoint =  %d", fpoint);
00259 
00260             if (XPnts[fpoint].cross) {  /* realy use to break line */
00261                 XPnts[fpoint].used = 1;
00262             }
00263 
00264             /* break or write last segment of broken line */
00265             if ((j == (Points->n_points - 1) && broken) ||
00266                 XPnts[fpoint].cross) {
00267                 Vect_reset_line(BPoints);
00268                 for (k = last; k <= j; k++) {
00269                     Vect_append_point(BPoints, Points->x[k], Points->y[k],
00270                                       Points->z[k]);
00271                 }
00272 
00273                 /* Result may collapse to one point */
00274                 Vect_line_prune(BPoints);
00275                 if (BPoints->n_points > 1) {
00276                     ret = Vect_write_line(Map, ltype, BPoints, Cats);
00277                     G_debug(3,
00278                             "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
00279                             ret, j, Points->n_points, BPoints->n_points);
00280                 }
00281 
00282                 if (!broken)
00283                     Vect_delete_line(Map, i);   /* not yet deleted */
00284 
00285                 last = j;
00286                 broken = 1;
00287                 nbreaks++;
00288             }
00289         }
00290         if (!broken && n_orig_points > Points->n_points) {      /* was pruned before -> rewrite */
00291             if (Points->n_points > 1) {
00292                 Vect_rewrite_line(Map, i, ltype, Points, Cats);
00293                 G_debug(3, "Line %d pruned, npoints = %d", i,
00294                         Points->n_points);
00295             }
00296             else {
00297                 Vect_delete_line(Map, i);
00298                 G_debug(3, "Line %d was deleted", i);
00299             }
00300         }
00301         else {
00302             G_debug(3, "Line %d was not changed", i);
00303         }
00304     }
00305 
00306     /* Write points on breaks */
00307     if (Err) {
00308         Vect_reset_cats(Cats);
00309         for (i = 1; i < npoints; i++) {
00310             if (XPnts[i].used) {
00311                 Vect_reset_line(Points);
00312                 Vect_append_point(Points, XPnts[i].x, XPnts[i].y, 0);
00313                 Vect_write_line(Err, GV_POINT, Points, Cats);
00314             }
00315         }
00316     }
00317 
00318     G_free(XPnts);
00319     RTreeDestroyNode(RTree);
00320 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines