00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdlib.h>
00018 #include <math.h>
00019 #include <grass/gis.h>
00020 #include <grass/Vect.h>
00021
00022 typedef struct {
00023 double x, y;
00024 double a1, a2;
00025 char cross;
00026 char used;
00027
00028
00029 } XPNT;
00030
00031 static int fpoint;
00032
00033 void srch (int id, int *arg)
00034 {
00035 fpoint = id;
00036 }
00037
00056 void
00057 Vect_break_polygons ( struct Map_info *Map, int type, struct Map_info *Err, FILE *msgout )
00058 {
00059 struct line_pnts *BPoints, *Points;
00060 struct line_cats *Cats;
00061 int i, j, k, ret, ltype, broken, last, nlines;
00062 int nbreaks;
00063 struct Node *RTree;
00064 int apoints, npoints, nallpoints, nmarks;
00065 XPNT *XPnts;
00066 struct Rect rect;
00067 double dx, dy, a1=0, a2=0;
00068 int closed, last_point, cross;
00069 int printed;
00070
00071 RTree = RTreeNewIndex();
00072
00073 BPoints = Vect_new_line_struct ();
00074 Points = Vect_new_line_struct ();
00075 Cats = Vect_new_cats_struct ();
00076
00077 nlines = Vect_get_num_lines (Map);
00078
00079 G_debug (3, "nlines = %d", nlines );
00080
00081
00082
00083 apoints = 0;
00084 nmarks = 0;
00085 npoints = 1;
00086 nallpoints = 0;
00087 XPnts = NULL;
00088
00089 printed = 0;
00090 for ( i = 1; i <= nlines; i++ ){
00091 G_debug (3, "i = %d", i);
00092 if ( !Vect_line_alive ( Map, i ) ) continue;
00093
00094 ltype = Vect_read_line (Map, Points, Cats, i);
00095 if ( !(ltype & type) ) continue;
00096
00097
00098
00099 Vect_line_prune ( Points );
00100
00101
00102
00103
00104 last_point = Points->n_points-1;
00105 if ( Points->x[0] == Points->x[last_point] && Points->y[0] == Points->y[last_point] )
00106 closed = 1;
00107 else
00108 closed = 0;
00109
00110 for ( j = 0; j < Points->n_points; j++ ){
00111 G_debug (3, "j = %d", j);
00112 nallpoints++;
00113
00114 if ( j == last_point && closed ) continue;
00115
00116
00117 rect.boundary[0] = Points->x[j]; rect.boundary[3] = Points->x[j];
00118 rect.boundary[1] = Points->y[j]; rect.boundary[4] = Points->y[j];
00119 rect.boundary[2] = 0; rect.boundary[5] = 0;
00120
00121
00122 fpoint = -1;
00123 RTreeSearch(RTree, &rect, (void *)srch, 0);
00124 G_debug (3, "fpoint = %d", fpoint);
00125
00126 if ( Points->n_points <= 2 || (!closed && (j == 0 || j == last_point) ) ) {
00127 cross = 1;
00128 } else {
00129 cross = 0;
00130 if ( j == 0 && closed ) {
00131 dx = Points->x[last_point] - Points->x[0]; dy = Points->y[last_point] - Points->y[0];
00132 a1 = atan2 ( dy, dx );
00133 dx = Points->x[1] - Points->x[0]; dy = Points->y[1] - Points->y[0];
00134 a2 = atan2 ( dy, dx );
00135 } else {
00136 dx = Points->x[j-1] - Points->x[j]; dy = Points->y[j-1] - Points->y[j];
00137 a1 = atan2 ( dy, dx );
00138 dx = Points->x[j+1] - Points->x[j]; dy = Points->y[j+1] - Points->y[j];
00139 a2 = atan2 ( dy, dx );
00140 }
00141 }
00142
00143 if ( fpoint > 0 ) {
00144 if ( XPnts[fpoint].cross == 1 ) continue;
00145
00146
00147 if ( cross ) {
00148 XPnts[fpoint].cross = 1;
00149 nmarks++;
00150 } else {
00151 G_debug (3, "a1 = %f xa1 = %f a2 = %f xa2 = %f", a1, XPnts[fpoint].a1,
00152 a2, XPnts[fpoint].a2);
00153 if ( ( a1 == XPnts[fpoint].a1 && a2 == XPnts[fpoint].a2 ) ||
00154 ( a1 == XPnts[fpoint].a2 && a2 == XPnts[fpoint].a1 ) )
00155 {
00156
00157 } else {
00158 XPnts[fpoint].cross = 1;
00159 nmarks++;
00160 }
00161 }
00162 } else {
00163
00164 RTreeInsertRect( &rect, npoints, &RTree, 0);
00165 if ( npoints >= apoints ) {
00166 apoints += 10000;
00167 XPnts = (XPNT *) G_realloc ( XPnts, (apoints + 1) * sizeof (XPNT) );
00168 }
00169 XPnts[npoints].x = Points->x[j];
00170 XPnts[npoints].y = Points->y[j];
00171 XPnts[npoints].used = 0;
00172 if ( j == 0 || j == (Points->n_points - 1) || Points->n_points < 3 ) {
00173 XPnts[npoints].a1 = 0;
00174 XPnts[npoints].a2 = 0;
00175 XPnts[npoints].cross = 1;
00176 nmarks++;
00177 } else {
00178 XPnts[npoints].a1 = a1;
00179 XPnts[npoints].a2 = a2;
00180 XPnts[npoints].cross = 0;
00181 }
00182
00183 npoints++;
00184 }
00185 }
00186 if ( msgout && printed > 5000 ) {
00187 fprintf (msgout, "\rRegistering points ... %d", npoints - 1);
00188 fflush ( msgout );
00189 printed = 0;
00190 }
00191 printed++;
00192 }
00193 if ( msgout ) {
00194 fprintf (msgout, "\rRegistering points ... %d", npoints - 1 );
00195 fprintf ( msgout, "\n" );
00196 fprintf ( msgout, "All points (vertices): %5d\n", nallpoints );
00197 fprintf ( msgout, "Registered points (unique coordinates): %5d\n", npoints - 1 );
00198 fprintf ( msgout, "Points marked for break: %5d\n", nmarks );
00199 }
00200
00201
00202
00203 nbreaks = 0;
00204 if ( msgout ) fprintf (msgout, "Breaks: %5d", nbreaks );
00205 printed = 0;
00206
00207
00208
00209 for ( i = 1; i <= nlines; i++ ){
00210 int n_orig_points;
00211
00212 G_debug (3, "i = %d", i);
00213 if ( !Vect_line_alive ( Map, i ) ) continue;
00214
00215 ltype = Vect_read_line (Map, Points, Cats, i);
00216 if ( !(ltype & type) ) continue;
00217 if ( !(ltype & GV_LINES ) ) continue;
00218
00219
00220 n_orig_points = Points->n_points;
00221 Vect_line_prune ( Points );
00222
00223 broken = 0; last = 0;
00224 G_debug (3, "n_points = %d", Points->n_points);
00225 for ( j = 1; j < Points->n_points; j++ ){
00226 G_debug (3, "j = %d", j);
00227 nallpoints++;
00228
00229
00230 rect.boundary[0] = Points->x[j]; rect.boundary[3] = Points->x[j];
00231 rect.boundary[1] = Points->y[j]; rect.boundary[4] = Points->y[j];
00232 rect.boundary[2] = 0; rect.boundary[5] = 0;
00233
00234 if ( Points->n_points <= 1 || (j == (Points->n_points-1) && !broken) ) break;
00235
00236
00237
00238 RTreeSearch(RTree, &rect, (void *)srch, 0);
00239 G_debug (3, "fpoint = %d", fpoint);
00240
00241 if ( XPnts[fpoint].cross ) {
00242 XPnts[fpoint].used = 1;
00243 }
00244
00245
00246 if ( (j == (Points->n_points-1) && broken ) || XPnts[fpoint].cross ) {
00247 Vect_reset_line ( BPoints );
00248 for ( k = last; k <= j; k++ ){
00249 Vect_append_point ( BPoints, Points->x[k], Points->y[k], Points->z[k] );
00250 }
00251
00252
00253 Vect_line_prune ( BPoints );
00254 if ( BPoints->n_points > 1 ) {
00255 ret = Vect_write_line ( Map, ltype, BPoints, Cats );
00256 G_debug (3, "Line %d written j = %d n_points(orig,pruned) = %d n_points(new) = %d",
00257 ret, j, Points->n_points, BPoints->n_points);
00258 }
00259
00260 if ( !broken ) Vect_delete_line (Map, i);
00261
00262 last = j;
00263 broken = 1;
00264 nbreaks++;
00265 }
00266
00267 if ( msgout && printed > 1000 ) {
00268 fprintf ( msgout, "\rBreaks: %5d", nbreaks );
00269 fflush ( msgout );
00270 printed = 0;
00271 }
00272 printed++;
00273 }
00274 if ( !broken && n_orig_points > Points->n_points ) {
00275 if ( Points->n_points > 1 ) {
00276 Vect_rewrite_line ( Map, i, ltype, Points, Cats );
00277 G_debug (3, "Line %d pruned, npoints = %d", i, Points->n_points );
00278 } else {
00279 Vect_delete_line (Map, i);
00280 G_debug (3, "Line %d was deleted", i);
00281 }
00282 } else {
00283 G_debug (3, "Line %d was not changed", i);
00284 }
00285 }
00286
00287 if ( msgout )
00288 fprintf ( msgout, "\rBreaks: %5d", nbreaks );
00289
00290
00291 if ( Err ) {
00292 Vect_reset_cats ( Cats );
00293 for ( i = 1; i < npoints; i++ ){
00294 if ( XPnts[i].used ) {
00295 Vect_reset_line ( Points );
00296 Vect_append_point ( Points, XPnts[i].x, XPnts[i].y, 0 );
00297 Vect_write_line ( Err, GV_POINT, Points, Cats );
00298 }
00299 }
00300 }
00301
00302 if ( msgout ) fprintf ( msgout, "\n" );
00303
00304 G_free ( XPnts );
00305 RTreeDestroyNode ( RTree);
00306 }
00307