Generic Image Library: numeric.hpp Source File

###  | #  |

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
21namespace boost { namespace gil {
22
34inline 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
63inline 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
91template <typename T = float, typename Allocator = std::allocator<T>>
92 inline auto generate_normalized_mean(std::size_t side_length)
93 -> detail::kernel_2d<T, Allocator>
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
111template <typename T = float, typename Allocator = std::allocator<T>>
112 inline auto generate_unnormalized_mean(std::size_t side_length)
113 -> detail::kernel_2d<T, Allocator>
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
131template <typename T = float, typename Allocator = std::allocator<T>>
132 inline auto generate_gaussian_kernel(std::size_t side_length, double sigma)
133 -> detail::kernel_2d<T, Allocator>
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 const middle = side_length / 2;
140 std::vector<T, Allocator> values(side_length * side_length);
141 T sum{0};
142 for (std::size_t y = 0; y < side_length; ++y)
143 {
144 for (std::size_t x = 0; x < side_length; ++x)
145 {
146 const auto delta_x = x - middle;
147 const auto delta_y = y - middle;
148 const auto power = static_cast<double>(delta_x * delta_x + delta_y * delta_y) / (2 * sigma * sigma);
149 const double nominator = std::exp(-power);
150 const auto value = static_cast<T>(nominator / denominator);
151 values[y * side_length + x] = value;
152 sum += value;
153 }
154 }
155
156 // normalize so that Gaussian kernel sums up to 1.
157 std::transform(values.begin(), values.end(), values.begin(), [&sum](const auto & v) { return v/sum; });
158
159 return detail::kernel_2d<T, Allocator>(values.begin(), values.size(), middle, middle);
160}
161
169template <typename T = float, typename Allocator = std::allocator<T>>
170 inline auto generate_dx_sobel(unsigned int degree = 1)
171 -> detail::kernel_2d<T, Allocator>
172{
173 switch (degree)
174 {
175 case 0:
176 {
177 return detail::get_identity_kernel<T, Allocator>();
178 }
179 case 1:
180 {
181 detail::kernel_2d<T, Allocator> result(3, 1, 1);
182 std::copy(detail::dx_sobel.begin(), detail::dx_sobel.end(), result.begin());
183 return result;
184 }
185 default:
186 throw std::logic_error("not supported yet");
187 }
188
189 //to not upset compiler
190 throw std::runtime_error("unreachable statement");
191}
192
200template <typename T = float, typename Allocator = std::allocator<T>>
201 inline auto generate_dx_scharr(unsigned int degree = 1)
202 -> detail::kernel_2d<T, Allocator>
203{
204 switch (degree)
205 {
206 case 0:
207 {
208 return detail::get_identity_kernel<T, Allocator>();
209 }
210 case 1:
211 {
212 detail::kernel_2d<T, Allocator> result(3, 1, 1);
213 std::copy(detail::dx_scharr.begin(), detail::dx_scharr.end(), result.begin());
214 return result;
215 }
216 default:
217 throw std::logic_error("not supported yet");
218 }
219
220 //to not upset compiler
221 throw std::runtime_error("unreachable statement");
222}
223
231template <typename T = float, typename Allocator = std::allocator<T>>
232 inline auto generate_dy_sobel(unsigned int degree = 1)
233 -> detail::kernel_2d<T, Allocator>
234{
235 switch (degree)
236 {
237 case 0:
238 {
239 return detail::get_identity_kernel<T, Allocator>();
240 }
241 case 1:
242 {
243 detail::kernel_2d<T, Allocator> result(3, 1, 1);
244 std::copy(detail::dy_sobel.begin(), detail::dy_sobel.end(), result.begin());
245 return result;
246 }
247 default:
248 throw std::logic_error("not supported yet");
249 }
250
251 //to not upset compiler
252 throw std::runtime_error("unreachable statement");
253}
254
262template <typename T = float, typename Allocator = std::allocator<T>>
263 inline auto generate_dy_scharr(unsigned int degree = 1)
264 -> detail::kernel_2d<T, Allocator>
265{
266 switch (degree)
267 {
268 case 0:
269 {
270 return detail::get_identity_kernel<T, Allocator>();
271 }
272 case 1:
273 {
274 detail::kernel_2d<T, Allocator> result(3, 1, 1);
275 std::copy(detail::dy_scharr.begin(), detail::dy_scharr.end(), result.begin());
276 return result;
277 }
278 default:
279 throw std::logic_error("not supported yet");
280 }
281
282 //to not upset compiler
283 throw std::runtime_error("unreachable statement");
284}
285
295template <typename GradientView, typename OutputView>
296 inline void compute_hessian_entries(
297 GradientView dx,
298 GradientView dy,
299 OutputView ddxx,
300 OutputView dxdy,
301 OutputView ddyy)
302{
303 auto sobel_x = generate_dx_sobel();
304 auto sobel_y = generate_dy_sobel();
305 detail::convolve_2d(dx, sobel_x, ddxx);
306 detail::convolve_2d(dx, sobel_y, dxdy);
307 detail::convolve_2d(dy, sobel_y, ddyy);
308}
309
310}} // namespace boost::gil
311
312#endif
boost::gil::detail::kernel_2d
variable-size kernel
Definition kernel.hpp:273
boost::gil::generate_gaussian_kernel
auto generate_gaussian_kernel(std::size_t side_length, double sigma) -> detail::kernel_2d< T, Allocator >
Generate Gaussian kernel.
Definition numeric.hpp:132
boost::gil::generate_dy_scharr
auto generate_dy_scharr(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generate Scharr operator in vertical direction.
Definition numeric.hpp:263
boost::gil::compute_hessian_entries
void compute_hessian_entries(GradientView dx, GradientView dy, OutputView ddxx, OutputView dxdy, OutputView ddyy)
Compute xy gradient, and second order x and y gradients.
Definition numeric.hpp:296
boost::gil::generate_dx_scharr
auto generate_dx_scharr(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generate Scharr operator in horizontal direction.
Definition numeric.hpp:201
boost::gil::generate_unnormalized_mean
auto generate_unnormalized_mean(std::size_t side_length) -> detail::kernel_2d< T, Allocator >
Generate kernel with all 1s.
Definition numeric.hpp:112
boost::gil::generate_dy_sobel
auto generate_dy_sobel(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generates Sobel operator in vertical direction.
Definition numeric.hpp:232
boost::gil::generate_dx_sobel
auto generate_dx_sobel(unsigned int degree=1) -> detail::kernel_2d< T, Allocator >
Generates Sobel operator in horizontal direction.
Definition numeric.hpp:170
boost::gil::lanczos
double lanczos(double x, std::ptrdiff_t a)
Lanczos response at point x.
Definition numeric.hpp:46
boost::gil::generate_normalized_mean
auto generate_normalized_mean(std::size_t side_length) -> detail::kernel_2d< T, Allocator >
Generate mean kernel.
Definition numeric.hpp:92
std::copy
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
boost
defined(BOOST_NO_CXX17_HDR_MEMORY_RESOURCE)
Definition algorithm.hpp:36