gil/image_processing/numeric.hpp

0001 //
0002 // Copyright 2019 Olzhas Zhumabek <anonymous.from.applecity@gmail.com>
0003 // Copyright 2021 Pranam Lashkari <plashkari628@gmail.com>
0004 //
0005 // Use, modification and distribution are subject to the Boost Software License,
0006 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0007 // http://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 #ifndef BOOST_GIL_IMAGE_PROCESSING_NUMERIC_HPP
0010 #define BOOST_GIL_IMAGE_PROCESSING_NUMERIC_HPP
0011
0012 #include <boost/gil/image_processing/kernel.hpp>
0013 #include <boost/gil/image_processing/convolve.hpp>
0014 #include <boost/gil/image_view.hpp>
0015 #include <boost/gil/typedefs.hpp>
0016 #include <boost/gil/detail/math.hpp>
0017 // fixes ambigious call to std::abs, https://stackoverflow.com/a/30084734/4593721
0018 #include <cstdlib>
0019 #include <cmath>
0020
0021 namespace boost { namespace gil {
0022
0023 /// \defgroup ImageProcessingMath
0024 /// \brief Math operations for IP algorithms
0025 ///
0026 /// This is mostly handful of mathemtical operations that are required by other
0027 /// image processing algorithms
0028 ///
0029 /// \brief Normalized cardinal sine
0030 /// \ingroup ImageProcessingMath
0031 ///
0032 /// normalized_sinc(x) = sin(pi * x) / (pi * x)
0033 ///
0034 inline double normalized_sinc(double x)
0035 {
0036     return std::sin(x * boost::gil::detail::pi) / (x * boost::gil::detail::pi);
0037 }
0038
0039 /// \brief Lanczos response at point x
0040 /// \ingroup ImageProcessingMath
0041 ///
0042 /// Lanczos response is defined as:
0043 /// x == 0: 1
0044 /// -a < x && x < a: 0
0045 /// otherwise: normalized_sinc(x) / normalized_sinc(x / a)
0046 inline double lanczos(double x, std::ptrdiff_t a)
0047 {
0048     // means == but <= avoids compiler warning
0049     if (0 <= x && x <= 0)
0050         return 1;
0051
0052     if (static_cast<double>(-a) < x && x < static_cast<double>(a))
0053         return normalized_sinc(x) / normalized_sinc(x / static_cast<double>(a));
0054
0055     return 0;
0056 }
0057
0058 #if BOOST_WORKAROUND(BOOST_MSVC, >= 1400)
0059 #pragma warning(push)
0060 #pragma warning(disable:4244) // 'argument': conversion from 'const Channel' to 'BaseChannelValue', possible loss of data
0061 #endif
0062
0063 inline void compute_tensor_entries(
0064     boost::gil::gray16s_view_t dx,
0065     boost::gil::gray16s_view_t dy,
0066     boost::gil::gray32f_view_t m11,
0067     boost::gil::gray32f_view_t m12_21,
0068     boost::gil::gray32f_view_t m22)
0069 {
0070     for (std::ptrdiff_t y = 0; y < dx.height(); ++y) {
0071         for (std::ptrdiff_t x = 0; x < dx.width(); ++x) {
0072             auto dx_value = dx(x, y);
0073             auto dy_value = dy(x, y);
0074             m11(x, y) = dx_value * dx_value;
0075             m12_21(x, y) = dx_value * dy_value;
0076             m22(x, y) = dy_value * dy_value;
0077         }
0078     }
0079 }
0080
0081 #if BOOST_WORKAROUND(BOOST_MSVC, >= 1400)
0082 #pragma warning(pop)
0083 #endif
0084
0085 /// \brief Generate mean kernel
0086 /// \ingroup ImageProcessingMath
0087 ///
0088 /// Fills supplied view with normalized mean
0089 /// in which all entries will be equal to
0090 /// \code 1 / (dst.size()) \endcode
0091 template <typename T = float, typename Allocator = std::allocator<T>>
0092 inline auto generate_normalized_mean(std::size_t side_length)
0093     -> detail::kernel_2d<T, Allocator>
0094 {
0095     if (side_length % 2 != 1)
0096         throw std::invalid_argument("kernel dimensions should be odd and equal");
0097     const float entry = 1.0f / static_cast<float>(side_length * side_length);
0098
0099     detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
0100     for (auto& cell: result) {
0101         cell = entry;
0102     }
0103
0104     return result;
0105 }
0106
0107 /// \brief Generate kernel with all 1s
0108 /// \ingroup ImageProcessingMath
0109 ///
0110 /// Fills supplied view with 1s (ones)
0111 template <typename T = float, typename Allocator = std::allocator<T>>
0112 inline auto generate_unnormalized_mean(std::size_t side_length)
0113     -> detail::kernel_2d<T, Allocator>
0114 {
0115     if (side_length % 2 != 1)
0116         throw std::invalid_argument("kernel dimensions should be odd and equal");
0117
0118     detail::kernel_2d<T, Allocator> result(side_length, side_length / 2, side_length / 2);
0119     for (auto& cell: result) {
0120         cell = 1.0f;
0121     }
0122
0123     return result;
0124 }
0125
0126 /// \brief Generate Gaussian kernel
0127 /// \ingroup ImageProcessingMath
0128 ///
0129 /// Fills supplied view with values taken from Gaussian distribution. See
0130 /// https://en.wikipedia.org/wiki/Gaussian_blur
0131 template <typename T = float, typename Allocator = std::allocator<T>>
0132 inline auto generate_gaussian_kernel(std::size_t side_length, double sigma)
0133     -> detail::kernel_2d<T, Allocator>
0134 {
0135     if (side_length % 2 != 1)
0136         throw std::invalid_argument("kernel dimensions should be odd and equal");
0137
0138     const double denominator = 2 * boost::gil::detail::pi * sigma * sigma;
0139     auto const middle = side_length / 2;
0140     std::vector<T, Allocator> values(side_length * side_length);
0141     T sum{0};
0142     for (std::size_t y = 0; y < side_length; ++y)
0143     {
0144         for (std::size_t x = 0; x < side_length; ++x)
0145         {
0146             const auto delta_x = x - middle;
0147             const auto delta_y = y - middle;
0148             const auto power = static_cast<double>(delta_x * delta_x + delta_y * delta_y) / (2 * sigma * sigma);
0149             const double nominator = std::exp(-power);
0150             const auto value = static_cast<T>(nominator / denominator);
0151             values[y * side_length + x] = value;
0152             sum += value;
0153         }
0154     }
0155
0156     // normalize so that Gaussian kernel sums up to 1.
0157     std::transform(values.begin(), values.end(), values.begin(), [&sum](const auto & v) { return v/sum; });
0158
0159     return detail::kernel_2d<T, Allocator>(values.begin(), values.size(), middle, middle);
0160 }
0161
0162 /// \brief Generates Sobel operator in horizontal direction
0163 /// \ingroup ImageProcessingMath
0164 ///
0165 /// Generates a kernel which will represent Sobel operator in
0166 /// horizontal direction of specified degree (no need to convolve multiple times
0167 /// to obtain the desired degree).
0168 /// https://www.researchgate.net/publication/239398674_An_Isotropic_3_3_Image_Gradient_Operator
0169 template <typename T = float, typename Allocator = std::allocator<T>>
0170 inline auto generate_dx_sobel(unsigned int degree = 1)
0171     -> detail::kernel_2d<T, Allocator>
0172 {
0173     switch (degree)
0174     {
0175         case 0:
0176         {
0177             return detail::get_identity_kernel<T, Allocator>();
0178         }
0179         case 1:
0180         {
0181             detail::kernel_2d<T, Allocator> result(3, 1, 1);
0182             std::copy(detail::dx_sobel.begin(), detail::dx_sobel.end(), result.begin());
0183             return result;
0184         }
0185         default:
0186             throw std::logic_error("not supported yet");
0187     }
0188
0189     //to not upset compiler
0190     throw std::runtime_error("unreachable statement");
0191 }
0192
0193 /// \brief Generate Scharr operator in horizontal direction
0194 /// \ingroup ImageProcessingMath
0195 ///
0196 /// Generates a kernel which will represent Scharr operator in
0197 /// horizontal direction of specified degree (no need to convolve multiple times
0198 /// to obtain the desired degree).
0199 /// https://www.researchgate.net/profile/Hanno_Scharr/publication/220955743_Optimal_Filters_for_Extended_Optical_Flow/links/004635151972eda98f000000/Optimal-Filters-for-Extended-Optical-Flow.pdf
0200 template <typename T = float, typename Allocator = std::allocator<T>>
0201 inline auto generate_dx_scharr(unsigned int degree = 1)
0202     -> detail::kernel_2d<T, Allocator>
0203 {
0204     switch (degree)
0205     {
0206         case 0:
0207         {
0208             return detail::get_identity_kernel<T, Allocator>();
0209         }
0210         case 1:
0211         {
0212             detail::kernel_2d<T, Allocator> result(3, 1, 1);
0213             std::copy(detail::dx_scharr.begin(), detail::dx_scharr.end(), result.begin());
0214             return result;
0215         }
0216         default:
0217             throw std::logic_error("not supported yet");
0218     }
0219
0220     //to not upset compiler
0221     throw std::runtime_error("unreachable statement");
0222 }
0223
0224 /// \brief Generates Sobel operator in vertical direction
0225 /// \ingroup ImageProcessingMath
0226 ///
0227 /// Generates a kernel which will represent Sobel operator in
0228 /// vertical direction of specified degree (no need to convolve multiple times
0229 /// to obtain the desired degree).
0230 /// https://www.researchgate.net/publication/239398674_An_Isotropic_3_3_Image_Gradient_Operator
0231 template <typename T = float, typename Allocator = std::allocator<T>>
0232 inline auto generate_dy_sobel(unsigned int degree = 1)
0233     -> detail::kernel_2d<T, Allocator>
0234 {
0235     switch (degree)
0236     {
0237         case 0:
0238         {
0239             return detail::get_identity_kernel<T, Allocator>();
0240         }
0241         case 1:
0242         {
0243             detail::kernel_2d<T, Allocator> result(3, 1, 1);
0244             std::copy(detail::dy_sobel.begin(), detail::dy_sobel.end(), result.begin());
0245             return result;
0246         }
0247         default:
0248             throw std::logic_error("not supported yet");
0249     }
0250
0251     //to not upset compiler
0252     throw std::runtime_error("unreachable statement");
0253 }
0254
0255 /// \brief Generate Scharr operator in vertical direction
0256 /// \ingroup ImageProcessingMath
0257 ///
0258 /// Generates a kernel which will represent Scharr operator in
0259 /// vertical direction of specified degree (no need to convolve multiple times
0260 /// to obtain the desired degree).
0261 /// https://www.researchgate.net/profile/Hanno_Scharr/publication/220955743_Optimal_Filters_for_Extended_Optical_Flow/links/004635151972eda98f000000/Optimal-Filters-for-Extended-Optical-Flow.pdf
0262 template <typename T = float, typename Allocator = std::allocator<T>>
0263 inline auto generate_dy_scharr(unsigned int degree = 1)
0264     -> detail::kernel_2d<T, Allocator>
0265 {
0266     switch (degree)
0267     {
0268         case 0:
0269         {
0270             return detail::get_identity_kernel<T, Allocator>();
0271         }
0272         case 1:
0273         {
0274             detail::kernel_2d<T, Allocator> result(3, 1, 1);
0275             std::copy(detail::dy_scharr.begin(), detail::dy_scharr.end(), result.begin());
0276             return result;
0277         }
0278         default:
0279             throw std::logic_error("not supported yet");
0280     }
0281
0282     //to not upset compiler
0283     throw std::runtime_error("unreachable statement");
0284 }
0285
0286 /// \brief Compute xy gradient, and second order x and y gradients
0287 /// \ingroup ImageProcessingMath
0288 ///
0289 /// Hessian matrix is defined as a matrix of partial derivates
0290 /// for 2d case, it is [[ddxx, dxdy], [dxdy, ddyy].
0291 /// d stands for derivative, and x or y stand for direction.
0292 /// For example, dx stands for derivative (gradient) in horizontal
0293 /// direction, and ddxx means second order derivative in horizon direction
0294 /// https://en.wikipedia.org/wiki/Hessian_matrix
0295 template <typename GradientView, typename OutputView>
0296 inline void compute_hessian_entries(
0297     GradientView dx,
0298     GradientView dy,
0299     OutputView ddxx,
0300     OutputView dxdy,
0301     OutputView ddyy)
0302 {
0303     auto sobel_x = generate_dx_sobel();
0304     auto sobel_y = generate_dy_sobel();
0305     detail::convolve_2d(dx, sobel_x, ddxx);
0306     detail::convolve_2d(dx, sobel_y, dxdy);
0307     detail::convolve_2d(dy, sobel_y, ddyy);
0308 }
0309
0310 }} // namespace boost::gil
0311
0312 #endif