00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <stdlib.h>
00018 #include <grass/gis.h>
00019 #include <grass/Vect.h>
00020
00034 int
00035 Vect_remove_small_areas ( struct Map_info *Map, double thresh, struct Map_info *Err, FILE *msgout,
00036 double *removed_area )
00037 {
00038 int area;
00039 int nremoved = 0;
00040 struct ilist *List;
00041 struct ilist *AList;
00042 struct line_pnts *Points;
00043 struct line_cats *Cats;
00044 double size_removed = 0.0;
00045
00046 List = Vect_new_list ();
00047 AList = Vect_new_list ();
00048 Points = Vect_new_line_struct ();
00049 Cats = Vect_new_cats_struct ();
00050
00051 if ( msgout ) fprintf (msgout, "Removed areas: %5d", nremoved );
00052
00053 for ( area = 1; area <= Vect_get_num_areas(Map); area++ ){
00054 int i, j, centroid, dissolve_neighbour;
00055 double length, size;
00056
00057 G_debug (3, "area = %d", area );
00058 if ( !Vect_area_alive ( Map, area ) ) continue;
00059
00060 size = Vect_get_area_area(Map,area);
00061 if ( size > thresh ) continue;
00062 size_removed += size;
00063
00064
00065
00066
00067 centroid = Vect_get_area_centroid ( Map, area );
00068 if ( centroid > 0 ) {
00069 if ( Err ) {
00070 Vect_read_line ( Map, Points, Cats, centroid );
00071 Vect_write_line ( Err, GV_CENTROID, Points, Cats );
00072 }
00073 Vect_delete_line (Map, centroid);
00074 }
00075
00076
00077
00078 Vect_get_area_boundaries ( Map, area, List );
00079
00080
00081 Vect_reset_list ( AList );
00082 for ( i = 0; i < List->n_values; i++ ) {
00083 int line, left, right, neighbour;
00084
00085 line = List->value[i];
00086
00087 if ( !Vect_line_alive(Map,abs(line)) )
00088 G_fatal_error ("Area is composed of dead boundary" );
00089
00090 Vect_get_line_areas ( Map, abs(line), &left, &right );
00091 if ( line > 0 ) neighbour = left;
00092 else neighbour = right;
00093
00094 G_debug (4, " line = %d left = %d right = %d neighbour = %d", line, left, right, neighbour );
00095
00096 Vect_list_append ( AList, neighbour );
00097 }
00098 G_debug (3, "num neighbours = %d", AList->n_values );
00099
00100
00101 dissolve_neighbour = 0;
00102 length = -1.0;
00103 for ( i = 0; i < AList->n_values; i++ ) {
00104 int neighbour1;
00105 double l = 0.0;
00106
00107 neighbour1 = AList->value[i];
00108 G_debug (4, " neighbour1 = %d", neighbour1 );
00109
00110 for ( j = 0; j < List->n_values; j++ ) {
00111 int line, left, right, neighbour2;
00112
00113 line = List->value[j];
00114 Vect_get_line_areas ( Map, abs(line), &left, &right );
00115 if ( line > 0 ) neighbour2= left;
00116 else neighbour2 = right;
00117
00118 if ( neighbour2 == neighbour1 ) {
00119 Vect_read_line ( Map, Points, NULL, abs(line) );
00120 l += Vect_line_length ( Points );
00121 }
00122 }
00123 if ( l > length ) {
00124 length = l;
00125 dissolve_neighbour = neighbour1;
00126 }
00127 }
00128
00129 G_debug (3, "dissolve_neighbour = %d", dissolve_neighbour );
00130
00131
00132 Vect_reset_list ( AList );
00133 for ( i = 0; i < List->n_values; i++ ) {
00134 int line, left, right, neighbour;
00135
00136 line = List->value[i];
00137 Vect_get_line_areas ( Map, abs(line), &left, &right );
00138 if ( line > 0 ) neighbour= left;
00139 else neighbour = right;
00140
00141 G_debug (3, " neighbour = %d", neighbour );
00142
00143 if ( neighbour == dissolve_neighbour ) {
00144 Vect_list_append ( AList, abs(line) );
00145 }
00146 }
00147
00148
00149 for ( i = 0; i < AList->n_values; i++ ) {
00150 int line;
00151
00152 line = AList->value[i];
00153
00154 if ( Err ) {
00155 Vect_read_line ( Map, Points, Cats, line );
00156 Vect_write_line ( Err, GV_BOUNDARY, Points, Cats );
00157 }
00158 Vect_delete_line ( Map, line );
00159 }
00160
00161 nremoved++;
00162 if ( msgout ) {
00163 fprintf (msgout, "\rRemoved areas: %5d", nremoved );
00164 fflush ( stderr );
00165 }
00166 }
00167 if ( msgout ) fprintf (stderr, "\n" );
00168
00169 if ( removed_area )
00170 *removed_area = size_removed;
00171
00172 return (nremoved);
00173 }
00174