hessian.hpp 2.8 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. #ifndef BOOST_GIL_IMAGE_PROCESSING_HESSIAN_HPP
  2. #define BOOST_GIL_IMAGE_PROCESSING_HESSIAN_HPP
  3. #include <boost/gil/image_view.hpp>
  4. #include <boost/gil/typedefs.hpp>
  5. #include <boost/gil/extension/numeric/kernel.hpp>
  6. #include <stdexcept>
  7. namespace boost { namespace gil {
  8. /// \brief Computes Hessian response
  9. ///
  10. /// Computes Hessian response based on computed entries of Hessian matrix, e.g. second order
  11. /// derivates in x and y, and derivatives in both x, y.
  12. /// d stands for derivative, and x or y stand for derivative direction. For example,
  13. /// ddxx means taking two derivatives (gradients) in horizontal direction.
  14. /// Weights change perception of surroinding pixels.
  15. /// Additional filtering is strongly advised.
  16. template <typename GradientView, typename T, typename Allocator, typename OutputView>
  17. inline void compute_hessian_responses(
  18. GradientView ddxx,
  19. GradientView dxdy,
  20. GradientView ddyy,
  21. const detail::kernel_2d<T, Allocator>& weights,
  22. OutputView dst)
  23. {
  24. if (ddxx.dimensions() != ddyy.dimensions()
  25. || ddyy.dimensions() != dxdy.dimensions()
  26. || dxdy.dimensions() != dst.dimensions()
  27. || weights.center_x() != weights.center_y())
  28. {
  29. throw std::invalid_argument("dimensions of views are not the same"
  30. " or weights don't have equal width and height"
  31. " or weights' dimensions are not odd");
  32. }
  33. // Use pixel type of output, as values will be written to output
  34. using pixel_t = typename std::remove_reference<decltype(std::declval<OutputView>()(0, 0))>::type;
  35. using channel_t = typename std::remove_reference
  36. <
  37. decltype(std::declval<pixel_t>().at(std::integral_constant<int, 0>{}))
  38. >::type;
  39. auto center = weights.center_y();
  40. for (auto y = center; y < dst.height() - center; ++y)
  41. {
  42. for (auto x = center; x < dst.width() - center; ++x)
  43. {
  44. auto ddxx_i = channel_t();
  45. auto ddyy_i = channel_t();
  46. auto dxdy_i = channel_t();
  47. for (typename OutputView::coord_t w_y = 0; w_y < weights.size(); ++w_y)
  48. {
  49. for (typename OutputView::coord_t w_x = 0; w_x < weights.size(); ++w_x)
  50. {
  51. ddxx_i += ddxx(x + w_x - center, y + w_y - center)
  52. .at(std::integral_constant<int, 0>{}) * weights.at(w_x, w_y);
  53. ddyy_i += ddyy(x + w_x - center, y + w_y - center)
  54. .at(std::integral_constant<int, 0>{}) * weights.at(w_x, w_y);
  55. dxdy_i += dxdy(x + w_x - center, y + w_y - center)
  56. .at(std::integral_constant<int, 0>{}) * weights.at(w_x, w_y);
  57. }
  58. }
  59. auto determinant = ddxx_i * ddyy_i - dxdy_i * dxdy_i;
  60. dst(x, y).at(std::integral_constant<int, 0>{}) = determinant;
  61. }
  62. }
  63. }
  64. }} // namespace boost::gil
  65. #endif