boost/gil/image_processing/numeric.hpp

   1 //
   2 // Copyright 2019 Olzhas Zhumabek <anonymous.from.applecity@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_NUMERIC_HPP
   9 #define BOOST_GIL_IMAGE_PROCESSING_NUMERIC_HPP
  10
  11 #include <boost/gil/extension/numeric/kernel.hpp>
  12 #include <boost/gil/extension/numeric/convolve.hpp>
  13 #include <boost/gil/image_view.hpp>
  14 #include <boost/gil/typedefs.hpp>
  15 #include <boost/gil/detail/math.hpp>
  16 // fixes ambigious call to std::abs, https://stackoverflow.com/a/30084734/4593721
  17 #include <cstdlib>
  18 #include <cmath>
  19
  20 namespace boost { namespace gil {
  21
  22 /// \defgroup ImageProcessingMath
  23 /// \brief Math operations for IP algorithms
  24 ///
  25 /// This is mostly handful of mathemtical operations that are required by other
  26 /// image processing algorithms
  27 ///
  28 /// \brief Normalized cardinal sine
  29 /// \ingroup ImageProcessingMath
  30 ///
  31 /// normalized_sinc(x) = sin(pi * x) / (pi * x)
  32 ///
  33 inline double normalized_sinc(double x)
  34 {
  35     return std::sin(x * boost::gil::pi) / (x * boost::gil::pi);
  36 }
  37
  38 /// \brief Lanczos response at point x
  39 /// \ingroup ImageProcessingMath
  40 ///
  41 /// Lanczos response is defined as:
  42 /// x == 0: 1
  43 /// -a < x && x < a: 0
  44 /// otherwise: normalized_sinc(x) / normalized_sinc(x / a)
  45 inline double lanczos(double x, std::ptrdiff_t a)
  46 {
  47     // means == but <= avoids compiler warning
  48     if (0 <= x && x <= 0)
  49         return 1;
  50
  51     if (-a < x && x < a)
  52         return normalized_sinc(x) / normalized_sinc(x / static_cast<double>(a));
  53
  54     return 0;
  55 }
  56
  57 inline void compute_tensor_entries(
  58     boost::gil::gray16s_view_t dx,
  59     boost::gil::gray16s_view_t dy,
  60     boost::gil::gray32f_view_t m11,
  61     boost::gil::gray32f_view_t m12_21,
  62     boost::gil::gray32f_view_t m22)
  63 {
  64     for (std::ptrdiff_t y = 0; y < dx.height(); ++y) {
  65         for (std::ptrdiff_t x = 0; x < dx.width(); ++x) {
  66             auto dx_value = dx(x, y);
  67             auto dy_value = dy(x, y);
  68             m11(x, y) = dx_value * dx_value;
  69             m12_21(x, y) = dx_value * dy_value;
  70             m22(x, y) = dy_value * dy_value;
  71         }
  72     }
  73 }
  74
  75 /// \brief Generate mean kernel
  76 /// \ingroup ImageProcessingMath
  77 ///
  78 /// Fills supplied view with normalized mean
  79 /// in which all entries will be equal to
  80 /// \code 1 / (dst.size()) \endcode
  81 template <typename T = float, typename Allocator = std::allocator<T>>
  82 inline detail::kernel_2d<T, Allocator> generate_normalized_mean(std::size_t side_length)
  83 {
  84     if (side_length % 2 != 1)
  85         throw std::invalid_argument("kernel dimensions should be odd and equal");
  86     const float entry = 1.0f / static_cast<float>(side_length * side_length);
  87
  88     detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
  89     for (auto& cell: result) {
  90         cell = entry;
  91     }
  92
  93     return result;
  94 }
  95
  96 /// \brief Generate kernel with all 1s
  97 /// \ingroup ImageProcessingMath
  98 ///
  99 /// Fills supplied view with 1s (ones)
 100 template <typename T = float, typename Allocator = std::allocator<T>>
 101 inline detail::kernel_2d<T, Allocator> generate_unnormalized_mean(std::size_t side_length)
 102 {
 103     if (side_length % 2 != 1)
 104         throw std::invalid_argument("kernel dimensions should be odd and equal");
 105
 106     detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
 107     for (auto& cell: result) {
 108         cell = 1.0f;
 109     }
 110
 111     return result;
 112 }
 113
 114 /// \brief Generate Gaussian kernel
 115 /// \ingroup ImageProcessingMath
 116 ///
 117 /// Fills supplied view with values taken from Gaussian distribution. See
 118 /// https://en.wikipedia.org/wiki/Gaussian_blur
 119 template <typename T = float, typename Allocator = std::allocator<T>>
 120 inline detail::kernel_2d<T, Allocator> generate_gaussian_kernel(std::size_t side_length, double sigma)
 121 {
 122     if (side_length % 2 != 1)
 123         throw std::invalid_argument("kernel dimensions should be odd and equal");
 124
 125
 126     const double denominator = 2 * boost::gil::pi * sigma * sigma;
 127     auto middle = side_length / 2;
 128     std::vector<T, Allocator> values(side_length * side_length);
 129     for (std::size_t y = 0; y < side_length; ++y)
 130     {
 131         for (std::size_t x = 0; x < side_length; ++x)
 132         {
 133             const auto delta_x = middle > x ? middle - x : x - middle;
 134             const auto delta_y = middle > y ? middle - y : y - middle;
 135             const double power = (delta_x * delta_x +  delta_y * delta_y) / (2 * sigma * sigma);
 136             const double nominator = std::exp(-power);
 137             const float value = static_cast<float>(nominator / denominator);
 138             values[y * side_length + x] = value;
 139         }
 140     }
 141
 142     return detail::kernel_2d<T, Allocator>(values.begin(), values.size(), middle, middle);
 143 }
 144
 145 /// \brief Generates Sobel operator in horizontal direction
 146 /// \ingroup ImageProcessingMath
 147 ///
 148 /// Generates a kernel which will represent Sobel operator in
 149 /// horizontal direction of specified degree (no need to convolve multiple times
 150 /// to obtain the desired degree).
 151 /// https://www.researchgate.net/publication/239398674_An_Isotropic_3_3_Image_Gradient_Operator
 152 template <typename T = float, typename Allocator = std::allocator<T>>
 153 inline detail::kernel_2d<T, Allocator> generate_dx_sobel(unsigned int degree = 1)
 154 {
 155     switch (degree)
 156     {
 157         case 0:
 158         {
 159             return get_identity_kernel<T, Allocator>();
 160         }
 161         case 1:
 162         {
 163             detail::kernel_2d<T, Allocator> result(3, 1, 1);
 164             std::copy(dx_sobel.begin(), dx_sobel.end(), result.begin());
 165             return result;
 166         }
 167         default:
 168             throw std::logic_error("not supported yet");
 169     }
 170
 171     //to not upset compiler
 172     throw std::runtime_error("unreachable statement");
 173 }
 174
 175 /// \brief Generate Scharr operator in horizontal direction
 176 /// \ingroup ImageProcessingMath
 177 ///
 178 /// Generates a kernel which will represent Scharr operator in
 179 /// horizontal direction of specified degree (no need to convolve multiple times
 180 /// to obtain the desired degree).
 181 /// https://www.researchgate.net/profile/Hanno_Scharr/publication/220955743_Optimal_Filters_for_Extended_Optical_Flow/links/004635151972eda98f000000/Optimal-Filters-for-Extended-Optical-Flow.pdf
 182 template <typename T = float, typename Allocator = std::allocator<T>>
 183 inline detail::kernel_2d<T, Allocator> generate_dx_scharr(unsigned int degree = 1)
 184 {
 185     switch (degree)
 186     {
 187         case 0:
 188         {
 189             return get_identity_kernel<T, Allocator>();
 190         }
 191         case 1:
 192         {
 193             detail::kernel_2d<T, Allocator> result(3, 1, 1);
 194             std::copy(dx_scharr.begin(), dx_scharr.end(), result.begin());
 195             return result;
 196         }
 197         default:
 198             throw std::logic_error("not supported yet");
 199     }
 200
 201     //to not upset compiler
 202     throw std::runtime_error("unreachable statement");
 203 }
 204
 205 /// \brief Generates Sobel operator in vertical direction
 206 /// \ingroup ImageProcessingMath
 207 ///
 208 /// Generates a kernel which will represent Sobel operator in
 209 /// vertical direction of specified degree (no need to convolve multiple times
 210 /// to obtain the desired degree).
 211 /// https://www.researchgate.net/publication/239398674_An_Isotropic_3_3_Image_Gradient_Operator
 212 template <typename T = float, typename Allocator = std::allocator<T>>
 213 inline detail::kernel_2d<T, Allocator> generate_dy_sobel(unsigned int degree = 1)
 214 {
 215     switch (degree)
 216     {
 217         case 0:
 218         {
 219             return get_identity_kernel<T, Allocator>();
 220         }
 221         case 1:
 222         {
 223             detail::kernel_2d<T, Allocator> result(3, 1, 1);
 224             std::copy(dy_sobel.begin(), dy_sobel.end(), result.begin());
 225             return result;
 226         }
 227         default:
 228             throw std::logic_error("not supported yet");
 229     }
 230
 231     //to not upset compiler
 232     throw std::runtime_error("unreachable statement");
 233 }
 234
 235 /// \brief Generate Scharr operator in vertical direction
 236 /// \ingroup ImageProcessingMath
 237 ///
 238 /// Generates a kernel which will represent Scharr operator in
 239 /// vertical direction of specified degree (no need to convolve multiple times
 240 /// to obtain the desired degree).
 241 /// https://www.researchgate.net/profile/Hanno_Scharr/publication/220955743_Optimal_Filters_for_Extended_Optical_Flow/links/004635151972eda98f000000/Optimal-Filters-for-Extended-Optical-Flow.pdf
 242 template <typename T = float, typename Allocator = std::allocator<T>>
 243 inline detail::kernel_2d<T, Allocator> generate_dy_scharr(unsigned int degree = 1)
 244 {
 245     switch (degree)
 246     {
 247         case 0:
 248         {
 249             return get_identity_kernel<T, Allocator>();
 250         }
 251         case 1:
 252         {
 253             detail::kernel_2d<T, Allocator> result(3, 1, 1);
 254             std::copy(dy_scharr.begin(), dy_scharr.end(), result.begin());
 255             return result;
 256         }
 257         default:
 258             throw std::logic_error("not supported yet");
 259     }
 260
 261     //to not upset compiler
 262     throw std::runtime_error("unreachable statement");
 263 }
 264
 265 /// \brief Compute xy gradient, and second order x and y gradients
 266 /// \ingroup ImageProcessingMath
 267 ///
 268 /// Hessian matrix is defined as a matrix of partial derivates
 269 /// for 2d case, it is [[ddxx, dxdy], [dxdy, ddyy].
 270 /// d stands for derivative, and x or y stand for direction.
 271 /// For example, dx stands for derivative (gradient) in horizontal
 272 /// direction, and ddxx means second order derivative in horizon direction
 273 /// https://en.wikipedia.org/wiki/Hessian_matrix
 274 template <typename GradientView, typename OutputView>
 275 inline void compute_hessian_entries(
 276     GradientView dx,
 277     GradientView dy,
 278     OutputView ddxx,
 279     OutputView dxdy,
 280     OutputView ddyy)
 281 {
 282     auto sobel_x = generate_dx_sobel();
 283     auto sobel_y = generate_dy_sobel();
 284     detail::convolve_2d(dx, sobel_x, ddxx);
 285     detail::convolve_2d(dx, sobel_y, dxdy);
 286     detail::convolve_2d(dy, sobel_y, ddyy);
 287 }
 288
 289 }} // namespace boost::gil
 290
 291 #endif