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 ///
normalized_sinc(double x)33 inline double normalized_sinc(double x)
34 {
35 return std::sin(x * boost::gil::detail::pi) / (x * boost::gil::detail::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)
lanczos(double x,std::ptrdiff_t 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 #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
compute_tensor_entries(boost::gil::gray16s_view_t dx,boost::gil::gray16s_view_t dy,boost::gil::gray32f_view_t m11,boost::gil::gray32f_view_t m12_21,boost::gil::gray32f_view_t m22)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
84 /// \brief Generate mean kernel
85 /// \ingroup ImageProcessingMath
86 ///
87 /// Fills supplied view with normalized mean
88 /// in which all entries will be equal to
89 /// \code 1 / (dst.size()) \endcode
90 template <typename T = float, typename Allocator = std::allocator<T>>
generate_normalized_mean(std::size_t side_length)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
105 /// \brief Generate kernel with all 1s
106 /// \ingroup ImageProcessingMath
107 ///
108 /// Fills supplied view with 1s (ones)
109 template <typename T = float, typename Allocator = std::allocator<T>>
generate_unnormalized_mean(std::size_t side_length)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
123 /// \brief Generate Gaussian kernel
124 /// \ingroup ImageProcessingMath
125 ///
126 /// Fills supplied view with values taken from Gaussian distribution. See
127 /// https://en.wikipedia.org/wiki/Gaussian_blur
128 template <typename T = float, typename Allocator = std::allocator<T>>
generate_gaussian_kernel(std::size_t side_length,double sigma)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
154 /// \brief Generates Sobel operator in horizontal direction
155 /// \ingroup ImageProcessingMath
156 ///
157 /// Generates a kernel which will represent Sobel operator in
158 /// horizontal direction of specified degree (no need to convolve multiple times
159 /// to obtain the desired degree).
160 /// https://www.researchgate.net/publication/239398674_An_Isotropic_3_3_Image_Gradient_Operator
161 template <typename T = float, typename Allocator = std::allocator<T>>
generate_dx_sobel(unsigned int degree=1)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
184 /// \brief Generate Scharr operator in horizontal direction
185 /// \ingroup ImageProcessingMath
186 ///
187 /// Generates a kernel which will represent Scharr operator in
188 /// horizontal direction of specified degree (no need to convolve multiple times
189 /// to obtain the desired degree).
190 /// https://www.researchgate.net/profile/Hanno_Scharr/publication/220955743_Optimal_Filters_for_Extended_Optical_Flow/links/004635151972eda98f000000/Optimal-Filters-for-Extended-Optical-Flow.pdf
191 template <typename T = float, typename Allocator = std::allocator<T>>
generate_dx_scharr(unsigned int degree=1)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
214 /// \brief Generates Sobel operator in vertical direction
215 /// \ingroup ImageProcessingMath
216 ///
217 /// Generates a kernel which will represent Sobel operator in
218 /// vertical direction of specified degree (no need to convolve multiple times
219 /// to obtain the desired degree).
220 /// https://www.researchgate.net/publication/239398674_An_Isotropic_3_3_Image_Gradient_Operator
221 template <typename T = float, typename Allocator = std::allocator<T>>
generate_dy_sobel(unsigned int degree=1)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
244 /// \brief Generate Scharr operator in vertical direction
245 /// \ingroup ImageProcessingMath
246 ///
247 /// Generates a kernel which will represent Scharr operator in
248 /// vertical direction of specified degree (no need to convolve multiple times
249 /// to obtain the desired degree).
250 /// https://www.researchgate.net/profile/Hanno_Scharr/publication/220955743_Optimal_Filters_for_Extended_Optical_Flow/links/004635151972eda98f000000/Optimal-Filters-for-Extended-Optical-Flow.pdf
251 template <typename T = float, typename Allocator = std::allocator<T>>
generate_dy_scharr(unsigned int degree=1)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
274 /// \brief Compute xy gradient, and second order x and y gradients
275 /// \ingroup ImageProcessingMath
276 ///
277 /// Hessian matrix is defined as a matrix of partial derivates
278 /// for 2d case, it is [[ddxx, dxdy], [dxdy, ddyy].
279 /// d stands for derivative, and x or y stand for direction.
280 /// For example, dx stands for derivative (gradient) in horizontal
281 /// direction, and ddxx means second order derivative in horizon direction
282 /// https://en.wikipedia.org/wiki/Hessian_matrix
283 template <typename GradientView, typename OutputView>
compute_hessian_entries(GradientView dx,GradientView dy,OutputView ddxx,OutputView dxdy,OutputView ddyy)284 inline void compute_hessian_entries(
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
301