boost/gil/image_processing/hessian.hpp

   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>
  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