diff --git a/example/adaptive_he.cpp b/example/adaptive_he.cpp new file mode 100644 index 0000000000..0610932c98 --- /dev/null +++ b/example/adaptive_he.cpp @@ -0,0 +1,25 @@ +// +// 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; + +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; +} diff --git a/include/boost/gil/image_processing/adaptive_histogram_equalization.hpp b/include/boost/gil/image_processing/adaptive_histogram_equalization.hpp new file mode 100644 index 0000000000..9dcbfba928 --- /dev/null +++ b/include/boost/gil/image_processing/adaptive_histogram_equalization.hpp @@ -0,0 +1,307 @@ +// +// 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_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP +#define BOOST_GIL_IMAGE_PROCESSING_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP + +#include +#include +#include +#include +#include + +#include +#include +#include + +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 +/// __.._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 +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 +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 +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> src_mask = {}) +{ + 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 destination views must have same color space"); + + using source_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(); + std::size_t pixel_max = std::numeric_limits::max(); + std::size_t pixel_min = std::numeric_limits::min(); + + // Find control points + + std::vector 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 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> prev_row(new_width / tile_width_x), + next_row((new_width / tile_width_x)); + std::vector> 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( + std::min(tile_width_x + j - sample_x1, bottom_right_x) - + (j - sample_x1), + 0), + std::max( + std::min((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( + static_cast(channel_view(i, j))); + } + else + { + dst_view(j - top_left_x, i - top_left_y) = + channel_convert(static_cast(numerator)); + } + } + } + } +} + +}} //namespace boost::gil + +#endif diff --git a/test/core/image_processing/adaptive_he.cpp b/test/core/image_processing/adaptive_he.cpp new file mode 100644 index 0000000000..df7068243f --- /dev/null +++ b/test/core/image_processing/adaptive_he.cpp @@ -0,0 +1,113 @@ +// +// 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 + +namespace gil = boost::gil; + +double epsilon = 1.0; + +std::uint8_t image_matrix[] = +{ + 1, 1, 1, 1, + 3, 3, 3, 3, + 5, 5, 5, 5, + 7, 7, 7, 7 +}; +gil::gray8c_view_t gray_view = gil::interleaved_view(4, 4, reinterpret_cast(image_matrix), 4); + +void check_actual_clip_limit() +{ + gil::histogram h; + for(std::size_t i = 0; i < 100; i++) + { + if (i % 40 == 0) + { + h(i) = 60; + } + else + { + h(i) = 5; + } + } + double limit = 0.01; + double value = gil::detail::actual_clip_limit(h, limit); + + long actual_limit = round(value * h.sum()), max_bin_val = 0; + double excess = 0; + for(std::size_t i = 0; i < 100; i++) + { + if (h(i) > actual_limit) + excess += actual_limit - h(i); + max_bin_val = std::max(max_bin_val, h(i)); + } + BOOST_TEST((abs(excess / h.size() + actual_limit) - limit * h.sum()) < epsilon); +} + +void check_clip_and_redistribute() +{ + gil::histogram h, h2; + for(std::size_t i = 0; i < 100; i++) + { + if (i % 50 == 0) + { + h(i) = 60; + } + else + { + h(i) = 5; + } + } + bool check = true; + double limit = 0.001; + gil::detail::clip_and_redistribute(h, h2, limit); + for(std::size_t i = 0; i < 100; i++) + { + check = check & (abs(limit * h.sum() - h2(i)) < epsilon); + } + BOOST_TEST(check); +} + +void check_non_overlapping_interpolated_clahe() +{ + { + gil::gray8_image_t img1(4, 4), img2(4, 4), img3(4, 4); + gil::histogram_equalization(gray_view, view(img2)); + gil::non_overlapping_interpolated_clahe(gray_view, view(img3), 8, 8, 1.0); + BOOST_TEST(gil::equal_pixels(view(img2), view(img3))); + } + { + gil::gray8_image_t img1(8, 8), img2(8, 8), img3(4, 4); + gil::copy_pixels(gray_view, gil::subimage_view(view(img1), 0, 0, 4, 4)); + gil::copy_pixels(gray_view, gil::subimage_view(view(img1), 0, 4, 4, 4)); + gil::copy_pixels(gray_view, gil::subimage_view(view(img1), 4, 0, 4, 4)); + gil::copy_pixels(gray_view, gil::subimage_view(view(img1), 4, 4, 4, 4)); + gil::histogram_equalization(gray_view, view(img3)); + gil::non_overlapping_interpolated_clahe(view(img1), view(img2), 8, 8, 1.0); + BOOST_TEST(gil::equal_pixels(gil::subimage_view(view(img2), 0, 0, 4, 4), view(img3))); + BOOST_TEST(gil::equal_pixels(gil::subimage_view(view(img2), 0, 4, 4, 4), view(img3))); + BOOST_TEST(gil::equal_pixels(gil::subimage_view(view(img2), 4, 0, 4, 4), view(img3))); + BOOST_TEST(gil::equal_pixels(gil::subimage_view(view(img2), 4, 4, 4, 4), view(img3))); + } +} + +int main() +{ + check_actual_clip_limit(); + check_clip_and_redistribute(); + check_non_overlapping_interpolated_clahe(); + + return boost::report_errors(); +}