Add histogram matching algorithm

doc/image_processing/contrast_enhancement/histogram_matching.rst

@@ -0,0 +1,87 @@
+.. _hm:
+
+##################
+Histogram Matching 
+##################
+
+Description
+-----------
+
+Histogram Matching is a technique to match the histograms of two images.
+
+One use case of this would be when two images of the same location have been taken 
+under the same local illumination but with different sensors, bringing out different 
+features in either image.
+
+The famous histogram equalization is a special case of this algorithm when the reference image
+is expected to have a uniform histogram.
+
+
+Algorithm
+---------
+
+#. Calculate the histogram corresponding to input image and reference image.
+#. If it is a multi channeled image (e.g. RGB), convert both to an independent color space
+ (like YCbCr, HSV etc.).
+#. Then calculate the cumulative histogram over the input image and reference image.
+#. Normalize both 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 cumulative histogram of input image is H(p\ :sub:`x`\) and of reference image is R(p\ :sub:`x'`\) 
+ p\ :sub:`x`\ & p\ :sub:`x'`\ in [0, 255], then apply the transformation 
+ p\ :sub:`x'`\ = R\ :sup:`-1`\ (H(p\ :sub:`x`\ ))
+
+**Explanation**
+
+Since we will be transforming the image to match a reference image, we match
+the cumulative histogram of the image to the cumulative histogram of the reference histogram.
+
+Hence,
+
+ => R(p\ :sub:`x'`\) = H(p\ :sub:`x`\ )
+
+ => p\ :sub:`x'`\ = R\ :sup:`-1`\ (H(p\ :sub:`x`\ ))
+
+Results
+-------
+The algorithm is applied on a few standard images. One of the transformations in shown below:
+
+**Original Image(left) & Reference Image(right)**
+
+.. figure:: matching.jpg
+ :width: 600px
+ :align: center
+ :height: 300px
+ :alt: Could not load image.
+ :figclass: align-center
+
+**Histogram matched Image**
+
+.. figure:: matching_out.jpg
+ :width: 300px
+ :align: center
+ :height: 300px
+ :alt: Could not load image.
+ :figclass: align-center
+
+
+Demo
+----
+
+Usage Syntax:
+
+ .. code-block:: cpp
+
+ gray8_image_t inp_img, ref_img;
+ read_image("your_image.png", inp_img, png_tag{});
+ read_image("your_ref_image.png", ref_img, png_tag{});
+ gray8_image_t dst_img(inp_img.dimensions());
+ histogram_matching(view(inp_img), view(ref_image), view(dst_img));
+
+ // To specify mask over input image
+
+ vector<vector<bool>> mask(inp_img.height(), vector<bool>(inp_img.width(), true));
+ histogram_matching(view(inp_img), view(ref_image), view(dst_img), true, mask);
+
+ .. tip:: Convert an RGB image to a channel independent color space
+ before trying the histogram matching algorithm.
+

example/histogram_matching.cpp

@@ -0,0 +1,45 @@
+//
+// 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)
+//
+
+#include <boost/gil.hpp>
+#include <boost/gil/extension/io/png.hpp>
+#include <boost/gil/image_processing/histogram_matching.hpp>
+
+#include <vector>
+
+using namespace boost::gil;
+
+std::vector<std::vector<bool>> get_mask(gray8_view_t const& mask)
+{
+ std::vector<std::vector<bool>> mask_vec(mask.height(), std::vector<bool>(mask.width(), 0));
+ for (std::size_t i = 0; i < mask.height(); i++)
+ {
+ for (std::size_t j = 0; j < mask.width(); j++)
+ {
+ mask_vec[i][j] = mask(j, i);
+ }
+ }
+ return mask_vec;
+}
+
+int main()
+{
+ gray8_image_t img;
+ read_image("test_adaptive.png", img, png_tag{});
+
+ gray8_image_t ref_img;
+ read_image("test_adaptive.png", ref_img, png_tag{});
+
+ gray8_image_t img_out(img.dimensions());
+
+ boost::gil::histogram_matching(view(img), view(ref_img), view(img_out));
+
+ write_view("histogram_gray_matching.png", view(img_out), png_tag{});
+
+ return 0;
+}

