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