303 lines
12 KiB
C++
303 lines
12 KiB
C++
//
|
|
// 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 (std::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::ptrdiff_t tile_width_x = 20,
|
|
std::ptrdiff_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();
|
|
|
|
// Find control points
|
|
|
|
std::vector<coord_t> sample_x;
|
|
coord_t sample_x1 = tile_width_x / 2;
|
|
coord_t sample_y1 = tile_width_y / 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
|