Implement Sobel and Scharr operators (#392)

* Implement Sobel and Scharr operators

This commit adds Sobel and Scharr
operators with support for 0th and 1st
degrees with other degrees planned for
later

* Migrate and fix Harris example

Generate Harris entries now uses
signed image view.
The Harris corner detector example
now uses the Scharr filter generator
and convolve_2d to reduce amount
of code needed.

* Fix and migrate Hessian example

The Hessian example now uses signed
image views and uses newly added kernel
generators to compute gradients

* Fix Harris and Hessian tests

The tests broke due to migration to
signed views in algorithms, but tests
were not adjusted

* Fix Jamfile for example/sobel_scharr.cpp

* Cosmetic changes

* Commented out fail tests

* Fixed pixel16 used in image16s

In Harris and Hessian tests, unsigned
pixel values was used to construct
signed image, which was causing
appveyor to error out.

* Reenable failing targets

* Unify kernel generator interface

This commit makes all kernel
generator functions to return kernel_2d
and adapts dependant threshold
function to use the new interface

* Migrate Hessian and Harris tests

Migrate Hessian and Harris tests to new
interface for kernel generators

* Migrate Harris and Hessian examples

Harris and Hessian examples now use
new interface for kernel generation

* Migrate simple_kernels tests

simple_kernels are now using kernel_2d
interface

* Add missing return

Normalized mean generation had missing
return at the end of the function

* Adapt code to namespace move

This commit reacts to kernel_2d,
convolve_2d being moved to
namespace detail

.gitignore

@@ -33,3 +33,7 @@ doc/warnings.txt
 # Conan
 conaninfo.txt
 graph_info.json
+
+# clangd
+compile_commands.json
+.clangd/

example/Jamfile

@@ -23,6 +23,7 @@ local sources =
     dynamic_image.cpp
     histogram.cpp
     harris.cpp
+    sobel_scharr.cpp
     ;
 
 local targets ;

example/harris.cpp

@@ -10,6 +10,7 @@
 #include <boost/gil/extension/io/png.hpp>
 #include <boost/gil/image_processing/numeric.hpp>
 #include <boost/gil/image_processing/harris.hpp>
+#include <boost/gil/extension/numeric/convolve.hpp>
 #include <vector>
 #include <functional>
 #include <set>
@@ -97,64 +98,6 @@ void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_
     }
 }
 
-void calculate_gradients(gil::gray8_view_t input,
-    gil::gray16_view_t x_gradient,
-    gil::gray16_view_t y_gradient)
-{
-    using x_coord_t = gil::gray16_view_t::x_coord_t;
-    using y_coord_t = gil::gray16_view_t::y_coord_t;
-    using pixel_t = std::remove_reference<decltype(x_gradient(0, 0))>::type;
-    using channel_t = typename std::remove_reference<
-                        decltype(
-                            std::declval<gil::gray16_pixel_t>().at(
-                                std::integral_constant<int, 0>{}
-                            )
-                        )
-                       >::type;
-
-    constexpr double x_kernel[3][3] =
-    {
-        {1, 0, -1},
-        {2, 0, -2},
-        {1, 0, -1}
-    };
-    constexpr double y_kernel[3][3] =
-    {
-        {1, 2, 1},
-        {0, 0, 0},
-        {-1, -2, -1}
-    };
-    constexpr auto chosen_channel = std::integral_constant<int, 0>{};
-    for (y_coord_t y = 1; y < input.height() - 1; ++y)
-    {
-        for (x_coord_t x = 1; x < input.width() - 1; ++x)
-        {
-            gil::gray16_pixel_t x_result;
-            boost::gil::static_transform(x_result, x_result,
-                [](channel_t) { return static_cast<channel_t>(0); });
-            gil::gray16_pixel_t y_result;
-            boost::gil::static_transform(y_result, y_result,
-                [](channel_t) { return static_cast<channel_t>(0); });
-            for (y_coord_t y_filter = 0; y_filter < 2; ++y_filter)
-            {
-                for (x_coord_t x_filter = 0; x_filter < 2; ++x_filter)
-                {
-                    auto adjusted_y = y + y_filter - 1;
-                    auto adjusted_x = x + x_filter - 1;
-                    x_result.at(std::integral_constant<int, 0>{}) +=
-                        input(adjusted_x, adjusted_y).at(chosen_channel)
-                        * x_kernel[y_filter][x_filter];
-                    y_result.at(std::integral_constant<int, 0>{}) +=
-                        input(adjusted_x, adjusted_y).at(chosen_channel)
-                        * y_kernel[y_filter][x_filter];
-                }
-            }
-            x_gradient(x, y) = static_cast<std::uint8_t>(x_result.at(chosen_channel));
-            y_gradient(x, y) = static_cast<std::uint8_t>(y_result.at(chosen_channel));
-        }
-    }
-}
-
 std::vector<gil::point_t> suppress(
     gil::gray32f_view_t harris_response,
     double harris_response_threshold)
