1 #include "vil_histogram_equalise.h"
2 //:
3 // \file
4 // \brief Apply histogram equalisation to given image
5 // \author Tim Cootes
6
7 #include <vil/algo/vil_histogram.h>
8 #include "vil/vil_math.h"
9
10
11 //: Apply histogram equalisation to given byte image
12 void
vil_histogram_equalise(vil_image_view<vxl_byte> & image)13 vil_histogram_equalise(vil_image_view<vxl_byte> & image)
14 {
15 std::vector<double> histo(256);
16 vil_histogram_byte(image, histo);
17
18 // Create cumulative frequency curve
19 double sum = 0.0;
20 for (unsigned i = 0; i < 256; ++i)
21 {
22 sum += histo[i];
23 histo[i] = sum;
24 }
25
26 // Parameters of mapping
27 int lo = 0;
28 // Find smallest value in image
29 while (histo[lo] == 0)
30 lo++;
31 double x0 = histo[lo];
32 double s = 255.1 / (sum - x0); // Smallest values get mapped to zero
33
34 std::vector<vxl_byte> lookup(256);
35 vxl_byte * lup = &lookup[0];
36 for (unsigned i = 0; i < 256; ++i)
37 {
38 lup[i] = vxl_byte(s * (histo[i] - x0));
39 }
40
41 unsigned ni = image.ni(), nj = image.nj(), np = image.nplanes();
42 std::ptrdiff_t istep = image.istep(), jstep = image.jstep(), pstep = image.planestep();
43 vxl_byte * plane = image.top_left_ptr();
44 for (unsigned p = 0; p < np; ++p, plane += pstep)
45 {
46 vxl_byte * row = plane;
47 for (unsigned j = 0; j < nj; ++j, row += jstep)
48 {
49 vxl_byte * pixel = row;
50 for (unsigned i = 0; i < ni; ++i, pixel += istep)
51 *pixel = lup[*pixel];
52 }
53 }
54 }
55
56 //: Apply histogram equalisation to given float image
57 void
vil_histogram_equalise(vil_image_view<float> & image)58 vil_histogram_equalise(vil_image_view<float> & image)
59 {
60 unsigned n_bins = 4000;
61 std::vector<double> histo(n_bins);
62
63 // Find smallest and largest values in image
64 float min_v, max_v;
65 vil_math_value_range(image, min_v, max_v);
66
67 // Generate histogram
68 vil_histogram(image, histo, min_v, max_v, n_bins);
69
70 // Create cumulative frequency curve
71 double sum = 0.0;
72 for (unsigned i = 0; i < n_bins; ++i)
73 {
74 sum += histo[i];
75 histo[i] = sum;
76 }
77
78 // Get mapping parameters and generate lookup tables
79 int lo = 0;
80 while (histo[lo] == 0)
81 lo++;
82 double x0 = histo[lo];
83 // To map bins to 256 output value range (smallest values gets mapped to zero)
84 double s_b2o = 255.1 / (sum - x0);
85 std::vector<unsigned> lookup_b2o(n_bins);
86 unsigned * lup_b2o = &lookup_b2o[0];
87 for (unsigned i = 0; i < n_bins; ++i)
88 {
89 lup_b2o[i] = unsigned(s_b2o * (histo[i] - x0));
90 }
91 // To map input values to bins
92 double s_i2b = double(n_bins - 1) / (double(max_v - min_v));
93 std::vector<unsigned> lookup_i2b(max_v - min_v + 1);
94 unsigned * lup_i2b = &lookup_i2b[0];
95 for (unsigned i = 0; i <= (max_v - min_v); ++i)
96 {
97 lup_i2b[i] = unsigned(s_i2b * i);
98 }
99
100 // Update image values
101 unsigned ni = image.ni(), nj = image.nj(), np = image.nplanes();
102 std::ptrdiff_t istep = image.istep(), jstep = image.jstep(), pstep = image.planestep();
103 float * plane = image.top_left_ptr();
104 for (unsigned p = 0; p < np; ++p, plane += pstep)
105 {
106 float * row = plane;
107 for (unsigned j = 0; j < nj; ++j, row += jstep)
108 {
109 float * pixel = row;
110 for (unsigned i = 0; i < ni; ++i, pixel += istep)
111 *pixel = lup_b2o[lup_i2b[(unsigned)(*pixel - min_v)]];
112 }
113 }
114 }
115