remove_areas.c

Go to the documentation of this file.
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         /* The area is smaller than the limit -> remove */
00069 
00070         /* Remove centroid */
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         /* Find the adjacent area with which the longest boundary is shared */
00081 
00082         Vect_get_area_boundaries(Map, area, List);
00083 
00084         /* Create a list of neighbour areas */
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)))       /* Should not happen */
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); /* this checks for duplicity */
00104         }
00105         G_debug(3, "num neighbours = %d", AList->n_values);
00106 
00107         /* Go through the list of neghours and find that with the longest boundary */
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         /* Make list of boundaries to be removed */
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         /* Remove boundaries */
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 }

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