@@ -209,8 +152,8 @@ int main(int argc, char* argv[])
         return -1;
     }
 
-    long int window_size = std::stoi(argv[2]);
-    double discrimnation_constant = std::stod(argv[3]);
+    std::size_t window_size = std::stoul(argv[2]);
+    double discrimnation_constant = std::stof(argv[3]);
     long harris_response_threshold = std::stol(argv[4]);
 
     gil::rgb8_image_t input_image;
@@ -222,12 +165,15 @@ int main(int argc, char* argv[])
     gil::gray8_image_t smoothed_image(grayscaled.dimensions());
     auto smoothed = gil::view(smoothed_image);
     apply_gaussian_blur(gil::view(grayscaled), smoothed);
-    gil::gray16_image_t x_gradient_image(grayscaled.dimensions());
-    gil::gray16_image_t y_gradient_image(grayscaled.dimensions());
+    gil::gray16s_image_t x_gradient_image(grayscaled.dimensions());
+    gil::gray16s_image_t y_gradient_image(grayscaled.dimensions());
 
     auto x_gradient = gil::view(x_gradient_image);
     auto y_gradient = gil::view(y_gradient_image);
-    calculate_gradients(smoothed, x_gradient, y_gradient);
+    auto scharr_x = gil::generate_dx_scharr();
+    gil::detail::convolve_2d(smoothed, scharr_x, x_gradient);
+    auto scharr_y = gil::generate_dy_scharr();
+    gil::detail::convolve_2d(smoothed, scharr_y, y_gradient);
 
     gil::gray32f_image_t m11(x_gradient.dimensions());
     gil::gray32f_image_t m12_21(x_gradient.dimensions());
@@ -241,13 +187,12 @@ int main(int argc, char* argv[])
     );
 
     gil::gray32f_image_t harris_response(x_gradient.dimensions());
-    gil::gray32f_image_t gaussian_kernel(gil::point_t(5, 5));
-    gil::generate_gaussian_kernel(gil::view(gaussian_kernel), 0.84089642);
+    auto gaussian_kernel = gil::generate_gaussian_kernel(window_size, 0.84089642);
     gil::compute_harris_responses(
         gil::view(m11),
         gil::view(m12_21),
         gil::view(m22),
-        gil::view(gaussian_kernel),
+        gaussian_kernel,
         discrimnation_constant,
         gil::view(harris_response)
     );

example/hessian.cpp

@@ -103,64 +103,6 @@ void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_
     }
 }
 
-void calculate_gradients(gil::gray8_view_t input,
-    gil::gray16_view_t x_gradient,
-    gil::gray16_view_t y_gradient)
-{
-    using x_coord_t = gil::gray16_view_t::x_coord_t;
-    using y_coord_t = gil::gray16_view_t::y_coord_t;
-    using pixel_t = std::remove_reference<decltype(x_gradient(0, 0))>::type;
-    using channel_t = typename std::remove_reference<
-                        decltype(
-                            std::declval<gil::gray16_pixel_t>().at(
-                                std::integral_constant<int, 0>{}
-                            )
-                        )
-                       >::type;
-
-    constexpr double x_kernel[3][3] =
-    {
-        {1, 0, -1},
-        {2, 0, -2},
-        {1, 0, -1}
-    };
-    constexpr double y_kernel[3][3] =
-    {
-        {1, 2, 1},
-        {0, 0, 0},
-        {-1, -2, -1}
-    };
-    constexpr auto chosen_channel = std::integral_constant<int, 0>{};
-    for (y_coord_t y = 1; y < input.height() - 1; ++y)
-    {
-        for (x_coord_t x = 1; x < input.width() - 1; ++x)
-        {
-            gil::gray16_pixel_t x_result;
-            boost::gil::static_transform(x_result, x_result,
-                [](channel_t) { return static_cast<channel_t>(0); });
-            gil::gray16_pixel_t y_result;
-            boost::gil::static_transform(y_result, y_result,
-                [](channel_t) { return static_cast<channel_t>(0); });
-            for (y_coord_t y_filter = 0; y_filter < 2; ++y_filter)
-            {
-                for (x_coord_t x_filter = 0; x_filter < 2; ++x_filter)
-                {
-                    auto adjusted_y = y + y_filter - 1;
-                    auto adjusted_x = x + x_filter - 1;
-                    x_result.at(std::integral_constant<int, 0>{}) +=
-                        input(adjusted_x, adjusted_y).at(chosen_channel)
-                        * x_kernel[y_filter][x_filter];
-                    y_result.at(std::integral_constant<int, 0>{}) +=
-                        input(adjusted_x, adjusted_y).at(chosen_channel)
-                        * y_kernel[y_filter][x_filter];
-                }
-            }
-            x_gradient(x, y) = static_cast<std::uint8_t>(x_result.at(chosen_channel));
-            y_gradient(x, y) = static_cast<std::uint8_t>(y_result.at(chosen_channel));
-        }
-    }
-}
-
 std::vector<gil::point_t> suppress(
     gil::gray32f_view_t harris_response,
     double harris_response_threshold)
