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_debug(3, "i =  %d", i);
00095         if (!Vect_line_alive(Map, i))
00096             continue;
00097 
00098         ltype = Vect_read_line(Map, Points, Cats, i);
00099         if (!(ltype & type))
00100             continue;
00101 
00102         /* This would be confused by duplicate coordinates (angle cannot be calculated) ->
00103          * prune line first */
00104         Vect_line_prune(Points);
00105 
00106         /* If first and last point are identical it is close polygon, we dont need to register last point
00107          * and we can calculate angle for first.
00108          * If first and last point are not identical we have to mark for break both */
00109         last_point = Points->n_points - 1;
00110         if (Points->x[0] == Points->x[last_point] &&
00111             Points->y[0] == Points->y[last_point])
00112             closed = 1;
00113         else
00114             closed = 0;
00115 
00116         for (j = 0; j < Points->n_points; j++) {
00117             G_debug(3, "j =  %d", j);
00118             nallpoints++;
00119 
00120             if (j == last_point && closed)
00121                 continue;       /* do not register last of close polygon */
00122 
00123             /* Box */
00124             rect.boundary[0] = Points->x[j];
00125             rect.boundary[3] = Points->x[j];
00126             rect.boundary[1] = Points->y[j];
00127             rect.boundary[4] = Points->y[j];
00128             rect.boundary[2] = 0;
00129             rect.boundary[5] = 0;
00130 
00131             /* Already in DB? */
00132             fpoint = -1;
00133             RTreeSearch(RTree, &rect, (void *)srch, 0);
00134             G_debug(3, "fpoint =  %d", fpoint);
00135 
00136             if (Points->n_points <= 2 ||
00137                 (!closed && (j == 0 || j == last_point))) {
00138                 cross = 1;      /* mark for cross in any case */
00139             }
00140             else {              /* calculate angles */
00141                 cross = 0;
00142                 if (j == 0 && closed) { /* closed polygon */
00143                     dx = Points->x[last_point] - Points->x[0];
00144                     dy = Points->y[last_point] - Points->y[0];
00145                     a1 = atan2(dy, dx);
00146                     dx = Points->x[1] - Points->x[0];
00147                     dy = Points->y[1] - Points->y[0];
00148                     a2 = atan2(dy, dx);
00149                 }
00150                 else {
00151                     dx = Points->x[j - 1] - Points->x[j];
00152                     dy = Points->y[j - 1] - Points->y[j];
00153                     a1 = atan2(dy, dx);
00154                     dx = Points->x[j + 1] - Points->x[j];
00155                     dy = Points->y[j + 1] - Points->y[j];
00156                     a2 = atan2(dy, dx);
00157                 }
00158             }
00159 
00160             if (fpoint > 0) {   /* Found */
00161                 if (XPnts[fpoint].cross == 1)
00162                     continue;   /* already marked */
00163 
00164                 /* Check angles */
00165                 if (cross) {
00166                     XPnts[fpoint].cross = 1;
00167                     nmarks++;
00168                 }
00169                 else {
00170                     G_debug(3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1,
00171                             XPnts[fpoint].a1, a2, XPnts[fpoint].a2);
00172                     if ((a1 == XPnts[fpoint].a1 && a2 == XPnts[fpoint].a2) || (a1 == XPnts[fpoint].a2 && a2 == XPnts[fpoint].a1)) {     /* identical */
00173 
00174                     }
00175                     else {
00176                         XPnts[fpoint].cross = 1;
00177                         nmarks++;
00178                     }
00179                 }
00180             }
00181             else {
00182                 /* Add to tree and to structure */
00183                 RTreeInsertRect(&rect, npoints, &RTree, 0);
00184                 if (npoints >= apoints) {
00185                     apoints += 10000;
00186                     XPnts =
00187                         (XPNT *) G_realloc(XPnts,
00188                                            (apoints + 1) * sizeof(XPNT));
00189                 }
00190                 XPnts[npoints].x = Points->x[j];
00191                 XPnts[npoints].y = Points->y[j];
00192                 XPnts[npoints].used = 0;
00193                 if (j == 0 || j == (Points->n_points - 1) ||
00194                     Points->n_points < 3) {
00195                     XPnts[npoints].a1 = 0;
00196                     XPnts[npoints].a2 = 0;
00197                     XPnts[npoints].cross = 1;
00198                     nmarks++;
00199                 }
00200                 else {
00201                     XPnts[npoints].a1 = a1;
00202                     XPnts[npoints].a2 = a2;
00203                     XPnts[npoints].cross = 0;
00204                 }
00205 
00206                 npoints++;
00207             }
00208         }
00209     }
00210 
00211     /* G_sleep (10); */
00212 
00213     nbreaks = 0;
00214 
00215     /* Second loop through lines (existing when loop is started, no need to process lines written again)
00216      * and break at points marked for break */
00217     for (i = 1; i <= nlines; i++) {
00218         int n_orig_points;
00219 
00220         G_debug(3, "i =  %d", i);
00221         if (!Vect_line_alive(Map, i))
00222             continue;
00223 
00224         ltype = Vect_read_line(Map, Points, Cats, i);
00225         if (!(ltype & type))
00226             continue;
00227         if (!(ltype & GV_LINES))
00228             continue;           /* Nonsense to break points */
00229 
00230         /* Duplicates would result in zero length lines -> prune line first */
00231         n_orig_points = Points->n_points;
00232         Vect_line_prune(Points);
00233 
00234         broken = 0;
00235         last = 0;
00236         G_debug(3, "n_points =  %d", Points->n_points);
00237         for (j = 1; j < Points->n_points; j++) {
00238             G_debug(3, "j =  %d", j);
00239             nallpoints++;
00240 
00241             /* Box */
00242             rect.boundary[0] = Points->x[j];
00243             rect.boundary[3] = Points->x[j];
00244             rect.boundary[1] = Points->y[j];
00245             rect.boundary[4] = Points->y[j];
00246             rect.boundary[2] = 0;
00247             rect.boundary[5] = 0;
00248 
00249             if (Points->n_points <= 1 ||
00250                 (j == (Points->n_points - 1) && !broken))
00251                 break;
00252             /* One point only or 
00253              * last point and line is not broken, do nothing */
00254 
00255             RTreeSearch(RTree, &rect, (void *)srch, 0);
00256             G_debug(3, "fpoint =  %d", fpoint);
00257 
00258             if (XPnts[fpoint].cross) {  /* realy use to break line */
00259                 XPnts[fpoint].used = 1;
00260             }
00261 
00262             /* break or write last segment of broken line */
00263             if ((j == (Points->n_points - 1) && broken) ||
00264                 XPnts[fpoint].cross) {
00265                 Vect_reset_line(BPoints);
00266                 for (k = last; k <= j; k++) {
00267                     Vect_append_point(BPoints, Points->x[k], Points->y[k],
00268                                       Points->z[k]);
00269                 }
00270 
00271                 /* Result may collapse to one point */
00272                 Vect_line_prune(BPoints);
00273                 if (BPoints->n_points > 1) {
00274                     ret = Vect_write_line(Map, ltype, BPoints, Cats);
00275                     G_debug(3,
00276                             "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
00277                             ret, j, Points->n_points, BPoints->n_points);
00278                 }
00279 
00280                 if (!broken)
00281                     Vect_delete_line(Map, i);   /* not yet deleted */
00282 
00283                 last = j;
00284                 broken = 1;
00285                 nbreaks++;
00286             }
00287         }
00288         if (!broken && n_orig_points > Points->n_points) {      /* was pruned before -> rewrite */
00289             if (Points->n_points > 1) {
00290                 Vect_rewrite_line(Map, i, ltype, Points, Cats);
00291                 G_debug(3, "Line %d pruned, npoints = %d", i,
00292                         Points->n_points);
00293             }
00294             else {
00295                 Vect_delete_line(Map, i);
00296                 G_debug(3, "Line %d was deleted", i);
00297             }
00298         }
00299         else {
00300             G_debug(3, "Line %d was not changed", i);
00301         }
00302     }
00303 
00304     /* Write points on breaks */
00305     if (Err) {
00306         Vect_reset_cats(Cats);
00307         for (i = 1; i < npoints; i++) {
00308             if (XPnts[i].used) {
00309                 Vect_reset_line(Points);
00310                 Vect_append_point(Points, XPnts[i].x, XPnts[i].y, 0);
00311                 Vect_write_line(Err, GV_POINT, Points, Cats);
00312             }
00313         }
00314     }
00315 
00316     G_free(XPnts);
00317     RTreeDestroyNode(RTree);
00318 }

Generated on Sat Oct 24 03:25:19 2009 for GRASS Programmer's Manual by  doxygen 1.6.1