include/boost/gil/image_processing/histogram_matching.hpp

@@ -0,0 +1,206 @@
+//
+// 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_MATCHING_HPP
+#define BOOST_GIL_IMAGE_PROCESSING_HISTOGRAM_MATCHING_HPP
+
+#include <boost/gil/algorithm.hpp>
+#include <boost/gil/histogram.hpp>
+#include <boost/gil/image.hpp>
+
+#include <algorithm>
+#include <cmath>
+#include <map>
+#include <vector>
+
+namespace boost { namespace gil {
+
+/////////////////////////////////////////
+/// Histogram Matching(HM)
+/////////////////////////////////////////
+/// \defgroup HM HM
+/// \brief Contains implementation and description of the algorithm used to compute
+/// global histogram matching of input images.
+///
+/// Algorithm :-
+/// 1. Calculate histogram A(pixel) of input image and G(pixel) of reference image.
+/// 2. Compute the normalized cumulative(CDF) histograms of A and G.
+/// 3. Match the histograms using transofrmation => CDF(A(px)) = CDF(G(px'))
+/// => px' = Inv-CDF (CDF(px))
+///
+
+/// \fn histogram_matching
+/// \ingroup HM
+/// \tparam SrcKeyType Key Type of input histogram
+/// @param src_hist INPUT Input source histogram
+/// @param ref_hist INPUT Input reference histogram
+/// \brief Overload for histogram matching algorithm, takes in a single source histogram &
+/// reference histogram and returns the color map used for histogram matching.
+///
+template <typename SrcKeyType, typename RefKeyType>
+std::map<SrcKeyType, SrcKeyType>
+ histogram_matching(histogram<SrcKeyType> const& src_hist, histogram<RefKeyType> const& ref_hist)
+{
+ histogram<SrcKeyType> dst_hist;
+ return histogram_matching(src_hist, ref_hist, dst_hist);
+}
+
+/// \overload histogram_matching
+/// \ingroup HM
+/// \tparam SrcKeyType Key Type of input histogram
+/// \tparam RefKeyType Key Type of reference histogram
+/// \tparam DstKeyType Key Type of output histogram
+/// @param src_hist INPUT source histogram
+/// @param ref_hist INPUT reference histogram
+/// @param dst_hist OUTPUT Output histogram
+/// \brief Overload for histogram matching algorithm, takes in source histogram, reference 
+/// histogram & destination histogram and returns the color map used for histogram
+/// matching as well as transforming the destination histogram.
+///
+template <typename SrcKeyType, typename RefKeyType, typename DstKeyType>
+std::map<SrcKeyType, DstKeyType> histogram_matching(
+ histogram<SrcKeyType> const& src_hist,
+ histogram<RefKeyType> const& ref_hist,
+ histogram<DstKeyType>& dst_hist)
+{
+ static_assert(
+ std::is_integral<SrcKeyType>::value &&
+ std::is_integral<RefKeyType>::value &&
+ std::is_integral<DstKeyType>::value,
+ "Source, Refernce or Destination histogram type is not appropriate.");
+
+ using value_t = typename histogram<SrcKeyType>::value_type;
+ dst_hist.clear();
+ double src_sum = src_hist.sum();
+ double ref_sum = ref_hist.sum();
+ auto cumltv_srchist = cumulative_histogram(src_hist);
+ auto cumltv_refhist = cumulative_histogram(ref_hist);
+ std::map<SrcKeyType, RefKeyType> inverse_mapping;
+
+ std::vector<typename histogram<RefKeyType>::key_type> src_keys, ref_keys;
+ src_keys = src_hist.sorted_keys();
+ ref_keys = ref_hist.sorted_keys();
+ std::ptrdiff_t start = ref_keys.size() - 1;
+ RefKeyType ref_max;
+ if (start >= 0)
+ ref_max = std::get<0>(ref_keys[start]);
+
+ for (std::ptrdiff_t j = src_keys.size() - 1; j >= 0; --j)
+ {
+ double src_val = (cumltv_srchist[src_keys[j]] * ref_sum) / src_sum;
+ while (cumltv_refhist[ref_keys[start]] > src_val && start > 0)
+ {
+ start--;
+ }
+ if (abs(cumltv_refhist[ref_keys[start]] - src_val) >
+ abs(cumltv_refhist(std::min<RefKeyType>(ref_max, std::get<0>(ref_keys[start + 1]))) -
+ src_val))
+ {
+ inverse_mapping[std::get<0>(src_keys[j])] = 
+ std::min<RefKeyType>(ref_max, std::get<0>(ref_keys[start + 1]));
+ }
+ else
+ {
+ inverse_mapping[std::get<0>(src_keys[j])] = std::get<0>(ref_keys[start]);
+ }
+ if (j == 0)
+ break;
+ }
+ std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
+ dst_hist[inverse_mapping[std::get<0>(v.first)]] += v.second;
+ });
+ return inverse_mapping;
+}
+
+/// \overload histogram_matching
+/// \ingroup HM
+/// @param src_view INPUT source image view
+/// @param ref_view INPUT Reference 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
+/// @param ref_mask INPUT Mask vector over reference image
+/// \brief Overload for histogram matching algorithm, takes in both source, reference & 
+/// destination image views and histogram matches the input image using the 
+/// reference image.
+///
+template <typename SrcView, typename ReferenceView, typename DstView>
+void histogram_matching(
+ SrcView const& src_view,
+ ReferenceView const& ref_view,
+ DstView const& dst_view,
+ std::size_t bin_width = 1,
+ bool mask = false,
+ std::vector<std::vector<bool>> src_mask = {},
+ std::vector<std::vector<bool>> ref_mask = {})
+{
+ gil_function_requires<ImageViewConcept<SrcView>>();
+ gil_function_requires<ImageViewConcept<ReferenceView>>();
+ gil_function_requires<MutableImageViewConcept<DstView>>();
+
+ static_assert(
+ color_spaces_are_compatible<
+ typename color_space_type<SrcView>::type,
+ typename color_space_type<ReferenceView>::type>::value,
+ "Source and reference view must have same color space");
+
+ static_assert(
+ color_spaces_are_compatible<
+ typename color_space_type<SrcView>::type,
+ typename color_space_type<DstView>::type>::value,
+ "Source and destination view must have same color space");
+
+ // Defining channel type
+ using source_channel_t = typename channel_type<SrcView>::type;
+ using ref_channel_t = typename channel_type<ReferenceView>::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();
+ source_channel_t src_pixel_min = std::numeric_limits<source_channel_t>::min();
+ source_channel_t src_pixel_max = std::numeric_limits<source_channel_t>::max();
+ ref_channel_t ref_pixel_min = std::numeric_limits<ref_channel_t>::min();
+ ref_channel_t ref_pixel_max = std::numeric_limits<ref_channel_t>::max();
+ dst_channel_t dst_pixel_min = std::numeric_limits<dst_channel_t>::min();
+ dst_channel_t dst_pixel_max = std::numeric_limits<dst_channel_t>::max();
+
+ for (std::size_t i = 0; i < channels; i++)
+ {
+ histogram<source_channel_t> src_histogram;
+ histogram<ref_channel_t> ref_histogram;
+ fill_histogram(
+ nth_channel_view(src_view, i), src_histogram, bin_width, false, false, mask, src_mask,
+ std::tuple<source_channel_t>(src_pixel_min),
+ std::tuple<source_channel_t>(src_pixel_max), true);
+ fill_histogram(
+ nth_channel_view(ref_view, i), ref_histogram, bin_width, false, false, mask, ref_mask,
+ std::tuple<ref_channel_t>(ref_pixel_min), std::tuple<ref_channel_t>(ref_pixel_max),
+ true);
+ auto inverse_mapping = histogram_matching(src_histogram, ref_histogram);
+ 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] = src_it[src_x][0];
+ else
+ dst_it[src_x][0] =
+ static_cast<dst_channel_t>(inverse_mapping[src_it[src_x][0]]);
+ }
+ }
+ }
+}
+
+}} //namespace boost::gil
+
+#endif