1
2 /**************************************************************
3 * Rast_histogram_eq (histo, map, min, max)
4 *
5 * struct Histogram *histo; histogram as returned by Rast_read_histogram()
6 * unsigned char **map; equalized category mapping
7 * CELL *min, *max; min,max category for map
8 *
9 * perform histogram equalization
10 * inputs are histo, output is map,min,max
11 ****************************************************************/
12 #include <grass/gis.h>
13 #include <grass/raster.h>
14
Rast_histogram_eq(const struct Histogram * histo,unsigned char ** map,CELL * min,CELL * max)15 void Rast_histogram_eq(const struct Histogram *histo,
16 unsigned char **map, CELL * min, CELL * max)
17 {
18 int i;
19 int x;
20 CELL cat, prev;
21 double total;
22 double sum;
23 double span;
24 int ncats;
25 long count;
26 unsigned char *xmap;
27 int len;
28 int first, last;
29
30 ncats = Rast_get_histogram_num(histo);
31 if (ncats == 1) {
32 *min = *max = Rast_get_histogram_cat(0, histo);
33 *map = xmap = (unsigned char *)G_malloc(1);
34 *xmap = 0;
35 return;
36 }
37 if ((*min = Rast_get_histogram_cat(first = 0, histo)) == 0)
38 *min = Rast_get_histogram_cat(++first, histo);
39 if ((*max = Rast_get_histogram_cat(last = ncats - 1, histo)) == 0)
40 *max = Rast_get_histogram_cat(--last, histo);
41 len = *max - *min + 1;
42 *map = xmap = (unsigned char *)G_malloc(len);
43
44 total = 0;
45 for (i = first; i <= last; i++) {
46 if (Rast_get_histogram_cat(i, histo) == 0)
47 continue;
48 count = Rast_get_histogram_count(i, histo);
49 if (count > 0)
50 total += count;
51 }
52 if (total <= 0) {
53 for (i = 0; i < len; i++)
54 xmap[i] = 0;
55 return;
56 }
57
58 span = total / 256;
59
60 sum = 0.0;
61 cat = *min - 1;
62 for (i = first; i <= last; i++) {
63 prev = cat + 1;
64 cat = Rast_get_histogram_cat(i, histo);
65 count = Rast_get_histogram_count(i, histo);
66 if (count < 0 || cat == 0)
67 count = 0;
68 x = (sum + (count / 2.0)) / span;
69 if (x < 0)
70 x = 0;
71 else if (x > 255)
72 x = 255;
73 sum += count;
74
75 while (prev++ <= cat)
76 *xmap++ = x;
77 }
78 }
79