threshold.hpp 16 KB

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