• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 >
threshold_impl(SrcView const & src_view,DstView const & dst_view,Operator const & threshold_op)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 checking 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 specifies 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>
threshold_binary(SrcView const & src_view,DstView const & dst_view,typename channel_type<DstView>::type threshold_value,typename channel_type<DstView>::type max_value,threshold_direction direction=threshold_direction::regular)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>
threshold_binary(SrcView const & src_view,DstView const & dst_view,typename channel_type<DstView>::type threshold_value,threshold_direction direction=threshold_direction::regular)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 performs 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 or equal to threshold_value will be set to threshold_value else no change
168 /// If mode is zero and direction is regular:
169 /// values less than or equal to 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>
threshold_truncate(SrcView const & src_view,DstView const & dst_view,typename channel_type<DstView>::type threshold_value,threshold_truncate_mode mode=threshold_truncate_mode::threshold,threshold_direction direction=threshold_direction::regular)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>
otsu_impl(SrcView const & src_view,DstView const & dst_view,threshold_direction direction)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>
threshold_optimal(SrcView const & src_view,DstView const & dst_view,threshold_optimal_value mode=threshold_optimal_value::otsu,threshold_direction direction=threshold_direction::regular)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 >
adaptive_impl(SrcView const & src_view,SrcView const & convolved_view,DstView const & dst_view,Operator const & threshold_op)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 checking 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>
threshold_adaptive(SrcView const & src_view,DstView const & dst_view,typename channel_type<DstView>::type max_value,std::size_t kernel_size,threshold_adaptive_method method=threshold_adaptive_method::mean,threshold_direction direction=threshold_direction::regular,typename channel_type<DstView>::type constant=0)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>
threshold_adaptive(SrcView const & src_view,DstView const & dst_view,std::size_t kernel_size,threshold_adaptive_method method=threshold_adaptive_method::mean,threshold_direction direction=threshold_direction::regular,int constant=0)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
465