Skip to content

Commit

Permalink
Add code for ahe algorithm (#516)
Browse files Browse the repository at this point in the history
  • Loading branch information
@codejaeger
codejaeger committed Jan 24, 2021
1 parent b82aed8 commit a68a95d
Show file tree
Hide file tree
Showing 3 changed files with 445 additions and 0 deletions.
25 changes: 25 additions & 0 deletions example/adaptive_he.cpp
View file
@@ -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;
}
307 changes: 307 additions & 0 deletions include/boost/gil/image_processing/adaptive_histogram_equalization.hpp
View file
@@ -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

0 comments on commit a68a95d

Please sign in to comment.