gil/image_processing/threshold.hpp

0001 //
0002 // Copyright 2019 Miral Shah <miralshah2211@gmail.com>
0003 // Copyright 2021 Pranam Lashkari <plashkari628@gmail.com>
0004 //
0005 // Use, modification and distribution are subject to the Boost Software License,
0006 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0007 // http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 #ifndef BOOST_GIL_IMAGE_PROCESSING_THRESHOLD_HPP
0010 #define BOOST_GIL_IMAGE_PROCESSING_THRESHOLD_HPP
0011
0012 #include <limits>
0013 #include <array>
0014 #include <type_traits>
0015 #include <cstddef>
0016 #include <algorithm>
0017 #include <vector>
0018 #include <cmath>
0019
0020 #include <boost/assert.hpp>
0021
0022 #include <boost/gil/image.hpp>
0023 #include <boost/gil/image_processing/kernel.hpp>
0024 #include <boost/gil/image_processing/convolve.hpp>
0025 #include <boost/gil/image_processing/numeric.hpp>
0026
0027 namespace boost { namespace gil {
0028
0029 namespace detail {
0030
0031 template
0032 <
0033     typename SourceChannelT,
0034     typename ResultChannelT,
0035     typename SrcView,
0036     typename DstView,
0037     typename Operator
0038 >
0039 void threshold_impl(SrcView const& src_view, DstView const& dst_view, Operator const& threshold_op)
0040 {
0041     gil_function_requires<ImageViewConcept<SrcView>>();
0042     gil_function_requires<MutableImageViewConcept<DstView>>();
0043     static_assert(color_spaces_are_compatible
0044     <
0045         typename color_space_type<SrcView>::type,
0046         typename color_space_type<DstView>::type
0047     >::value, "Source and destination views must have pixels with the same color space");
0048
0049     //iterate over the image checking each pixel value for the threshold
0050     for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
0051     {
0052         typename SrcView::x_iterator src_it = src_view.row_begin(y);
0053         typename DstView::x_iterator dst_it = dst_view.row_begin(y);
0054
0055         for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
0056         {
0057             static_transform(src_it[x], dst_it[x], threshold_op);
0058         }
0059     }
0060 }
0061
0062 } //namespace boost::gil::detail
0063
0064 /// \addtogroup ImageProcessing
0065 /// @{
0066 ///
0067 /// \brief Direction of image segmentation.
0068 /// The direction specifies which pixels are considered as corresponding to object
0069 /// and which pixels correspond to background.
0070 enum class threshold_direction
0071 {
0072     regular, ///< Consider values greater than threshold value
0073     inverse  ///< Consider values less than or equal to threshold value
0074 };
0075
0076 /// \ingroup ImageProcessing
0077 /// \brief Method of optimal threshold value calculation.
0078 enum class threshold_optimal_value
0079 {
0080     otsu        ///< \todo TODO
0081 };
0082
0083 /// \ingroup ImageProcessing
0084 /// \brief TODO
0085 enum class threshold_truncate_mode
0086 {
0087     threshold,  ///< \todo TODO
0088     zero        ///< \todo TODO
0089 };
0090
0091 enum class threshold_adaptive_method
0092 {
0093     mean,
0094     gaussian
0095 };
0096
0097 /// \ingroup ImageProcessing
0098 /// \brief Applies fixed threshold to each pixel of image view.
0099 /// Performs image binarization by thresholding channel value of each
0100 /// pixel of given image view.
0101 /// \param src_view - TODO
0102 /// \param dst_view - TODO
0103 /// \param threshold_value - TODO
0104 /// \param max_value - TODO
0105 /// \param threshold_direction - if regular, values greater than threshold_value are
0106 /// set to max_value else set to 0; if inverse, values greater than threshold_value are
0107 /// set to 0 else set to max_value.
0108 template <typename SrcView, typename DstView>
0109 void threshold_binary(
0110     SrcView const& src_view,
0111     DstView const& dst_view,
0112     typename channel_type<DstView>::type threshold_value,
0113     typename channel_type<DstView>::type max_value,
0114     threshold_direction direction = threshold_direction::regular
0115 )
0116 {
0117     //deciding output channel type and creating functor
0118     using source_channel_t = typename channel_type<SrcView>::type;
0119     using result_channel_t = typename channel_type<DstView>::type;
0120
0121     if (direction == threshold_direction::regular)
0122     {
0123         detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
0124             [threshold_value, max_value](source_channel_t px) -> result_channel_t {
0125                 return px > threshold_value ? max_value : 0;
0126             });
0127     }
0128     else
0129     {
0130         detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
0131             [threshold_value, max_value](source_channel_t px) -> result_channel_t {
0132                 return px > threshold_value ? 0 : max_value;
0133             });
0134     }
0135 }
0136
0137 /// \ingroup ImageProcessing
0138 /// \brief Applies fixed threshold to each pixel of image view.
0139 /// Performs image binarization by thresholding channel value of each
0140 /// pixel of given image view.
0141 /// This variant of threshold_binary automatically deduces maximum value for each channel
0142 /// of pixel based on channel type.
0143 /// If direction is regular, values greater than threshold_value will be set to maximum
0144 /// numeric limit of channel else 0.
0145 /// If direction is inverse, values greater than threshold_value will be set to 0 else maximum
0146 /// numeric limit of channel.
0147 template <typename SrcView, typename DstView>
0148 void threshold_binary(
0149     SrcView const& src_view,
0150     DstView const& dst_view,
0151     typename channel_type<DstView>::type threshold_value,
0152     threshold_direction direction = threshold_direction::regular
0153 )
0154 {
0155     //deciding output channel type and creating functor
0156     using result_channel_t = typename channel_type<DstView>::type;
0157
0158     result_channel_t max_value = (std::numeric_limits<result_channel_t>::max)();
0159     threshold_binary(src_view, dst_view, threshold_value, max_value, direction);
0160 }
0161
0162 /// \ingroup ImageProcessing
0163 /// \brief Applies truncating threshold to each pixel of image view.
0164 /// Takes an image view and performs truncating threshold operation on each chennel.
0165 /// If mode is threshold and direction is regular:
0166 /// values greater than threshold_value will be set to threshold_value else no change
0167 /// If mode is threshold and direction is inverse:
0168 /// values less than or equal to threshold_value will be set to threshold_value else no change
0169 /// If mode is zero and direction is regular:
0170 /// values less than or equal to threshold_value will be set to 0 else no change
0171 /// If mode is zero and direction is inverse:
0172 /// values more than threshold_value will be set to 0 else no change
0173 template <typename SrcView, typename DstView>
0174 void threshold_truncate(
0175     SrcView const& src_view,
0176     DstView const& dst_view,
0177     typename channel_type<DstView>::type threshold_value,
0178     threshold_truncate_mode mode = threshold_truncate_mode::threshold,
0179     threshold_direction direction = threshold_direction::regular
0180 )
0181 {
0182     //deciding output channel type and creating functor
0183     using source_channel_t = typename channel_type<SrcView>::type;
0184     using result_channel_t = typename channel_type<DstView>::type;
0185
0186     std::function<result_channel_t(source_channel_t)> threshold_logic;
0187
0188     if (mode == threshold_truncate_mode::threshold)
0189     {
0190         if (direction == threshold_direction::regular)
0191         {
0192             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
0193                 [threshold_value](source_channel_t px) -> result_channel_t {
0194                     return px > threshold_value ? threshold_value : px;
0195                 });
0196         }
0197         else
0198         {
0199             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
0200                 [threshold_value](source_channel_t px) -> result_channel_t {
0201                     return px > threshold_value ? px : threshold_value;
0202                 });
0203         }
0204     }
0205     else
0206     {
0207         if (direction == threshold_direction::regular)
0208         {
0209             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
0210                 [threshold_value](source_channel_t px) -> result_channel_t {
0211                     return px > threshold_value ? px : 0;
0212                 });
0213         }
0214         else
0215         {
0216             detail::threshold_impl<source_channel_t, result_channel_t>(src_view, dst_view,
0217                 [threshold_value](source_channel_t px) -> result_channel_t {
0218                     return px > threshold_value ? 0 : px;
0219                 });
0220         }
0221     }
0222 }
0223
0224 namespace detail{
0225
0226 template <typename SrcView, typename DstView>
0227 void otsu_impl(SrcView const& src_view, DstView const& dst_view, threshold_direction direction)
0228 {
0229     //deciding output channel type and creating functor
0230     using source_channel_t = typename channel_type<SrcView>::type;
0231
0232     std::array<std::size_t, 256> histogram{};
0233     //initial value of min is set to maximum possible value to compare histogram data
0234     //initial value of max is set to minimum possible value to compare histogram data
0235     auto min = (std::numeric_limits<source_channel_t>::max)(),
0236         max = (std::numeric_limits<source_channel_t>::min)();
0237
0238     if (sizeof(source_channel_t) > 1 || std::is_signed<source_channel_t>::value)
0239     {
0240         //iterate over the image to find the min and max pixel values
0241         for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
0242         {
0243             typename SrcView::x_iterator src_it = src_view.row_begin(y);
0244             for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
0245             {
0246                 if (src_it[x] < min) min = src_it[x];
0247                 if (src_it[x] > min) min = src_it[x];
0248             }
0249         }
0250
0251         //making histogram
0252         for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
0253         {
0254             typename SrcView::x_iterator src_it = src_view.row_begin(y);
0255
0256             for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
0257             {
0258                 histogram[((src_it[x] - min) * 255) / (max - min)]++;
0259             }
0260         }
0261     }
0262     else
0263     {
0264         //making histogram
0265         for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
0266         {
0267             typename SrcView::x_iterator src_it = src_view.row_begin(y);
0268
0269             for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
0270             {
0271                 histogram[src_it[x]]++;
0272             }
0273         }
0274     }
0275
0276     //histData = histogram data
0277     //sum = total (background + foreground)
0278     //sumB = sum background
0279     //wB = weight background
0280     //wf = weight foreground
0281     //varMax = tracking the maximum known value of between class variance
0282     //mB = mu background
0283     //mF = mu foreground
0284     //varBeetween = between class variance
0285     //http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html
0286     //https://www.ipol.im/pub/art/2016/158/
0287     std::ptrdiff_t total_pixel = src_view.height() * src_view.width();
0288     std::ptrdiff_t sum_total = 0, sum_back = 0;
0289     std::size_t weight_back = 0, weight_fore = 0, threshold = 0;
0290     double var_max = 0, mean_back, mean_fore, var_intra_class;
0291
0292     for (std::size_t t = 0; t < 256; t++)
0293     {
0294         sum_total += t * histogram[t];
0295     }
0296
0297     for (int t = 0; t < 256; t++)
0298     {
0299         weight_back += histogram[t];               // Weight Background
0300         if (weight_back == 0) continue;
0301
0302         weight_fore = total_pixel - weight_back;          // Weight Foreground
0303         if (weight_fore == 0) break;
0304
0305         sum_back += t * histogram[t];
0306
0307         mean_back = sum_back / weight_back;            // Mean Background
0308         mean_fore = (sum_total - sum_back) / weight_fore;    // Mean Foreground
0309
0310         // Calculate Between Class Variance
0311         var_intra_class = weight_back * weight_fore * (mean_back - mean_fore) * (mean_back - mean_fore);
0312
0313         // Check if new maximum found
0314         if (var_intra_class > var_max) {
0315             var_max = var_intra_class;
0316             threshold = t;
0317         }
0318     }
0319     if (sizeof(source_channel_t) > 1 && std::is_unsigned<source_channel_t>::value)
0320     {
0321         threshold_binary(src_view, dst_view, (threshold * (max - min) / 255) + min, direction);
0322     }
0323     else {
0324         threshold_binary(src_view, dst_view, threshold, direction);
0325     }
0326 }
0327 } //namespace detail
0328
0329 template <typename SrcView, typename DstView>
0330 void threshold_optimal
0331 (
0332     SrcView const& src_view,
0333     DstView const& dst_view,
0334     threshold_optimal_value mode = threshold_optimal_value::otsu,
0335     threshold_direction direction = threshold_direction::regular
0336 )
0337 {
0338     if (mode == threshold_optimal_value::otsu)
0339     {
0340         for (std::size_t i = 0; i < src_view.num_channels(); i++)
0341         {
0342             detail::otsu_impl
0343                 (nth_channel_view(src_view, i), nth_channel_view(dst_view, i), direction);
0344         }
0345     }
0346 }
0347
0348 namespace detail {
0349
0350 template
0351 <
0352     typename SourceChannelT,
0353     typename ResultChannelT,
0354     typename SrcView,
0355     typename DstView,
0356     typename Operator
0357 >
0358 void adaptive_impl
0359 (
0360     SrcView const& src_view,
0361     SrcView const& convolved_view,
0362     DstView const& dst_view,
0363     Operator const& threshold_op
0364 )
0365 {
0366     //template argument validation
0367     gil_function_requires<ImageViewConcept<SrcView>>();
0368     gil_function_requires<MutableImageViewConcept<DstView>>();
0369
0370     static_assert(color_spaces_are_compatible
0371     <
0372         typename color_space_type<SrcView>::type,
0373         typename color_space_type<DstView>::type
0374     >::value, "Source and destination views must have pixels with the same color space");
0375
0376     //iterate over the image checking each pixel value for the threshold
0377     for (std::ptrdiff_t y = 0; y < src_view.height(); y++)
0378     {
0379         typename SrcView::x_iterator src_it = src_view.row_begin(y);
0380         typename SrcView::x_iterator convolved_it = convolved_view.row_begin(y);
0381         typename DstView::x_iterator dst_it = dst_view.row_begin(y);
0382
0383         for (std::ptrdiff_t x = 0; x < src_view.width(); x++)
0384         {
0385             static_transform(src_it[x], convolved_it[x], dst_it[x], threshold_op);
0386         }
0387     }
0388 }
0389 } //namespace boost::gil::detail
0390
0391 template <typename SrcView, typename DstView>
0392 void threshold_adaptive
0393 (
0394     SrcView const& src_view,
0395     DstView const& dst_view,
0396     typename channel_type<DstView>::type max_value,
0397     std::size_t kernel_size,
0398     threshold_adaptive_method method = threshold_adaptive_method::mean,
0399     threshold_direction direction = threshold_direction::regular,
0400     typename channel_type<DstView>::type constant = 0
0401 )
0402 {
0403     BOOST_ASSERT_MSG((kernel_size % 2 != 0), "Kernel size must be an odd number");
0404
0405     typedef typename channel_type<SrcView>::type source_channel_t;
0406     typedef typename channel_type<DstView>::type result_channel_t;
0407
0408     image<typename SrcView::value_type> temp_img(src_view.width(), src_view.height());
0409     typename image<typename SrcView::value_type>::view_t temp_view = view(temp_img);
0410     SrcView temp_conv(temp_view);
0411
0412     if (method == threshold_adaptive_method::mean)
0413     {
0414         std::vector<float> mean_kernel_values(kernel_size, 1.0f/kernel_size);
0415         kernel_1d<float> kernel(mean_kernel_values.begin(), kernel_size, kernel_size/2);
0416
0417         detail::convolve_1d
0418         <
0419             pixel<float, typename SrcView::value_type::layout_t>
0420         >(src_view, kernel, temp_view);
0421     }
0422     else if (method == threshold_adaptive_method::gaussian)
0423     {
0424         detail::kernel_2d<float> kernel = generate_gaussian_kernel(kernel_size, 1.0);
0425         convolve_2d(src_view, kernel, temp_view);
0426     }
0427
0428     if (direction == threshold_direction::regular)
0429     {
0430         detail::adaptive_impl<source_channel_t, result_channel_t>(src_view, temp_conv, dst_view,
0431             [max_value, constant](source_channel_t px, source_channel_t threshold) -> result_channel_t
0432         { return px > (threshold - constant) ? max_value : 0; });
0433     }
0434     else
0435     {
0436         detail::adaptive_impl<source_channel_t, result_channel_t>(src_view, temp_conv, dst_view,
0437             [max_value, constant](source_channel_t px, source_channel_t threshold) -> result_channel_t
0438         { return px > (threshold - constant) ? 0 : max_value; });
0439     }
0440 }
0441
0442 template <typename SrcView, typename DstView>
0443 void threshold_adaptive
0444 (
0445     SrcView const& src_view,
0446     DstView const& dst_view,
0447     std::size_t kernel_size,
0448     threshold_adaptive_method method = threshold_adaptive_method::mean,
0449     threshold_direction direction = threshold_direction::regular,
0450     int constant = 0
0451 )
0452 {
0453     //deciding output channel type and creating functor
0454     typedef typename channel_type<DstView>::type result_channel_t;
0455
0456     result_channel_t max_value = (std::numeric_limits<result_channel_t>::max)();
0457
0458     threshold_adaptive(src_view, dst_view, max_value, kernel_size, method, direction, constant);
0459 }
0460
0461 /// @}
0462
0463 }} //namespace boost::gil
0464
0465 #endif //BOOST_GIL_IMAGE_PROCESSING_THRESHOLD_HPP