Generic Image Library: adaptive_histogram_equalization.hpp Source File

###  | #  |

adaptive_histogram_equalization.hpp
1//
2// Copyright 2020 Debabrata Mandal <mandaldebabrata123@gmail.com>
3//
4// Use, modification and distribution are subject to the Boost Software License,
5// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
6// http://www.boost.org/LICENSE_1_0.txt)
7//
8
9#ifndef BOOST_GIL_IMAGE_PROCESSING_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP
10#define BOOST_GIL_IMAGE_PROCESSING_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP
11
12#include <boost/gil/algorithm.hpp>
13#include <boost/gil/histogram.hpp>
14#include <boost/gil/image.hpp>
15#include <boost/gil/image_processing/histogram_equalization.hpp>
16#include <boost/gil/image_view_factory.hpp>
17
18#include <cmath>
19#include <map>
20#include <vector>
21
22namespace boost { namespace gil {
23
36
37namespace detail {
38
41
48template <typename SrcHist>
49double actual_clip_limit(SrcHist const& src_hist, double cliplimit = 0.03)
50{
51 double epsilon = 1.0;
52 using value_t = typename SrcHist::value_type;
53 double sum = src_hist.sum();
54 std::size_t num_bins = src_hist.size();
55
56 cliplimit = sum * cliplimit;
57 long low = 0, high = cliplimit, middle = low;
58 while (high - low >= 1)
59 {
60 middle = (low + high + 1) >> 1;
61 long excess = 0;
62 std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
63 if (v.second > middle)
64 excess += v.second - middle;
65 });
66 if (std::abs(excess - (cliplimit - middle) * num_bins) < epsilon)
67 break;
68 else if (excess > (cliplimit - middle) * num_bins)
69 high = middle - 1;
70 else
71 low = middle + 1;
72 }
73 return middle / sum;
74}
75
83template <typename SrcHist, typename DstHist>
84void clip_and_redistribute(SrcHist const& src_hist, DstHist& dst_hist, double clip_limit = 0.03)
85{
86 using value_t = typename SrcHist::value_type;
87 double sum = src_hist.sum();
88 double actual_clip_value = detail::actual_clip_limit(src_hist, clip_limit);
89 // double actual_clip_value = clip_limit;
90 long actual_clip_limit = actual_clip_value * sum;
91 double excess = 0;
92 std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
93 if (v.second > actual_clip_limit)
94 excess += v.second - actual_clip_limit;
95 });
96 std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
97 if (v.second >= actual_clip_limit)
98 dst_hist[dst_hist.key_from_tuple(v.first)] = clip_limit * sum;
99 else
100 dst_hist[dst_hist.key_from_tuple(v.first)] = v.second + excess / src_hist.size();
101 });
102 long rem = long(excess) % src_hist.size();
103 if (rem == 0)
104 return;
105 long period = round(src_hist.size() / rem);
106 std::size_t index = 0;
107 while (rem)
108 {
109 if (dst_hist(index) >= clip_limit * sum)
110 {
111 index = (index + 1) % src_hist.size();
112 }
113 dst_hist(index)++;
114 rem--;
115 index = (index + period) % src_hist.size();
116 }
117}
118
119} // namespace detail
120
136template <typename SrcView, typename DstView>
137void non_overlapping_interpolated_clahe(
138 SrcView const& src_view,
139 DstView const& dst_view,
140 std::ptrdiff_t tile_width_x = 20,
141 std::ptrdiff_t tile_width_y = 20,
142 double clip_limit = 0.03,
143 std::size_t bin_width = 1.0,
144 bool mask = false,
145 std::vector<std::vector<bool>> src_mask = {})
146{
147 gil_function_requires<ImageViewConcept<SrcView>>();
148 gil_function_requires<MutableImageViewConcept<DstView>>();
149
150 static_assert(
151 color_spaces_are_compatible<
152 typename color_space_type<SrcView>::type,
153 typename color_space_type<DstView>::type>::value,
154 "Source and destination views must have same color space");
155
156 using source_channel_t = typename channel_type<SrcView>::type;
157 using dst_channel_t = typename channel_type<DstView>::type;
158 using coord_t = typename SrcView::x_coord_t;
159
160 std::size_t const channels = num_channels<SrcView>::value;
161 coord_t const width = src_view.width();
162 coord_t const height = src_view.height();
163
164 // Find control points
165
166 std::vector<coord_t> sample_x;
167 coord_t sample_x1 = tile_width_x / 2;
168 coord_t sample_y1 = tile_width_y / 2;
169
170 auto extend_left = tile_width_x;
171 auto extend_top = tile_width_y;
172 auto extend_right = (tile_width_x - width % tile_width_x) % tile_width_x + tile_width_x;
173 auto extend_bottom = (tile_width_y - height % tile_width_y) % tile_width_y + tile_width_y;
174
175 auto new_width = width + extend_left + extend_right;
176 auto new_height = height + extend_top + extend_bottom;
177
178 image<typename SrcView::value_type> padded_img(new_width, new_height);
179
180 auto top_left_x = tile_width_x;
181 auto top_left_y = tile_width_y;
182 auto bottom_right_x = tile_width_x + width;
183 auto bottom_right_y = tile_width_y + height;
184
185 copy_pixels(src_view, subimage_view(view(padded_img), top_left_x, top_left_y, width, height));
186
187 for (std::size_t k = 0; k < channels; k++)
188 {
189 std::vector<histogram<source_channel_t>> prev_row(new_width / tile_width_x),
190 next_row((new_width / tile_width_x));
191 std::vector<std::map<source_channel_t, source_channel_t>> prev_map(
192 new_width / tile_width_x),
193 next_map((new_width / tile_width_x));
194
195 coord_t prev = 0, next = 1;
196 auto channel_view = nth_channel_view(view(padded_img), k);
197
198 for (std::ptrdiff_t i = top_left_y; i < bottom_right_y; ++i)
199 {
200 if ((i - sample_y1) / tile_width_y >= next || i == top_left_y)
201 {
202 if (i != top_left_y)
203 {
204 prev = next;
205 next++;
206 }
207 prev_row = next_row;
208 prev_map = next_map;
209 for (std::ptrdiff_t j = sample_x1; j < new_width; j += tile_width_x)
210 {
211 auto img_view = subimage_view(
212 channel_view, j - sample_x1, next * tile_width_y,
213 std::max<int>(
214 std::min<int>(tile_width_x + j - sample_x1, bottom_right_x) -
215 (j - sample_x1),
216 0),
217 std::max<int>(
218 std::min<int>((next + 1) * tile_width_y, bottom_right_y) -
219 next * tile_width_y,
220 0));
221
222 fill_histogram(
223 img_view, next_row[(j - sample_x1) / tile_width_x], bin_width, false,
224 false);
225
226 detail::clip_and_redistribute(
227 next_row[(j - sample_x1) / tile_width_x],
228 next_row[(j - sample_x1) / tile_width_x], clip_limit);
229
230 next_map[(j - sample_x1) / tile_width_x] =
231 histogram_equalization(next_row[(j - sample_x1) / tile_width_x]);
232 }
233 }
234 bool prev_row_mask = 1, next_row_mask = 1;
235 if (prev == 0)
236 prev_row_mask = false;
237 else if (next + 1 == new_height / tile_width_y)
238 next_row_mask = false;
239 for (std::ptrdiff_t j = top_left_x; j < bottom_right_x; ++j)
240 {
241 bool prev_col_mask = true, next_col_mask = true;
242 if ((j - sample_x1) / tile_width_x == 0)
243 prev_col_mask = false;
244 else if ((j - sample_x1) / tile_width_x + 1 == new_width / tile_width_x - 1)
245 next_col_mask = false;
246
247 // Bilinear interpolation
248 point_t top_left(
249 (j - sample_x1) / tile_width_x * tile_width_x + sample_x1,
250 prev * tile_width_y + sample_y1);
251 point_t top_right(top_left.x + tile_width_x, top_left.y);
252 point_t bottom_left(top_left.x, top_left.y + tile_width_y);
253 point_t bottom_right(top_left.x + tile_width_x, top_left.y + tile_width_y);
254
255 long double x_diff = top_right.x - top_left.x;
256 long double y_diff = bottom_left.y - top_left.y;
257
258 long double x1 = (j - top_left.x) / x_diff;
259 long double x2 = (top_right.x - j) / x_diff;
260 long double y1 = (i - top_left.y) / y_diff;
261 long double y2 = (bottom_left.y - i) / y_diff;
262
263 if (prev_row_mask == 0)
264 y1 = 1;
265 else if (next_row_mask == 0)
266 y2 = 1;
267 if (prev_col_mask == 0)
268 x1 = 1;
269 else if (next_col_mask == 0)
270 x2 = 1;
271
272 long double numerator =
273 ((prev_row_mask & prev_col_mask) * x2 *
274 prev_map[(top_left.x - sample_x1) / tile_width_x][channel_view(j, i)] +
275 (prev_row_mask & next_col_mask) * x1 *
276 prev_map[(top_right.x - sample_x1) / tile_width_x][channel_view(j, i)]) *
277 y2 +
278 ((next_row_mask & prev_col_mask) * x2 *
279 next_map[(bottom_left.x - sample_x1) / tile_width_x][channel_view(j, i)] +
280 (next_row_mask & next_col_mask) * x1 *
281 next_map[(bottom_right.x - sample_x1) / tile_width_x][channel_view(j, i)]) *
282 y1;
283
284 if (mask && !src_mask[i - top_left_y][j - top_left_x])
285 {
286 dst_view(j - top_left_x, i - top_left_y) =
287 channel_convert<dst_channel_t>(
288 static_cast<source_channel_t>(channel_view(i, j)));
289 }
290 else
291 {
292 dst_view(j - top_left_x, i - top_left_y) =
293 channel_convert<dst_channel_t>(static_cast<source_channel_t>(numerator));
294 }
295 }
296 }
297 }
298}
299
300}} //namespace boost::gil
301
302#endif
boost::gil::channel_convert
auto channel_convert(SrcChannel const &src) -> typename channel_traits< DstChannel >::value_type
Converting from one channel type to another.
Definition channel_algorithm.hpp:461
boost::gil::view
auto view(image< Pixel, IsPlanar, Alloc > &img) -> typename image< Pixel, IsPlanar, Alloc >::view_t const &
Returns the non-constant-pixel view of an image.
Definition image.hpp:565
boost::gil::copy_pixels
BOOST_FORCEINLINE void copy_pixels(const View1 &src, const View2 &dst)
std::copy for image views
Definition algorithm.hpp:292
boost
defined(BOOST_NO_CXX17_HDR_MEMORY_RESOURCE)
Definition algorithm.hpp:36