remove_areas.c

Go to the documentation of this file.
00001 /* **************************************************************
00002  * 
00003  * MODULE:       vector library
00004  * 
00005  * AUTHOR(S):    Radim Blazek
00006  *               
00007  * PURPOSE:      Clean areas
00008  *               
00009  * COPYRIGHT:    (C) 2001 by the GRASS Development Team
00010  *
00011  *               This program is free software under the 
00012  *               GNU General Public License (>=v2). 
00013  *               Read the file COPYING that comes with GRASS
00014  *               for details.
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         /* The area is smaller than the limit -> remove */
00065 
00066         /* Remove centroid */
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         /* Find the adjacent area with which the longest boundary is shared */
00077         
00078         Vect_get_area_boundaries ( Map, area, List );
00079 
00080         /* Create a list of neighbour areas */
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)) ) /* Should not happen */
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 ); /* this checks for duplicity */
00097         }
00098         G_debug (3, "num neighbours = %d", AList->n_values );
00099 
00100         /* Go through the list of neghours and find that with the longest boundary */
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         /* Make list of boundaries to be removed */
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         /* Remove boundaries */
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 

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5