Boost GIL


numeric.hpp
1 //
2 // Copyright 2019 Olzhas Zhumabek <anonymous.from.applecity@gmail.com>
3 // Copyright 2021 Pranam Lashkari <plashkari628@gmail.com>
4 //
5 // Use, modification and distribution are subject to the Boost Software License,
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
8 //
9 #ifndef BOOST_GIL_IMAGE_PROCESSING_NUMERIC_HPP
10 #define BOOST_GIL_IMAGE_PROCESSING_NUMERIC_HPP
11 
12 #include <boost/gil/image_processing/kernel.hpp>
13 #include <boost/gil/image_processing/convolve.hpp>
14 #include <boost/gil/image_view.hpp>
15 #include <boost/gil/typedefs.hpp>
16 #include <boost/gil/detail/math.hpp>
17 // fixes ambigious call to std::abs, https://stackoverflow.com/a/30084734/4593721
18 #include <cstdlib>
19 #include <cmath>
20 
21 namespace boost { namespace gil {
22 
34 inline double normalized_sinc(double x)
35 {
36  return std::sin(x * boost::gil::detail::pi) / (x * boost::gil::detail::pi);
37 }
38 
46 inline double lanczos(double x, std::ptrdiff_t a)
47 {
48  // means == but <= avoids compiler warning
49  if (0 <= x && x <= 0)
50  return 1;
51 
52  if (static_cast<double>(-a) < x && x < static_cast<double>(a))
53  return normalized_sinc(x) / normalized_sinc(x / static_cast<double>(a));
54 
55  return 0;
56 }
57 
58 #if BOOST_WORKAROUND(BOOST_MSVC, >= 1400)
59 #pragma warning(push)
60 #pragma warning(disable:4244) // 'argument': conversion from 'const Channel' to 'BaseChannelValue', possible loss of data
61 #endif
62 
63 inline void compute_tensor_entries(
64  boost::gil::gray16s_view_t dx,
65  boost::gil::gray16s_view_t dy,
66  boost::gil::gray32f_view_t m11,
67  boost::gil::gray32f_view_t m12_21,
68  boost::gil::gray32f_view_t m22)
69 {
70  for (std::ptrdiff_t y = 0; y < dx.height(); ++y) {
71  for (std::ptrdiff_t x = 0; x < dx.width(); ++x) {
72  auto dx_value = dx(x, y);
73  auto dy_value = dy(x, y);
74  m11(x, y) = dx_value * dx_value;
75  m12_21(x, y) = dx_value * dy_value;
76  m22(x, y) = dy_value * dy_value;
77  }
78  }
79 }
80 
81 #if BOOST_WORKAROUND(BOOST_MSVC, >= 1400)
82 #pragma warning(pop)
83 #endif
84 
91 template <typename T = float, typename Allocator = std::allocator<T>>
92 inline auto generate_normalized_mean(std::size_t side_length)
94 {
95  if (side_length % 2 != 1)
96  throw std::invalid_argument("kernel dimensions should be odd and equal");
97  const float entry = 1.0f / static_cast<float>(side_length * side_length);
98 
99  detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
100  for (auto& cell: result) {
101  cell = entry;
102  }
103 
104  return result;
105 }
106 
111 template <typename T = float, typename Allocator = std::allocator<T>>
112 inline auto generate_unnormalized_mean(std::size_t side_length)
114 {
115  if (side_length % 2 != 1)
116  throw std::invalid_argument("kernel dimensions should be odd and equal");
117 
118  detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
119  for (auto& cell: result) {
120  cell = 1.0f;
121  }
122 
123  return result;
124 }
125 
131 template <typename T = float, typename Allocator = std::allocator<T>>
132 inline auto generate_gaussian_kernel(std::size_t side_length, double sigma)
134 {
135  if (side_length % 2 != 1)
136  throw std::invalid_argument("kernel dimensions should be odd and equal");
137 
138  const double denominator = 2 * boost::gil::detail::pi * sigma * sigma;
139  auto middle = side_length / 2;
140  std::vector<T, Allocator> values(side_length * side_length);
141  for (std::size_t y = 0; y < side_length; ++y)
142  {
143  for (std::size_t x = 0; x < side_length; ++x)
144  {
145  const auto delta_x = middle > x ? middle - x : x - middle;
146  const auto delta_y = middle > y ? middle - y : y - middle;
147  const double power = (delta_x * delta_x + delta_y * delta_y) / (2 * sigma * sigma);
148  const double nominator = std::exp(-power);
149  const float value = static_cast<float>(nominator / denominator);
150  values[y * side_length + x] = value;
151  }
152  }
153 
154  return detail::kernel_2d<T, Allocator>(values.begin(), values.size(), middle, middle);
155 }
156 
164 template <typename T = float, typename Allocator = std::allocator<T>>
165 inline auto generate_dx_sobel(unsigned int degree = 1)
167 {
168  switch (degree)
169  {
170  case 0:
171  {
172  return detail::get_identity_kernel<T, Allocator>();
173  }
174  case 1:
175  {
176  detail::kernel_2d<T, Allocator> result(3, 1, 1);
177  std::copy(detail::dx_sobel.begin(), detail::dx_sobel.end(), result.begin());
178  return result;
179  }
180  default:
181  throw std::logic_error("not supported yet");
182  }
183 
184  //to not upset compiler
185  throw std::runtime_error("unreachable statement");
186 }
187 
195 template <typename T = float, typename Allocator = std::allocator<T>>
196 inline auto generate_dx_scharr(unsigned int degree = 1)
198 {
199  switch (degree)
200  {
201  case 0:
202  {
203  return detail::get_identity_kernel<T, Allocator>();
204  }
205  case 1:
206  {
207  detail::kernel_2d<T, Allocator> result(3, 1, 1);
208  std::copy(detail::dx_scharr.begin(), detail::dx_scharr.end(), result.begin());
209  return result;
210  }
211  default:
212  throw std::logic_error("not supported yet");
213  }
214 
215  //to not upset compiler
216  throw std::runtime_error("unreachable statement");
217 }
218 
226 template <typename T = float, typename Allocator = std::allocator<T>>
227 inline auto generate_dy_sobel(unsigned int degree = 1)
229 {
230  switch (degree)
231  {
232  case 0:
233  {
234  return detail::get_identity_kernel<T, Allocator>();
235  }
236  case 1:
237  {
238  detail::kernel_2d<T, Allocator> result(3, 1, 1);
239  std::copy(detail::dy_sobel.begin(), detail::dy_sobel.end(), result.begin());
240  return result;
241  }
242  default:
243  throw std::logic_error("not supported yet");
244  }
245 
246  //to not upset compiler
247  throw std::runtime_error("unreachable statement");
248 }
249 
257 template <typename T = float, typename Allocator = std::allocator<T>>
258 inline auto generate_dy_scharr(unsigned int degree = 1)
260 {
261  switch (degree)
262  {
263  case 0:
264  {
265  return detail::get_identity_kernel<T, Allocator>();
266  }
267  case 1:
268  {
269  detail::kernel_2d<T, Allocator> result(3, 1, 1);
270  std::copy(detail::dy_scharr.begin(), detail::dy_scharr.end(), result.begin());
271  return result;
272  }
273  default:
274  throw std::logic_error("not supported yet");
275  }
276 
277  //to not upset compiler
278  throw std::runtime_error("unreachable statement");
279 }
280 
290 template <typename GradientView, typename OutputView>
292  GradientView dx,
293  GradientView dy,
294  OutputView ddxx,
295  OutputView dxdy,
296  OutputView ddyy)
297 {
298  auto sobel_x = generate_dx_sobel();
299  auto sobel_y = generate_dy_sobel();
300  detail::convolve_2d(dx, sobel_x, ddxx);
301  detail::convolve_2d(dx, sobel_y, dxdy);
302  detail::convolve_2d(dy, sobel_y, ddyy);
303 }
304 
305 }} // namespace boost::gil
306 
307 #endif
auto generate_dx_sobel(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generates Sobel operator in horizontal directionGenerates a kernel which will represent Sobel operato...
Definition: numeric.hpp:165
defined(BOOST_NO_CXX17_HDR_MEMORY_RESOURCE)
Definition: algorithm.hpp:36
auto generate_normalized_mean(std::size_t side_length) -> detail::kernel_2d< T, Allocator >
Generate mean kernelFills supplied view with normalized mean in which all entries will be equal to.
Definition: numeric.hpp:92
auto generate_dx_scharr(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generate Scharr operator in horizontal directionGenerates a kernel which will represent Scharr operat...
Definition: numeric.hpp:196
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:291
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:46
auto generate_unnormalized_mean(std::size_t side_length) -> detail::kernel_2d< T, Allocator >
Generate kernel with all 1sFills supplied view with 1s (ones)
Definition: numeric.hpp:112
auto generate_dy_sobel(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generates Sobel operator in vertical directionGenerates a kernel which will represent Sobel operator ...
Definition: numeric.hpp:227
auto generate_dy_scharr(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generate Scharr operator in vertical directionGenerates a kernel which will represent Scharr operator...
Definition: numeric.hpp:258
auto generate_gaussian_kernel(std::size_t side_length, double sigma) -> detail::kernel_2d< T, Allocator >
Generate Gaussian kernelFills supplied view with values taken from Gaussian distribution....
Definition: numeric.hpp:132
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:145
variable-size kernel
Definition: kernel.hpp:272