GRASS Programmer's Manual
6.4.1(2011)
|
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 }