Re-adding histogram equalization

boostorg · Mar 1, 2020 · c1e15a3 · c1e15a3
1 parent f3d3862
commit c1e15a3
Show file tree

Hide file tree

Showing 3 changed files with 281 additions and 0 deletions.
diff --git a/example/histogram_equalization.cpp b/example/histogram_equalization.cpp
@@ -0,0 +1,18 @@
+#include <boost/gil.hpp>
+#include <boost/gil/extension/io/png.hpp>
+#include <boost/gil/image_processing/histogram_equalization.hpp>
+
+using namespace boost::gil;
+
+int main()
+{
+    gray8_image_t img;
+    read_image("test-hist-equaliz.png", img, png_tag{});
+    gray8_image_t img_out(img.dimensions());
+
+    boost::gil::histogram_equalization(view(img),view(img_out));
+
+    write_view("out-histogram-equalized-image.png", view(img_out), png_tag{});
+
+    return 0;
+}
diff --git a/include/boost/gil/image_processing/histogram_equalization.hpp b/include/boost/gil/image_processing/histogram_equalization.hpp
@@ -0,0 +1,160 @@
+//
+// Copyright 2019 Debabrata Mandal <[email protected]>
+//
+// Use, modification and distribution are subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+#ifndef BOOST_GIL_IMAGE_PROCESSING_HISTOGRAM_EQUALIZATION_HPP
+#define BOOST_GIL_IMAGE_PROCESSING_HISTOGRAM_EQUALIZATION_HPP
+
+#include <boost/gil/image.hpp>
+
+#include <array>
+
+namespace boost { namespace gil {
+
+namespace detail {
+
+template <typename SrcView,std::size_t channels>
+void histogram_impl(SrcView const& src_view,
+    std::array<std::array<std::size_t, 256>, channels>& histogram,
+    std::array<typename channel_type<SrcView>::type, channels>& min,
+    std::array<typename channel_type<SrcView>::type, channels>& max)
+{
+    using source_channel_t = typename channel_type<SrcView>::type;
+
+    //Checking channel sizes i.e. for channel sizes more than 1 byte or signed values,
+    //need to find max and min pixel intensities and scale the values to [0,255] appropriately
+    if(sizeof(source_channel_t) > 1 || std::is_signed<source_channel_t>::value)
+    {
+        //Per channel min and max values found independently
+        for(std::ptrdiff_t src_y = 0; src_y < src_view.height(); ++src_y)
+        {
+            typename SrcView::x_iterator src_it = src_view.row_begin(src_y);
+
+            for(std::ptrdiff_t src_x = 0; src_x < src_view.width(); ++src_x)
+            {
+                std::ptrdiff_t c = 0;
+                static_for_each(src_it[src_x],[&c,&min](source_channel_t px) {
+                    min[c] = px < min[c] ? px : min[c];
+                    c++;
+                });
+                c = 0;
+                static_for_each(src_it[src_x],[&c,&max](source_channel_t px) {
+                    max[c] = px > max[c] ? px : max[c];
+                    c++;
+                });
+            }
+        }
+
+        //Histogram calculation after scaling intensities to lie in [0,255]
+        for(std::ptrdiff_t src_y = 0; src_y < src_view.height(); ++src_y)
+        {
+            typename SrcView::x_iterator src_it = src_view.row_begin(src_y);
+
+            for(std::ptrdiff_t src_x = 0; src_x < src_view.width(); ++src_x)
+            {
+                std::ptrdiff_t c = 0;
+                static_for_each(src_it[src_x],[&c,&histogram,&min,&max](source_channel_t px) {
+                    histogram[c][((px - min[c]) * 255) / (max[c] - min[c])]++;
+                    c++;
+                });
+            }
+        }
+    }
+    else
+    {
+        //Default intitialzation when unsigned 8 bit pixel depths
+        min.fill(0);
+        max.fill(255);
+        //Histogram calculation for default case
+        for(std::ptrdiff_t src_y = 0; src_y < src_view.height(); ++src_y)
+        {
+            typename SrcView::x_iterator src_it = src_view.row_begin(src_y);
+
+            for(std::ptrdiff_t src_x = 0; src_x < src_view.width(); ++src_x)
+            {
+                std::ptrdiff_t c = 0;
+                static_for_each(src_it[src_x],[&c,&histogram](source_channel_t px) {
+                    histogram[c][px]++;
+                    c++;
+                });
+            }
+        }
+    }
+}
+
+template <std::size_t channels>
+void histogram_cdf_impl(
+    std::array<std::array<std::size_t, 256>, channels>& hist_cdf)
+{
+    //Compute per channel Cumulative Density Function 
+    for(std::size_t c = 0; c < channels; ++c)
+    {
+        for(std::size_t i = 1; i < 256; ++i)
+        {
+            hist_cdf[c][i] += hist_cdf[c][i-1];
+        }
+    }
+}
+
+}//namespace boost::gil::detail
+
+template <typename SrcView, typename DstView>
+void histogram_equalization(SrcView const& src_view, DstView& dst_view)
+{
+    // gil_function_requires<ImageViewConcept<SrcView>>();
+    // gil_function_requires<MutableImageViewConcept<DstView>>();
+    static_assert(color_spaces_are_compatible
+    <
+        typename color_space_type<SrcView>::type,
+        typename color_space_type<DstView>::type
+    >::value, "Source and destination views must have same color space");
+
+    //Deciding channel types 
+    using source_channel_t = typename channel_type<SrcView>::type;
+    using x_coord_t = typename SrcView::x_coord_t;
+    using y_coord_t = typename SrcView::y_coord_t;
+
+    x_coord_t const width = src_view.width();
+    y_coord_t const height = src_view.height();
+    std::size_t const channels = num_channels<SrcView>::value;
+
+    //Initializations for min and max arrays to be filled while finding per channel
+    //extremas
+    std::array<source_channel_t, channels>  min{
+                std::numeric_limits<source_channel_t>::max()},
+                                            max{
+                std::numeric_limits<source_channel_t>::min()};
+
+    std::array<std::array<std::size_t, 256>, channels> histogram{};
+
+    //Passing channel size as template parameter
+    detail::histogram_impl<SrcView,channels>(src_view, histogram, min, max);
+
+    detail::histogram_cdf_impl<channels>(histogram);
+
+    std::size_t num_pixels = (height * width);
+
+    //Using the histogram equalization function , if pdf A is to be equalized to the
+    //uniform pdf A' for which CDF(A') = A' , then CDF(A') = CDF(A) => A' = CDF(A)
+    //hence the pixel transform , px => histogram[px](for that channel).
+    for(std::ptrdiff_t src_y = 0; src_y < height; ++src_y)
+    {
+        for(std::ptrdiff_t src_x = 0; src_x < width; ++src_x)
+        {
+            typename SrcView::x_iterator src_it = src_view.row_begin(src_y);
+            typename SrcView::x_iterator dst_it = dst_view.row_begin(src_y);
+
+            for(std::size_t c = 0; c < channels; ++c)
+            {
+                dst_it[src_x][c] = (histogram[c][src_it[src_x][c]] * (max[c] - min[c])) / num_pixels;
+            }
+        }
+    }
+}
+
+}} //namespace boost::gil
+
+#endif // !BOOST_GIL_IMAGE_PROCESSING_HISTOGRAM_EQUALIZATION_HPP
diff --git a/test/core/image_processing/histogram_equalization.cpp b/test/core/image_processing/histogram_equalization.cpp
@@ -0,0 +1,103 @@
+//
+// Copyright 2019 Debabrata Mandal <[email protected]>
+//
+// Use, modification and distribution are subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+//
+#include <boost/gil/image_processing/histogram_equalization.hpp>
+#include <boost/core/lightweight_test.hpp>
+
+#include <vector>
+
+namespace gil = boost::gil;
+
+boost::gil::gray8_image_t original(5, 5);
+boost::gil::gray8_image_t processed(5, 5), expected(5, 5);
+std::vector<std::vector<int> > test1_random{
+    {  1,   10,  10,  10,  10},
+    {  20,  25,  25,  55,  20},
+    {  0,   55,  55,  55,  20},
+    {  20,  255, 255, 255, 0},
+    {  100, 100, 100, 10,  0}};
+std::vector<std::vector<int> > expected_test1{
+    {  40,  91,  91,  91,  91},
+    {  132, 153, 153, 193, 132},
+    {  30,  193, 193, 193, 132},
+    {  132, 255, 255, 255, 30},
+    {  224, 224, 224, 91,  30}};
+
+std::vector<std::vector<int> > test2_uniform{
+    {  0,   10,  20,  30,  40},
+    {  50,  60,  70,  80,  90},
+    {  100, 110, 120, 130, 140},
+    {  150, 160, 170, 180, 190},
+    {  200, 210, 220, 230, 240}};
+std::vector<std::vector<int> > expected_test2{
+    {  10,  20,  30,  40,  51},
+    {  61,  71,  81,  91,  102},
+    {  112, 122, 132, 142, 153},
+    {  163, 173, 183, 193, 204},
+    {  214, 224, 234, 244, 255}};
+
+std::vector<std::vector<int> > test3_2peaks{
+    {  0,   0,   0,   0,   10},
+    {  40,  43,  44,  46,  50},
+    {  55,  56,  44,  46,  44},
+    {  200, 201, 202, 203, 200},
+    {  201, 202, 201, 201, 22}};
+std::vector<std::vector<int> > expected_test3{
+    {  40,  40,  40,  40,  51},
+    {  71,  81,  112, 132, 142},
+    {  153, 163, 112, 132, 112},
+    {  183, 224, 244, 255, 183},
+    {  224, 244, 224, 224, 61}};
+
+void vector_to_gray_image(boost::gil::gray8_image_t& img,
+    std::vector<std::vector<int> >& grid)
+{
+    for(std::ptrdiff_t y = 0; y < grid.size(); ++y)
+    {
+        for(std::ptrdiff_t x = 0; x < grid[0].size(); ++x)
+        {
+            boost::gil::view(img)(x,y) = boost::gil::gray8_pixel_t(grid[y][x]);      
+        }
+    }
+}
+
+void test_random_image()
+{
+    vector_to_gray_image(original,test1_random);
+    vector_to_gray_image(expected,expected_test1);
+    histogram_equalization(boost::gil::view(original),boost::gil::view(processed));
+    BOOST_TEST(boost::gil::equal_pixels(boost::gil::view(processed), boost::gil::view(expected)));
+}
+
+void test_uniform_image()
+{
+    vector_to_gray_image(original,test2_uniform);
+    vector_to_gray_image(expected,expected_test2);
+    histogram_equalization(boost::gil::view(original),boost::gil::view(processed));
+    BOOST_TEST(boost::gil::equal_pixels(boost::gil::view(processed), boost::gil::view(expected)));
+}
+
+void test_double_peaked_image()
+{
+    vector_to_gray_image(original,test3_2peaks);
+    vector_to_gray_image(expected,expected_test3);
+    histogram_equalization(boost::gil::view(original),boost::gil::view(processed));
+    BOOST_TEST(boost::gil::equal_pixels(boost::gil::view(processed), boost::gil::view(expected)));
+}
+
+int main()
+{
+    //Basic tests for grayscale histogram_equalization
+    test_random_image();
+    test_uniform_image();
+    test_double_peaked_image();
+
+    //Basic test for multichanneled image
+    // test_independent_channel_image();
+
+    return boost::report_errors();
+}