Boost GIL


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 
33 inline double normalized_sinc(double x)
34 {
35  return std::sin(x * boost::gil::detail::pi) / (x * boost::gil::detail::pi);
36 }
37 
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 #if BOOST_WORKAROUND(BOOST_MSVC, >= 1400)
58 #pragma warning(push)
59 #pragma warning(disable:4244) // 'argument': conversion from 'const Channel' to 'BaseChannelValue', possible loss of data
60 #endif
61 
62 inline void compute_tensor_entries(
63  boost::gil::gray16s_view_t dx,
64  boost::gil::gray16s_view_t dy,
65  boost::gil::gray32f_view_t m11,
66  boost::gil::gray32f_view_t m12_21,
67  boost::gil::gray32f_view_t m22)
68 {
69  for (std::ptrdiff_t y = 0; y < dx.height(); ++y) {
70  for (std::ptrdiff_t x = 0; x < dx.width(); ++x) {
71  auto dx_value = dx(x, y);
72  auto dy_value = dy(x, y);
73  m11(x, y) = dx_value * dx_value;
74  m12_21(x, y) = dx_value * dy_value;
75  m22(x, y) = dy_value * dy_value;
76  }
77  }
78 }
79 
80 #if BOOST_WORKAROUND(BOOST_MSVC, >= 1400)
81 #pragma warning(pop)
82 #endif
83 
90 template <typename T = float, typename Allocator = std::allocator<T>>
91 inline detail::kernel_2d<T, Allocator> generate_normalized_mean(std::size_t side_length)
92 {
93  if (side_length % 2 != 1)
94  throw std::invalid_argument("kernel dimensions should be odd and equal");
95  const float entry = 1.0f / static_cast<float>(side_length * side_length);
96 
97  detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
98  for (auto& cell: result) {
99  cell = entry;
100  }
101 
102  return result;
103 }
104 
109 template <typename T = float, typename Allocator = std::allocator<T>>
110 inline detail::kernel_2d<T, Allocator> generate_unnormalized_mean(std::size_t side_length)
111 {
112  if (side_length % 2 != 1)
113  throw std::invalid_argument("kernel dimensions should be odd and equal");
114 
115  detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
116  for (auto& cell: result) {
117  cell = 1.0f;
118  }
119 
120  return result;
121 }
122 
128 template <typename T = float, typename Allocator = std::allocator<T>>
129 inline detail::kernel_2d<T, Allocator> generate_gaussian_kernel(std::size_t side_length, double sigma)
130 {
131  if (side_length % 2 != 1)
132  throw std::invalid_argument("kernel dimensions should be odd and equal");
133 
134 
135  const double denominator = 2 * boost::gil::detail::pi * sigma * sigma;
136  auto middle = side_length / 2;
137  std::vector<T, Allocator> values(side_length * side_length);
138  for (std::size_t y = 0; y < side_length; ++y)
139  {
140  for (std::size_t x = 0; x < side_length; ++x)
141  {
142  const auto delta_x = middle > x ? middle - x : x - middle;
143  const auto delta_y = middle > y ? middle - y : y - middle;
144  const double power = (delta_x * delta_x + delta_y * delta_y) / (2 * sigma * sigma);
145  const double nominator = std::exp(-power);
146  const float value = static_cast<float>(nominator / denominator);
147  values[y * side_length + x] = value;
148  }
149  }
150 
151  return detail::kernel_2d<T, Allocator>(values.begin(), values.size(), middle, middle);
152 }
153 
161 template <typename T = float, typename Allocator = std::allocator<T>>
162 inline detail::kernel_2d<T, Allocator> generate_dx_sobel(unsigned int degree = 1)
163 {
164  switch (degree)
165  {
166  case 0:
167  {
168  return detail::get_identity_kernel<T, Allocator>();
169  }
170  case 1:
171  {
172  detail::kernel_2d<T, Allocator> result(3, 1, 1);
173  std::copy(detail::dx_sobel.begin(), detail::dx_sobel.end(), result.begin());
174  return result;
175  }
176  default:
177  throw std::logic_error("not supported yet");
178  }
179 
180  //to not upset compiler
181  throw std::runtime_error("unreachable statement");
182 }
183 
191 template <typename T = float, typename Allocator = std::allocator<T>>
192 inline detail::kernel_2d<T, Allocator> generate_dx_scharr(unsigned int degree = 1)
193 {
194  switch (degree)
195  {
196  case 0:
197  {
198  return detail::get_identity_kernel<T, Allocator>();
199  }
200  case 1:
201  {
202  detail::kernel_2d<T, Allocator> result(3, 1, 1);
203  std::copy(detail::dx_scharr.begin(), detail::dx_scharr.end(), result.begin());
204  return result;
205  }
206  default:
207  throw std::logic_error("not supported yet");
208  }
209 
210  //to not upset compiler
211  throw std::runtime_error("unreachable statement");
212 }
213 
221 template <typename T = float, typename Allocator = std::allocator<T>>
222 inline detail::kernel_2d<T, Allocator> generate_dy_sobel(unsigned int degree = 1)
223 {
224  switch (degree)
225  {
226  case 0:
227  {
228  return detail::get_identity_kernel<T, Allocator>();
229  }
230  case 1:
231  {
232  detail::kernel_2d<T, Allocator> result(3, 1, 1);
233  std::copy(detail::dy_sobel.begin(), detail::dy_sobel.end(), result.begin());
234  return result;
235  }
236  default:
237  throw std::logic_error("not supported yet");
238  }
239 
240  //to not upset compiler
241  throw std::runtime_error("unreachable statement");
242 }
243 
251 template <typename T = float, typename Allocator = std::allocator<T>>
252 inline detail::kernel_2d<T, Allocator> generate_dy_scharr(unsigned int degree = 1)
253 {
254  switch (degree)
255  {
256  case 0:
257  {
258  return detail::get_identity_kernel<T, Allocator>();
259  }
260  case 1:
261  {
262  detail::kernel_2d<T, Allocator> result(3, 1, 1);
263  std::copy(detail::dy_scharr.begin(), detail::dy_scharr.end(), result.begin());
264  return result;
265  }
266  default:
267  throw std::logic_error("not supported yet");
268  }
269 
270  //to not upset compiler
271  throw std::runtime_error("unreachable statement");
272 }
273 
283 template <typename GradientView, typename OutputView>
285  GradientView dx,
286  GradientView dy,
287  OutputView ddxx,
288  OutputView dxdy,
289  OutputView ddyy)
290 {
291  auto sobel_x = generate_dx_sobel();
292  auto sobel_y = generate_dy_sobel();
293  detail::convolve_2d(dx, sobel_x, ddxx);
294  detail::convolve_2d(dx, sobel_y, dxdy);
295  detail::convolve_2d(dy, sobel_y, ddyy);
296 }
297 
298 }} // namespace boost::gil
299 
300 #endif
detail::kernel_2d< T, Allocator > generate_dy_sobel(unsigned int degree=1)
Generates Sobel operator in vertical directionGenerates a kernel which will represent Sobel operator ...
Definition: numeric.hpp:222
void compute_hessian_entries(GradientView dx, GradientView dy, OutputView ddxx, OutputView dxdy, OutputView ddyy)
Compute xy gradient, and second order x and y gradientsHessian matrix is defined as a matrix of parti...
Definition: numeric.hpp:284
double lanczos(double x, std::ptrdiff_t a)
Lanczos response at point xLanczos response is defined as: x == 0: 1 -a < x && x < a: 0 otherwise: no...
Definition: numeric.hpp:45
BOOST_FORCEINLINE auto copy(boost::gil::pixel< T, CS > *first, boost::gil::pixel< T, CS > *last, boost::gil::pixel< T, CS > *dst) -> boost::gil::pixel< T, CS > *
Copy when both src and dst are interleaved and of the same type can be just memmove.
Definition: algorithm.hpp:139
detail::kernel_2d< T, Allocator > generate_dx_scharr(unsigned int degree=1)
Generate Scharr operator in horizontal directionGenerates a kernel which will represent Scharr operat...
Definition: numeric.hpp:192
detail::kernel_2d< T, Allocator > generate_dx_sobel(unsigned int degree=1)
Generates Sobel operator in horizontal directionGenerates a kernel which will represent Sobel operato...
Definition: numeric.hpp:162
detail::kernel_2d< T, Allocator > generate_unnormalized_mean(std::size_t side_length)
Generate kernel with all 1sFills supplied view with 1s (ones)
Definition: numeric.hpp:110
detail::kernel_2d< T, Allocator > generate_dy_scharr(unsigned int degree=1)
Generate Scharr operator in vertical directionGenerates a kernel which will represent Scharr operator...
Definition: numeric.hpp:252
detail::kernel_2d< T, Allocator > generate_gaussian_kernel(std::size_t side_length, double sigma)
Generate Gaussian kernelFills supplied view with values taken from Gaussian distribution....
Definition: numeric.hpp:129
detail::kernel_2d< T, Allocator > generate_normalized_mean(std::size_t side_length)
Generate mean kernelFills supplied view with normalized mean in which all entries will be equal to.
Definition: numeric.hpp:91