GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /************************************************************** 00003 * G_histogram_eq (histo, map, min, max) 00004 * 00005 * struct Histogram *histo; histogram as returned by G_read_histogram() 00006 * unsigned char **map; equalized category mapping 00007 * CELL *min, *max; min,max category for map 00008 * 00009 * perform histogram equalization 00010 * inputs are histo, output is map,min,max 00011 ****************************************************************/ 00012 #include <grass/gis.h> 00013 int G_histogram_eq(const struct Histogram *histo, 00014 unsigned char **map, CELL * min, CELL * max) 00015 { 00016 int i; 00017 int x; 00018 CELL cat, prev; 00019 double total; 00020 double sum; 00021 double span; 00022 int ncats; 00023 long count; 00024 unsigned char *xmap; 00025 int len; 00026 int first, last; 00027 00028 ncats = G_get_histogram_num(histo); 00029 if (ncats == 1) { 00030 *min = *max = G_get_histogram_cat(0, histo); 00031 *map = xmap = (unsigned char *)G_malloc(1); 00032 *xmap = 0; 00033 return 0; 00034 } 00035 if ((*min = G_get_histogram_cat(first = 0, histo)) == 0) 00036 *min = G_get_histogram_cat(++first, histo); 00037 if ((*max = G_get_histogram_cat(last = ncats - 1, histo)) == 0) 00038 *max = G_get_histogram_cat(--last, histo); 00039 len = *max - *min + 1; 00040 *map = xmap = (unsigned char *)G_malloc(len); 00041 00042 total = 0; 00043 for (i = first; i <= last; i++) { 00044 if (G_get_histogram_cat(i, histo) == 0) 00045 continue; 00046 count = G_get_histogram_count(i, histo); 00047 if (count > 0) 00048 total += count; 00049 } 00050 if (total <= 0) { 00051 for (i = 0; i < len; i++) 00052 xmap[i] = 0; 00053 return 0; 00054 } 00055 00056 span = total / 256; 00057 00058 sum = 0.0; 00059 cat = *min - 1; 00060 for (i = first; i <= last; i++) { 00061 prev = cat + 1; 00062 cat = G_get_histogram_cat(i, histo); 00063 count = G_get_histogram_count(i, histo); 00064 if (count < 0 || cat == 0) 00065 count = 0; 00066 x = (sum + (count / 2.0)) / span; 00067 if (x < 0) 00068 x = 0; 00069 else if (x > 255) 00070 x = 255; 00071 sum += count; 00072 00073 while (prev++ <= cat) 00074 *xmap++ = x; 00075 } 00076 00077 return 0; 00078 }