Added histogram equalization for images

boostorg · Feb 19, 2020 · f9eeea4 · f9eeea4
1 parent 2ff00f6
commit f9eeea4
Show file tree

Hide file tree

Showing 3 changed files with 267 additions and 0 deletions.
diff --git a/example/histogram_equalization.cpp b/example/histogram_equalization.cpp
@@ -0,0 +1,19 @@
+#include <boost/gil.hpp>
+#include <boost/gil/extension/io/png.hpp>
+#include <boost/gil/image_view_factory.hpp>
+#include <boost/gil/image_processing/histogram_equalization.hpp>
+
+using namespace boost::gil;
+
+int main()
+{
+    gray8_image_t img;
+    read_image("test_he.png", img, png_tag{});
+    gray8_image_t img_out(img.dimensions());
+
+    boost::gil::histogram_equalization(view(img),view(img_out));
+
+    write_view("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,162 @@
+//
+// 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_HE_HPP
+#define BOOST_GIL_IMAGE_PROCESSING_HE_HPP
+
+#include <array>
+
+#include <boost/gil.hpp>
+#include <boost/assert.hpp>
+
+namespace boost { namespace gil {
+
+template <typename SrcView,std::size_t channels>
+void compute_histogram(
+    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 compute_histogram_cdf(
+    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];
+        }
+    }
+}
+
+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 result_channel_t = typename channel_type<DstView>::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
+    compute_histogram<SrcView,channels>(src_view, histogram, min, max);
+
+    compute_histogram_cdf<channels>(histogram);
+
+    source_channel_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::ptrdiff_t c = 0; c < channels; ++c)
+            {
+                //Automatic rounding off signiicant digits after decimal by clipping
+                dst_it[src_x][c] = (histogram[c][src_it[src_x][c]] * (max[c] - min[c]) ) / num_pixels ;
+            }
+        }
+    }
+}
+
+}}
+
+#endif // !BOOST_GIL_IMAGE_PROCESSING_HE_HPP
diff --git a/test/core/image_processing/histogram_equalization.cpp b/test/core/image_processing/histogram_equalization.cpp
@@ -0,0 +1,86 @@
+#include <vector>
+
+#include <boost/gil/image_processing/histogram_equalization.hpp>
+#include <boost/core/lightweight_test.hpp>
+
+using namespace boost::gil;
+
+int main()
+{
+    //Basic tests for histogram_equalization
+    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}};
+    for(std::ptrdiff_t y=0; y<expected_test1.size(); ++y)
+    {
+        for(std::ptrdiff_t x=0; x<expected_test1[0].size(); ++x)
+        {
+            boost::gil::view(original)(x,y) = boost::gil::gray8_pixel_t(test1_random[y][x]);
+            boost::gil::view(expected)(x,y) = boost::gil::gray8_pixel_t(expected_test1[y][x]);
+        }
+    }
+    //start 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)));
+    //end test1
+    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}};    
+    for(std::ptrdiff_t y=0; y<expected_test2.size(); ++y)
+    {
+        for(std::ptrdiff_t x=0; x<expected_test2[0].size(); ++x)
+        {
+            boost::gil::view(original)(x,y) = boost::gil::gray8_pixel_t(test2_uniform[y][x]);      
+            boost::gil::view(expected)(x,y) = boost::gil::gray8_pixel_t(expected_test2[y][x]);
+        }
+    }
+    //start 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)));
+    //end test2
+    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}};    
+    for(std::ptrdiff_t y=0; y<expected_test2.size(); ++y)
+    {
+        for(std::ptrdiff_t x=0; x<expected_test2[0].size(); ++x)
+        {
+            boost::gil::view(original)(x,y) = boost::gil::gray8_pixel_t(test3_2peaks[y][x]);      
+            boost::gil::view(expected)(x,y) = boost::gil::gray8_pixel_t(expected_test3[y][x]);
+        }
+    }
+    //start 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)));
+    //end test3
+    return boost::report_errors();
+}