From 2b40c82076e4532403eaa8b2d94f7b7dcd883e37 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 3 Aug 2019 22:52:33 +0000 Subject: [PATCH 01/11] Hessian plain determinant This commit partially implements Hessian corner detector, but only uses determinant as Hessian response --- example/hessian.cpp | 252 ++++++++++++++++++ .../boost/gil/image_processing/hessian.hpp | 77 ++++++ .../boost/gil/image_processing/numeric.hpp | 70 +++++ 3 files changed, 399 insertions(+) create mode 100644 example/hessian.cpp create mode 100644 include/boost/gil/image_processing/hessian.hpp diff --git a/example/hessian.cpp b/example/hessian.cpp new file mode 100644 index 0000000000..ae8c8d70ec --- /dev/null +++ b/example/hessian.cpp @@ -0,0 +1,252 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace gil = boost::gil; + +// some images might produce artifacts +// when converted to grayscale, +// which was previously observed on +// canny edge detector for test input +// used for this example +gil::gray8_image_t to_grayscale(gil::rgb8_view_t original) +{ + gil::gray8_image_t output_image(original.dimensions()); + auto output = gil::view(output_image); + constexpr double max_channel_intensity = (std::numeric_limits::max)(); + for (long int y = 0; y < original.height(); ++y) { + for (long int x = 0; x < original.width(); ++x) { + // scale the values into range [0, 1] and calculate linear intensity + double red_intensity = original(x, y).at(std::integral_constant{}) + / max_channel_intensity; + double green_intensity = original(x, y).at(std::integral_constant{}) + / max_channel_intensity; + double blue_intensity = original(x, y).at(std::integral_constant{}) + / max_channel_intensity; + auto linear_luminosity = 0.2126 * red_intensity + + 0.7152 * green_intensity + + 0.0722 * blue_intensity; + + // perform gamma adjustment + double gamma_compressed_luminosity = 0; + if (linear_luminosity < 0.0031308) { + gamma_compressed_luminosity = linear_luminosity * 12.92; + } else { + gamma_compressed_luminosity = 1.055 * std::pow(linear_luminosity, 1 / 2.4) - 0.055; + } + + // since now it is scaled, descale it back + output(x, y) = gamma_compressed_luminosity * max_channel_intensity; + } + } + + return output_image; +} + +void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_view) +{ + constexpr static auto filterHeight = 5ull; + constexpr static auto filterWidth = 5ull; + constexpr static double filter[filterHeight][filterWidth] = + { + 2, 4, 6, 4, 2, + 4, 9, 12, 9, 4, + 5, 12, 15, 12, 5, + 4, 9, 12, 9, 4, + 2, 4, 5, 4, 2, + }; + constexpr double factor = 1.0 / 159; + constexpr double bias = 0.0; + + const auto height = input_view.height(); + const auto width = input_view.width(); + for (long x = 0; x < width; ++x) { + for (long y = 0; y < height; ++y) { + double intensity = 0.0; + for (size_t filter_y = 0; filter_y < filterHeight; ++filter_y) { + for (size_t filter_x = 0; filter_x < filterWidth; ++filter_x) { + int image_x = x - filterWidth / 2 + filter_x; + int image_y = y - filterHeight / 2 + filter_y; + if (image_x >= input_view.width() || image_x < 0 + || image_y >= input_view.height() || image_y < 0) { + continue; + } + auto& pixel = input_view(image_x, image_y); + intensity += pixel.at(std::integral_constant{}) + * filter[filter_y][filter_x]; + } + } + auto& pixel = output_view(gil::point_t(x, y)); + pixel = (std::min)((std::max)(int(factor * intensity + bias), 0), 255); + } + + } +} + +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::type; + using channel_t = typename std::remove_reference< + decltype( + std::declval().at( + std::integral_constant{} + ) + ) + >::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{}; + 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(0); }); + gil::gray16_pixel_t y_result; + boost::gil::static_transform(y_result, y_result, + [](channel_t) { return static_cast(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{}) += + input(adjusted_x, adjusted_y).at(chosen_channel) + * x_kernel[y_filter][x_filter]; + y_result.at(std::integral_constant{}) += + input(adjusted_x, adjusted_y).at(chosen_channel) + * y_kernel[y_filter][x_filter]; + } + } + x_gradient(x, y) = static_cast(x_result.at(chosen_channel)); + y_gradient(x, y) = static_cast(y_result.at(chosen_channel)); + } + } +} + +std::vector suppress( + gil::gray32f_view_t harris_response, + double harris_response_threshold) +{ + std::vector corner_points; + for (gil::gray32f_view_t::coord_t y = 1; y < harris_response.height() - 1; ++y) + { + for (gil::gray32f_view_t::coord_t x = 1; x < harris_response.width() - 1; ++x) + { + auto value = [](gil::gray32f_pixel_t pixel) { + return pixel.at(std::integral_constant{}); + }; + double values[9] = { + value(harris_response(x - 1, y - 1)), + value(harris_response(x, y - 1)), + value(harris_response(x + 1, y - 1)), + value(harris_response(x - 1, y)), + value(harris_response(x, y)), + value(harris_response(x + 1, y)), + value(harris_response(x - 1, y + 1)), + value(harris_response(x, y + 1)), + value(harris_response(x + 1, y + 1)) + }; + + auto maxima = *std::max_element( + values, + values + 9, + [](double lhs, double rhs) + { + return lhs < rhs; + } + ); + + if (maxima == value(harris_response(x, y)) + && std::count(values, values + 9, maxima) == 1 + && maxima >= harris_response_threshold) + { + corner_points.emplace_back(x, y); + } + } + } + + return corner_points; +} + +int main(int argc, char* argv[]) { + if (argc != 6) { + std::cout << "usage: " << argv[0] << " " + " \n"; + return -1; + } + + long int window_size = std::stoi(argv[2]); + long hessian_determinant_threshold = std::stol(argv[3]); + double discrimination_constant = std::stod(argv[4]); + + gil::rgb8_image_t input_image; + + gil::read_image(argv[1], input_image, gil::png_tag{}); + + auto input_view = gil::view(input_image); + auto grayscaled = to_grayscale(input_view); + 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()); + + auto x_gradient = gil::view(x_gradient_image); + auto y_gradient = gil::view(y_gradient_image); + calculate_gradients(smoothed, x_gradient, y_gradient); + + gil::gray32f_image_t m11(x_gradient.dimensions()); + gil::gray32f_image_t m12_21(x_gradient.dimensions()); + gil::gray32f_image_t m22(x_gradient.dimensions()); + gil::compute_hessian_entries( + x_gradient, + y_gradient, + gil::view(m11), + gil::view(m12_21), + gil::view(m22) + ); + + gil::gray32f_image_t hessian_determinant(x_gradient.dimensions()); + gil::gray32f_image_t gaussian_kernel(gil::point_t(5, 5)); + gil::generate_gaussian_kernel(gil::view(gaussian_kernel), 0.84089642); + gil::compute_hessian_determinant( + gil::view(m11), + gil::view(m12_21), + gil::view(m22), + gil::view(gaussian_kernel), + discrimination_constant, + gil::view(hessian_determinant) + ); + + auto corner_points = suppress(gil::view(hessian_determinant), hessian_determinant_threshold); + for (auto point: corner_points) { + input_view(point) = gil::rgb8_pixel_t(0, 0, 0); + input_view(point).at(std::integral_constant{}) = 255; + } + gil::write_view(argv[5], input_view, gil::png_tag{}); +} diff --git a/include/boost/gil/image_processing/hessian.hpp b/include/boost/gil/image_processing/hessian.hpp new file mode 100644 index 0000000000..d737fac8cb --- /dev/null +++ b/include/boost/gil/image_processing/hessian.hpp @@ -0,0 +1,77 @@ +#ifndef BOOST_GIL_IMAGE_PROCESSING_HESSIAN_HPP +#define BOOST_GIL_IMAGE_PROCESSING_HESSIAN_HPP + +#include +#include +#include + +namespace boost { namespace gil { + +template +inline void compute_hessian_determinant( + GradientView ddxx, + GradientView dxdy, + GradientView ddyy, + Weights weights, + double discrimination_constant, + OutputView dst) +{ + if (ddxx.dimensions() != ddyy.dimensions() + || ddyy.dimensions() != dxdy.dimensions() + || dxdy.dimensions() != dst.dimensions() + || weights.width() != weights.height() + || weights.width() % 2 != 1) + { + throw std::invalid_argument("dimensions of views are not the same" + " or weights don't have equal width and height" + " or weights' dimensions are not odd"); + } + // Use pixel type of output, as values will be written to output + using pixel_t = typename std::remove_reference()(0, 0))>::type; + + using channel_t = typename std::remove_reference + < + decltype(std::declval().at(std::integral_constant{})) + >::type; + + auto half_size = weights.width() / 2; + for (auto y = half_size; y < dst.height() - half_size; ++y) + { + for (auto x = half_size; x < dst.width() - half_size; ++x) + { + auto ddxx_i = channel_t(); + auto ddyy_i = channel_t(); + auto dxdy_i = channel_t(); + constexpr std::integral_constant chosen_channel{}; + // for (typename OutputView::coord_t w_y = 0; w_y < weights.height(); ++w_y) + // { + // for (typename OutputView::coord_t w_x = 0; w_x < weights.width(); ++w_x) + // { + // ddxx_i += ddxx(x + w_x - half_size, y + w_y - half_size) + // .at(std::integral_constant{}) + // // * weights(w_x, w_y).at(std::integral_constant{}) + // ; + // ddyy_i += ddyy(x + w_x - half_size, y + w_y - half_size) + // .at(std::integral_constant{}) + // // * weights(w_x, w_y).at(std::integral_constant{}) + // ; + // dxdy_i += dxdy(x + w_x - half_size, y + w_y - half_size) + // .at(std::integral_constant{}) + // // * weights(w_x, w_y).at(std::integral_constant{}) + // ; + // } + // } + + ddxx_i = ddxx(x, y).at(chosen_channel); + ddyy_i = ddyy(x, y).at(chosen_channel); + dxdy_i = dxdy(x, y).at(chosen_channel); + + auto determinant = ddxx_i * ddyy_i - dxdy_i * dxdy_i; + dst(x, y).at(std::integral_constant{}) = determinant; + } + } +} + +}} // namespace boost::gil + +#endif diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index 912de8584f..1305680ed1 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -50,6 +50,7 @@ inline double lanczos(double x, std::ptrdiff_t a) return 0; } +<<<<<<< HEAD inline void compute_tensor_entries( boost::gil::gray16_view_t dx, boost::gil::gray16_view_t dy, @@ -68,6 +69,75 @@ inline void compute_tensor_entries( } } +template +inline void compute_hessian_entries( + GradientView dx, + GradientView dy, + OutputView ddxx, + OutputView dxdy, + OutputView ddyy) +{ + using x_coord_t = typename OutputView::x_coord_t; + using y_coord_t = typename OutputView::y_coord_t; + using pixel_t = typename std::remove_reference::type; + using channel_t = typename std::remove_reference< + decltype( + std::declval().at( + std::integral_constant{} + ) + ) + >::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{}; + for (y_coord_t y = 1; y < ddxx.height() - 1; ++y) + { + for (x_coord_t x = 1; x < ddxx.width() - 1; ++x) + { + pixel_t ddxx_i; + boost::gil::static_transform(ddxx_i, ddxx_i, + [](channel_t) { return static_cast(0); }); + pixel_t dxdy_i; + boost::gil::static_transform(dxdy_i, dxdy_i, + [](channel_t) { return static_cast(0); }); + pixel_t ddyy_i; + boost::gil::static_transform(ddyy_i, ddyy_i, + [](channel_t) { return static_cast(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; + ddxx_i.at(std::integral_constant{}) += + dx(adjusted_x, adjusted_y).at(chosen_channel) + * x_kernel[y_filter][x_filter]; + dxdy_i.at(std::integral_constant{}) += + dx(adjusted_x, adjusted_y).at(chosen_channel) + * y_kernel[y_filter][x_filter]; + ddyy_i.at(std::integral_constant{}) += + dy(adjusted_x, adjusted_y).at(chosen_channel) + * y_kernel[y_filter][x_filter]; + } + } + ddxx(x, y) = ddxx_i; + dxdy(x, y) = dxdy_i; + ddyy(x, y) = ddyy_i; + } + } +} + /// \brief Generate mean kernel /// \ingroup ImageProcessingMath /// From 238d124a7e187614688bcf6ced4ceb9c989b1ebb Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 3 Aug 2019 22:57:20 +0000 Subject: [PATCH 02/11] Implement full Hessian corner detector This commit complements last one by summing in a window and applying det - k * trace * trace formula to final Hessian response --- .../boost/gil/image_processing/hessian.hpp | 45 ++++++++++--------- 1 file changed, 23 insertions(+), 22 deletions(-) diff --git a/include/boost/gil/image_processing/hessian.hpp b/include/boost/gil/image_processing/hessian.hpp index d737fac8cb..ae7522730a 100644 --- a/include/boost/gil/image_processing/hessian.hpp +++ b/include/boost/gil/image_processing/hessian.hpp @@ -43,31 +43,32 @@ inline void compute_hessian_determinant( auto ddyy_i = channel_t(); auto dxdy_i = channel_t(); constexpr std::integral_constant chosen_channel{}; - // for (typename OutputView::coord_t w_y = 0; w_y < weights.height(); ++w_y) - // { - // for (typename OutputView::coord_t w_x = 0; w_x < weights.width(); ++w_x) - // { - // ddxx_i += ddxx(x + w_x - half_size, y + w_y - half_size) - // .at(std::integral_constant{}) - // // * weights(w_x, w_y).at(std::integral_constant{}) - // ; - // ddyy_i += ddyy(x + w_x - half_size, y + w_y - half_size) - // .at(std::integral_constant{}) - // // * weights(w_x, w_y).at(std::integral_constant{}) - // ; - // dxdy_i += dxdy(x + w_x - half_size, y + w_y - half_size) - // .at(std::integral_constant{}) - // // * weights(w_x, w_y).at(std::integral_constant{}) - // ; - // } - // } + for (typename OutputView::coord_t w_y = 0; w_y < weights.height(); ++w_y) + { + for (typename OutputView::coord_t w_x = 0; w_x < weights.width(); ++w_x) + { + ddxx_i += ddxx(x + w_x - half_size, y + w_y - half_size) + .at(std::integral_constant{}) + * weights(w_x, w_y).at(std::integral_constant{}) + ; + ddyy_i += ddyy(x + w_x - half_size, y + w_y - half_size) + .at(std::integral_constant{}) + * weights(w_x, w_y).at(std::integral_constant{}) + ; + dxdy_i += dxdy(x + w_x - half_size, y + w_y - half_size) + .at(std::integral_constant{}) + * weights(w_x, w_y).at(std::integral_constant{}) + ; + } + } - ddxx_i = ddxx(x, y).at(chosen_channel); - ddyy_i = ddyy(x, y).at(chosen_channel); - dxdy_i = dxdy(x, y).at(chosen_channel); + // ddxx_i = ddxx(x, y).at(chosen_channel); + // ddyy_i = ddyy(x, y).at(chosen_channel); + // dxdy_i = dxdy(x, y).at(chosen_channel); auto determinant = ddxx_i * ddyy_i - dxdy_i * dxdy_i; - dst(x, y).at(std::integral_constant{}) = determinant; + auto trace = ddyy_i + ddxx_i; + dst(x, y).at(std::integral_constant{}) = determinant - discrimination_constant * trace * trace; } } } From e040b7410e04e9b02e6b1ef9ca5bd6220aabc2c7 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 3 Aug 2019 23:30:48 +0000 Subject: [PATCH 03/11] Add docs and make code align with docs This commit adds docs to new functions, and makes function and variable names align with docs to not confuse readers --- example/hessian.cpp | 8 ++++---- include/boost/gil/image_processing/hessian.hpp | 13 ++++++++----- include/boost/gil/image_processing/numeric.hpp | 7 ++++++- 3 files changed, 18 insertions(+), 10 deletions(-) diff --git a/example/hessian.cpp b/example/hessian.cpp index ae8c8d70ec..f65bf66b28 100644 --- a/example/hessian.cpp +++ b/example/hessian.cpp @@ -231,19 +231,19 @@ int main(int argc, char* argv[]) { gil::view(m22) ); - gil::gray32f_image_t hessian_determinant(x_gradient.dimensions()); + 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); - gil::compute_hessian_determinant( + gil::compute_hessian_response( gil::view(m11), gil::view(m12_21), gil::view(m22), gil::view(gaussian_kernel), discrimination_constant, - gil::view(hessian_determinant) + gil::view(hessian_response) ); - auto corner_points = suppress(gil::view(hessian_determinant), hessian_determinant_threshold); + auto corner_points = suppress(gil::view(hessian_response), hessian_determinant_threshold); for (auto point: corner_points) { input_view(point) = gil::rgb8_pixel_t(0, 0, 0); input_view(point).at(std::integral_constant{}) = 255; diff --git a/include/boost/gil/image_processing/hessian.hpp b/include/boost/gil/image_processing/hessian.hpp index ae7522730a..dca7c536a2 100644 --- a/include/boost/gil/image_processing/hessian.hpp +++ b/include/boost/gil/image_processing/hessian.hpp @@ -7,8 +7,15 @@ namespace boost { namespace gil { +/// \brief Computes Hessian response +/// +/// Computes Hessian response based on computed entries of Hessian matrix, e.g. second order +/// derivates in x and y, and derivatives in both x, y. +/// Weights change perception of surroinding pixels. +/// Discrimination constant is used to differentiate between +/// edges and corners, though additional filtering is strongly advised. template -inline void compute_hessian_determinant( +inline void compute_hessian_response( GradientView ddxx, GradientView dxdy, GradientView ddyy, @@ -62,10 +69,6 @@ inline void compute_hessian_determinant( } } - // ddxx_i = ddxx(x, y).at(chosen_channel); - // ddyy_i = ddyy(x, y).at(chosen_channel); - // dxdy_i = dxdy(x, y).at(chosen_channel); - auto determinant = ddxx_i * ddyy_i - dxdy_i * dxdy_i; auto trace = ddyy_i + ddxx_i; dst(x, y).at(std::integral_constant{}) = determinant - discrimination_constant * trace * trace; diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index 1305680ed1..2df3c5c39e 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -50,7 +50,6 @@ inline double lanczos(double x, std::ptrdiff_t a) return 0; } -<<<<<<< HEAD inline void compute_tensor_entries( boost::gil::gray16_view_t dx, boost::gil::gray16_view_t dy, @@ -69,6 +68,12 @@ inline void compute_tensor_entries( } } +/// \brief Compute xy gradient, and second order x and y gradients +/// \ingroup ImageProcessingMath +/// +/// Hessian matrix is defined as a matrix of partial derivates +/// for 2d case, it is [[ddxx, dxdy], [dxdy, ddyy]. +/// https://en.wikipedia.org/wiki/Hessian_matrix template inline void compute_hessian_entries( GradientView dx, From 688167dee753a95e8fe4f8a5f9a4c046c9a096ba Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 10 Aug 2019 00:49:44 +0000 Subject: [PATCH 04/11] Use determinant as response function A-KAZE uses only determinant in it's response, and since for now Hessian is only a mean to advance A-KAZE implementation, response function is adjusted to use only determinant --- example/hessian.cpp | 8 +++----- include/boost/gil/image_processing/hessian.hpp | 17 +++++------------ 2 files changed, 8 insertions(+), 17 deletions(-) diff --git a/example/hessian.cpp b/example/hessian.cpp index f65bf66b28..f42a7e519e 100644 --- a/example/hessian.cpp +++ b/example/hessian.cpp @@ -194,15 +194,14 @@ std::vector suppress( } int main(int argc, char* argv[]) { - if (argc != 6) { + if (argc != 5) { std::cout << "usage: " << argv[0] << " " - " \n"; + " \n"; return -1; } long int window_size = std::stoi(argv[2]); long hessian_determinant_threshold = std::stol(argv[3]); - double discrimination_constant = std::stod(argv[4]); gil::rgb8_image_t input_image; @@ -239,7 +238,6 @@ int main(int argc, char* argv[]) { gil::view(m12_21), gil::view(m22), gil::view(gaussian_kernel), - discrimination_constant, gil::view(hessian_response) ); @@ -248,5 +246,5 @@ int main(int argc, char* argv[]) { input_view(point) = gil::rgb8_pixel_t(0, 0, 0); input_view(point).at(std::integral_constant{}) = 255; } - gil::write_view(argv[5], input_view, gil::png_tag{}); + gil::write_view(argv[4], input_view, gil::png_tag{}); } diff --git a/include/boost/gil/image_processing/hessian.hpp b/include/boost/gil/image_processing/hessian.hpp index dca7c536a2..2dcbb02f55 100644 --- a/include/boost/gil/image_processing/hessian.hpp +++ b/include/boost/gil/image_processing/hessian.hpp @@ -12,15 +12,13 @@ namespace boost { namespace gil { /// Computes Hessian response based on computed entries of Hessian matrix, e.g. second order /// derivates in x and y, and derivatives in both x, y. /// Weights change perception of surroinding pixels. -/// Discrimination constant is used to differentiate between -/// edges and corners, though additional filtering is strongly advised. +/// Additional filtering is strongly advised. template inline void compute_hessian_response( GradientView ddxx, GradientView dxdy, GradientView ddyy, Weights weights, - double discrimination_constant, OutputView dst) { if (ddxx.dimensions() != ddyy.dimensions() @@ -56,22 +54,17 @@ inline void compute_hessian_response( { ddxx_i += ddxx(x + w_x - half_size, y + w_y - half_size) .at(std::integral_constant{}) - * weights(w_x, w_y).at(std::integral_constant{}) - ; + * weights(w_x, w_y).at(std::integral_constant{}); ddyy_i += ddyy(x + w_x - half_size, y + w_y - half_size) .at(std::integral_constant{}) - * weights(w_x, w_y).at(std::integral_constant{}) - ; + * weights(w_x, w_y).at(std::integral_constant{}); dxdy_i += dxdy(x + w_x - half_size, y + w_y - half_size) .at(std::integral_constant{}) - * weights(w_x, w_y).at(std::integral_constant{}) - ; + * weights(w_x, w_y).at(std::integral_constant{}); } } - auto determinant = ddxx_i * ddyy_i - dxdy_i * dxdy_i; - auto trace = ddyy_i + ddxx_i; - dst(x, y).at(std::integral_constant{}) = determinant - discrimination_constant * trace * trace; + dst(x, y).at(std::integral_constant{}) = determinant; } } } From 6f439d9a0200210b9564803bd9ecb7bd8ece3077 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 10 Aug 2019 01:15:54 +0000 Subject: [PATCH 05/11] Create simple test for Hessian detector This commit adds an s at the end of function name to make it uniform with another detector, and adds a simple test for sanity check --- .../boost/gil/image_processing/hessian.hpp | 2 +- test/core/image_processing/CMakeLists.txt | 3 +- test/core/image_processing/Jamfile | 1 + test/core/image_processing/hessian.cpp | 72 +++++++++++++++++++ 4 files changed, 76 insertions(+), 2 deletions(-) create mode 100644 test/core/image_processing/hessian.cpp diff --git a/include/boost/gil/image_processing/hessian.hpp b/include/boost/gil/image_processing/hessian.hpp index 2dcbb02f55..ebc407c604 100644 --- a/include/boost/gil/image_processing/hessian.hpp +++ b/include/boost/gil/image_processing/hessian.hpp @@ -14,7 +14,7 @@ namespace boost { namespace gil { /// Weights change perception of surroinding pixels. /// Additional filtering is strongly advised. template -inline void compute_hessian_response( +inline void compute_hessian_responses( GradientView ddxx, GradientView dxdy, GradientView ddyy, diff --git a/test/core/image_processing/CMakeLists.txt b/test/core/image_processing/CMakeLists.txt index 25eda96536..9e75ae4d2e 100644 --- a/test/core/image_processing/CMakeLists.txt +++ b/test/core/image_processing/CMakeLists.txt @@ -30,7 +30,8 @@ endforeach() foreach(_name lanczos_scaling simple_kernels - harris) + harris + hessian) set(_test t_core_image_processing_${_name}) set(_target test_core_image_processing_${_name}) diff --git a/test/core/image_processing/Jamfile b/test/core/image_processing/Jamfile index 19e2348c62..f1c32fc0ac 100644 --- a/test/core/image_processing/Jamfile +++ b/test/core/image_processing/Jamfile @@ -20,3 +20,4 @@ run threshold_otsu.cpp ; run lanczos_scaling.cpp ; run simple_kernels.cpp ; run harris.cpp ; +run hessian.cpp ; diff --git a/test/core/image_processing/hessian.cpp b/test/core/image_processing/hessian.cpp new file mode 100644 index 0000000000..dd96e409f8 --- /dev/null +++ b/test/core/image_processing/hessian.cpp @@ -0,0 +1,72 @@ +// +// Copyright 2019 Olzhas Zhumabek +// +// Use, modification and distribution are subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// +#include +#include +#include +#include +#include + +namespace gil = boost::gil; + +bool are_equal(gil::gray32f_view_t expected, gil::gray32f_view_t actual) { + if (expected.dimensions() != actual.dimensions()) + return false; + + for (long int y = 0; y < expected.height(); ++y) + { + for (long int x = 0; x < expected.width(); ++x) + { + if (expected(x, y) != actual(x, y)) + { + return false; + } + } + } + + return true; +} + +void test_blank_image() +{ + const gil::point_t dimensions(20, 20); + gil::gray16_image_t dx(dimensions, gil::gray16_pixel_t(0), 0); + gil::gray16_image_t dy(dimensions, gil::gray16_pixel_t(0), 0); + + gil::gray32f_image_t m11(dimensions); + gil::gray32f_image_t m12_21(dimensions); + gil::gray32f_image_t m22(dimensions); + gil::gray32f_image_t expected(dimensions, gil::gray32f_pixel_t(0), 0); + gil::compute_hessian_entries( + gil::view(dx), + gil::view(dy), + gil::view(m11), + gil::view(m12_21), + gil::view(m22) + ); + BOOST_TEST(are_equal(gil::view(expected), gil::view(m11))); + BOOST_TEST(are_equal(gil::view(expected), gil::view(m12_21))); + BOOST_TEST(are_equal(gil::view(expected), gil::view(m22))); + + gil::gray32f_image_t hessian_response(dimensions, gil::gray32f_pixel_t(0), 0); + gil::gray32f_image_t unnormalized_mean(gil::point_t(5, 5)); + gil::generate_unnormalized_mean(gil::view(unnormalized_mean)); + gil::compute_hessian_responses( + gil::view(m11), + gil::view(m12_21), + gil::view(m22), + gil::view(unnormalized_mean), + gil::view(hessian_response) + ); + BOOST_TEST(are_equal(gil::view(expected), gil::view(hessian_response))); +} + +int main(int argc, char* argv[]) +{ + test_blank_image(); + return boost::report_errors(); +} From f96df0daeea74ef6b4ef5e211d8e76ab75333e13 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 12 Aug 2019 03:06:24 +0000 Subject: [PATCH 06/11] Improve documentations for d params dx and similar naming seems to be confusing, improved documentation to explain the naming convention --- include/boost/gil/image_processing/hessian.hpp | 2 ++ include/boost/gil/image_processing/numeric.hpp | 3 +++ 2 files changed, 5 insertions(+) diff --git a/include/boost/gil/image_processing/hessian.hpp b/include/boost/gil/image_processing/hessian.hpp index ebc407c604..5d36599f3d 100644 --- a/include/boost/gil/image_processing/hessian.hpp +++ b/include/boost/gil/image_processing/hessian.hpp @@ -11,6 +11,8 @@ namespace boost { namespace gil { /// /// Computes Hessian response based on computed entries of Hessian matrix, e.g. second order /// derivates in x and y, and derivatives in both x, y. +/// d stands for derivative, and x or y stand for derivative direction. For example, +/// ddxx means taking two derivatives (gradients) in horizontal direction. /// Weights change perception of surroinding pixels. /// Additional filtering is strongly advised. template diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index 2df3c5c39e..92c0b7df3e 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -73,6 +73,9 @@ inline void compute_tensor_entries( /// /// Hessian matrix is defined as a matrix of partial derivates /// for 2d case, it is [[ddxx, dxdy], [dxdy, ddyy]. +/// d stands for derivative, and x or y stand for direction. +/// For example, dx stands for derivative (gradient) in horizontal +/// direction, and ddxx means second order derivative in horizon direction /// https://en.wikipedia.org/wiki/Hessian_matrix template inline void compute_hessian_entries( From 8683b6345d2de7d8b82e35779673c1ec3991eef1 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 12 Aug 2019 03:10:38 +0000 Subject: [PATCH 07/11] Address minor comments about style --- example/hessian.cpp | 24 ++++++++++++------- .../boost/gil/image_processing/numeric.hpp | 6 ++--- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/example/hessian.cpp b/example/hessian.cpp index f42a7e519e..c5095d8e8e 100644 --- a/example/hessian.cpp +++ b/example/hessian.cpp @@ -21,8 +21,10 @@ gil::gray8_image_t to_grayscale(gil::rgb8_view_t original) gil::gray8_image_t output_image(original.dimensions()); auto output = gil::view(output_image); constexpr double max_channel_intensity = (std::numeric_limits::max)(); - for (long int y = 0; y < original.height(); ++y) { - for (long int x = 0; x < original.width(); ++x) { + for (long int y = 0; y < original.height(); ++y) + { + for (long int x = 0; x < original.width(); ++x) + { // scale the values into range [0, 1] and calculate linear intensity double red_intensity = original(x, y).at(std::integral_constant{}) / max_channel_intensity; @@ -36,9 +38,11 @@ gil::gray8_image_t to_grayscale(gil::rgb8_view_t original) // perform gamma adjustment double gamma_compressed_luminosity = 0; - if (linear_luminosity < 0.0031308) { + if (linear_luminosity < 0.0031308) + { gamma_compressed_luminosity = linear_luminosity * 12.92; - } else { + } else + { gamma_compressed_luminosity = 1.055 * std::pow(linear_luminosity, 1 / 2.4) - 0.055; } @@ -70,12 +74,15 @@ void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_ for (long x = 0; x < width; ++x) { for (long y = 0; y < height; ++y) { double intensity = 0.0; - for (size_t filter_y = 0; filter_y < filterHeight; ++filter_y) { - for (size_t filter_x = 0; filter_x < filterWidth; ++filter_x) { + for (size_t filter_y = 0; filter_y < filterHeight; ++filter_y) + { + for (size_t filter_x = 0; filter_x < filterWidth; ++filter_x) + { int image_x = x - filterWidth / 2 + filter_x; int image_y = y - filterHeight / 2 + filter_y; if (image_x >= input_view.width() || image_x < 0 - || image_y >= input_view.height() || image_y < 0) { + || image_y >= input_view.height() || image_y < 0) + { continue; } auto& pixel = input_view(image_x, image_y); @@ -194,7 +201,8 @@ std::vector suppress( } int main(int argc, char* argv[]) { - if (argc != 5) { + if (argc != 5) + { std::cout << "usage: " << argv[0] << " " " \n"; return -1; diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index 92c0b7df3e..d0ae0ac020 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -114,13 +114,13 @@ inline void compute_hessian_entries( for (x_coord_t x = 1; x < ddxx.width() - 1; ++x) { pixel_t ddxx_i; - boost::gil::static_transform(ddxx_i, ddxx_i, + static_transform(ddxx_i, ddxx_i, [](channel_t) { return static_cast(0); }); pixel_t dxdy_i; - boost::gil::static_transform(dxdy_i, dxdy_i, + static_transform(dxdy_i, dxdy_i, [](channel_t) { return static_cast(0); }); pixel_t ddyy_i; - boost::gil::static_transform(ddyy_i, ddyy_i, + static_transform(ddyy_i, ddyy_i, [](channel_t) { return static_cast(0); }); for (y_coord_t y_filter = 0; y_filter < 2; ++y_filter) { From 8ee318aae3effe2f18e50dc3ea9096f0d4bee591 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 12 Aug 2019 03:15:46 +0000 Subject: [PATCH 08/11] Address type based issues Mostly changes to constness and integral types, with small cosmetic changes mixed in --- example/hessian.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/example/hessian.cpp b/example/hessian.cpp index c5095d8e8e..5e385941fb 100644 --- a/example/hessian.cpp +++ b/example/hessian.cpp @@ -56,9 +56,9 @@ gil::gray8_image_t to_grayscale(gil::rgb8_view_t original) void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_view) { - constexpr static auto filterHeight = 5ull; - constexpr static auto filterWidth = 5ull; - constexpr static double filter[filterHeight][filterWidth] = + constexpr static auto filter_height = 5ull; + constexpr static auto filter_width = 5ull; + constexpr static double filter[filter_height][filter_width] = { 2, 4, 6, 4, 2, 4, 9, 12, 9, 4, @@ -71,21 +71,21 @@ void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_ const auto height = input_view.height(); const auto width = input_view.width(); - for (long x = 0; x < width; ++x) { - for (long y = 0; y < height; ++y) { + for (std::ptrdiff_t x = 0; x < width; ++x) { + for (std::ptrdiff_t y = 0; y < height; ++y) { double intensity = 0.0; - for (size_t filter_y = 0; filter_y < filterHeight; ++filter_y) + for (std::ptrdiff_t filter_y = 0; filter_y < filter_height; ++filter_y) { - for (size_t filter_x = 0; filter_x < filterWidth; ++filter_x) + for (std::ptrdiff_t filter_x = 0; filter_x < filter_width; ++filter_x) { - int image_x = x - filterWidth / 2 + filter_x; - int image_y = y - filterHeight / 2 + filter_y; + int image_x = x - filter_width / 2 + filter_x; + int image_y = y - filter_height / 2 + filter_y; if (image_x >= input_view.width() || image_x < 0 || image_y >= input_view.height() || image_y < 0) { continue; } - auto& pixel = input_view(image_x, image_y); + const auto& pixel = input_view(image_x, image_y); intensity += pixel.at(std::integral_constant{}) * filter[filter_y][filter_x]; } From 75856655783349c158630db31e77a42707c057ce Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 12 Aug 2019 03:23:17 +0000 Subject: [PATCH 09/11] Fix typo and address review comment Fixes typo in call for Hessian and addresses a review comment about replacing multiple exact indexing calls to one with a reference --- example/hessian.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/example/hessian.cpp b/example/hessian.cpp index 5e385941fb..826f1fb675 100644 --- a/example/hessian.cpp +++ b/example/hessian.cpp @@ -26,11 +26,12 @@ gil::gray8_image_t to_grayscale(gil::rgb8_view_t original) for (long int x = 0; x < original.width(); ++x) { // scale the values into range [0, 1] and calculate linear intensity - double red_intensity = original(x, y).at(std::integral_constant{}) + auto& p = original(x, y); + double red_intensity = p.at(std::integral_constant{}) / max_channel_intensity; - double green_intensity = original(x, y).at(std::integral_constant{}) + double green_intensity = p.at(std::integral_constant{}) / max_channel_intensity; - double blue_intensity = original(x, y).at(std::integral_constant{}) + double blue_intensity = p.at(std::integral_constant{}) / max_channel_intensity; auto linear_luminosity = 0.2126 * red_intensity + 0.7152 * green_intensity @@ -241,7 +242,7 @@ 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); - gil::compute_hessian_response( + gil::compute_hessian_responses( gil::view(m11), gil::view(m12_21), gil::view(m22), From 9370e3d49cb5d66ec8de62a9f81531475195d4ec Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 12 Aug 2019 03:27:41 +0000 Subject: [PATCH 10/11] Reorder includes in Hessian example --- example/hessian.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/hessian.cpp b/example/hessian.cpp index 826f1fb675..cc752f2b66 100644 --- a/example/hessian.cpp +++ b/example/hessian.cpp @@ -1,8 +1,8 @@ #include #include -#include #include #include +#include #include #include #include From b232d89d5b258ac6ddb6fea26afe4214751f4910 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 20 Aug 2019 06:08:13 +0300 Subject: [PATCH 11/11] Address review comments Add literature reference to luminosity computation, and perform some cosmetic changes. --- example/hessian.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/example/hessian.cpp b/example/hessian.cpp index cc752f2b66..aca22f4ce6 100644 --- a/example/hessian.cpp +++ b/example/hessian.cpp @@ -15,7 +15,10 @@ namespace gil = boost::gil; // when converted to grayscale, // which was previously observed on // canny edge detector for test input -// used for this example +// used for this example. +// the algorithm here follows sRGB gamma definition +// taken from here (luminance calculation): +// https://en.wikipedia.org/wiki/Grayscale gil::gray8_image_t to_grayscale(gil::rgb8_view_t original) { gil::gray8_image_t output_image(original.dimensions()); @@ -72,8 +75,10 @@ void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_ const auto height = input_view.height(); const auto width = input_view.width(); - for (std::ptrdiff_t x = 0; x < width; ++x) { - for (std::ptrdiff_t y = 0; y < height; ++y) { + for (std::ptrdiff_t x = 0; x < width; ++x) + { + for (std::ptrdiff_t y = 0; y < height; ++y) + { double intensity = 0.0; for (std::ptrdiff_t filter_y = 0; filter_y < filter_height; ++filter_y) { @@ -81,8 +86,8 @@ void apply_gaussian_blur(gil::gray8_view_t input_view, gil::gray8_view_t output_ { int image_x = x - filter_width / 2 + filter_x; int image_y = y - filter_height / 2 + filter_y; - if (image_x >= input_view.width() || image_x < 0 - || image_y >= input_view.height() || image_y < 0) + if (image_x >= input_view.width() || image_x < 0 || + image_y >= input_view.height() || image_y < 0) { continue; }