Add histogram equalization feature

doc/image_processing/contrast_enhancement/histogram_equalization.rst

@@ -0,0 +1,100 @@
+.. _he:
+
+######################
+Histogram Equalization 
+######################
+
+Description
+-----------
+
+Histogram equalization also known as histogram flattening, is a non-linear image enhancement
+algorithm that follows the idea that not only should an image cover the entire grayscale space
+but also be uniformly distributed over that range. 
+
+An ideal image would be the one having a flat histogram.
+
+Although care should be taken before applying a non-linear transformation on the image
+histogram, there are good mathematical reasons why a flat histogram is the desired goal.
+
+A simple scenario would be an image with pixels concentrated in an interval, in which case
+histogram equalization transforms pixels to achieve a flat histogram image. Thus enhancing
+the image contrast.
+
+.. figure:: he_chart.png
+    :width: 200px
+    :align: center
+    :height: 100px
+    :alt: Could not load image.
+    :figclass: align-center
+
+    Pixels concentrated in an interval spread out.
+
+Algorithm
+---------
+
+#. First calculate the histogram corresponding to input image.
+#. If it is a multi channeled image (e.g. RGB), convert it to a independent color space
+   (like YCbCr, HSV etc.).
+#. Then calculate the cumulative histogram over the input image.
+#. Normalize the histogram to bring bin values between 0-1. For multi-channeled images 
+   normalize each channel independently (by the number of pixels in image).
+#. If the histogram of image is H(p\ :sub:`x`\)  p\ :sub:`x`\  in  [0, 255], then apply 
+   the transformation p\ :sub:`x'`\  = H(p\ :sub:`x`\),  p\ :sub:`x'`\  is pixel in output
+   image.
+
+**Explanation**
+
+Since we will be transforming the image to match a flat histogram, we match
+the cumulative histogram of the image to the cumulative histogram of a flat histogram.
+
+Cumulative histogram of flat image is  H(p\ :sub:`x'`\)  =  p\ :sub:`x'` .
+
+Hence,
+
+        =>  H(p\ :sub:`x'`\)  =  H(p\ :sub:`x`\)
+
+        =>  p\ :sub:`x'`\     =  H(p\ :sub:`x`\)
+
+Results
+-------
+The algorithm is applied on a few standard images. One of the transformations in shown below:
+
+**Grayscale Image**
+
+.. figure:: barbara.jpg
+    :width: 512px
+    :align: center
+    :height: 256px
+    :alt: Could not load image.
+    :figclass: align-center
+
+**RGB**
+
+.. figure:: church.jpg
+    :width: 900px
+    :align: center
+    :height: 300px
+    :alt: Could not load image.
+    :figclass: align-center
+
+
+Demo
+----
+
+Usage Syntax:
+
+    .. code-block:: cpp
+
+            gray8_image_t inp_img;
+            read_image("your_image.png", inp_img, png_tag{});
+            gray8_image_t dst_img(inp_img.dimensions());
+            histogram_equalization(view(inp_img), view(dst_img));
+
+            // To specify mask over input image
+
+            vector<vector<bool>> mask(inp_img.height(), vector<bool>(inp_img.width(), true));
+            histogram_equalization(view(inp_img), view(dst_img), true, mask);
+
+    .. tip::  Convert an RGB image to a channel independent color space
+              before trying the histogram equalization algorithm.
+

example/histogram_equalization.cpp

@@ -0,0 +1,20 @@
+#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_adaptive.png", img, png_tag{});
+    gray8_image_t img_out(img.dimensions());
+
+    // Consider changing image to independent color space, e.g. cmyk
+    boost::gil::histogram_equalization(view(img),view(img_out));
+
+    write_view("histogram_gray_equalized.png", view(img_out), png_tag{});
+
+    return 0;
+}

include/boost/gil/image_processing/histogram_equalization.hpp

