diff --git a/doc/image_processing/contrast_enhancement/histogram_matching.rst b/doc/image_processing/contrast_enhancement/histogram_matching.rst new file mode 100644 index 0000000000..be7883a0c9 --- /dev/null +++ b/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> mask(inp_img.height(), vector(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. + diff --git a/doc/image_processing/contrast_enhancement/matching.jpg b/doc/image_processing/contrast_enhancement/matching.jpg new file mode 100644 index 0000000000..d03e07a040 Binary files /dev/null and b/doc/image_processing/contrast_enhancement/matching.jpg differ diff --git a/doc/image_processing/contrast_enhancement/matching_out.jpg b/doc/image_processing/contrast_enhancement/matching_out.jpg new file mode 100644 index 0000000000..869043accf Binary files /dev/null and b/doc/image_processing/contrast_enhancement/matching_out.jpg differ diff --git a/example/histogram_matching.cpp b/example/histogram_matching.cpp new file mode 100644 index 0000000000..494dcfd81b --- /dev/null +++ b/example/histogram_matching.cpp @@ -0,0 +1,45 @@ +// +// Copyright 2020 Debabrata Mandal +// +// 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 +#include +#include + +#include + +using namespace boost::gil; + +std::vector> get_mask(gray8_view_t const& mask) +{ + std::vector> mask_vec(mask.height(), std::vector(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; +} diff --git a/include/boost/gil/image_processing/histogram_matching.hpp b/include/boost/gil/image_processing/histogram_matching.hpp new file mode 100644 index 0000000000..ea49a5c153 --- /dev/null +++ b/include/boost/gil/image_processing/histogram_matching.hpp @@ -0,0 +1,206 @@ +// +// Copyright 2020 Debabrata Mandal +// +// 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 +#include +#include + +#include +#include +#include +#include + +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 +std::map + histogram_matching(histogram const& src_hist, histogram const& ref_hist) +{ + histogram 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 +std::map histogram_matching( + histogram const& src_hist, + histogram const& ref_hist, + histogram& dst_hist) +{ + static_assert( + std::is_integral::value && + std::is_integral::value && + std::is_integral::value, + "Source, Refernce or Destination histogram type is not appropriate."); + + using value_t = typename histogram::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 inverse_mapping; + + std::vector::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(ref_max, std::get<0>(ref_keys[start + 1]))) - + src_val)) + { + inverse_mapping[std::get<0>(src_keys[j])] = + std::min(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 +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> src_mask = {}, + std::vector> ref_mask = {}) +{ + gil_function_requires>(); + gil_function_requires>(); + gil_function_requires>(); + + static_assert( + color_spaces_are_compatible< + typename color_space_type::type, + typename color_space_type::type>::value, + "Source and reference view must have same color space"); + + static_assert( + color_spaces_are_compatible< + typename color_space_type::type, + typename color_space_type::type>::value, + "Source and destination view must have same color space"); + + // Defining channel type + using source_channel_t = typename channel_type::type; + using ref_channel_t = typename channel_type::type; + using dst_channel_t = typename channel_type::type; + using coord_t = typename SrcView::x_coord_t; + + std::size_t const channels = num_channels::value; + coord_t const width = src_view.width(); + coord_t const height = src_view.height(); + source_channel_t src_pixel_min = std::numeric_limits::min(); + source_channel_t src_pixel_max = std::numeric_limits::max(); + ref_channel_t ref_pixel_min = std::numeric_limits::min(); + ref_channel_t ref_pixel_max = std::numeric_limits::max(); + dst_channel_t dst_pixel_min = std::numeric_limits::min(); + dst_channel_t dst_pixel_max = std::numeric_limits::max(); + + for (std::size_t i = 0; i < channels; i++) + { + histogram src_histogram; + histogram ref_histogram; + fill_histogram( + nth_channel_view(src_view, i), src_histogram, bin_width, false, false, mask, src_mask, + std::tuple(src_pixel_min), + std::tuple(src_pixel_max), true); + fill_histogram( + nth_channel_view(ref_view, i), ref_histogram, bin_width, false, false, mask, ref_mask, + std::tuple(ref_pixel_min), std::tuple(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(inverse_mapping[src_it[src_x][0]]); + } + } + } +} + +}} //namespace boost::gil + +#endif diff --git a/test/core/image_processing/histogram_matching.cpp b/test/core/image_processing/histogram_matching.cpp new file mode 100644 index 0000000000..56672e0494 --- /dev/null +++ b/test/core/image_processing/histogram_matching.cpp @@ -0,0 +1,137 @@ +// +// Copyright 2020 Debabrata Mandal +// +// 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 +#include +#include +#include +#include +#include + +#include +#include + +const int a = 5; +const double epsilon = 0.000001; // Decided by the value 5/255 i.e. an error of 5 px in 255 px +boost::gil::gray8_image_t original(a, a), reference(a, a); +boost::gil::gray8_image_t processed(a, a), processed2(a, a); +std::vector > 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 > test1_reference{ + { 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 > 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 > test2_reference{ + { 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 > test3_equal_image{ + { 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 > test3_reference{ + { 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}}; + +void vector_to_gray_image(boost::gil::gray8_image_t& img, + std::vector >& grid) +{ + for(std::ptrdiff_t y=0; y +bool equal_histograms(SrcView const& v1, SrcView const& v2, double threshold = epsilon) +{ + double sum=0.0; + using channel_t = typename boost::gil::channel_type::type; + boost::gil::histogram h1, h2; + using value_t = typename boost::gil::histogram::value_type; + channel_t max_p = std::numeric_limits::max(); + channel_t min_p = std::numeric_limits::min(); + long int num_pixels = v1.width() * v1.height(); + + boost::gil::fill_histogram(v1, h1, 1, false, false); + boost::gil::fill_histogram(v2, h2, 1, false, false); + auto ch1 = boost::gil::cumulative_histogram(h1); + auto ch2 = boost::gil::cumulative_histogram(h2); + std::for_each(ch1.begin(), ch1.end(), [&](value_t const& v) { + sum+=abs(v.second-ch1[v.first]); + }); + return ( abs(sum) / (ch1.size() * (max_p - min_p)) < threshold ); +} + +void test_random_image() +{ + vector_to_gray_image(original,test1_random); + vector_to_gray_image(reference,test1_reference); + histogram_matching(boost::gil::const_view(original),boost::gil::const_view(reference),boost::gil::view(processed)); + BOOST_TEST(equal_histograms(boost::gil::view(processed), boost::gil::view(reference))); + + histogram_matching(boost::gil::const_view(processed),boost::gil::const_view(reference),boost::gil::view(processed2)); + BOOST_TEST(equal_histograms(boost::gil::view(processed), boost::gil::view(processed2))); +} + +void test_uniform_image() +{ + vector_to_gray_image(original,test2_uniform); + vector_to_gray_image(reference,test2_reference); + histogram_matching(boost::gil::const_view(original),boost::gil::const_view(reference),boost::gil::view(processed)); + BOOST_TEST(equal_histograms(boost::gil::view(processed), boost::gil::view(reference))); + + histogram_matching(boost::gil::const_view(processed),boost::gil::const_view(reference),boost::gil::view(processed2)); + BOOST_TEST(equal_histograms(boost::gil::view(processed), boost::gil::view(processed2))); +} + +void test_equal_image() +{ + vector_to_gray_image(original,test3_equal_image); + vector_to_gray_image(reference,test3_reference); + histogram_matching(boost::gil::const_view(original),boost::gil::const_view(reference),boost::gil::view(processed)); + BOOST_TEST(equal_histograms(boost::gil::view(processed), boost::gil::view(reference))); + + histogram_matching(boost::gil::const_view(processed),boost::gil::const_view(reference),boost::gil::view(processed2)); + BOOST_TEST(equal_histograms(boost::gil::view(processed), boost::gil::view(processed2))); +} + +int main() +{ + //Basic tests for grayscale histogram_equalization + test_random_image(); + test_uniform_image(); + test_equal_image(); + + return boost::report_errors(); +}