123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302 |
- #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 {
- namespace detail {
- 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;
- }
- 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);
-
- 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();
- }
- }
- }
- 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();
-
- 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;
-
- 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));
- }
- }
- }
- }
- }
- }}
- #endif
|