@@ -0,0 +1,150 @@
+//
+// Copyright 2020 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/histogram.hpp>
+#include <boost/gil/image.hpp>
+
+#include <cmath>
+#include <map>
+#include <vector>
+
+namespace boost { namespace gil {
+
+
+/////////////////////////////////////////
+/// Histogram Equalization(HE)
+/////////////////////////////////////////
+/// \defgroup HE HE
+/// \brief Contains implementation and description of the algorithm used to compute
+///        global histogram equalization of input images.
+///
+///        Algorithm :-
+///        1. If histogram A is to be equalized compute the cumulative histogram of A.
+///        2. Let CFD(A) refer to the cumulative histogram of A
+///        3. For a uniform histogram A', CDF(A') = A'
+///        4. We need to transfrom A to A' such that
+///        5. CDF(A') = CDF(A) => A' = CDF(A)
+///        6. Hence the pixel transform , px => histogram_of_ith_channel[px].
+///
+
+/// \fn histogram_equalization
+/// \ingroup HE
+/// \tparam SrcKeyType Key Type of input histogram
+/// @param src_hist INPUT Input source histogram
+/// \brief Overload for histogram equalization algorithm, takes in a single source histogram 
+///        and returns the color map used for histogram equalization.
+///
+template <typename SrcKeyType>
+std::map<SrcKeyType, SrcKeyType> histogram_equalization(histogram<SrcKeyType> const& src_hist)
+{
+    histogram<SrcKeyType> dst_hist;
+    return histogram_equalization(src_hist, dst_hist);
+}
+
+/// \overload histogram_equalization
+/// \ingroup HE
+/// \tparam SrcKeyType Key Type of input histogram
+/// \tparam DstKeyType Key Type of output histogram
+/// @param src_hist INPUT source histogram
+/// @param dst_hist OUTPUT Output histogram
+/// \brief Overload for histogram equalization algorithm, takes in both source histogram &
+///        destination histogram and returns the color map used for histogram equalization
+///        as well as transforming the destination histogram.
+///
+template <typename SrcKeyType, typename DstKeyType>
+std::map<SrcKeyType, DstKeyType>
+    histogram_equalization(histogram<SrcKeyType> const& src_hist, histogram<DstKeyType>& dst_hist)
+{
+    static_assert(
+        std::is_integral<SrcKeyType>::value &&
+        std::is_integral<DstKeyType>::value,
+        "Source and destination histogram types are not appropriate");
+
+    using value_t = typename histogram<SrcKeyType>::value_type;
+    dst_hist.clear();
+    double sum          = src_hist.sum();
+    SrcKeyType min_key  = std::numeric_limits<DstKeyType>::min();
+    SrcKeyType max_key  = std::numeric_limits<DstKeyType>::max();
+    auto cumltv_srchist = cumulative_histogram(src_hist);
+    std::map<SrcKeyType, DstKeyType> color_map;
+    std::for_each(cumltv_srchist.begin(), cumltv_srchist.end(), [&](value_t const& v) {
+        DstKeyType trnsfrmd_key =
+            static_cast<DstKeyType>((v.second * (max_key - min_key)) / sum + min_key);
+        color_map[std::get<0>(v.first)] = trnsfrmd_key;
+    });
+    std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
+        dst_hist[color_map[std::get<0>(v.first)]] += v.second;
+    });
+    return color_map;
+}
+
+/// \overload histogram_equalization
+/// \ingroup HE
+/// @param src_view  INPUT source image view
+/// @param dst_view  OUTPUT Output image view
+/// @param bin_width INPUT Histogram bin width
+/// @param mask      INPUT Specify is mask is to be used
+/// @param src_mask  INPUT Mask vector over input image
+/// \brief Overload for histogram equalization algorithm, takes in both source & destination
+///        image views and histogram equalizes the input image.
+///
+template <typename SrcView, typename DstView>
+void histogram_equalization(
+    SrcView const& src_view,
+    DstView const& dst_view,
+    std::size_t bin_width = 1,
+    bool mask = false,
+    std::vector<std::vector<bool>> src_mask = {})
+{
+    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");
+
+    // Defining channel type
+    using source_channel_t = typename channel_type<SrcView>::type;
+    using dst_channel_t    = typename channel_type<DstView>::type;
+    using coord_t          = typename SrcView::x_coord_t;
+
+    std::size_t const channels = num_channels<SrcView>::value;
+    coord_t const width        = src_view.width();
+    coord_t const height       = src_view.height();
+    std::size_t pixel_max      = std::numeric_limits<dst_channel_t>::max();
+    std::size_t pixel_min      = std::numeric_limits<dst_channel_t>::min();
+
+    for (std::size_t i = 0; i < channels; i++)
+    {
+        histogram<source_channel_t> h;
+        fill_histogram(nth_channel_view(src_view, i), h, bin_width, false, false, mask, src_mask);
+        h.normalize();
+        auto h2 = cumulative_histogram(h);
+        for (std::ptrdiff_t src_y = 0; src_y < height; ++src_y)
+        {
+            auto src_it = nth_channel_view(src_view, i).row_begin(src_y);
+            auto dst_it = nth_channel_view(dst_view, i).row_begin(src_y);
+            for (std::ptrdiff_t src_x = 0; src_x < width; ++src_x)
+            {
+                if (mask && !src_mask[src_y][src_x])
+                    dst_it[src_x][0] = channel_convert<dst_channel_t>(src_it[src_x][0]);
+                else
+                    dst_it[src_x][0] = static_cast<dst_channel_t>(
+                        h2[src_it[src_x][0]] * (pixel_max - pixel_min) + pixel_min);
+            }
+        }
+    }
+}
+
+}}  //namespace boost::gil
+
+#endif