Add code for ahe algorithm (boostorg#516)

example/adaptive_he.cpp

@@ -0,0 +1,25 @@
+//
+// Copyright 2020 Debabrata Mandal <mandaldebabrata123@gmail.com>
+//
+// 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/extension/numeric/algorithm.hpp>
+#include <boost/gil/image_processing/adaptive_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());
+
+    boost::gil::non_overlapping_interpolated_clahe(view(img), view(img_out));
+    write_view("out-adaptive.png", view(img_out), png_tag{});
+    return 0;
+}

include/boost/gil/image_processing/adaptive_histogram_equalization.hpp

@@ -0,0 +1,307 @@
+//
+// Copyright 2020 Debabrata Mandal <mandaldebabrata123@gmail.com>
+//
+// 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_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP
+#define BOOST_GIL_IMAGE_PROCESSING_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP
+
+#include <boost/gil/algorithm.hpp>
+#include <boost/gil/histogram.hpp>
+#include <boost/gil/image.hpp>
+#include <boost/gil/image_processing/histogram_equalization.hpp>
+#include <boost/gil/image_view_factory.hpp>
+
+#include <cmath>
+#include <map>
+#include <vector>
+
+namespace boost { namespace gil {
+
+/////////////////////////////////////////
+/// Adaptive Histogram Equalization(AHE)
+/////////////////////////////////////////
+/// \defgroup AHE AHE
+/// \brief Contains implementation and description of the algorithm used to compute
+///        adaptive histogram equalization of input images. Naming for the AHE functions
+///        are done in the following way 
+///             <feature-1>_<feature-2>_.._<feature-n>ahe
+///        For example, for AHE done using local (non-overlapping) tiles/blocks and 
+///        final output interpolated among tiles , it is called
+///             non_overlapping_interpolated_clahe
+///
+
+namespace detail {
+
+/// \defgroup AHE-helpers AHE-helpers
+/// \brief AHE helper functions
+
+/// \fn double actual_clip_limit
+/// \ingroup AHE-helpers
+/// \brief Computes the actual clip limit given a clip limit value using binary search.
+///        Reference -  Adaptive Histogram Equalization and Its Variations
+///                     (http://www.cs.unc.edu/techreports/86-013.pdf, Pg - 15)
+///
+template <typename SrcHist>
+double actual_clip_limit(SrcHist const& src_hist, double cliplimit = 0.03)
+{
+    double epsilon       = 1.0;
+    using value_t        = typename SrcHist::value_type;
+    double sum           = src_hist.sum();
+    std::size_t num_bins = src_hist.size();
+
+    cliplimit = sum * cliplimit;
+    long low = 0, high = cliplimit, middle = low;
+    while (high - low >= 1)
+    {
+        middle      = (low + high + 1) >> 1;
+        long excess = 0;
+        std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
+            if (v.second > middle)
+                excess += v.second - middle;
+        });
+        if (abs(excess - (cliplimit - middle) * num_bins) < epsilon)
+            break;
+        else if (excess > (cliplimit - middle) * num_bins)
+            high = middle - 1;
+        else
+            low = middle + 1;
+    }
+    return middle / sum;
+}
+
+/// \fn void clip_and_redistribute
+/// \ingroup AHE-helpers
+/// \brief Clips and redistributes excess pixels based on the actual clip limit value 
+///        obtained from the other helper function actual_clip_limit
+///        Reference - Graphic Gems 4, Pg. 474
+///        (http://cas.xav.free.fr/Graphics%20Gems%204%20-%20Paul%20S.%20Heckbert.pdf)
+/// 
+template <typename SrcHist, typename DstHist>
+void clip_and_redistribute(SrcHist const& src_hist, DstHist& dst_hist, double clip_limit = 0.03)
+{
+    using value_t            = typename SrcHist::value_type;
+    double sum               = src_hist.sum();
+    double actual_clip_value = detail::actual_clip_limit(src_hist, clip_limit);
+    // double actual_clip_value = clip_limit;
+    long actual_clip_limit = actual_clip_value * sum;
+    double excess          = 0;
+    std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
+        if (v.second > actual_clip_limit)
+            excess += v.second - actual_clip_limit;
+    });
+    std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
+        if (v.second >= actual_clip_limit)
+            dst_hist[dst_hist.key_from_tuple(v.first)] = clip_limit * sum;
+        else
+            dst_hist[dst_hist.key_from_tuple(v.first)] = v.second + excess / src_hist.size();
+    });
+    long rem = long(excess) % src_hist.size();
+    if (rem == 0)
+        return;
+    long period       = round(src_hist.size() / rem);
+    std::size_t index = 0;
+    while (rem)
+    {
+        if (dst_hist(index) >= clip_limit * sum)
+        {
+            index = (index + 1) % src_hist.size();
+        }
+        dst_hist(index)++;
+        rem--;
+        index = (index + period) % src_hist.size();
+    }
+}
+
+}  // namespace detail
+
+
+/// \fn void non_overlapping_interpolated_clahe
+/// \ingroup AHE
+/// @param src_view      Input   Source image view
+/// @param dst_view      Output  Output image view
+/// @param tile_width_x  Input   Tile width along x-axis to apply HE
+/// @param tile_width_y  Input   Tile width along x-axis to apply HE
+/// @param clip_limit    Input   Clipping limit to be applied
+/// @param bin_width     Input   Bin widths for histogram
+/// @param mask          Input   Specify if mask is to be used
+/// @param src_mask      Input   Mask on input image to ignore specified pixels
+/// \brief Performs local histogram equalization on tiles of size (tile_width_x, tile_width_y)
+///        Then uses the clip limit to redistribute excess pixels above the limit uniformly to
+///        other bins. The clip limit is specified as a fraction i.e. a bin's value is clipped 
+///        if bin_value >= clip_limit * (Total number of pixels in the tile) 
+///
+template <typename SrcView, typename DstView>
+void non_overlapping_interpolated_clahe(
+    SrcView const& src_view,
+    DstView const& dst_view,
+    std::size_t tile_width_x                = 20,
+    std::size_t tile_width_y                = 20,
+    double clip_limit                       = 0.03,
+    std::size_t bin_width                   = 1.0,
+    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");
+
+    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();
+
+    // Find control points
+
+    std::vector<coord_t> sample_x;
+    coord_t sample_x1 = tile_width_x / 2;
+    coord_t sample_x2 = (tile_width_x + 1) / 2;
+    coord_t sample_y1 = tile_width_y / 2;
+    coord_t sample_y2 = (tile_width_y + 1) / 2;
+
+    auto extend_left   = tile_width_x;
+    auto extend_top    = tile_width_y;
+    auto extend_right  = (tile_width_x - width % tile_width_x) % tile_width_x + tile_width_x;
+    auto extend_bottom = (tile_width_y - height % tile_width_y) % tile_width_y + tile_width_y;
+
+    auto new_width  = width + extend_left + extend_right;
+    auto new_height = height + extend_top + extend_bottom;
+
+    image<typename SrcView::value_type> padded_img(new_width, new_height);
+
+    auto top_left_x     = tile_width_x;
+    auto top_left_y     = tile_width_y;
+    auto bottom_right_x = tile_width_x + width;
+    auto bottom_right_y = tile_width_y + height;
+
+    copy_pixels(src_view, subimage_view(view(padded_img), top_left_x, top_left_y, width, height));
+
+    for (std::size_t k = 0; k < channels; k++)
+    {
+        std::vector<histogram<source_channel_t>> prev_row(new_width / tile_width_x),
+            next_row((new_width / tile_width_x));
+        std::vector<std::map<source_channel_t, source_channel_t>> prev_map(
+            new_width / tile_width_x),
+            next_map((new_width / tile_width_x));
+
+        coord_t prev = 0, next = 1;
+        auto channel_view = nth_channel_view(view(padded_img), k);
+
+        for (std::ptrdiff_t i = top_left_y; i < bottom_right_y; ++i)
+        {
+            if ((i - sample_y1) / tile_width_y >= next || i == top_left_y)
+            {
+                if (i != top_left_y)
+                {
+                    prev = next;
+                    next++;
+                }
+                prev_row = next_row;
+                prev_map = next_map;
+                for (std::ptrdiff_t j = sample_x1; j < new_width; j += tile_width_x)
+                {
+                    auto img_view = subimage_view(
+                        channel_view, j - sample_x1, next * tile_width_y,
+                        std::max<int>(
+                            std::min<int>(tile_width_x + j - sample_x1, bottom_right_x) -
+                                (j - sample_x1),
+                            0),
+                        std::max<int>(
+                            std::min<int>((next + 1) * tile_width_y, bottom_right_y) -
+                                next * tile_width_y,
+                            0));
+
+                    fill_histogram(
+                        img_view, next_row[(j - sample_x1) / tile_width_x], bin_width, false,
+                        false);
+
+                    detail::clip_and_redistribute(
+                        next_row[(j - sample_x1) / tile_width_x],
+                        next_row[(j - sample_x1) / tile_width_x], clip_limit);
+
+                    next_map[(j - sample_x1) / tile_width_x] =
+                        histogram_equalization(next_row[(j - sample_x1) / tile_width_x]);
+                }
+            }
+            bool prev_row_mask = 1, next_row_mask = 1;
+            if (prev == 0)
+                prev_row_mask = false;
+            else if (next + 1 == new_height / tile_width_y)
+                next_row_mask = false;
+            for (std::ptrdiff_t j = top_left_x; j < bottom_right_x; ++j)
+            {
+                bool prev_col_mask = true, next_col_mask = true;
+                if ((j - sample_x1) / tile_width_x == 0)
+                    prev_col_mask = false;
+                else if ((j - sample_x1) / tile_width_x + 1 == new_width / tile_width_x - 1)
+                    next_col_mask = false;
+
+                // Bilinear interpolation
+                point_t top_left(
+                    (j - sample_x1) / tile_width_x * tile_width_x + sample_x1,
+                                    prev * tile_width_y + sample_y1);
+                point_t top_right(top_left.x + tile_width_x, top_left.y);
+                point_t bottom_left(top_left.x, top_left.y + tile_width_y);
+                point_t bottom_right(top_left.x + tile_width_x, top_left.y + tile_width_y);
+
+                long double x_diff = top_right.x - top_left.x;
+                long double y_diff = bottom_left.y - top_left.y;
+
+                long double x1 = (j - top_left.x) / x_diff;
+                long double x2 = (top_right.x - j) / x_diff;
+                long double y1 = (i - top_left.y) / y_diff;
+                long double y2 = (bottom_left.y - i) / y_diff;
+
+                if (prev_row_mask == 0)
+                    y1 = 1;
+                else if (next_row_mask == 0)
+                    y2 = 1;
+                if (prev_col_mask == 0)
+                    x1 = 1;
+                else if (next_col_mask == 0)
+                    x2 = 1;
+
+                long double numerator =
+                    ((prev_row_mask & prev_col_mask) * x2 *
+                         prev_map[(top_left.x - sample_x1) / tile_width_x][channel_view(j, i)] +
+                     (prev_row_mask & next_col_mask) * x1 *
+                         prev_map[(top_right.x - sample_x1) / tile_width_x][channel_view(j, i)]) *
+                        y2 +
+                    ((next_row_mask & prev_col_mask) * x2 *
+                         next_map[(bottom_left.x - sample_x1) / tile_width_x][channel_view(j, i)] +
+                     (next_row_mask & next_col_mask) * x1 *
+                         next_map[(bottom_right.x - sample_x1) / tile_width_x][channel_view(j, i)]) *
+                        y1;
+
+                if (mask && !src_mask[i - top_left_y][j - top_left_x])
+                {
+                    dst_view(j - top_left_x, i - top_left_y) = 
+                        channel_convert<dst_channel_t>(
+                            static_cast<source_channel_t>(channel_view(i, j)));
+                }
+                else
+                {
+                    dst_view(j - top_left_x, i - top_left_y) = 
+                        channel_convert<dst_channel_t>(static_cast<source_channel_t>(numerator));
+                }
+            }
+        }
+    }
+}
+
+}}  //namespace boost::gil
+
+#endif