@@ -214,7 +156,7 @@ int main(int argc, char* argv[]) {
         return -1;
     }
 
-    long int window_size = std::stoi(argv[2]);
+    std::size_t window_size = std::stoul(argv[2]);
     long hessian_determinant_threshold = std::stol(argv[3]);
 
     gil::rgb8_image_t input_image;
@@ -226,12 +168,15 @@ int main(int argc, char* argv[]) {
     gil::gray8_image_t smoothed_image(grayscaled.dimensions());
     auto smoothed = gil::view(smoothed_image);
     apply_gaussian_blur(gil::view(grayscaled), smoothed);
-    gil::gray16_image_t x_gradient_image(grayscaled.dimensions());
-    gil::gray16_image_t y_gradient_image(grayscaled.dimensions());
+    gil::gray16s_image_t x_gradient_image(grayscaled.dimensions());
+    gil::gray16s_image_t y_gradient_image(grayscaled.dimensions());
 
     auto x_gradient = gil::view(x_gradient_image);
     auto y_gradient = gil::view(y_gradient_image);
-    calculate_gradients(smoothed, x_gradient, y_gradient);
+    auto scharr_x = gil::generate_dx_scharr();
+    gil::detail::convolve_2d(smoothed, scharr_x, x_gradient);
+    auto scharr_y = gil::generate_dy_scharr();
+    gil::detail::convolve_2d(smoothed, scharr_y, y_gradient);
 
     gil::gray32f_image_t m11(x_gradient.dimensions());
     gil::gray32f_image_t m12_21(x_gradient.dimensions());
@@ -245,13 +190,12 @@ int main(int argc, char* argv[]) {
     );
 
     gil::gray32f_image_t hessian_response(x_gradient.dimensions());
-    gil::gray32f_image_t gaussian_kernel(gil::point_t(5, 5));
-    gil::generate_gaussian_kernel(gil::view(gaussian_kernel), 0.84089642);
+    auto gaussian_kernel = gil::generate_gaussian_kernel(window_size, 0.84089642);
     gil::compute_hessian_responses(
         gil::view(m11),
         gil::view(m12_21),
         gil::view(m22),
-        gil::view(gaussian_kernel),
+        gaussian_kernel,
         gil::view(hessian_response)
     );
 

example/sobel_scharr.cpp

@@ -0,0 +1,45 @@
+#include <boost/gil/typedefs.hpp>
+#include <boost/gil/image_processing/numeric.hpp>
+#include <boost/gil/extension/io/png.hpp>
+#include <boost/gil/extension/numeric/convolve.hpp>
+#include <string>
+#include <iostream>
+
+namespace gil = boost::gil;
+
+int main(int argc, char* argv[])
+{
+    if (argc != 5)
+    {
+        std::cerr << "usage: " << argv[0] << ": <input.png> <sobel|scharr> <output-x.png> <output-y.png>\n";
+        return -1;
+    }
+
+    gil::gray8_image_t input_image;
+    gil::read_image(argv[1], input_image, gil::png_tag{});
+    auto input = gil::view(input_image);
+    auto filter_type = std::string(argv[2]);
+
+    gil::gray16_image_t dx_image(input_image.dimensions());
+    auto dx = gil::view(dx_image);
+    gil::gray16_image_t dy_image(input_image.dimensions());
+    auto dy = gil::view(dy_image);
+    if (filter_type == "sobel")
+    {
+        gil::convolve_2d(input, gil::generate_dx_sobel(1), dx);
+        gil::convolve_2d(input, gil::generate_dy_sobel(1), dy);
+    }
+    else if (filter_type == "scharr")
+    {
+        gil::convolve_2d(input, gil::generate_dx_scharr(1), dx);
+        gil::convolve_2d(input, gil::generate_dy_scharr(1), dy);
+    }
+    else
+    {
+        std::cerr << "unrecognized gradient filter type. Must be either sobel or scharr\n";
+        return -1;
+    }
+
+    gil::write_view(argv[3], dx, gil::png_tag{});
+    gil::write_view(argv[4], dy, gil::png_tag{});
+}

include/boost/gil/detail/math.hpp

@@ -8,10 +8,27 @@
 #ifndef BOOST_GIL_IMAGE_PROCESSING_DETAIL_MATH_HPP
 #define BOOST_GIL_IMAGE_PROCESSING_DETAIL_MATH_HPP
 
+#include <array>
+#include <boost/gil/extension/numeric/kernel.hpp>
+
 namespace boost { namespace gil {
 
 static constexpr double pi = 3.14159265358979323846;
 
+static constexpr std::array<float, 9> dx_sobel = {-1, 0, 1, -2, 0, 2, -1, 0, 1};
+static constexpr std::array<float, 9> dx_scharr = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
+
+static constexpr std::array<float, 9> dy_sobel = {1, 2, 1, 0, 0, 0, -1, -2, -1};
+static constexpr std::array<float, 9> dy_scharr = {1, 1, 1, 0, 0, 0, -1, -1, -1};
+
+template <typename T, typename Allocator>
+inline detail::kernel_2d<T, Allocator> get_identity_kernel()
+{
+    detail::kernel_2d<T, Allocator> kernel(1, 0, 0);
+    kernel[0] = 1;
+    return kernel;
+}
+
 }} // namespace boost::gil
 
 #endif

include/boost/gil/image_processing/harris.hpp

@@ -10,6 +10,7 @@
 
 #include <boost/gil/image_view.hpp>
 #include <boost/gil/typedefs.hpp>
+#include <boost/gil/extension/numeric/kernel.hpp>
 
 namespace boost { namespace gil {
 /// \defgroup CornerDetectionAlgorithms
@@ -29,11 +30,12 @@ namespace boost { namespace gil {
 /// to compute sum of corresponding entries. k is a discrimination
 /// constant against edges (usually in range 0.04 to 0.06).
 /// harris_response is an out parameter that will contain the Harris responses.
+template <typename T, typename Allocator>
 void compute_harris_responses(
     boost::gil::gray32f_view_t m11,
     boost::gil::gray32f_view_t m12_21,
     boost::gil::gray32f_view_t m22,
-    boost::gil::gray32f_view_t weights,
+    boost::gil::detail::kernel_2d<T, Allocator> weights,
     float k,
     boost::gil::gray32f_view_t harris_response)
 {
@@ -42,7 +44,7 @@ void compute_harris_responses(
             " tensor from the same image");
     }
 
-    auto const window_length = weights.dimensions().x;
+    auto const window_length = weights.size();
     auto const width = m11.width();
     auto const height = m11.height();
     auto const half_length = window_length / 2;
@@ -61,14 +63,11 @@ void compute_harris_responses(
                     x_kernel < window_length;
                     ++x_kernel) {
                     ddxx += m11(x + x_kernel - half_length, y + y_kernel - half_length)
-                        .at(std::integral_constant<int, 0>{})
-                        * weights(x_kernel, y_kernel).at(std::integral_constant<int, 0>{});
+                        .at(std::integral_constant<int, 0>{}) * weights.at(x_kernel, y_kernel);
                     dxdy += m12_21(x + x_kernel - half_length, y + y_kernel - half_length)
-                        .at(std::integral_constant<int, 0>{})
-                        * weights(x_kernel, y_kernel).at(std::integral_constant<int, 0>{});
+                        .at(std::integral_constant<int, 0>{}) * weights.at(x_kernel, y_kernel);
                     ddyy += m22(x + x_kernel - half_length, y + y_kernel - half_length)
-                        .at(std::integral_constant<int, 0>{})
-                        * weights(x_kernel, y_kernel).at(std::integral_constant<int, 0>{});
+                        .at(std::integral_constant<int, 0>{}) * weights.at(x_kernel, y_kernel);
                 }
             }
             auto det = (ddxx * ddyy) - dxdy * dxdy;