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