ceph/src/boost/boost/gil/image_processing/threshold.hpp

   1 //
   2 // Copyright 2019 Miral Shah <miralshah2211@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 #ifndef BOOST_GIL_IMAGE_PROCESSING_THRESHOLD_HPP
   9 #define BOOST_GIL_IMAGE_PROCESSING_THRESHOLD_HPP
  10
  11 #include <limits>
  12 #include <array>
  13 #include <type_traits>
  14 #include <cstddef>
  15 #include <algorithm>
  16 #include <vector>
  17 #include <cmath>
  18
  19 #include <boost/assert.hpp>
  20
  21 #include <boost/gil/image.hpp>
  22 #include <boost/gil/extension/numeric/kernel.hpp>
  23 #include <boost/gil/extension/numeric/convolve.hpp>
  24 #include <boost/gil/image_processing/numeric.hpp>
  25
  26 namespace boost { namespace gil {
  27
  28 namespace detail {
  29
  30 template
  31 <
  32     typename SourceChannelT,
  33     typename ResultChannelT,
  34     typename SrcView,
  35     typename DstView,
  36     typename Operator
  37 >
  38 void threshold_impl(SrcView const& src_view, DstView const& dst_view, Operator const& threshold_op)
  39 {
  40     gil_function_requires<ImageViewConcept<SrcView>>();
  41     gil_function_requires<MutableImageViewConcept<DstView>>();
  42     static_assert(color_spaces_are_compatible
  43     <
  44         typename color_space_type<SrcView>::type,
  45         typename color_space_type<DstView>::type
  46     >::value, "Source and destination views must have pixels with the same color space");
  47
  48     //iterate over the image chaecking each pixel value for the threshold
  49     for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
  50     {
  51         typename SrcView::x_iterator src_it = src_view.row_begin(y);
  52         typename DstView::x_iterator dst_it = dst_view.row_begin(y);
  53
  54         for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
  55         {
  56             static_transform(src_it[x], dst_it[x], threshold_op);
  57         }
  58     }
  59 }
  60
  61 } //namespace boost::gil::detail
  62
  63 /// \addtogroup ImageProcessing
  64 /// @{
  65 ///
  66 /// \brief Direction of image segmentation.
  67 /// The direction specifieds which pixels are considered as corresponding to object
  68 /// and which pixels correspond to background.
  69 enum class threshold_direction
  70 {
  71     regular, ///< Consider values greater than threshold value
  72     inverse  ///< Consider values less than or equal to threshold value
  73 };
  74
  75 /// \ingroup ImageProcessing
  76 /// \brief Method of optimal threshold value calculation.
  77 enum class threshold_optimal_value
  78 {
  79     otsu        ///< \todo TODO
  80 };
  81
  82 /// \ingroup ImageProcessing
  83 /// \brief TODO
  84 enum class threshold_truncate_mode
  85 {
  86     threshold,  ///< \todo TODO
  87     zero        ///< \todo TODO
  88 };
  89
  90 enum class threshold_adaptive_method
  91 {
  92     mean,
  93     gaussian
  94 };
  95
  96 /// \ingroup ImageProcessing
  97 /// \brief Applies fixed threshold to each pixel of image view.
  98 /// Performs image binarization by thresholding channel value of each
  99 /// pixel of given image view.
 100 /// \param src_view - TODO
 101 /// \param dst_view - TODO
 102 /// \param threshold_value - TODO
 103 /// \param max_value - TODO
 104 /// \param threshold_direction - if regular, values greater than threshold_value are
 105 /// set to max_value else set to 0; if inverse, values greater than threshold_value are
 106 /// set to 0 else set to max_value.
 107 template <typename SrcView, typename DstView>
 108 void threshold_binary(
 109     SrcView const& src_view,
 110     DstView const& dst_view,
 111     typename channel_type<DstView>::type threshold_value,
 112     typename channel_type<DstView>::type max_value,
 113     threshold_direction direction = threshold_direction::regular
 114 )
 115 {
 116     //deciding output channel type and creating functor
 117     using source_channel_t = typename channel_type<SrcView>::type;
 118     using result_channel_t = typename channel_type<DstView>::type;
 119
 120     if (direction == threshold_direction::regular)
 121     {
 122         detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
 123             [threshold_value, max_value](source_channel_t px) -> result_channel_t {
 124                 return px > threshold_value ? max_value : 0;
 125             });
 126     }
 127     else
 128     {
 129         detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
 130             [threshold_value, max_value](source_channel_t px) -> result_channel_t {
 131                 return px > threshold_value ? 0 : max_value;
 132             });
 133     }
 134 }
 135
 136 /// \ingroup ImageProcessing
 137 /// \brief Applies fixed threshold to each pixel of image view.
 138 /// Performs image binarization by thresholding channel value of each
 139 /// pixel of given image view.
 140 /// This variant of threshold_binary automatically deduces maximum value for each channel
 141 /// of pixel based on channel type.
 142 /// If direction is regular, values greater than threshold_value will be set to maximum
 143 /// numeric limit of channel else 0.
 144 /// If direction is inverse, values greater than threshold_value will be set to 0 else maximum
 145 /// numeric limit of channel.
 146 template <typename SrcView, typename DstView>
 147 void threshold_binary(
 148     SrcView const& src_view,
 149     DstView const& dst_view,
 150     typename channel_type<DstView>::type threshold_value,
 151     threshold_direction direction = threshold_direction::regular
 152 )
 153 {
 154     //deciding output channel type and creating functor
 155     using result_channel_t = typename channel_type<DstView>::type;
 156
 157     result_channel_t max_value = (std::numeric_limits<result_channel_t>::max)();
 158     threshold_binary(src_view, dst_view, threshold_value, max_value, direction);
 159 }
 160
 161 /// \ingroup ImageProcessing
 162 /// \brief Applies truncating threshold to each pixel of image view.
 163 /// Takes an image view and performes truncating threshold operation on each chennel.
 164 /// If mode is threshold and direction is regular:
 165 /// values greater than threshold_value will be set to threshold_value else no change
 166 /// If mode is threshold and direction is inverse:
 167 /// values less than threshold_value will be set to threshold_value else no change
 168 /// If mode is zero and direction is regular:
 169 /// values less than threshold_value will be set to 0 else no change
 170 /// If mode is zero and direction is inverse:
 171 /// values more than threshold_value will be set to 0 else no change
 172 template <typename SrcView, typename DstView>
 173 void threshold_truncate(
 174     SrcView const& src_view,
 175     DstView const& dst_view,
 176     typename channel_type<DstView>::type threshold_value,
 177     threshold_truncate_mode mode = threshold_truncate_mode::threshold,
 178     threshold_direction direction = threshold_direction::regular
 179 )
 180 {
 181     //deciding output channel type and creating functor
 182     using source_channel_t = typename channel_type<SrcView>::type;
 183     using result_channel_t = typename channel_type<DstView>::type;
 184
 185     std::function<result_channel_t(source_channel_t)> threshold_logic;
 186
 187     if (mode == threshold_truncate_mode::threshold)
 188     {
 189         if (direction == threshold_direction::regular)
 190         {
 191             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
 192                 [threshold_value](source_channel_t px) -> result_channel_t {
 193                     return px > threshold_value ? threshold_value : px;
 194                 });
 195         }
 196         else
 197         {
 198             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
 199                 [threshold_value](source_channel_t px) -> result_channel_t {
 200                     return px > threshold_value ? px : threshold_value;
 201                 });
 202         }
 203     }
 204     else
 205     {
 206         if (direction == threshold_direction::regular)
 207         {
 208             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
 209                 [threshold_value](source_channel_t px) -> result_channel_t {
 210                     return px > threshold_value ? px : 0;
 211                 });
 212         }
 213         else
 214         {
 215             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
 216                 [threshold_value](source_channel_t px) -> result_channel_t {
 217                     return px > threshold_value ? 0 : px;
 218                 });
 219         }
 220     }
 221 }
 222
 223 namespace detail{
 224
 225 template <typename SrcView, typename DstView>
 226 void otsu_impl(SrcView const& src_view, DstView const& dst_view, threshold_direction direction)
 227 {
 228     //deciding output channel type and creating functor
 229     using source_channel_t = typename channel_type<SrcView>::type;
 230
 231     std::array<std::size_t, 256> histogram{};
 232     //initial value of min is set to maximum possible value to compare histogram data
 233     //initial value of max is set to minimum possible value to compare histogram data
 234     auto min = (std::numeric_limits<source_channel_t>::max)(),
 235         max = (std::numeric_limits<source_channel_t>::min)();
 236
 237     if (sizeof(source_channel_t) > 1 || std::is_signed<source_channel_t>::value)
 238     {
 239         //iterate over the image to find the min and max pixel values
 240         for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
 241         {
 242             typename SrcView::x_iterator src_it = src_view.row_begin(y);
 243             for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
 244             {
 245                 if (src_it[x] < min) min = src_it[x];
 246                 if (src_it[x] > min) min = src_it[x];
 247             }
 248         }
 249
 250         //making histogram
 251         for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
 252         {
 253             typename SrcView::x_iterator src_it = src_view.row_begin(y);
 254
 255             for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
 256             {
 257                 histogram[((src_it[x] - min) * 255) / (max - min)]++;
 258             }
 259         }
 260     }
 261     else
 262     {
 263         //making histogram
 264         for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
 265         {
 266             typename SrcView::x_iterator src_it = src_view.row_begin(y);
 267
 268             for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
 269             {
 270                 histogram[src_it[x]]++;
 271             }
 272         }
 273     }
 274
 275     //histData = histogram data
 276     //sum = total (background + foreground)
 277     //sumB = sum background
 278     //wB = weight background
 279     //wf = weight foreground
 280     //varMax = tracking the maximum known value of between class variance
 281     //mB = mu background
 282     //mF = mu foreground
 283     //varBeetween = between class variance
 284     //http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html
 285     //https://www.ipol.im/pub/art/2016/158/
 286     std::ptrdiff_t total_pixel = src_view.height() * src_view.width();
 287     std::ptrdiff_t sum_total = 0, sum_back = 0;
 288     std::size_t weight_back = 0, weight_fore = 0, threshold = 0;
 289     double var_max = 0, mean_back, mean_fore, var_intra_class;
 290
 291     for (std::size_t t = 0; t < 256; t++)
 292     {
 293         sum_total += t * histogram[t];
 294     }
 295
 296     for (int t = 0; t < 256; t++)
 297     {
 298         weight_back += histogram[t];               // Weight Background
 299         if (weight_back == 0) continue;
 300
 301         weight_fore = total_pixel - weight_back;          // Weight Foreground
 302         if (weight_fore == 0) break;
 303
 304         sum_back += t * histogram[t];
 305
 306         mean_back = sum_back / weight_back;            // Mean Background
 307         mean_fore = (sum_total - sum_back) / weight_fore;    // Mean Foreground
 308
 309         // Calculate Between Class Variance
 310         var_intra_class = weight_back * weight_fore * (mean_back - mean_fore) * (mean_back - mean_fore);
 311
 312         // Check if new maximum found
 313         if (var_intra_class > var_max) {
 314             var_max = var_intra_class;
 315             threshold = t;
 316         }
 317     }
 318     if (sizeof(source_channel_t) > 1 && std::is_unsigned<source_channel_t>::value)
 319     {
 320         threshold_binary(src_view, dst_view, (threshold * (max - min) / 255) + min, direction);
 321     }
 322     else {
 323         threshold_binary(src_view, dst_view, threshold, direction);
 324     }
 325 }
 326 } //namespace detail
 327
 328 template <typename SrcView, typename DstView>
 329 void threshold_optimal
 330 (
 331     SrcView const& src_view,
 332     DstView const& dst_view,
 333     threshold_optimal_value mode = threshold_optimal_value::otsu,
 334     threshold_direction direction = threshold_direction::regular
 335 )
 336 {
 337     if (mode == threshold_optimal_value::otsu)
 338     {
 339         for (std::size_t i = 0; i < src_view.num_channels(); i++)
 340         {
 341             detail::otsu_impl
 342                 (nth_channel_view(src_view, i), nth_channel_view(dst_view, i), direction);
 343         }
 344     }
 345 }
 346
 347 namespace detail {
 348
 349 template
 350 <
 351     typename SourceChannelT,
 352     typename ResultChannelT,
 353     typename SrcView,
 354     typename DstView,
 355     typename Operator
 356 >
 357 void adaptive_impl
 358 (
 359     SrcView const& src_view,
 360     SrcView const& convolved_view,
 361     DstView const& dst_view,
 362     Operator const& threshold_op
 363 )
 364 {
 365     //template argument validation
 366     gil_function_requires<ImageViewConcept<SrcView>>();
 367     gil_function_requires<MutableImageViewConcept<DstView>>();
 368
 369     static_assert(color_spaces_are_compatible
 370     <
 371         typename color_space_type<SrcView>::type,
 372         typename color_space_type<DstView>::type
 373     >::value, "Source and destination views must have pixels with the same color space");
 374
 375     //iterate over the image chaecking each pixel value for the threshold
 376     for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
 377     {
 378         typename SrcView::x_iterator src_it = src_view.row_begin(y);
 379         typename SrcView::x_iterator convolved_it = convolved_view.row_begin(y);
 380         typename DstView::x_iterator dst_it = dst_view.row_begin(y);
 381
 382         for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
 383         {
 384             static_transform(src_it[x], convolved_it[x], dst_it[x], threshold_op);
 385         }
 386     }
 387 }
 388 } //namespace boost::gil::detail
 389
 390 template <typename SrcView, typename DstView>
 391 void threshold_adaptive
 392 (
 393     SrcView const& src_view,
 394     DstView const& dst_view,
 395     typename channel_type<DstView>::type max_value,
 396     std::size_t kernel_size,
 397     threshold_adaptive_method method = threshold_adaptive_method::mean,
 398     threshold_direction direction = threshold_direction::regular,
 399     typename channel_type<DstView>::type constant = 0
 400 )
 401 {
 402     BOOST_ASSERT_MSG((kernel_size % 2 != 0), "Kernel size must be an odd number");
 403
 404     typedef typename channel_type<SrcView>::type source_channel_t;
 405     typedef typename channel_type<DstView>::type result_channel_t;
 406
 407     image<typename SrcView::value_type> temp_img(src_view.width(), src_view.height());
 408     typename image<typename SrcView::value_type>::view_t temp_view = view(temp_img);
 409     SrcView temp_conv(temp_view);
 410
 411     if (method == threshold_adaptive_method::mean)
 412     {
 413         std::vector<float> mean_kernel_values(kernel_size, 1.0f/kernel_size);
 414         kernel_1d<float> kernel(mean_kernel_values.begin(), kernel_size, kernel_size/2);
 415
 416         detail::convolve_1d
 417         <
 418             pixel<float, typename SrcView::value_type::layout_t>
 419         >(src_view, kernel, temp_view);
 420     }
 421     else if (method == threshold_adaptive_method::gaussian)
 422     {
 423         detail::kernel_2d<float> kernel = generate_gaussian_kernel(kernel_size, 1.0);
 424         convolve_2d(src_view, kernel, temp_view);
 425     }
 426
 427     if (direction == threshold_direction::regular)
 428     {
 429         detail::adaptive_impl<source_channel_t, result_channel_t>(src_view, temp_conv, dst_view,
 430             [max_value, constant](source_channel_t px, source_channel_t threshold) -> result_channel_t
 431         { return px > (threshold - constant) ? max_value : 0; });
 432     }
 433     else
 434     {
 435         detail::adaptive_impl<source_channel_t, result_channel_t>(src_view, temp_conv, dst_view,
 436             [max_value, constant](source_channel_t px, source_channel_t threshold) -> result_channel_t
 437         { return px > (threshold - constant) ? 0 : max_value; });
 438     }
 439 }
 440
 441 template <typename SrcView, typename DstView>
 442 void threshold_adaptive
 443 (
 444     SrcView const& src_view,
 445     DstView const& dst_view,
 446     std::size_t kernel_size,
 447     threshold_adaptive_method method = threshold_adaptive_method::mean,
 448     threshold_direction direction = threshold_direction::regular,
 449     int constant = 0
 450 )
 451 {
 452     //deciding output channel type and creating functor
 453     typedef typename channel_type<DstView>::type result_channel_t;
 454
 455     result_channel_t max_value = (std::numeric_limits<result_channel_t>::max)();
 456
 457     threshold_adaptive(src_view, dst_view, max_value, kernel_size, method, direction, constant);
 458 }
 459
 460 /// @}
 461
 462 }} //namespace boost::gil
 463
 464 #endif //BOOST_GIL_IMAGE_PROCESSING_THRESHOLD_HPP