gil/image_processing/adaptive_histogram_equalization.hpp

0001 //
0002 // Copyright 2020 Debabrata Mandal <mandaldebabrata123@gmail.com>
0003 //
0004 // Use, modification and distribution are subject to the Boost Software License,
0005 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0006 // http://www.boost.org/LICENSE_1_0.txt)
0007 //
0008
0009 #ifndef BOOST_GIL_IMAGE_PROCESSING_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP
0010 #define BOOST_GIL_IMAGE_PROCESSING_ADAPTIVE_HISTOGRAM_EQUALIZATION_HPP
0011
0012 #include <boost/gil/algorithm.hpp>
0013 #include <boost/gil/histogram.hpp>
0014 #include <boost/gil/image.hpp>
0015 #include <boost/gil/image_processing/histogram_equalization.hpp>
0016 #include <boost/gil/image_view_factory.hpp>
0017
0018 #include <cmath>
0019 #include <map>
0020 #include <vector>
0021
0022 namespace boost { namespace gil {
0023
0024 /////////////////////////////////////////
0025 /// Adaptive Histogram Equalization(AHE)
0026 /////////////////////////////////////////
0027 /// \defgroup AHE AHE
0028 /// \brief Contains implementation and description of the algorithm used to compute
0029 ///        adaptive histogram equalization of input images. Naming for the AHE functions
0030 ///        are done in the following way
0031 ///             <feature-1>_<feature-2>_.._<feature-n>ahe
0032 ///        For example, for AHE done using local (non-overlapping) tiles/blocks and
0033 ///        final output interpolated among tiles , it is called
0034 ///             non_overlapping_interpolated_clahe
0035 ///
0036
0037 namespace detail {
0038
0039 /// \defgroup AHE-helpers AHE-helpers
0040 /// \brief AHE helper functions
0041
0042 /// \fn double actual_clip_limit
0043 /// \ingroup AHE-helpers
0044 /// \brief Computes the actual clip limit given a clip limit value using binary search.
0045 ///        Reference -  Adaptive Histogram Equalization and Its Variations
0046 ///                     (http://www.cs.unc.edu/techreports/86-013.pdf, Pg - 15)
0047 ///
0048 template <typename SrcHist>
0049 double actual_clip_limit(SrcHist const& src_hist, double cliplimit = 0.03)
0050 {
0051     double epsilon       = 1.0;
0052     using value_t        = typename SrcHist::value_type;
0053     double sum           = src_hist.sum();
0054     std::size_t num_bins = src_hist.size();
0055
0056     cliplimit = sum * cliplimit;
0057     long low = 0, high = cliplimit, middle = low;
0058     while (high - low >= 1)
0059     {
0060         middle      = (low + high + 1) >> 1;
0061         long excess = 0;
0062         std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
0063             if (v.second > middle)
0064                 excess += v.second - middle;
0065         });
0066         if (std::abs(excess - (cliplimit - middle) * num_bins) < epsilon)
0067             break;
0068         else if (excess > (cliplimit - middle) * num_bins)
0069             high = middle - 1;
0070         else
0071             low = middle + 1;
0072     }
0073     return middle / sum;
0074 }
0075
0076 /// \fn void clip_and_redistribute
0077 /// \ingroup AHE-helpers
0078 /// \brief Clips and redistributes excess pixels based on the actual clip limit value
0079 ///        obtained from the other helper function actual_clip_limit
0080 ///        Reference - Graphic Gems 4, Pg. 474
0081 ///        (http://cas.xav.free.fr/Graphics%20Gems%204%20-%20Paul%20S.%20Heckbert.pdf)
0082 ///
0083 template <typename SrcHist, typename DstHist>
0084 void clip_and_redistribute(SrcHist const& src_hist, DstHist& dst_hist, double clip_limit = 0.03)
0085 {
0086     using value_t            = typename SrcHist::value_type;
0087     double sum               = src_hist.sum();
0088     double actual_clip_value = detail::actual_clip_limit(src_hist, clip_limit);
0089     // double actual_clip_value = clip_limit;
0090     long actual_clip_limit = actual_clip_value * sum;
0091     double excess          = 0;
0092     std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
0093         if (v.second > actual_clip_limit)
0094             excess += v.second - actual_clip_limit;
0095     });
0096     std::for_each(src_hist.begin(), src_hist.end(), [&](value_t const& v) {
0097         if (v.second >= actual_clip_limit)
0098             dst_hist[dst_hist.key_from_tuple(v.first)] = clip_limit * sum;
0099         else
0100             dst_hist[dst_hist.key_from_tuple(v.first)] = v.second + excess / src_hist.size();
0101     });
0102     long rem = long(excess) % src_hist.size();
0103     if (rem == 0)
0104         return;
0105     long period       = round(src_hist.size() / rem);
0106     std::size_t index = 0;
0107     while (rem)
0108     {
0109         if (dst_hist(index) >= clip_limit * sum)
0110         {
0111             index = (index + 1) % src_hist.size();
0112         }
0113         dst_hist(index)++;
0114         rem--;
0115         index = (index + period) % src_hist.size();
0116     }
0117 }
0118
0119 } // namespace detail
0120
0121 /// \fn void non_overlapping_interpolated_clahe
0122 /// \ingroup AHE
0123 /// @param src_view      Input   Source image view
0124 /// @param dst_view      Output  Output image view
0125 /// @param tile_width_x  Input   Tile width along x-axis to apply HE
0126 /// @param tile_width_y  Input   Tile width along x-axis to apply HE
0127 /// @param clip_limit    Input   Clipping limit to be applied
0128 /// @param bin_width     Input   Bin widths for histogram
0129 /// @param mask          Input   Specify if mask is to be used
0130 /// @param src_mask      Input   Mask on input image to ignore specified pixels
0131 /// \brief Performs local histogram equalization on tiles of size (tile_width_x, tile_width_y)
0132 ///        Then uses the clip limit to redistribute excess pixels above the limit uniformly to
0133 ///        other bins. The clip limit is specified as a fraction i.e. a bin's value is clipped
0134 ///        if bin_value >= clip_limit * (Total number of pixels in the tile)
0135 ///
0136 template <typename SrcView, typename DstView>
0137 void non_overlapping_interpolated_clahe(
0138     SrcView const& src_view,
0139     DstView const& dst_view,
0140     std::ptrdiff_t tile_width_x             = 20,
0141     std::ptrdiff_t tile_width_y             = 20,
0142     double clip_limit                       = 0.03,
0143     std::size_t bin_width                   = 1.0,
0144     bool mask                               = false,
0145     std::vector<std::vector<bool>> src_mask = {})
0146 {
0147     gil_function_requires<ImageViewConcept<SrcView>>();
0148     gil_function_requires<MutableImageViewConcept<DstView>>();
0149
0150     static_assert(
0151         color_spaces_are_compatible<
0152             typename color_space_type<SrcView>::type,
0153             typename color_space_type<DstView>::type>::value,
0154         "Source and destination views must have same color space");
0155
0156     using source_channel_t = typename channel_type<SrcView>::type;
0157     using dst_channel_t    = typename channel_type<DstView>::type;
0158     using coord_t          = typename SrcView::x_coord_t;
0159
0160     std::size_t const channels = num_channels<SrcView>::value;
0161     coord_t const width        = src_view.width();
0162     coord_t const height       = src_view.height();
0163
0164     // Find control points
0165
0166     std::vector<coord_t> sample_x;
0167     coord_t sample_x1 = tile_width_x / 2;
0168     coord_t sample_y1 = tile_width_y / 2;
0169
0170     auto extend_left   = tile_width_x;
0171     auto extend_top    = tile_width_y;
0172     auto extend_right  = (tile_width_x - width % tile_width_x) % tile_width_x + tile_width_x;
0173     auto extend_bottom = (tile_width_y - height % tile_width_y) % tile_width_y + tile_width_y;
0174
0175     auto new_width  = width + extend_left + extend_right;
0176     auto new_height = height + extend_top + extend_bottom;
0177
0178     image<typename SrcView::value_type> padded_img(new_width, new_height);
0179
0180     auto top_left_x     = tile_width_x;
0181     auto top_left_y     = tile_width_y;
0182     auto bottom_right_x = tile_width_x + width;
0183     auto bottom_right_y = tile_width_y + height;
0184
0185     copy_pixels(src_view, subimage_view(view(padded_img), top_left_x, top_left_y, width, height));
0186
0187     for (std::size_t k = 0; k < channels; k++)
0188     {
0189         std::vector<histogram<source_channel_t>> prev_row(new_width / tile_width_x),
0190             next_row((new_width / tile_width_x));
0191         std::vector<std::map<source_channel_t, source_channel_t>> prev_map(
0192             new_width / tile_width_x),
0193             next_map((new_width / tile_width_x));
0194
0195         coord_t prev = 0, next = 1;
0196         auto channel_view = nth_channel_view(view(padded_img), k);
0197
0198         for (std::ptrdiff_t i = top_left_y; i < bottom_right_y; ++i)
0199         {
0200             if ((i - sample_y1) / tile_width_y >= next || i == top_left_y)
0201             {
0202                 if (i != top_left_y)
0203                 {
0204                     prev = next;
0205                     next++;
0206                 }
0207                 prev_row = next_row;
0208                 prev_map = next_map;
0209                 for (std::ptrdiff_t j = sample_x1; j < new_width; j += tile_width_x)
0210                 {
0211                     auto img_view = subimage_view(
0212                         channel_view, j - sample_x1, next * tile_width_y,
0213                         std::max<int>(
0214                             std::min<int>(tile_width_x + j - sample_x1, bottom_right_x) -
0215                                 (j - sample_x1),
0216                             0),
0217                         std::max<int>(
0218                             std::min<int>((next + 1) * tile_width_y, bottom_right_y) -
0219                                 next * tile_width_y,
0220                             0));
0221
0222                     fill_histogram(
0223                         img_view, next_row[(j - sample_x1) / tile_width_x], bin_width, false,
0224                         false);
0225
0226                     detail::clip_and_redistribute(
0227                         next_row[(j - sample_x1) / tile_width_x],
0228                         next_row[(j - sample_x1) / tile_width_x], clip_limit);
0229
0230                     next_map[(j - sample_x1) / tile_width_x] =
0231                         histogram_equalization(next_row[(j - sample_x1) / tile_width_x]);
0232                 }
0233             }
0234             bool prev_row_mask = 1, next_row_mask = 1;
0235             if (prev == 0)
0236                 prev_row_mask = false;
0237             else if (next + 1 == new_height / tile_width_y)
0238                 next_row_mask = false;
0239             for (std::ptrdiff_t j = top_left_x; j < bottom_right_x; ++j)
0240             {
0241                 bool prev_col_mask = true, next_col_mask = true;
0242                 if ((j - sample_x1) / tile_width_x == 0)
0243                     prev_col_mask = false;
0244                 else if ((j - sample_x1) / tile_width_x + 1 == new_width / tile_width_x - 1)
0245                     next_col_mask = false;
0246
0247                 // Bilinear interpolation
0248                 point_t top_left(
0249                     (j - sample_x1) / tile_width_x * tile_width_x + sample_x1,
0250                                     prev * tile_width_y + sample_y1);
0251                 point_t top_right(top_left.x + tile_width_x, top_left.y);
0252                 point_t bottom_left(top_left.x, top_left.y + tile_width_y);
0253                 point_t bottom_right(top_left.x + tile_width_x, top_left.y + tile_width_y);
0254
0255                 long double x_diff = top_right.x - top_left.x;
0256                 long double y_diff = bottom_left.y - top_left.y;
0257
0258                 long double x1 = (j - top_left.x) / x_diff;
0259                 long double x2 = (top_right.x - j) / x_diff;
0260                 long double y1 = (i - top_left.y) / y_diff;
0261                 long double y2 = (bottom_left.y - i) / y_diff;
0262
0263                 if (prev_row_mask == 0)
0264                     y1 = 1;
0265                 else if (next_row_mask == 0)
0266                     y2 = 1;
0267                 if (prev_col_mask == 0)
0268                     x1 = 1;
0269                 else if (next_col_mask == 0)
0270                     x2 = 1;
0271
0272                 long double numerator =
0273                     ((prev_row_mask & prev_col_mask) * x2 *
0274                          prev_map[(top_left.x - sample_x1) / tile_width_x][channel_view(j, i)] +
0275                      (prev_row_mask & next_col_mask) * x1 *
0276                          prev_map[(top_right.x - sample_x1) / tile_width_x][channel_view(j, i)]) *
0277                         y2 +
0278                     ((next_row_mask & prev_col_mask) * x2 *
0279                          next_map[(bottom_left.x - sample_x1) / tile_width_x][channel_view(j, i)] +
0280                      (next_row_mask & next_col_mask) * x1 *
0281                          next_map[(bottom_right.x - sample_x1) / tile_width_x][channel_view(j, i)]) *
0282                         y1;
0283
0284                 if (mask && !src_mask[i - top_left_y][j - top_left_x])
0285                 {
0286                     dst_view(j - top_left_x, i - top_left_y) =
0287                         channel_convert<dst_channel_t>(
0288                             static_cast<source_channel_t>(channel_view(i, j)));
0289                 }
0290                 else
0291                 {
0292                     dst_view(j - top_left_x, i - top_left_y) =
0293                         channel_convert<dst_channel_t>(static_cast<source_channel_t>(numerator));
0294                 }
0295             }
0296         }
0297     }
0298 }
0299
0300 }}  //namespace boost::gil
0301
0302 #endif