From 5a89ea301ad4015c94bfbaa971ad5293859cddae Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 2 Jun 2020 22:55:33 +0600 Subject: [PATCH 01/48] Starting point --- example/anisotropic_diffusion.cpp | 160 ++++++++++++++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 example/anisotropic_diffusion.cpp diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp new file mode 100644 index 0000000000..355427ade9 --- /dev/null +++ b/example/anisotropic_diffusion.cpp @@ -0,0 +1,160 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace gil = boost::gil; + +using gray_channel = std::integral_constant; + +namespace boost{namespace gil{ +using gray64f_view_t = gil::image_view>>; +}} + +enum direction { + north = 0, + south = 1, + west = 2, + east = 3, + north_east = 4, + south_east = 5, + south_west = 6, + north_west = 7 +}; + +void compute_nabla(gil::gray64f_view_t view, const std::vector& nabla) { + for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) + { + for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) + { + nabla[north](x, y) = view(x, y - 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); + nabla[south](x, y) = view(x, y + 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); + nabla[west](x, y) = view(x - 1, y).at(gray_channel{}) - view(x, y).at(gray_channel{}); + nabla[east](x, y) = view(x + 1, y).at(gray_channel{}) - view(x, y).at(gray_channel{}); + + nabla[north_east](x, y) = view(x + 1, y - 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); + nabla[south_east](x, y) = view(x + 1, y + 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); + nabla[south_west](x, y) = view(x - 1, y + 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); + nabla[north_west](x, y) = view(x - 1, y - 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); + } + } +} + +void calculate_diffusvity(std::vector nablas, double kappa, const std::vector diffusivities) +{ + BOOST_ASSERT(nablas.size() == diffusivities.size()); + for (std::size_t i = 0; i < nablas.size(); ++i) { + gil::transform_pixels(nablas[i], diffusivities[i], [kappa](gil::gray32f_pixel_t p) + { + float value = p.at(gray_channel{}) / kappa; + float result = std::exp(-value * value); + return result; + }); + } +} + +template +void print_images(const std::vector& views, const std::string& prefix = "./nabla") +{ + unsigned int counter = 0; + for (const auto& view: views) + { + gil::write_view(prefix + std::to_string(counter) + ".png", gil::color_converted_view(view), gil::png_tag{}); + ++counter; + } +} + +void anisotropic_diffusion(gil::gray64f_view_t input, unsigned int num_iter, double delta_t, double kappa, gil::gray64f_view_t output) +{ + gil::copy_pixels(input, output); + + for (unsigned int i = 0; i < num_iter; ++i) { + std::vector nabla_images(8, gil::gray32f_image_t(input.dimensions())); + std::vector nabla; + std::transform(nabla_images.begin(), nabla_images.end(), + std::back_inserter(nabla), [](gil::gray32f_image_t& i) + { + return gil::view(i); + }); + + std::vector diffusivity_images(8, gil::gray32f_image_t(input.dimensions())); + std::vector diffusivity; + std::transform(diffusivity_images.begin(), diffusivity_images.end(), + std::back_inserter(diffusivity), [](gil::gray32f_image_t& i) + { + return gil::view(i); + }); + + compute_nabla(output, nabla); + calculate_diffusvity(nabla, kappa, diffusivity); + + // print_images(nabla); + // print_images(diffusivity, "diffusivity"); + + float half = float(1.0f / 2); + for (std::ptrdiff_t y = 0; y < output.height(); ++y) + { + for (std::ptrdiff_t x = 0; x < output.width(); ++ x) + { + float delta = delta_t * ( + diffusivity[north](x, y).at(gray_channel{}) * nabla[north](x, y).at(gray_channel{}) + diffusivity[south](x, y).at(gray_channel{}) * nabla[south](x, y).at(gray_channel{}) + + diffusivity[west](x, y).at(gray_channel{}) * nabla[west](x, y).at(gray_channel{}) + diffusivity[east](x, y).at(gray_channel{}) * nabla[east](x, y).at(gray_channel{}) + + half * diffusivity[north_east](x, y).at(gray_channel{}) * nabla[north_east](x, y).at(gray_channel{}) + half * diffusivity[south_east](x, y).at(gray_channel{}) * nabla[south_east](x, y).at(gray_channel{}) + + half * diffusivity[south_west](x, y).at(gray_channel{}) * nabla[south_west](x, y).at(gray_channel{}) + half * diffusivity[north_west](x, y).at(gray_channel{}) * nabla[north_west](x, y).at(gray_channel{}) + ); + output(x, y) = output(x, y).at(gray_channel{}) + delta; + } + } + } +} + +#include + +int main(int argc, char* argv[]) +{ + if (argc != 5) + { + std::cerr << "usage: " << argv[0] << " \n"; + return -1; + } + std::string input_path = argv[1]; + std::string output_path = argv[2]; + + unsigned int iteration_count = static_cast(std::stoul(argv[3])); + unsigned int kappa = static_cast(std::stoul(argv[4])); + + gil::gray8_image_t input; + gil::read_image(input_path, input, gil::png_tag{}); + auto input_view = gil::view(input); + + gil::gray64f_image_t gray(input.dimensions()); + auto gray_view = gil::view(gray); + + gil::transform_pixels(input_view, gray_view, + [](const gil::gray8_pixel_t& p) + { + return p.at(gray_channel{}); + }); + gil::gray64f_image_t output(gray.dimensions()); + auto output_view = gil::view(output); + + anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); + + gil::gray8_image_t true_output(output.dimensions()); + gil::transform_pixels(output_view, gil::view(true_output), [](gil::gray64f_pixel_t p) + { + return p.at(gray_channel{}); + }); + + gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); +} \ No newline at end of file From e9282ac8c233858f00fb0e734ff5705c965c9764 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Wed, 3 Jun 2020 00:41:14 +0600 Subject: [PATCH 02/48] Working example version --- example/anisotropic_diffusion.cpp | 87 +++++++++++++++++++------------ 1 file changed, 55 insertions(+), 32 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index 355427ade9..9266f968c0 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -20,7 +21,7 @@ using gray_channel = std::integral_constant; namespace boost{namespace gil{ using gray64f_view_t = gil::image_view>>; }} - +namespace boost{ namespace gil{ enum direction { north = 0, south = 1, @@ -32,33 +33,48 @@ enum direction { north_west = 7 }; -void compute_nabla(gil::gray64f_view_t view, const std::vector& nabla) { +template +void compute_nabla(InputView view, const std::vector& nabla) { + constexpr auto input_num_channels = num_channels{}; + static_assert(num_channels{} == input_num_channels); for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) { for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) { - nabla[north](x, y) = view(x, y - 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); - nabla[south](x, y) = view(x, y + 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); - nabla[west](x, y) = view(x - 1, y).at(gray_channel{}) - view(x, y).at(gray_channel{}); - nabla[east](x, y) = view(x + 1, y).at(gray_channel{}) - view(x, y).at(gray_channel{}); - - nabla[north_east](x, y) = view(x + 1, y - 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); - nabla[south_east](x, y) = view(x + 1, y + 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); - nabla[south_west](x, y) = view(x - 1, y + 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); - nabla[north_west](x, y) = view(x - 1, y - 1).at(gray_channel{}) - view(x, y).at(gray_channel{}); + for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; ++channel_index) + { + nabla[north](x, y) = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[south](x, y) = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[west](x, y) = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; + nabla[east](x, y) = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; + + nabla[north_east](x, y) = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[south_east](x, y) = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[south_west](x, y) = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[north_west](x, y) = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; + } } } } -void calculate_diffusvity(std::vector nablas, double kappa, const std::vector diffusivities) +template +void calculate_diffusvity(std::vector nablas, double kappa, const std::vector diffusivities) { + using pixel_type = typename View::value_type; + using channel_type = typename channel_type::type; BOOST_ASSERT(nablas.size() == diffusivities.size()); for (std::size_t i = 0; i < nablas.size(); ++i) { - gil::transform_pixels(nablas[i], diffusivities[i], [kappa](gil::gray32f_pixel_t p) + gil::transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) { - float value = p.at(gray_channel{}) / kappa; - float result = std::exp(-value * value); - return result; + auto op = [kappa](channel_type value) + { + value /= kappa; + float result = std::exp(-value * value); + return result; + }; + pixel_type result_pixel; + static_transform(p, result_pixel, op); + return result_pixel; }); } } @@ -74,25 +90,28 @@ void print_images(const std::vector& views, const std::string& prefix } } -void anisotropic_diffusion(gil::gray64f_view_t input, unsigned int num_iter, double delta_t, double kappa, gil::gray64f_view_t output) +template +void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, OutputView output) { gil::copy_pixels(input, output); + using element_type = typename OutputView::value_type; + using computation_image_type = image; for (unsigned int i = 0; i < num_iter; ++i) { - std::vector nabla_images(8, gil::gray32f_image_t(input.dimensions())); - std::vector nabla; + std::vector nabla_images(8, computation_image_type(input.dimensions())); + std::vector nabla; std::transform(nabla_images.begin(), nabla_images.end(), - std::back_inserter(nabla), [](gil::gray32f_image_t& i) + std::back_inserter(nabla), [](computation_image_type& img) { - return gil::view(i); + return gil::view(img); }); - std::vector diffusivity_images(8, gil::gray32f_image_t(input.dimensions())); - std::vector diffusivity; + std::vector diffusivity_images(8, computation_image_type(input.dimensions())); + std::vector diffusivity; std::transform(diffusivity_images.begin(), diffusivity_images.end(), - std::back_inserter(diffusivity), [](gil::gray32f_image_t& i) + std::back_inserter(diffusivity), [](computation_image_type& img) { - return gil::view(i); + return gil::view(img); }); compute_nabla(output, nabla); @@ -102,21 +121,25 @@ void anisotropic_diffusion(gil::gray64f_view_t input, unsigned int num_iter, dou // print_images(diffusivity, "diffusivity"); float half = float(1.0f / 2); + constexpr auto channel_count = num_channels{}; for (std::ptrdiff_t y = 0; y < output.height(); ++y) { for (std::ptrdiff_t x = 0; x < output.width(); ++ x) { - float delta = delta_t * ( - diffusivity[north](x, y).at(gray_channel{}) * nabla[north](x, y).at(gray_channel{}) + diffusivity[south](x, y).at(gray_channel{}) * nabla[south](x, y).at(gray_channel{}) - + diffusivity[west](x, y).at(gray_channel{}) * nabla[west](x, y).at(gray_channel{}) + diffusivity[east](x, y).at(gray_channel{}) * nabla[east](x, y).at(gray_channel{}) - + half * diffusivity[north_east](x, y).at(gray_channel{}) * nabla[north_east](x, y).at(gray_channel{}) + half * diffusivity[south_east](x, y).at(gray_channel{}) * nabla[south_east](x, y).at(gray_channel{}) - + half * diffusivity[south_west](x, y).at(gray_channel{}) * nabla[south_west](x, y).at(gray_channel{}) + half * diffusivity[north_west](x, y).at(gray_channel{}) * nabla[north_west](x, y).at(gray_channel{}) - ); - output(x, y) = output(x, y).at(gray_channel{}) + delta; + for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; ++channel_index) { + float delta = delta_t * ( + diffusivity[north](x, y)[channel_index] * nabla[north](x, y)[channel_index] + diffusivity[south](x, y)[channel_index] * nabla[south](x, y)[channel_index] + + diffusivity[west](x, y)[channel_index] * nabla[west](x, y)[channel_index] + diffusivity[east](x, y)[channel_index] * nabla[east](x, y)[channel_index] + + half * diffusivity[north_east](x, y)[channel_index] * nabla[north_east](x, y)[channel_index] + half * diffusivity[south_east](x, y)[channel_index] * nabla[south_east](x, y)[channel_index] + + half * diffusivity[south_west](x, y)[channel_index] * nabla[south_west](x, y)[channel_index] + half * diffusivity[north_west](x, y)[channel_index] * nabla[north_west](x, y)[channel_index] + ); + output(x, y)[channel_index] = output(x, y)[channel_index] + delta; + } } } } } +}} #include From 0c9bec250e25275ef9e1159cd496e18847a91e07 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Wed, 3 Jun 2020 00:53:39 +0600 Subject: [PATCH 03/48] Stylistic changes --- example/anisotropic_diffusion.cpp | 48 +++++++++++-------------------- 1 file changed, 16 insertions(+), 32 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index 9266f968c0..81c82a83ca 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -16,13 +16,11 @@ namespace gil = boost::gil; -using gray_channel = std::integral_constant; - namespace boost{namespace gil{ using gray64f_view_t = gil::image_view>>; }} namespace boost{ namespace gil{ -enum direction { +enum class direction: std::size_t { north = 0, south = 1, west = 2, @@ -43,15 +41,15 @@ void compute_nabla(InputView view, const std::vector& nabla) { { for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; ++channel_index) { - nabla[north](x, y) = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[south](x, y) = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[west](x, y) = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; - nabla[east](x, y) = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; - - nabla[north_east](x, y) = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[south_east](x, y) = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[south_west](x, y) = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[north_west](x, y) = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::north](x, y) = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south](x, y) = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::west](x, y) = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::east](x, y) = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; + + nabla[(std::size_t)direction::north_east](x, y) = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_east](x, y) = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_west](x, y) = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::north_west](x, y) = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; } } } @@ -79,17 +77,6 @@ void calculate_diffusvity(std::vector nablas, double kappa, const std::vec } } -template -void print_images(const std::vector& views, const std::string& prefix = "./nabla") -{ - unsigned int counter = 0; - for (const auto& view: views) - { - gil::write_view(prefix + std::to_string(counter) + ".png", gil::color_converted_view(view), gil::png_tag{}); - ++counter; - } -} - template void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, OutputView output) { @@ -117,9 +104,6 @@ void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_ compute_nabla(output, nabla); calculate_diffusvity(nabla, kappa, diffusivity); - // print_images(nabla); - // print_images(diffusivity, "diffusivity"); - float half = float(1.0f / 2); constexpr auto channel_count = num_channels{}; for (std::ptrdiff_t y = 0; y < output.height(); ++y) @@ -128,10 +112,10 @@ void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_ { for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; ++channel_index) { float delta = delta_t * ( - diffusivity[north](x, y)[channel_index] * nabla[north](x, y)[channel_index] + diffusivity[south](x, y)[channel_index] * nabla[south](x, y)[channel_index] - + diffusivity[west](x, y)[channel_index] * nabla[west](x, y)[channel_index] + diffusivity[east](x, y)[channel_index] * nabla[east](x, y)[channel_index] - + half * diffusivity[north_east](x, y)[channel_index] * nabla[north_east](x, y)[channel_index] + half * diffusivity[south_east](x, y)[channel_index] * nabla[south_east](x, y)[channel_index] - + half * diffusivity[south_west](x, y)[channel_index] * nabla[south_west](x, y)[channel_index] + half * diffusivity[north_west](x, y)[channel_index] * nabla[north_west](x, y)[channel_index] + diffusivity[(std::size_t)direction::north](x, y)[channel_index] * nabla[(std::size_t)direction::north](x, y)[channel_index] + diffusivity[(std::size_t)direction::south](x, y)[channel_index] * nabla[(std::size_t)direction::south](x, y)[channel_index] + + diffusivity[(std::size_t)direction::west](x, y)[channel_index] * nabla[(std::size_t)direction::west](x, y)[channel_index] + diffusivity[(std::size_t)direction::east](x, y)[channel_index] * nabla[(std::size_t)direction::east](x, y)[channel_index] + + half * diffusivity[(std::size_t)direction::north_east](x, y)[channel_index] * nabla[(std::size_t)direction::north_east](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::south_east](x, y)[channel_index] * nabla[(std::size_t)direction::south_east](x, y)[channel_index] + + half * diffusivity[(std::size_t)direction::south_west](x, y)[channel_index] * nabla[(std::size_t)direction::south_west](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::north_west](x, y)[channel_index] * nabla[(std::size_t)direction::north_west](x, y)[channel_index] ); output(x, y)[channel_index] = output(x, y)[channel_index] + delta; } @@ -166,7 +150,7 @@ int main(int argc, char* argv[]) gil::transform_pixels(input_view, gray_view, [](const gil::gray8_pixel_t& p) { - return p.at(gray_channel{}); + return p[0]; }); gil::gray64f_image_t output(gray.dimensions()); auto output_view = gil::view(output); @@ -176,7 +160,7 @@ int main(int argc, char* argv[]) gil::gray8_image_t true_output(output.dimensions()); gil::transform_pixels(output_view, gil::view(true_output), [](gil::gray64f_pixel_t p) { - return p.at(gray_channel{}); + return p[0]; }); gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); From 902fbb02c99e5ef12be9080385f4ebd7290c43fe Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Wed, 3 Jun 2020 01:04:32 +0600 Subject: [PATCH 04/48] Reduce amount of warnings to zero The compilation of example will still produce some warnings, but they are unrelated to the new code --- example/anisotropic_diffusion.cpp | 114 +----------------- .../boost/gil/image_processing/numeric.hpp | 107 ++++++++++++++++ 2 files changed, 110 insertions(+), 111 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index 81c82a83ca..afba7ede3e 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -2,12 +2,9 @@ #include #include #include -#include -#include -#include #include #include -#include +#include #include #include @@ -19,111 +16,6 @@ namespace gil = boost::gil; namespace boost{namespace gil{ using gray64f_view_t = gil::image_view>>; }} -namespace boost{ namespace gil{ -enum class direction: std::size_t { - north = 0, - south = 1, - west = 2, - east = 3, - north_east = 4, - south_east = 5, - south_west = 6, - north_west = 7 -}; - -template -void compute_nabla(InputView view, const std::vector& nabla) { - constexpr auto input_num_channels = num_channels{}; - static_assert(num_channels{} == input_num_channels); - for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) - { - for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) - { - for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; ++channel_index) - { - nabla[(std::size_t)direction::north](x, y) = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south](x, y) = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::west](x, y) = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::east](x, y) = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; - - nabla[(std::size_t)direction::north_east](x, y) = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_east](x, y) = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_west](x, y) = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::north_west](x, y) = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; - } - } - } -} - -template -void calculate_diffusvity(std::vector nablas, double kappa, const std::vector diffusivities) -{ - using pixel_type = typename View::value_type; - using channel_type = typename channel_type::type; - BOOST_ASSERT(nablas.size() == diffusivities.size()); - for (std::size_t i = 0; i < nablas.size(); ++i) { - gil::transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) - { - auto op = [kappa](channel_type value) - { - value /= kappa; - float result = std::exp(-value * value); - return result; - }; - pixel_type result_pixel; - static_transform(p, result_pixel, op); - return result_pixel; - }); - } -} - -template -void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, OutputView output) -{ - gil::copy_pixels(input, output); - using element_type = typename OutputView::value_type; - using computation_image_type = image; - - for (unsigned int i = 0; i < num_iter; ++i) { - std::vector nabla_images(8, computation_image_type(input.dimensions())); - std::vector nabla; - std::transform(nabla_images.begin(), nabla_images.end(), - std::back_inserter(nabla), [](computation_image_type& img) - { - return gil::view(img); - }); - - std::vector diffusivity_images(8, computation_image_type(input.dimensions())); - std::vector diffusivity; - std::transform(diffusivity_images.begin(), diffusivity_images.end(), - std::back_inserter(diffusivity), [](computation_image_type& img) - { - return gil::view(img); - }); - - compute_nabla(output, nabla); - calculate_diffusvity(nabla, kappa, diffusivity); - - float half = float(1.0f / 2); - constexpr auto channel_count = num_channels{}; - for (std::ptrdiff_t y = 0; y < output.height(); ++y) - { - for (std::ptrdiff_t x = 0; x < output.width(); ++ x) - { - for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; ++channel_index) { - float delta = delta_t * ( - diffusivity[(std::size_t)direction::north](x, y)[channel_index] * nabla[(std::size_t)direction::north](x, y)[channel_index] + diffusivity[(std::size_t)direction::south](x, y)[channel_index] * nabla[(std::size_t)direction::south](x, y)[channel_index] - + diffusivity[(std::size_t)direction::west](x, y)[channel_index] * nabla[(std::size_t)direction::west](x, y)[channel_index] + diffusivity[(std::size_t)direction::east](x, y)[channel_index] * nabla[(std::size_t)direction::east](x, y)[channel_index] - + half * diffusivity[(std::size_t)direction::north_east](x, y)[channel_index] * nabla[(std::size_t)direction::north_east](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::south_east](x, y)[channel_index] * nabla[(std::size_t)direction::south_east](x, y)[channel_index] - + half * diffusivity[(std::size_t)direction::south_west](x, y)[channel_index] * nabla[(std::size_t)direction::south_west](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::north_west](x, y)[channel_index] * nabla[(std::size_t)direction::north_west](x, y)[channel_index] - ); - output(x, y)[channel_index] = output(x, y)[channel_index] + delta; - } - } - } - } -} -}} #include @@ -155,12 +47,12 @@ int main(int argc, char* argv[]) gil::gray64f_image_t output(gray.dimensions()); auto output_view = gil::view(output); - anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); + gil::anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); gil::gray8_image_t true_output(output.dimensions()); gil::transform_pixels(output_view, gil::view(true_output), [](gil::gray64f_pixel_t p) { - return p[0]; + return static_cast(p[0]); }); gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index 2109c555f4..d91ac38e36 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -295,6 +295,113 @@ inline void compute_hessian_entries( detail::convolve_2d(dy, sobel_y, ddyy); } +namespace detail { +enum class direction: std::size_t { + north = 0, + south = 1, + west = 2, + east = 3, + north_east = 4, + south_east = 5, + south_west = 6, + north_west = 7 +}; + +template +void compute_nabla(InputView view, const std::vector& nabla) { + constexpr std::ptrdiff_t input_num_channels = num_channels{}; + static_assert(num_channels{} == input_num_channels, "input and output views must have the same amount of channels"); + for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) + { + for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) + { + for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; ++channel_index) + { + nabla[(std::size_t)direction::north](x, y) = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south](x, y) = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::west](x, y) = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::east](x, y) = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; + + nabla[(std::size_t)direction::north_east](x, y) = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_east](x, y) = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_west](x, y) = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::north_west](x, y) = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; + } + } + } +} + +template +void calculate_diffusvity(std::vector nablas, double kappa, const std::vector diffusivities) +{ + using pixel_type = typename View::value_type; + using channel_type = typename channel_type::type; + BOOST_ASSERT(nablas.size() == diffusivities.size()); + for (std::size_t i = 0; i < nablas.size(); ++i) { + gil::transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) + { + auto op = [kappa](channel_type value) + { + value /= kappa; + auto result = std::exp(-value * value); + return result; + }; + pixel_type result_pixel; + static_transform(p, result_pixel, op); + return result_pixel; + }); + } +} +} // namespace boost::gil::detail + +template +void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, OutputView output) +{ + gil::copy_pixels(input, output); + using element_type = typename OutputView::value_type; + using computation_image_type = image; + + for (unsigned int i = 0; i < num_iter; ++i) { + std::vector nabla_images(8, computation_image_type(input.dimensions())); + std::vector nabla; + std::transform(nabla_images.begin(), nabla_images.end(), + std::back_inserter(nabla), [](computation_image_type& img) + { + return gil::view(img); + }); + + std::vector diffusivity_images(8, computation_image_type(input.dimensions())); + std::vector diffusivity; + std::transform(diffusivity_images.begin(), diffusivity_images.end(), + std::back_inserter(diffusivity), [](computation_image_type& img) + { + return gil::view(img); + }); + + detail::compute_nabla(output, nabla); + detail::calculate_diffusvity(nabla, kappa, diffusivity); + + float half = float(1.0f / 2); + constexpr std::ptrdiff_t channel_count = num_channels{}; + for (std::ptrdiff_t y = 0; y < output.height(); ++y) + { + for (std::ptrdiff_t x = 0; x < output.width(); ++ x) + { + for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; ++channel_index) { + using detail::direction; + auto delta = delta_t * ( + diffusivity[(std::size_t)direction::north](x, y)[channel_index] * nabla[(std::size_t)direction::north](x, y)[channel_index] + diffusivity[(std::size_t)direction::south](x, y)[channel_index] * nabla[(std::size_t)direction::south](x, y)[channel_index] + + diffusivity[(std::size_t)direction::west](x, y)[channel_index] * nabla[(std::size_t)direction::west](x, y)[channel_index] + diffusivity[(std::size_t)direction::east](x, y)[channel_index] * nabla[(std::size_t)direction::east](x, y)[channel_index] + + half * diffusivity[(std::size_t)direction::north_east](x, y)[channel_index] * nabla[(std::size_t)direction::north_east](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::south_east](x, y)[channel_index] * nabla[(std::size_t)direction::south_east](x, y)[channel_index] + + half * diffusivity[(std::size_t)direction::south_west](x, y)[channel_index] * nabla[(std::size_t)direction::south_west](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::north_west](x, y)[channel_index] * nabla[(std::size_t)direction::north_west](x, y)[channel_index] + ); + output(x, y)[channel_index] = output(x, y)[channel_index] + delta; + } + } + } + } +} + }} // namespace boost::gil #endif From 3fcffe375589f96a5a86cb6eb9101d44ff73843b Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Wed, 3 Jun 2020 01:32:10 +0600 Subject: [PATCH 05/48] Add rgb version and squash a bug There was a bug in nabla computation that relied on gray layout --- example/anisotropic_diffusion.cpp | 78 ++++++++++++++----- .../boost/gil/image_processing/numeric.hpp | 18 ++--- 2 files changed, 68 insertions(+), 28 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index afba7ede3e..1ff4f6cff3 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -13,30 +13,15 @@ namespace gil = boost::gil; -namespace boost{namespace gil{ -using gray64f_view_t = gil::image_view>>; -}} - #include -int main(int argc, char* argv[]) +void gray_version(const std::string& input_path, const std::string& output_path, std::uint64_t iteration_count, unsigned int kappa) { - if (argc != 5) - { - std::cerr << "usage: " << argv[0] << " \n"; - return -1; - } - std::string input_path = argv[1]; - std::string output_path = argv[2]; - - unsigned int iteration_count = static_cast(std::stoul(argv[3])); - unsigned int kappa = static_cast(std::stoul(argv[4])); - gil::gray8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); auto input_view = gil::view(input); - gil::gray64f_image_t gray(input.dimensions()); + gil::gray32f_image_t gray(input.dimensions()); auto gray_view = gil::view(gray); gil::transform_pixels(input_view, gray_view, @@ -44,16 +29,71 @@ int main(int argc, char* argv[]) { return p[0]; }); - gil::gray64f_image_t output(gray.dimensions()); + gil::gray32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); gil::anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); gil::gray8_image_t true_output(output.dimensions()); - gil::transform_pixels(output_view, gil::view(true_output), [](gil::gray64f_pixel_t p) + gil::transform_pixels(output_view, gil::view(true_output), [](gil::gray32f_pixel_t p) { return static_cast(p[0]); }); gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); +} + +void rgb_version(const std::string& input_path, const std::string& output_path, std::uint64_t iteration_count, unsigned int kappa) +{ + gil::rgb8_image_t input; + gil::read_image(input_path, input, gil::png_tag{}); + auto input_view = gil::view(input); + + gil::rgb32f_image_t gray(input.dimensions()); + auto gray_view = gil::view(gray); + + gil::transform_pixels(input_view, gray_view, + [](const gil::rgb8_pixel_t& p) + { + return gil::rgb32f_pixel_t(p[0], p[1], p[2]); + }); + gil::rgb32f_image_t output(gray.dimensions()); + auto output_view = gil::view(output); + + gil::anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); + + gil::rgb8_image_t true_output(output.dimensions()); + gil::transform_pixels(output_view, gil::view(true_output), [](gil::rgb32f_pixel_t p) + { + return gil::rgb8_pixel_t(static_cast(p[0]), + static_cast(p[1]), + static_cast(p[2])); + }); + + gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); +} + +int main(int argc, char* argv[]) +{ + if (argc != 6) + { + std::cerr << "usage: " << argv[0] << " \n"; + return -1; + } + std::string input_path = argv[1]; + std::string output_path = argv[2]; + std::string colorspace = argv[3]; + + unsigned int iteration_count = static_cast(std::stoul(argv[4])); + unsigned int kappa = static_cast(std::stoul(argv[5])); + if (colorspace == "gray") + { + gray_version(input_path, output_path, iteration_count, kappa); + } else if (colorspace == "rgb") + { + rgb_version(input_path, output_path, iteration_count, kappa); + } else + { + std::cerr << "unknown colorspace option passed (did you type gray with A?)\n"; + } } \ No newline at end of file diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index d91ac38e36..03efdf84c1 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -317,15 +317,15 @@ void compute_nabla(InputView view, const std::vector& nabla) { { for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; ++channel_index) { - nabla[(std::size_t)direction::north](x, y) = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south](x, y) = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::west](x, y) = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::east](x, y) = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; - - nabla[(std::size_t)direction::north_east](x, y) = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_east](x, y) = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_west](x, y) = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::north_west](x, y) = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::north](x, y)[channel_index] = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south](x, y)[channel_index] = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::west](x, y)[channel_index] = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::east](x, y)[channel_index] = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; + + nabla[(std::size_t)direction::north_east](x, y)[channel_index] = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_east](x, y)[channel_index] = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_west](x, y)[channel_index] = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::north_west](x, y)[channel_index] = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; } } } From a27c941a588618852e2336baab9e5c44656affd0 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Wed, 3 Jun 2020 01:47:52 +0600 Subject: [PATCH 06/48] Add docs --- example/anisotropic_diffusion.cpp | 4 ++-- include/boost/gil/image_processing/numeric.hpp | 15 +++++++++++---- 2 files changed, 13 insertions(+), 6 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index 1ff4f6cff3..cf935fa93f 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -15,7 +15,7 @@ namespace gil = boost::gil; #include -void gray_version(const std::string& input_path, const std::string& output_path, std::uint64_t iteration_count, unsigned int kappa) +void gray_version(const std::string& input_path, const std::string& output_path, unsigned int iteration_count, unsigned int kappa) { gil::gray8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -43,7 +43,7 @@ void gray_version(const std::string& input_path, const std::string& output_path, gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); } -void rgb_version(const std::string& input_path, const std::string& output_path, std::uint64_t iteration_count, unsigned int kappa) +void rgb_version(const std::string& input_path, const std::string& output_path, unsigned int iteration_count, unsigned int kappa) { gil::rgb8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index 03efdf84c1..15c94f803e 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -342,7 +342,7 @@ void calculate_diffusvity(std::vector nablas, double kappa, const std::vec { auto op = [kappa](channel_type value) { - value /= kappa; + value /= static_cast(kappa); auto result = std::exp(-value * value); return result; }; @@ -354,12 +354,19 @@ void calculate_diffusvity(std::vector nablas, double kappa, const std::vec } } // namespace boost::gil::detail +/// \brief Performs diffusion according to Perona-Malik equation +/// +/// Anisotropic diffusion is a smoothing algorithm that respects +/// edge boundaries and can work as an edge detector if suitable +/// iteration count is set and grayscale image view is used +/// as an input template void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, OutputView output) { gil::copy_pixels(input, output); using element_type = typename OutputView::value_type; using computation_image_type = image; + using channel_type = typename channel_type::type; for (unsigned int i = 0; i < num_iter; ++i) { std::vector nabla_images(8, computation_image_type(input.dimensions())); @@ -381,7 +388,7 @@ void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_ detail::compute_nabla(output, nabla); detail::calculate_diffusvity(nabla, kappa, diffusivity); - float half = float(1.0f / 2); + channel_type half = channel_type(1.0f / 2); constexpr std::ptrdiff_t channel_count = num_channels{}; for (std::ptrdiff_t y = 0; y < output.height(); ++y) { @@ -389,12 +396,12 @@ void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_ { for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; ++channel_index) { using detail::direction; - auto delta = delta_t * ( + channel_type delta = static_cast(delta_t * ( diffusivity[(std::size_t)direction::north](x, y)[channel_index] * nabla[(std::size_t)direction::north](x, y)[channel_index] + diffusivity[(std::size_t)direction::south](x, y)[channel_index] * nabla[(std::size_t)direction::south](x, y)[channel_index] + diffusivity[(std::size_t)direction::west](x, y)[channel_index] * nabla[(std::size_t)direction::west](x, y)[channel_index] + diffusivity[(std::size_t)direction::east](x, y)[channel_index] * nabla[(std::size_t)direction::east](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::north_east](x, y)[channel_index] * nabla[(std::size_t)direction::north_east](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::south_east](x, y)[channel_index] * nabla[(std::size_t)direction::south_east](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::south_west](x, y)[channel_index] * nabla[(std::size_t)direction::south_west](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::north_west](x, y)[channel_index] * nabla[(std::size_t)direction::north_west](x, y)[channel_index] - ); + )); output(x, y)[channel_index] = output(x, y)[channel_index] + delta; } } From bb7d83a0b828aae95400e82238b3b96f946e7105 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 6 Jun 2020 02:42:14 +0600 Subject: [PATCH 07/48] Stylistical changes --- example/anisotropic_diffusion.cpp | 39 ++-- .../boost/gil/image_processing/diffusion.hpp | 175 ++++++++++++++++++ .../boost/gil/image_processing/numeric.hpp | 114 ------------ 3 files changed, 201 insertions(+), 127 deletions(-) create mode 100644 include/boost/gil/image_processing/diffusion.hpp diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index cf935fa93f..1d2447eb57 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -1,21 +1,30 @@ +// +// Copyright 2020 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 -#include +#include #include -#include +#include #include +#include #include -#include #include +#include namespace gil = boost::gil; -#include - -void gray_version(const std::string& input_path, const std::string& output_path, unsigned int iteration_count, unsigned int kappa) +void gray_version(std::string const& input_path, + std::string const& output_path, + unsigned int iteration_count, + unsigned int kappa) { gil::gray8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -43,7 +52,10 @@ void gray_version(const std::string& input_path, const std::string& output_path, gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); } -void rgb_version(const std::string& input_path, const std::string& output_path, unsigned int iteration_count, unsigned int kappa) +void rgb_version(const std::string& input_path, + const std::string& output_path, + unsigned int iteration_count, + unsigned int kappa) { gil::rgb8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -70,14 +82,15 @@ void rgb_version(const std::string& input_path, const std::string& output_path, static_cast(p[2])); }); - gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); + gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); } int main(int argc, char* argv[]) { if (argc != 6) { - std::cerr << "usage: " << argv[0] << " \n"; + std::cerr << "usage: " << argv[0] << " " + " \n"; return -1; } std::string input_path = argv[1]; @@ -85,15 +98,15 @@ int main(int argc, char* argv[]) std::string colorspace = argv[3]; unsigned int iteration_count = static_cast(std::stoul(argv[4])); - unsigned int kappa = static_cast(std::stoul(argv[5])); + unsigned int kappa = static_cast(std::stoul(argv[5])); if (colorspace == "gray") { gray_version(input_path, output_path, iteration_count, kappa); } else if (colorspace == "rgb") { rgb_version(input_path, output_path, iteration_count, kappa); - } else + } else { std::cerr << "unknown colorspace option passed (did you type gray with A?)\n"; } -} \ No newline at end of file +} diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp new file mode 100644 index 0000000000..06a86a94fa --- /dev/null +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -0,0 +1,175 @@ +// +// Copyright 2020 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 boost{ namespace gil{ +namespace detail { +enum class direction: std::size_t +{ + north = 0, + south = 1, + west = 2, + east = 3, + north_east = 4, + south_east = 5, + south_west = 6, + north_west = 7 +}; + +// nabla is directed difference in pixels, +// nabla north for example represents intensity +// change in the north direction (up), and nabla south +// represents intensity change in south direction (down). +template +void compute_nabla(InputView view, const std::vector& nabla) +{ + constexpr std::ptrdiff_t input_num_channels = num_channels{}; + static_assert(num_channels{} == input_num_channels, + "input and output views must have the same amount of channels"); + for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) + { + for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) + { + for (std::ptrdiff_t channel_index = 0; + channel_index < input_num_channels; + ++channel_index) + { + nabla[(std::size_t)direction::north](x, y)[channel_index] = + view(x, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south](x, y)[channel_index] = + view(x, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::west](x, y)[channel_index] = + view(x - 1, y)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::east](x, y)[channel_index] = + view(x + 1, y)[channel_index] - view(x, y)[channel_index]; + + nabla[(std::size_t)direction::north_east](x, y)[channel_index] = + view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_east](x, y)[channel_index] = + view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::south_west](x, y)[channel_index] = + view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; + nabla[(std::size_t)direction::north_west](x, y)[channel_index] = + view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; + } + } + } +} + +template +void calculate_diffusvity(std::vector nablas, + double kappa, + std::vector const& diffusivities) +{ + using pixel_type = typename View::value_type; + using channel_type = typename channel_type::type; + for (std::size_t i = 0; i < nablas.size(); ++i) { + gil::transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) + { + auto op = [kappa](channel_type value) + { + value /= static_cast(kappa); + auto result = std::exp(-value * value); + return result; + }; + pixel_type result_pixel; + static_transform(p, result_pixel, op); + return result_pixel; + }); + } +} +} // namespace boost::gil::detail + +/// \brief Performs diffusion according to Perona-Malik equation +/// +/// WARNING: Output channel type must be floating point, +/// otherwise there will be loss in accuracy which most +/// probably will lead to incorrect results (input will be unchanged). +/// Anisotropic diffusion is a smoothing algorithm that respects +/// edge boundaries and can work as an edge detector if suitable +/// iteration count is set and grayscale image view is used +/// as an input +template +void anisotropic_diffusion(InputView input, + unsigned int num_iter, + double delta_t, + double kappa, + OutputView output) +{ + gil::copy_pixels(input, output); + using element_type = typename OutputView::value_type; + using computation_image_type = image; + using channel_type = typename channel_type::type; + + for (unsigned int i = 0; i < num_iter; ++i) { + std::vector nabla_images(8, + computation_image_type(input.dimensions())); + std::vector nabla; + std::transform(nabla_images.begin(), nabla_images.end(), + std::back_inserter(nabla), [](computation_image_type& img) + { + return gil::view(img); + }); + + std::vector diffusivity_images(8, + computation_image_type(input.dimensions())); + std::vector diffusivity; + std::transform(diffusivity_images.begin(), diffusivity_images.end(), + std::back_inserter(diffusivity), [](computation_image_type& img) + { + return gil::view(img); + }); + + detail::compute_nabla(output, nabla); + detail::calculate_diffusvity(nabla, kappa, diffusivity); + + channel_type half = channel_type(1.0f / 2); + constexpr std::ptrdiff_t channel_count = num_channels{}; + for (std::ptrdiff_t y = 0; y < output.height(); ++y) + { + for (std::ptrdiff_t x = 0; x < output.width(); ++ x) + { + for (std::ptrdiff_t channel_index = 0; + channel_index < channel_count; + ++channel_index) + { + using detail::direction; + channel_type delta = static_cast(delta_t * ( + diffusivity[(std::size_t)direction::north](x, y)[channel_index] + * nabla[(std::size_t)direction::north](x, y)[channel_index] + + diffusivity[(std::size_t)direction::south](x, y)[channel_index] + * nabla[(std::size_t)direction::south](x, y)[channel_index] + + diffusivity[(std::size_t)direction::west](x, y)[channel_index] + * nabla[(std::size_t)direction::west](x, y)[channel_index] + + diffusivity[(std::size_t)direction::east](x, y)[channel_index] + * nabla[(std::size_t)direction::east](x, y)[channel_index] + + half + * diffusivity[(std::size_t)direction::north_east](x, y)[channel_index] + * nabla[(std::size_t)direction::north_east](x, y)[channel_index] + + half + * diffusivity[(std::size_t)direction::south_east](x, y)[channel_index] + * nabla[(std::size_t)direction::south_east](x, y)[channel_index] + + half + * diffusivity[(std::size_t)direction::south_west](x, y)[channel_index] + * nabla[(std::size_t)direction::south_west](x, y)[channel_index] + + half + * diffusivity[(std::size_t)direction::north_west](x, y)[channel_index] + * nabla[(std::size_t)direction::north_west](x, y)[channel_index] + )); + output(x, y)[channel_index] = output(x, y)[channel_index] + delta; + } + } + } + } +} + +}} // namespace boost::gil diff --git a/include/boost/gil/image_processing/numeric.hpp b/include/boost/gil/image_processing/numeric.hpp index 15c94f803e..2109c555f4 100644 --- a/include/boost/gil/image_processing/numeric.hpp +++ b/include/boost/gil/image_processing/numeric.hpp @@ -295,120 +295,6 @@ inline void compute_hessian_entries( detail::convolve_2d(dy, sobel_y, ddyy); } -namespace detail { -enum class direction: std::size_t { - north = 0, - south = 1, - west = 2, - east = 3, - north_east = 4, - south_east = 5, - south_west = 6, - north_west = 7 -}; - -template -void compute_nabla(InputView view, const std::vector& nabla) { - constexpr std::ptrdiff_t input_num_channels = num_channels{}; - static_assert(num_channels{} == input_num_channels, "input and output views must have the same amount of channels"); - for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) - { - for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) - { - for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; ++channel_index) - { - nabla[(std::size_t)direction::north](x, y)[channel_index] = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south](x, y)[channel_index] = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::west](x, y)[channel_index] = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::east](x, y)[channel_index] = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; - - nabla[(std::size_t)direction::north_east](x, y)[channel_index] = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_east](x, y)[channel_index] = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_west](x, y)[channel_index] = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::north_west](x, y)[channel_index] = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; - } - } - } -} - -template -void calculate_diffusvity(std::vector nablas, double kappa, const std::vector diffusivities) -{ - using pixel_type = typename View::value_type; - using channel_type = typename channel_type::type; - BOOST_ASSERT(nablas.size() == diffusivities.size()); - for (std::size_t i = 0; i < nablas.size(); ++i) { - gil::transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) - { - auto op = [kappa](channel_type value) - { - value /= static_cast(kappa); - auto result = std::exp(-value * value); - return result; - }; - pixel_type result_pixel; - static_transform(p, result_pixel, op); - return result_pixel; - }); - } -} -} // namespace boost::gil::detail - -/// \brief Performs diffusion according to Perona-Malik equation -/// -/// Anisotropic diffusion is a smoothing algorithm that respects -/// edge boundaries and can work as an edge detector if suitable -/// iteration count is set and grayscale image view is used -/// as an input -template -void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, OutputView output) -{ - gil::copy_pixels(input, output); - using element_type = typename OutputView::value_type; - using computation_image_type = image; - using channel_type = typename channel_type::type; - - for (unsigned int i = 0; i < num_iter; ++i) { - std::vector nabla_images(8, computation_image_type(input.dimensions())); - std::vector nabla; - std::transform(nabla_images.begin(), nabla_images.end(), - std::back_inserter(nabla), [](computation_image_type& img) - { - return gil::view(img); - }); - - std::vector diffusivity_images(8, computation_image_type(input.dimensions())); - std::vector diffusivity; - std::transform(diffusivity_images.begin(), diffusivity_images.end(), - std::back_inserter(diffusivity), [](computation_image_type& img) - { - return gil::view(img); - }); - - detail::compute_nabla(output, nabla); - detail::calculate_diffusvity(nabla, kappa, diffusivity); - - channel_type half = channel_type(1.0f / 2); - constexpr std::ptrdiff_t channel_count = num_channels{}; - for (std::ptrdiff_t y = 0; y < output.height(); ++y) - { - for (std::ptrdiff_t x = 0; x < output.width(); ++ x) - { - for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; ++channel_index) { - using detail::direction; - channel_type delta = static_cast(delta_t * ( - diffusivity[(std::size_t)direction::north](x, y)[channel_index] * nabla[(std::size_t)direction::north](x, y)[channel_index] + diffusivity[(std::size_t)direction::south](x, y)[channel_index] * nabla[(std::size_t)direction::south](x, y)[channel_index] - + diffusivity[(std::size_t)direction::west](x, y)[channel_index] * nabla[(std::size_t)direction::west](x, y)[channel_index] + diffusivity[(std::size_t)direction::east](x, y)[channel_index] * nabla[(std::size_t)direction::east](x, y)[channel_index] - + half * diffusivity[(std::size_t)direction::north_east](x, y)[channel_index] * nabla[(std::size_t)direction::north_east](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::south_east](x, y)[channel_index] * nabla[(std::size_t)direction::south_east](x, y)[channel_index] - + half * diffusivity[(std::size_t)direction::south_west](x, y)[channel_index] * nabla[(std::size_t)direction::south_west](x, y)[channel_index] + half * diffusivity[(std::size_t)direction::north_west](x, y)[channel_index] * nabla[(std::size_t)direction::north_west](x, y)[channel_index] - )); - output(x, y)[channel_index] = output(x, y)[channel_index] + delta; - } - } - } - } -} - }} // namespace boost::gil #endif From dc278a4b037c0d59784bab44c38c8958511820e6 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 6 Jun 2020 03:00:06 +0600 Subject: [PATCH 08/48] More stylistical changes --- include/boost/gil/image_processing/diffusion.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 06a86a94fa..b1c886c376 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -29,8 +29,8 @@ enum class direction: std::size_t // nabla north for example represents intensity // change in the north direction (up), and nabla south // represents intensity change in south direction (down). -template -void compute_nabla(InputView view, const std::vector& nabla) +template +void compute_nabla(InputView view, std::vector const& nabla) { constexpr std::ptrdiff_t input_num_channels = num_channels{}; static_assert(num_channels{} == input_num_channels, @@ -110,7 +110,8 @@ void anisotropic_diffusion(InputView input, using computation_image_type = image; using channel_type = typename channel_type::type; - for (unsigned int i = 0; i < num_iter; ++i) { + for (unsigned int i = 0; i < num_iter; ++i) + { std::vector nabla_images(8, computation_image_type(input.dimensions())); std::vector nabla; From f128ef87b1fa9d0d78fb5da155e239c21839472a Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 6 Jun 2020 03:06:02 +0600 Subject: [PATCH 09/48] Add missing include --- example/anisotropic_diffusion.cpp | 1 - include/boost/gil/image_processing/diffusion.hpp | 2 ++ 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index 1d2447eb57..fe56ffc891 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -17,7 +17,6 @@ #include #include #include -#include namespace gil = boost::gil; diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index b1c886c376..87da589e0f 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -11,6 +11,8 @@ #include #include +#include + namespace boost{ namespace gil{ namespace detail { enum class direction: std::size_t From 0c486e8fd003ad436bd2606b2d5fba0c441acfa4 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 6 Jun 2020 04:53:38 +0600 Subject: [PATCH 10/48] Dealth with the direction enum --- .../boost/gil/image_processing/diffusion.hpp | 55 ++++++++++--------- 1 file changed, 28 insertions(+), 27 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 87da589e0f..92bb422be9 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -15,8 +15,8 @@ namespace boost{ namespace gil{ namespace detail { -enum class direction: std::size_t -{ +namespace direction { +enum { north = 0, south = 1, west = 2, @@ -26,6 +26,7 @@ enum class direction: std::size_t south_west = 6, north_west = 7 }; +} // nabla is directed difference in pixels, // nabla north for example represents intensity @@ -45,22 +46,23 @@ void compute_nabla(InputView view, std::vector const& nabla) channel_index < input_num_channels; ++channel_index) { - nabla[(std::size_t)direction::north](x, y)[channel_index] = + using namespace direction; + nabla[north](x, y)[channel_index] = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south](x, y)[channel_index] = + nabla[south](x, y)[channel_index] = view(x, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::west](x, y)[channel_index] = + nabla[west](x, y)[channel_index] = view(x - 1, y)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::east](x, y)[channel_index] = + nabla[east](x, y)[channel_index] = view(x + 1, y)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::north_east](x, y)[channel_index] = + nabla[north_east](x, y)[channel_index] = view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_east](x, y)[channel_index] = + nabla[south_east](x, y)[channel_index] = view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::south_west](x, y)[channel_index] = + nabla[south_west](x, y)[channel_index] = view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[(std::size_t)direction::north_west](x, y)[channel_index] = + nabla[north_west](x, y)[channel_index] = view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; } } @@ -145,28 +147,27 @@ void anisotropic_diffusion(InputView input, channel_index < channel_count; ++channel_index) { - using detail::direction; + using namespace detail::direction; channel_type delta = static_cast(delta_t * ( - diffusivity[(std::size_t)direction::north](x, y)[channel_index] - * nabla[(std::size_t)direction::north](x, y)[channel_index] - + diffusivity[(std::size_t)direction::south](x, y)[channel_index] - * nabla[(std::size_t)direction::south](x, y)[channel_index] - + diffusivity[(std::size_t)direction::west](x, y)[channel_index] - * nabla[(std::size_t)direction::west](x, y)[channel_index] - + diffusivity[(std::size_t)direction::east](x, y)[channel_index] - * nabla[(std::size_t)direction::east](x, y)[channel_index] + diffusivity[north](x, y)[channel_index] * nabla[north](x, y)[channel_index] + + diffusivity[south](x, y)[channel_index] + * nabla[south](x, y)[channel_index] + + diffusivity[west](x, y)[channel_index] + * nabla[west](x, y)[channel_index] + + diffusivity[east](x, y)[channel_index] + * nabla[east](x, y)[channel_index] + half - * diffusivity[(std::size_t)direction::north_east](x, y)[channel_index] - * nabla[(std::size_t)direction::north_east](x, y)[channel_index] + * diffusivity[north_east](x, y)[channel_index] + * nabla[north_east](x, y)[channel_index] + half - * diffusivity[(std::size_t)direction::south_east](x, y)[channel_index] - * nabla[(std::size_t)direction::south_east](x, y)[channel_index] + * diffusivity[south_east](x, y)[channel_index] + * nabla[south_east](x, y)[channel_index] + half - * diffusivity[(std::size_t)direction::south_west](x, y)[channel_index] - * nabla[(std::size_t)direction::south_west](x, y)[channel_index] + * diffusivity[south_west](x, y)[channel_index] + * nabla[south_west](x, y)[channel_index] + half - * diffusivity[(std::size_t)direction::north_west](x, y)[channel_index] - * nabla[(std::size_t)direction::north_west](x, y)[channel_index] + * diffusivity[north_west](x, y)[channel_index] + * nabla[north_west](x, y)[channel_index] )); output(x, y)[channel_index] = output(x, y)[channel_index] + delta; } From 65ac1512dfa1068839e2ea5f46ff59c532983256 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 15 Jun 2020 01:49:18 +0600 Subject: [PATCH 11/48] All cases version --- .../boost/gil/image_processing/diffusion.hpp | 241 +++++++++++++----- 1 file changed, 180 insertions(+), 61 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 92bb422be9..3a8ecafd63 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -13,10 +13,11 @@ #include -namespace boost{ namespace gil{ +namespace boost { namespace gil { namespace detail { namespace direction { -enum { +enum +{ north = 0, south = 1, west = 2, @@ -33,20 +34,154 @@ enum { // change in the north direction (up), and nabla south // represents intensity change in south direction (down). template -void compute_nabla(InputView view, std::vector const& nabla) +void compute_nabla(InputView view, std::vector const &nabla) { constexpr std::ptrdiff_t input_num_channels = num_channels{}; + using pixel_type = typename OutputView::value_type; static_assert(num_channels{} == input_num_channels, - "input and output views must have the same amount of channels"); + "input and output views must have the same amount of channels"); + const auto width = view.width(); + const auto height = view.height(); + using namespace direction; + const auto zero_pixel = []() { + pixel_type pixel; + for (std::ptrdiff_t i = 0; i < input_num_channels; ++i) + { + pixel[i] = 0; + } + + return pixel; + }(); + auto minus = [](pixel_type lhs, pixel_type rhs) { + for (std::ptrdiff_t i = 0; i < input_num_channels; ++i) + { + lhs[i] -= rhs[i]; + } + + return lhs; + }; + // NW corner + { + const std::ptrdiff_t x = 0; + const std::ptrdiff_t y = 0; + const auto current = view(x, y); + nabla[north](x, y) = zero_pixel; + nabla[south](x, y) = minus(view(x, y + 1), current); + nabla[west](x, y) = zero_pixel; + nabla[east](x, y) = minus(view(x + 1, y), current); + nabla[north_east](x, y) = zero_pixel; + nabla[south_east](x, y) = minus(view(x + 1, y + 1), current); + nabla[south_west](x, y) = zero_pixel; + nabla[north_west](x, y) = zero_pixel; + } + + // NE corner + { + const std::ptrdiff_t x = width - 1; + const std::ptrdiff_t y = 0; + const auto current = view(x, y); + nabla[north](x, y) = zero_pixel; + nabla[south](x, y) = minus(view(x, y + 1), current); + nabla[west](x, y) = minus(view(x - 1, y), current); + nabla[east](x, y) = zero_pixel; + nabla[north_east](x, y) = zero_pixel; + nabla[south_east](x, y) = zero_pixel; + nabla[south_west](x, y) = minus(view(x - 1, y + 1), current); + nabla[north_west](x, y) = zero_pixel; + } + + // SW corner + { + const std::ptrdiff_t x = 0; + const std::ptrdiff_t y = height - 1; + const auto current = view(x, y); + nabla[north](x, y) = minus(view(x, y - 1), current); + nabla[south](x, y) = zero_pixel; + nabla[west](x, y) = zero_pixel; + nabla[east](x, y) = minus(view(x + 1, y), current); + nabla[north_east](x, y) = minus(view(x + 1, y - 1), current); + nabla[south_east](x, y) = zero_pixel; + nabla[south_west](x, y) = zero_pixel; + nabla[north_west](x, y) = zero_pixel; + } + + // SE corner + { + const std::ptrdiff_t x = width - 1; + const std::ptrdiff_t y = height - 1; + const auto current = view(x, y); + nabla[north](x, y) = minus(view(x, y - 1), current); + nabla[south](x, y) = zero_pixel; + nabla[west](x, y) = minus(view(x - 1, y), current); + nabla[east](x, y) = zero_pixel; + nabla[north_east](x, y) = zero_pixel; + nabla[south_east](x, y) = zero_pixel; + nabla[south_west](x, y) = zero_pixel; + nabla[north_west](x, y) = minus(view(x - 1, y - 1), current); + } + + // first and last rows + { + const std::ptrdiff_t y1 = 0; + const std::ptrdiff_t y2 = height - 1; + for (std::ptrdiff_t x = 1; x < width - 1; ++x) + { + const auto current1 = view(x, y1); + const auto current2 = view(x, y2); + + nabla[north](x, y1) = zero_pixel; + nabla[north](x, y2) = minus(view(x, y2 - 1), current2); + nabla[south](x, y1) = minus(view(x, y1 + 1), current1); + nabla[south](x, y2) = zero_pixel; + nabla[west](x, y1) = minus(view(x - 1, y1), current1); + nabla[west](x, y2) = minus(view(x - 1, y2), current2); + nabla[east](x, y1) = minus(view(x + 1, y1), current1); + nabla[east](x, y2) = minus(view(x + 1, y2), current2); + + nabla[north_east](x, y1) = zero_pixel; + nabla[north_east](x, y2) = minus(view(x + 1, y2 - 1), current2); + nabla[south_east](x, y1) = minus(view(x + 1, y1 + 1), current1); + nabla[south_east](x, y2) = zero_pixel; + nabla[south_west](x, y1) = minus(view(x - 1, y1 + 1), current1); + nabla[north_west](x, y1) = zero_pixel; + } + } + + // first and last columns + { + const std::ptrdiff_t x1 = 0; + const std::ptrdiff_t x2 = width - 1; + for (std::ptrdiff_t y = 1; y < height - 1; ++y) + { + const auto current1 = view(x1, y); + const auto current2 = view(x2, y); + nabla[north](x1, y) = minus(view(x1, y - 1), current1); + nabla[north](x2, y) = minus(view(x2, y - 1), current2); + nabla[south](x1, y) = minus(view(x1, y + 1), current1); + nabla[south](x2, y) = minus(view(x2, y + 1), current2); + nabla[west](x1, y) = zero_pixel; + nabla[west](x2, y) = minus(view(x2 - 1, y), current2); + nabla[east](x1, y) = minus(view(x1 + 1, y), current1); + nabla[east](x2, y) = zero_pixel; + + nabla[north_east](x1, y) = minus(view(x1 + 1, y - 1), current1); + nabla[north_west](x2, y) = zero_pixel; + nabla[south_east](x1, y) = minus(view(x1 + 1, y + 1), current1); + nabla[south_east](x2, y) = zero_pixel; + nabla[south_west](x1, y) = zero_pixel; + nabla[south_west](x2, y) = minus(view(x2 - 1, y + 1), current2); + nabla[north_west](x1, y) = zero_pixel; + nabla[north_west](x2, y) = minus(view(x2 - 1, y - 1), current2); + } + } + for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) { for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) { - for (std::ptrdiff_t channel_index = 0; - channel_index < input_num_channels; - ++channel_index) + for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; + ++channel_index) { - using namespace direction; nabla[north](x, y)[channel_index] = view(x, y - 1)[channel_index] - view(x, y)[channel_index]; nabla[south](x, y)[channel_index] = @@ -70,17 +205,15 @@ void compute_nabla(InputView view, std::vector const& nabla) } template -void calculate_diffusvity(std::vector nablas, - double kappa, - std::vector const& diffusivities) +void calculate_diffusvity(std::vector nablas, double kappa, + std::vector const &diffusivities) { using pixel_type = typename View::value_type; using channel_type = typename channel_type::type; - for (std::size_t i = 0; i < nablas.size(); ++i) { - gil::transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) - { - auto op = [kappa](channel_type value) - { + for (std::size_t i = 0; i < nablas.size(); ++i) + { + transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) { + auto op = [kappa](channel_type value) { value /= static_cast(kappa); auto result = std::exp(-value * value); return result; @@ -91,7 +224,7 @@ void calculate_diffusvity(std::vector nablas, }); } } -} // namespace boost::gil::detail +} // namespace detail /// \brief Performs diffusion according to Perona-Malik equation /// @@ -103,36 +236,28 @@ void calculate_diffusvity(std::vector nablas, /// iteration count is set and grayscale image view is used /// as an input template -void anisotropic_diffusion(InputView input, - unsigned int num_iter, - double delta_t, - double kappa, - OutputView output) +void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, + OutputView output) { - gil::copy_pixels(input, output); + copy_pixels(input, output); using element_type = typename OutputView::value_type; using computation_image_type = image; using channel_type = typename channel_type::type; for (unsigned int i = 0; i < num_iter; ++i) { - std::vector nabla_images(8, - computation_image_type(input.dimensions())); + std::vector nabla_images( + 8, computation_image_type(input.dimensions())); std::vector nabla; - std::transform(nabla_images.begin(), nabla_images.end(), - std::back_inserter(nabla), [](computation_image_type& img) - { - return gil::view(img); - }); + std::transform(nabla_images.begin(), nabla_images.end(), std::back_inserter(nabla), + [](computation_image_type &img) { return gil::view(img); }); - std::vector diffusivity_images(8, - computation_image_type(input.dimensions())); + std::vector diffusivity_images( + 8, computation_image_type(input.dimensions())); std::vector diffusivity; std::transform(diffusivity_images.begin(), diffusivity_images.end(), - std::back_inserter(diffusivity), [](computation_image_type& img) - { - return gil::view(img); - }); + std::back_inserter(diffusivity), + [](computation_image_type &img) { return gil::view(img); }); detail::compute_nabla(output, nabla); detail::calculate_diffusvity(nabla, kappa, diffusivity); @@ -141,34 +266,28 @@ void anisotropic_diffusion(InputView input, constexpr std::ptrdiff_t channel_count = num_channels{}; for (std::ptrdiff_t y = 0; y < output.height(); ++y) { - for (std::ptrdiff_t x = 0; x < output.width(); ++ x) + for (std::ptrdiff_t x = 0; x < output.width(); ++x) { - for (std::ptrdiff_t channel_index = 0; - channel_index < channel_count; - ++channel_index) + for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; + ++channel_index) { using namespace detail::direction; - channel_type delta = static_cast(delta_t * ( - diffusivity[north](x, y)[channel_index] * nabla[north](x, y)[channel_index] - + diffusivity[south](x, y)[channel_index] - * nabla[south](x, y)[channel_index] - + diffusivity[west](x, y)[channel_index] - * nabla[west](x, y)[channel_index] - + diffusivity[east](x, y)[channel_index] - * nabla[east](x, y)[channel_index] - + half - * diffusivity[north_east](x, y)[channel_index] - * nabla[north_east](x, y)[channel_index] - + half - * diffusivity[south_east](x, y)[channel_index] - * nabla[south_east](x, y)[channel_index] - + half - * diffusivity[south_west](x, y)[channel_index] - * nabla[south_west](x, y)[channel_index] - + half - * diffusivity[north_west](x, y)[channel_index] - * nabla[north_west](x, y)[channel_index] - )); + channel_type delta = static_cast( + delta_t * + (diffusivity[north](x, y)[channel_index] * + nabla[north](x, y)[channel_index] + + diffusivity[south](x, y)[channel_index] * + nabla[south](x, y)[channel_index] + + diffusivity[west](x, y)[channel_index] * nabla[west](x, y)[channel_index] + + diffusivity[east](x, y)[channel_index] * nabla[east](x, y)[channel_index] + + half * diffusivity[north_east](x, y)[channel_index] * + nabla[north_east](x, y)[channel_index] + + half * diffusivity[south_east](x, y)[channel_index] * + nabla[south_east](x, y)[channel_index] + + half * diffusivity[south_west](x, y)[channel_index] * + nabla[south_west](x, y)[channel_index] + + half * diffusivity[north_west](x, y)[channel_index] * + nabla[north_west](x, y)[channel_index])); output(x, y)[channel_index] = output(x, y)[channel_index] + delta; } } From de226761e5ea197b8bd8984d589084a6d423b860 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 15 Jun 2020 02:02:08 +0600 Subject: [PATCH 12/48] Modify example to show the properties --- example/anisotropic_diffusion.cpp | 77 ++++++++++++++++++------------- 1 file changed, 46 insertions(+), 31 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index fe56ffc891..6b25ed6fcc 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -8,8 +8,8 @@ #include #include #include -#include #include +#include #include #include @@ -20,10 +20,8 @@ namespace gil = boost::gil; -void gray_version(std::string const& input_path, - std::string const& output_path, - unsigned int iteration_count, - unsigned int kappa) +void gray_version(std::string const &input_path, std::string const &output_path, + unsigned int iteration_count, unsigned int kappa) { gil::gray8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -32,29 +30,29 @@ void gray_version(std::string const& input_path, gil::gray32f_image_t gray(input.dimensions()); auto gray_view = gil::view(gray); - gil::transform_pixels(input_view, gray_view, - [](const gil::gray8_pixel_t& p) - { - return p[0]; - }); + gil::transform_pixels(input_view, gray_view, [](const gil::gray8_pixel_t &p) { return p[0]; }); + double sum_before = 0; + gil::for_each_pixel(gray_view, [&sum_before](gil::gray32f_pixel_t p) { sum_before += p[0]; }); gil::gray32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); gil::anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); + double sum_after = 0; + gil::for_each_pixel(output_view, [&sum_after](gil::gray32f_pixel_t p) { sum_after += p[0]; }); gil::gray8_image_t true_output(output.dimensions()); - gil::transform_pixels(output_view, gil::view(true_output), [](gil::gray32f_pixel_t p) - { - return static_cast(p[0]); - }); + gil::transform_pixels(output_view, gil::view(true_output), + [](gil::gray32f_pixel_t p) { return static_cast(p[0]); }); gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); + + std::cout << "sum of intensity before diffusion: " << sum_before << '\n' + << "sum of intensity after diffusion: " << sum_after << '\n' + << "difference: " << sum_after - sum_before << '\n'; } -void rgb_version(const std::string& input_path, - const std::string& output_path, - unsigned int iteration_count, - unsigned int kappa) +void rgb_version(const std::string &input_path, const std::string &output_path, + unsigned int iteration_count, unsigned int kappa) { gil::rgb8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -63,33 +61,48 @@ void rgb_version(const std::string& input_path, gil::rgb32f_image_t gray(input.dimensions()); auto gray_view = gil::view(gray); - gil::transform_pixels(input_view, gray_view, - [](const gil::rgb8_pixel_t& p) - { + gil::transform_pixels(input_view, gray_view, [](const gil::rgb8_pixel_t &p) { return gil::rgb32f_pixel_t(p[0], p[1], p[2]); }); + double sum_before[3] = {}; + gil::for_each_pixel(gray_view, [&sum_before](gil::rgb32f_pixel_t p) { + sum_before[0] += p[0]; + sum_before[1] += p[1]; + sum_before[2] += p[2]; + }); gil::rgb32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); gil::anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); + double sum_after[3] = {}; + gil::for_each_pixel(gray_view, [&sum_after](gil::rgb32f_pixel_t p) { + sum_after[0] += p[0]; + sum_after[1] += p[1]; + sum_after[2] += p[2]; + }); gil::rgb8_image_t true_output(output.dimensions()); - gil::transform_pixels(output_view, gil::view(true_output), [](gil::rgb32f_pixel_t p) - { - return gil::rgb8_pixel_t(static_cast(p[0]), - static_cast(p[1]), - static_cast(p[2])); + gil::transform_pixels(output_view, gil::view(true_output), [](gil::rgb32f_pixel_t p) { + return gil::rgb8_pixel_t(static_cast(p[0]), static_cast(p[1]), + static_cast(p[2])); }); gil::write_view(output_path, gil::view(true_output), gil::png_tag{}); + std::cout << "sum of intensity before diffusion: (" << sum_before[0] << ", " << sum_before[1] + << ", " << sum_before[2] << ")\n" + << "sum of intensity after diffusion: (" << sum_after[0] << ", " << sum_after[1] + << ", " << sum_after[2] << ")\n" + << "difference: (" << sum_after[0] - sum_before[0] << ", " + << sum_after[1] - sum_before[1] << ", " << sum_after[2] - sum_before[2] << ")\n"; } -int main(int argc, char* argv[]) +int main(int argc, char *argv[]) { if (argc != 6) { - std::cerr << "usage: " << argv[0] << " " - " \n"; + std::cerr << "usage: " << argv[0] + << " " + " \n"; return -1; } std::string input_path = argv[1]; @@ -101,10 +114,12 @@ int main(int argc, char* argv[]) if (colorspace == "gray") { gray_version(input_path, output_path, iteration_count, kappa); - } else if (colorspace == "rgb") + } + else if (colorspace == "rgb") { rgb_version(input_path, output_path, iteration_count, kappa); - } else + } + else { std::cerr << "unknown colorspace option passed (did you type gray with A?)\n"; } From 3b3f4871e640c8e67e97e14dfc73eb67a5c537fc Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 16 Jun 2020 00:24:06 +0600 Subject: [PATCH 13/48] Improve documentation and static_assert The output type must be floating point, thus a check was added to make sure it is the case. Though I had to add specific cases for float32_t as std::is_floating_point does not consider it a floating point type --- .../boost/gil/image_processing/diffusion.hpp | 65 ++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 3a8ecafd63..b4b77bea0c 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -60,7 +60,24 @@ void compute_nabla(InputView view, std::vector const &nabla) return lhs; }; + + // north: x, y - 1 + // south: x, y + 1 + // east: x + 1, y + // west: x - 1, y + // NW: x - 1, y - 1 + // NE: x + 1, y - 1 + // SE: x + 1, y + 1 + // SW: x - 1, y + 1 + // NW corner + /* + XOOOO + OOOOO + OOOOO + OOOOO + OOOOO + */ { const std::ptrdiff_t x = 0; const std::ptrdiff_t y = 0; @@ -76,6 +93,13 @@ void compute_nabla(InputView view, std::vector const &nabla) } // NE corner + /* + OOOOX + OOOOO + OOOOO + OOOOO + OOOOO + */ { const std::ptrdiff_t x = width - 1; const std::ptrdiff_t y = 0; @@ -91,6 +115,13 @@ void compute_nabla(InputView view, std::vector const &nabla) } // SW corner + /* + OOOOO + OOOOO + OOOOO + OOOOO + XOOOO + */ { const std::ptrdiff_t x = 0; const std::ptrdiff_t y = height - 1; @@ -106,6 +137,13 @@ void compute_nabla(InputView view, std::vector const &nabla) } // SE corner + /* + OOOOO + OOOOO + OOOOO + OOOOO + OOOOX + */ { const std::ptrdiff_t x = width - 1; const std::ptrdiff_t y = height - 1; @@ -121,6 +159,13 @@ void compute_nabla(InputView view, std::vector const &nabla) } // first and last rows + /* + OXXXO + OOOOO + OOOOO + OOOOO + OXXXO + */ { const std::ptrdiff_t y1 = 0; const std::ptrdiff_t y2 = height - 1; @@ -148,6 +193,13 @@ void compute_nabla(InputView view, std::vector const &nabla) } // first and last columns + /* + OOOOO + XOOOX + XOOOX + XOOOX + OOOOO + */ { const std::ptrdiff_t x1 = 0; const std::ptrdiff_t x2 = width - 1; @@ -174,7 +226,14 @@ void compute_nabla(InputView view, std::vector const &nabla) nabla[north_west](x2, y) = minus(view(x2 - 1, y - 1), current2); } } - + // middle case + /* + OOOOO + OXXXO + OXXXO + OXXXO + OOOOO + */ for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) { for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) @@ -204,6 +263,7 @@ void compute_nabla(InputView view, std::vector const &nabla) } } +// formula is c = exp(-(nabla/kappa)^2) template void calculate_diffusvity(std::vector nablas, double kappa, std::vector const &diffusivities) @@ -243,6 +303,9 @@ void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_ using element_type = typename OutputView::value_type; using computation_image_type = image; using channel_type = typename channel_type::type; + static_assert( + std::is_floating_point::value || std::is_same::value, + "Output view must be floating point, otherwise accuracy loss will nullify diffusion"); for (unsigned int i = 0; i < num_iter; ++i) { From b1354387c69741a1cacecb11d93698f6a20339e0 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 16 Jun 2020 01:06:42 +0600 Subject: [PATCH 14/48] Add heat conservation test --- test/core/image_processing/CMakeLists.txt | 3 +- .../anisotropic_diffusion.cpp | 66 +++++++++++++++++++ 2 files changed, 68 insertions(+), 1 deletion(-) create mode 100644 test/core/image_processing/anisotropic_diffusion.cpp diff --git a/test/core/image_processing/CMakeLists.txt b/test/core/image_processing/CMakeLists.txt index 5ecd590629..2ad32f4f28 100644 --- a/test/core/image_processing/CMakeLists.txt +++ b/test/core/image_processing/CMakeLists.txt @@ -34,7 +34,8 @@ foreach(_name hessian box_filter median_filter - sobel_scharr) + sobel_scharr + anisotropic_diffusion) set(_test t_core_image_processing_${_name}) set(_target test_core_image_processing_${_name}) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp new file mode 100644 index 0000000000..e7321f2e90 --- /dev/null +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -0,0 +1,66 @@ +// +// Copyright 2020 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 + +#include + +#include +#include + +namespace gil = boost::gil; + +template void heat_conservation_test() +{ + std::mt19937 twister; + std::uniform_int_distribution dist; + + ImageType image(32, 32); + auto view = gil::view(image); + constexpr std::ptrdiff_t num_channels = gil::num_channels::value; + double before_diffusion[num_channels] = {0}; + for (auto &pixel : view) + { + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + pixel[channel_index] = dist(twister); + before_diffusion[channel_index] += pixel[channel_index]; + } + } + + gil::gray32f_image_t output(32, 32); + auto output_view = gil::view(output); + gil::anisotropic_diffusion(view, 30, 1 / 7.0, 30, output_view); + double after_diffusion[num_channels] = {0}; + for (const auto &pixel : output_view) + { + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + after_diffusion[channel_index] += pixel[channel_index]; + } + } + + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + const auto percentage = + std::abs(after_diffusion[channel_index] - before_diffusion[channel_index]) / + after_diffusion[channel_index] * 100.0; + BOOST_TEST(percentage < 10.0); + } +} + +int main() +{ + heat_conservation_test(); + heat_conservation_test(); + + return boost::report_errors(); +} From b3a5540036debf36de824ce3a3f0cafb17a57444 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 16 Jun 2020 01:36:05 +0600 Subject: [PATCH 15/48] Improve testing procedure, raise bar Now conservation is only within 10% --- .../anisotropic_diffusion.cpp | 24 +++++++++++++++---- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index e7321f2e90..4ba703fb62 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -18,9 +18,10 @@ namespace gil = boost::gil; -template void heat_conservation_test() +template +void heat_conservation_test(std::uint32_t seed) { - std::mt19937 twister; + std::mt19937 twister(seed); std::uniform_int_distribution dist; ImageType image(32, 32); @@ -38,7 +39,7 @@ template void heat_conservation_test() gil::gray32f_image_t output(32, 32); auto output_view = gil::view(output); - gil::anisotropic_diffusion(view, 30, 1 / 7.0, 30, output_view); + gil::anisotropic_diffusion(view, 10, 1 / 7.0, 30, output_view); double after_diffusion[num_channels] = {0}; for (const auto &pixel : output_view) { @@ -53,14 +54,27 @@ template void heat_conservation_test() const auto percentage = std::abs(after_diffusion[channel_index] - before_diffusion[channel_index]) / after_diffusion[channel_index] * 100.0; + std::cout << percentage << ' '; BOOST_TEST(percentage < 10.0); } + std::cout << '\n'; +} + +template +void convergence_to_mean_test() +{ } int main() { - heat_conservation_test(); - heat_conservation_test(); + for (std::uint32_t seed = 0; seed < 1000; ++seed) + { + std::cout << "seed: " << seed << '\n'; + std::cout << "gray8 test:\n"; + heat_conservation_test(seed); + std::cout << "rgb8 test:\n"; + heat_conservation_test(seed); + } return boost::report_errors(); } From 38d74158931a17bfe26ab4f322d8aaed9ca99f85 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 16 Jun 2020 01:56:34 +0600 Subject: [PATCH 16/48] Add convergence to mean test The test seems to be failing though, so revisiting the correctness of the algorithm is needed --- .../anisotropic_diffusion.cpp | 85 ++++++++++++++++++- 1 file changed, 83 insertions(+), 2 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 4ba703fb62..a6f1dc6c46 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -61,8 +61,81 @@ void heat_conservation_test(std::uint32_t seed) } template -void convergence_to_mean_test() +void convergence_to_mean_test(std::uint32_t seed) { + std::mt19937 twister(seed); + std::uniform_int_distribution dist; + + const std::size_t size = 32; + ImageType image(size, size); + auto view = gil::view(image); + constexpr std::ptrdiff_t num_channels = gil::num_channels::value; + double mean_before_diffusion[num_channels] = {0}; + double before_diffusion[num_channels] = {0}; + for (auto &pixel : view) + { + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + pixel[channel_index] = dist(twister); + mean_before_diffusion[channel_index] += pixel[channel_index]; + } + } + + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + mean_before_diffusion[channel_index] /= size * size; + } + + for (auto &pixel : view) + { + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + double difference = pixel[channel_index] - mean_before_diffusion[channel_index]; + before_diffusion[channel_index] += difference * difference; + } + } + + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + before_diffusion[channel_index] /= size * size; + before_diffusion[channel_index] = std::sqrt(before_diffusion[channel_index]); + } + + gil::gray32f_image_t output(32, 32); + auto output_view = gil::view(output); + gil::anisotropic_diffusion(view, 10, 1 / 7.0, 30, output_view); + + double mean_after_diffusion[num_channels] = {0}; + double after_diffusion[num_channels] = {0}; + + for (auto &pixel : output_view) + { + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + mean_after_diffusion[channel_index] += pixel[channel_index]; + } + } + + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + mean_after_diffusion[channel_index] /= size * size; + } + + for (auto &pixel : view) + { + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + double difference = pixel[channel_index] - mean_after_diffusion[channel_index]; + after_diffusion[channel_index] += difference * difference; + } + } + + for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + { + after_diffusion[channel_index] /= size * size; + after_diffusion[channel_index] = std::sqrt(after_diffusion[channel_index]); + BOOST_TEST(before_diffusion[channel_index] > after_diffusion[channel_index]); + } } int main() @@ -70,10 +143,18 @@ int main() for (std::uint32_t seed = 0; seed < 1000; ++seed) { std::cout << "seed: " << seed << '\n'; - std::cout << "gray8 test:\n"; + std::cout << "conservation of heat test:\n" + << "gray8 test:\n"; heat_conservation_test(seed); std::cout << "rgb8 test:\n"; heat_conservation_test(seed); + + std::cout << "convergence to mean test:\n" + << "gray8 test:\n"; + + convergence_to_mean_test(seed); + std::cout << "rgb8 test:\n"; + convergence_to_mean_test(seed); } return boost::report_errors(); From a5a8f098bc96a50e7122ae704707448a0e2410b5 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sun, 21 Jun 2020 01:22:12 +0600 Subject: [PATCH 17/48] 4 way version --- example/anisotropic_diffusion.cpp | 32 +- .../boost/gil/image_processing/diffusion.hpp | 399 ++++-------------- .../anisotropic_diffusion.cpp | 4 +- 3 files changed, 100 insertions(+), 335 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index 6b25ed6fcc..cfee1a05cf 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -20,8 +20,8 @@ namespace gil = boost::gil; -void gray_version(std::string const &input_path, std::string const &output_path, - unsigned int iteration_count, unsigned int kappa) +void gray_version(std::string const& input_path, std::string const& output_path, + unsigned int iteration_count, float kappa, float delta_t) { gil::gray8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -30,13 +30,13 @@ void gray_version(std::string const &input_path, std::string const &output_path, gil::gray32f_image_t gray(input.dimensions()); auto gray_view = gil::view(gray); - gil::transform_pixels(input_view, gray_view, [](const gil::gray8_pixel_t &p) { return p[0]; }); + gil::transform_pixels(input_view, gray_view, [](const gil::gray8_pixel_t& p) { return p[0]; }); double sum_before = 0; gil::for_each_pixel(gray_view, [&sum_before](gil::gray32f_pixel_t p) { sum_before += p[0]; }); gil::gray32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); - gil::anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); + gil::anisotropic_diffusion(gray_view, output_view, iteration_count, delta_t, kappa); double sum_after = 0; gil::for_each_pixel(output_view, [&sum_after](gil::gray32f_pixel_t p) { sum_after += p[0]; }); @@ -51,8 +51,8 @@ void gray_version(std::string const &input_path, std::string const &output_path, << "difference: " << sum_after - sum_before << '\n'; } -void rgb_version(const std::string &input_path, const std::string &output_path, - unsigned int iteration_count, unsigned int kappa) +void rgb_version(const std::string& input_path, const std::string& output_path, + unsigned int iteration_count, float kappa, float delta_t) { gil::rgb8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -61,7 +61,7 @@ void rgb_version(const std::string &input_path, const std::string &output_path, gil::rgb32f_image_t gray(input.dimensions()); auto gray_view = gil::view(gray); - gil::transform_pixels(input_view, gray_view, [](const gil::rgb8_pixel_t &p) { + gil::transform_pixels(input_view, gray_view, [](const gil::rgb8_pixel_t& p) { return gil::rgb32f_pixel_t(p[0], p[1], p[2]); }); double sum_before[3] = {}; @@ -73,9 +73,9 @@ void rgb_version(const std::string &input_path, const std::string &output_path, gil::rgb32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); - gil::anisotropic_diffusion(gray_view, iteration_count, 1 / 7., kappa, output_view); + gil::anisotropic_diffusion(gray_view, output_view, iteration_count, delta_t, kappa); double sum_after[3] = {}; - gil::for_each_pixel(gray_view, [&sum_after](gil::rgb32f_pixel_t p) { + gil::for_each_pixel(output_view, [&sum_after](gil::rgb32f_pixel_t p) { sum_after[0] += p[0]; sum_after[1] += p[1]; sum_after[2] += p[2]; @@ -96,13 +96,14 @@ void rgb_version(const std::string &input_path, const std::string &output_path, << sum_after[1] - sum_before[1] << ", " << sum_after[2] - sum_before[2] << ")\n"; } -int main(int argc, char *argv[]) +int main(int argc, char* argv[]) { - if (argc != 6) + if (argc != 7) { std::cerr << "usage: " << argv[0] << " " - " \n"; + " " + "\n"; return -1; } std::string input_path = argv[1]; @@ -110,14 +111,15 @@ int main(int argc, char *argv[]) std::string colorspace = argv[3]; unsigned int iteration_count = static_cast(std::stoul(argv[4])); - unsigned int kappa = static_cast(std::stoul(argv[5])); + float kappa = std::stof(argv[5]); + float delta_t = std::stof(argv[6]); if (colorspace == "gray") { - gray_version(input_path, output_path, iteration_count, kappa); + gray_version(input_path, output_path, iteration_count, kappa, delta_t); } else if (colorspace == "rgb") { - rgb_version(input_path, output_path, iteration_count, kappa); + rgb_version(input_path, output_path, iteration_count, kappa, delta_t); } else { diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index b4b77bea0c..e82a91ee49 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -6,356 +6,119 @@ // http://www.boost.org/LICENSE_1_0.txt) // #include +#include #include #include +#include #include +#include #include - +#include +#include #include namespace boost { namespace gil { -namespace detail { -namespace direction { -enum -{ - north = 0, - south = 1, - west = 2, - east = 3, - north_east = 4, - south_east = 5, - south_west = 6, - north_west = 7 -}; -} - -// nabla is directed difference in pixels, -// nabla north for example represents intensity -// change in the north direction (up), and nabla south -// represents intensity change in south direction (down). +/// \brief Performs diffusion according to Perona-Malik equation +/// +/// WARNING: Output channel type must be floating point, +/// otherwise there will be loss in accuracy which most +/// probably will lead to incorrect results (input will be unchanged). +/// Anisotropic diffusion is a smoothing algorithm that respects +/// edge boundaries and can work as an edge detector if suitable +/// iteration count is set and grayscale image view is used +/// as an input template -void compute_nabla(InputView view, std::vector const &nabla) +void anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, + typename channel_type::type delta_t, + typename channel_type::type kappa) { - constexpr std::ptrdiff_t input_num_channels = num_channels{}; + using input_pixel_type = typename InputView::value_type; using pixel_type = typename OutputView::value_type; - static_assert(num_channels{} == input_num_channels, - "input and output views must have the same amount of channels"); - const auto width = view.width(); - const auto height = view.height(); - using namespace direction; + using computation_image = image; + using channel_type = typename channel_type::type; const auto zero_pixel = []() { pixel_type pixel; - for (std::ptrdiff_t i = 0; i < input_num_channels; ++i) + for (std::ptrdiff_t i = 0; i < num_channels{}; ++i) { pixel[i] = 0; } return pixel; }(); - auto minus = [](pixel_type lhs, pixel_type rhs) { - for (std::ptrdiff_t i = 0; i < input_num_channels; ++i) + const auto width = input.width(); + const auto height = input.height(); + const point_t dims(width, height); + computation_image result_image(width + 2, height + 2, zero_pixel); + auto result = view(result_image); + auto result_region = subimage_view(result, 1, 1, width, height); + transform_pixels(input, result_region, [](const input_pixel_type& pixel) { + pixel_type converted; + for (std::ptrdiff_t i = 0; i < num_channels{}; ++i) { - lhs[i] -= rhs[i]; + converted[i] = pixel[i]; } + return converted; + }); - return lhs; + const auto minus = [](const pixel_type& lhs, const pixel_type& rhs) { + pixel_type result_pixel; + static_transform(lhs, rhs, result_pixel, std::minus{}); + return result_pixel; + }; + const auto multiply = [](const pixel_type& lhs, const pixel_type& rhs) { + pixel_type result_pixel; + static_transform(lhs, rhs, result_pixel, std::multiplies{}); + return result_pixel; + }; + const auto plus_delta_t = [delta_t](const pixel_type& lhs, const pixel_type& rhs) { + pixel_type result_pixel; + static_transform(lhs, rhs, result_pixel, [delta_t](channel_type ilhs, channel_type irhs) { + // std::cout << ilhs << ' ' << irhs << '\n'; + return (ilhs + irhs) * delta_t; + }); + return result_pixel; }; - // north: x, y - 1 - // south: x, y + 1 - // east: x + 1, y - // west: x - 1, y - // NW: x - 1, y - 1 - // NE: x + 1, y - 1 - // SE: x + 1, y + 1 - // SW: x - 1, y + 1 - - // NW corner - /* - XOOOO - OOOOO - OOOOO - OOOOO - OOOOO - */ - { - const std::ptrdiff_t x = 0; - const std::ptrdiff_t y = 0; - const auto current = view(x, y); - nabla[north](x, y) = zero_pixel; - nabla[south](x, y) = minus(view(x, y + 1), current); - nabla[west](x, y) = zero_pixel; - nabla[east](x, y) = minus(view(x + 1, y), current); - nabla[north_east](x, y) = zero_pixel; - nabla[south_east](x, y) = minus(view(x + 1, y + 1), current); - nabla[south_west](x, y) = zero_pixel; - nabla[north_west](x, y) = zero_pixel; - } - - // NE corner - /* - OOOOX - OOOOO - OOOOO - OOOOO - OOOOO - */ - { - const std::ptrdiff_t x = width - 1; - const std::ptrdiff_t y = 0; - const auto current = view(x, y); - nabla[north](x, y) = zero_pixel; - nabla[south](x, y) = minus(view(x, y + 1), current); - nabla[west](x, y) = minus(view(x - 1, y), current); - nabla[east](x, y) = zero_pixel; - nabla[north_east](x, y) = zero_pixel; - nabla[south_east](x, y) = zero_pixel; - nabla[south_west](x, y) = minus(view(x - 1, y + 1), current); - nabla[north_west](x, y) = zero_pixel; - } - - // SW corner - /* - OOOOO - OOOOO - OOOOO - OOOOO - XOOOO - */ - { - const std::ptrdiff_t x = 0; - const std::ptrdiff_t y = height - 1; - const auto current = view(x, y); - nabla[north](x, y) = minus(view(x, y - 1), current); - nabla[south](x, y) = zero_pixel; - nabla[west](x, y) = zero_pixel; - nabla[east](x, y) = minus(view(x + 1, y), current); - nabla[north_east](x, y) = minus(view(x + 1, y - 1), current); - nabla[south_east](x, y) = zero_pixel; - nabla[south_west](x, y) = zero_pixel; - nabla[north_west](x, y) = zero_pixel; - } - - // SE corner - /* - OOOOO - OOOOO - OOOOO - OOOOO - OOOOX - */ - { - const std::ptrdiff_t x = width - 1; - const std::ptrdiff_t y = height - 1; - const auto current = view(x, y); - nabla[north](x, y) = minus(view(x, y - 1), current); - nabla[south](x, y) = zero_pixel; - nabla[west](x, y) = minus(view(x - 1, y), current); - nabla[east](x, y) = zero_pixel; - nabla[north_east](x, y) = zero_pixel; - nabla[south_east](x, y) = zero_pixel; - nabla[south_west](x, y) = zero_pixel; - nabla[north_west](x, y) = minus(view(x - 1, y - 1), current); - } - - // first and last rows - /* - OXXXO - OOOOO - OOOOO - OOOOO - OXXXO - */ - { - const std::ptrdiff_t y1 = 0; - const std::ptrdiff_t y2 = height - 1; - for (std::ptrdiff_t x = 1; x < width - 1; ++x) - { - const auto current1 = view(x, y1); - const auto current2 = view(x, y2); - - nabla[north](x, y1) = zero_pixel; - nabla[north](x, y2) = minus(view(x, y2 - 1), current2); - nabla[south](x, y1) = minus(view(x, y1 + 1), current1); - nabla[south](x, y2) = zero_pixel; - nabla[west](x, y1) = minus(view(x - 1, y1), current1); - nabla[west](x, y2) = minus(view(x - 1, y2), current2); - nabla[east](x, y1) = minus(view(x + 1, y1), current1); - nabla[east](x, y2) = minus(view(x + 1, y2), current2); - - nabla[north_east](x, y1) = zero_pixel; - nabla[north_east](x, y2) = minus(view(x + 1, y2 - 1), current2); - nabla[south_east](x, y1) = minus(view(x + 1, y1 + 1), current1); - nabla[south_east](x, y2) = zero_pixel; - nabla[south_west](x, y1) = minus(view(x - 1, y1 + 1), current1); - nabla[north_west](x, y1) = zero_pixel; - } - } - - // first and last columns - /* - OOOOO - XOOOX - XOOOX - XOOOX - OOOOO - */ - { - const std::ptrdiff_t x1 = 0; - const std::ptrdiff_t x2 = width - 1; - for (std::ptrdiff_t y = 1; y < height - 1; ++y) - { - const auto current1 = view(x1, y); - const auto current2 = view(x2, y); - nabla[north](x1, y) = minus(view(x1, y - 1), current1); - nabla[north](x2, y) = minus(view(x2, y - 1), current2); - nabla[south](x1, y) = minus(view(x1, y + 1), current1); - nabla[south](x2, y) = minus(view(x2, y + 1), current2); - nabla[west](x1, y) = zero_pixel; - nabla[west](x2, y) = minus(view(x2 - 1, y), current2); - nabla[east](x1, y) = minus(view(x1 + 1, y), current1); - nabla[east](x2, y) = zero_pixel; - - nabla[north_east](x1, y) = minus(view(x1 + 1, y - 1), current1); - nabla[north_west](x2, y) = zero_pixel; - nabla[south_east](x1, y) = minus(view(x1 + 1, y + 1), current1); - nabla[south_east](x2, y) = zero_pixel; - nabla[south_west](x1, y) = zero_pixel; - nabla[south_west](x2, y) = minus(view(x2 - 1, y + 1), current2); - nabla[north_west](x1, y) = zero_pixel; - nabla[north_west](x2, y) = minus(view(x2 - 1, y - 1), current2); - } - } - // middle case - /* - OOOOO - OXXXO - OXXXO - OXXXO - OOOOO - */ - for (std::ptrdiff_t y = 1; y < view.height() - 1; ++y) - { - for (std::ptrdiff_t x = 1; x < view.width() - 1; ++x) - { - for (std::ptrdiff_t channel_index = 0; channel_index < input_num_channels; - ++channel_index) - { - nabla[north](x, y)[channel_index] = - view(x, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[south](x, y)[channel_index] = - view(x, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[west](x, y)[channel_index] = - view(x - 1, y)[channel_index] - view(x, y)[channel_index]; - nabla[east](x, y)[channel_index] = - view(x + 1, y)[channel_index] - view(x, y)[channel_index]; - - nabla[north_east](x, y)[channel_index] = - view(x + 1, y - 1)[channel_index] - view(x, y)[channel_index]; - nabla[south_east](x, y)[channel_index] = - view(x + 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[south_west](x, y)[channel_index] = - view(x - 1, y + 1)[channel_index] - view(x, y)[channel_index]; - nabla[north_west](x, y)[channel_index] = - view(x - 1, y - 1)[channel_index] - view(x, y)[channel_index]; - } - } - } -} - -// formula is c = exp(-(nabla/kappa)^2) -template -void calculate_diffusvity(std::vector nablas, double kappa, - std::vector const &diffusivities) -{ - using pixel_type = typename View::value_type; - using channel_type = typename channel_type::type; - for (std::size_t i = 0; i < nablas.size(); ++i) - { - transform_pixels(nablas[i], diffusivities[i], [kappa](pixel_type p) { - auto op = [kappa](channel_type value) { - value /= static_cast(kappa); - auto result = std::exp(-value * value); - return result; - }; - pixel_type result_pixel; - static_transform(p, result_pixel, op); - return result_pixel; + const auto g = [kappa](const pixel_type& input) { + pixel_type result_pixel(input); + static_transform(input, result_pixel, [kappa](channel_type value) { + value /= kappa; + return std::exp(-value * value); }); - } -} -} // namespace detail -/// \brief Performs diffusion according to Perona-Malik equation -/// -/// WARNING: Output channel type must be floating point, -/// otherwise there will be loss in accuracy which most -/// probably will lead to incorrect results (input will be unchanged). -/// Anisotropic diffusion is a smoothing algorithm that respects -/// edge boundaries and can work as an edge detector if suitable -/// iteration count is set and grayscale image view is used -/// as an input -template -void anisotropic_diffusion(InputView input, unsigned int num_iter, double delta_t, double kappa, - OutputView output) -{ - copy_pixels(input, output); - using element_type = typename OutputView::value_type; - using computation_image_type = image; - using channel_type = typename channel_type::type; - static_assert( - std::is_floating_point::value || std::is_same::value, - "Output view must be floating point, otherwise accuracy loss will nullify diffusion"); + return result_pixel; + }; - for (unsigned int i = 0; i < num_iter; ++i) + for (unsigned int iteration = 0; iteration < num_iter; ++iteration) { - std::vector nabla_images( - 8, computation_image_type(input.dimensions())); - std::vector nabla; - std::transform(nabla_images.begin(), nabla_images.end(), std::back_inserter(nabla), - [](computation_image_type &img) { return gil::view(img); }); - - std::vector diffusivity_images( - 8, computation_image_type(input.dimensions())); - std::vector diffusivity; - std::transform(diffusivity_images.begin(), diffusivity_images.end(), - std::back_inserter(diffusivity), - [](computation_image_type &img) { return gil::view(img); }); - - detail::compute_nabla(output, nabla); - detail::calculate_diffusvity(nabla, kappa, diffusivity); - - channel_type half = channel_type(1.0f / 2); - constexpr std::ptrdiff_t channel_count = num_channels{}; - for (std::ptrdiff_t y = 0; y < output.height(); ++y) + for (std::ptrdiff_t relative_y = 0; relative_y < height; ++relative_y) { - for (std::ptrdiff_t x = 0; x < output.width(); ++x) + for (std::ptrdiff_t relative_x = 0; relative_x < width; ++relative_x) { - for (std::ptrdiff_t channel_index = 0; channel_index < channel_count; - ++channel_index) - { - using namespace detail::direction; - channel_type delta = static_cast( - delta_t * - (diffusivity[north](x, y)[channel_index] * - nabla[north](x, y)[channel_index] + - diffusivity[south](x, y)[channel_index] * - nabla[south](x, y)[channel_index] + - diffusivity[west](x, y)[channel_index] * nabla[west](x, y)[channel_index] + - diffusivity[east](x, y)[channel_index] * nabla[east](x, y)[channel_index] + - half * diffusivity[north_east](x, y)[channel_index] * - nabla[north_east](x, y)[channel_index] + - half * diffusivity[south_east](x, y)[channel_index] * - nabla[south_east](x, y)[channel_index] + - half * diffusivity[south_west](x, y)[channel_index] * - nabla[south_west](x, y)[channel_index] + - half * diffusivity[north_west](x, y)[channel_index] * - nabla[north_west](x, y)[channel_index])); - output(x, y)[channel_index] = output(x, y)[channel_index] + delta; - } + const auto x = relative_x + 1; + const auto y = relative_y + 1; + const auto current = result(x, y); + std::array nabla{ + minus(result(x, y - 1), current), // north + minus(result(x, y + 1), current), // south + minus(result(x - 1, y), current), // west + minus(result(x + 1, y), current) // east + }; + std::array diffusivity{g(nabla[0]), g(nabla[1]), g(nabla[2]), + g(nabla[3])}; + std::array product; + std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), product.begin(), + multiply); + auto sum = + std::accumulate(product.begin(), product.end(), zero_pixel, plus_delta_t); + pixel_type new_pixel; + static_transform(result(x, y), sum, new_pixel, std::plus{}); + result(x, y) = new_pixel; } } } + + copy_pixels(result_region, output); } }} // namespace boost::gil diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index a6f1dc6c46..af0e9c1a7c 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -39,7 +39,7 @@ void heat_conservation_test(std::uint32_t seed) gil::gray32f_image_t output(32, 32); auto output_view = gil::view(output); - gil::anisotropic_diffusion(view, 10, 1 / 7.0, 30, output_view); + gil::anisotropic_diffusion(view, output_view, 10, 1 / 7.0f, 30); double after_diffusion[num_channels] = {0}; for (const auto &pixel : output_view) { @@ -103,7 +103,7 @@ void convergence_to_mean_test(std::uint32_t seed) gil::gray32f_image_t output(32, 32); auto output_view = gil::view(output); - gil::anisotropic_diffusion(view, 10, 1 / 7.0, 30, output_view); + gil::anisotropic_diffusion(view, output_view, 10, 1 / 7.0f, 30); double mean_after_diffusion[num_channels] = {0}; double after_diffusion[num_channels] = {0}; From 5991e21f0b264f3d48529846dafb156a6ffe91ab Mon Sep 17 00:00:00 2001 From: Olzhas Date: Sun, 21 Jun 2020 16:55:19 +0600 Subject: [PATCH 18/48] Reduce warnings and fix bug Previously algorithm did not use a buffer image and updated pixels on the fly, by using those that were already modified, which is incorrect. Also reduced number of warnings in the new files to zero --- .../boost/gil/image_processing/diffusion.hpp | 23 +++++++++---------- .../anisotropic_diffusion.cpp | 12 +++++----- 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index e82a91ee49..d2f2ddefe2 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -38,7 +38,7 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns using channel_type = typename channel_type::type; const auto zero_pixel = []() { pixel_type pixel; - for (std::ptrdiff_t i = 0; i < num_channels{}; ++i) + for (std::size_t i = 0; i < num_channels{}; ++i) { pixel[i] = 0; } @@ -50,10 +50,11 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns const point_t dims(width, height); computation_image result_image(width + 2, height + 2, zero_pixel); auto result = view(result_image); - auto result_region = subimage_view(result, 1, 1, width, height); - transform_pixels(input, result_region, [](const input_pixel_type& pixel) { + computation_image scratch_result_image(width + 2, height + 2, zero_pixel); + auto scratch_result = view(scratch_result_image); + transform_pixels(input, subimage_view(result, 1, 1, width, height), [](const input_pixel_type& pixel) { pixel_type converted; - for (std::ptrdiff_t i = 0; i < num_channels{}; ++i) + for (std::size_t i = 0; i < num_channels{}; ++i) { converted[i] = pixel[i]; } @@ -73,15 +74,14 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns const auto plus_delta_t = [delta_t](const pixel_type& lhs, const pixel_type& rhs) { pixel_type result_pixel; static_transform(lhs, rhs, result_pixel, [delta_t](channel_type ilhs, channel_type irhs) { - // std::cout << ilhs << ' ' << irhs << '\n'; return (ilhs + irhs) * delta_t; }); return result_pixel; }; - const auto g = [kappa](const pixel_type& input) { - pixel_type result_pixel(input); - static_transform(input, result_pixel, [kappa](channel_type value) { + const auto g = [kappa](const pixel_type& input_pixel) { + pixel_type result_pixel(input_pixel); + static_transform(input_pixel, result_pixel, [kappa](channel_type value) { value /= kappa; return std::exp(-value * value); }); @@ -111,14 +111,13 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns multiply); auto sum = std::accumulate(product.begin(), product.end(), zero_pixel, plus_delta_t); - pixel_type new_pixel; - static_transform(result(x, y), sum, new_pixel, std::plus{}); - result(x, y) = new_pixel; + static_transform(result(x, y), sum, scratch_result(x, y), std::plus{}); } } + std::swap(result, scratch_result); } - copy_pixels(result_region, output); + copy_pixels(subimage_view(result, 1, 1, width, height), output); } }} // namespace boost::gil diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index af0e9c1a7c..06e558166b 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -30,7 +30,7 @@ void heat_conservation_test(std::uint32_t seed) double before_diffusion[num_channels] = {0}; for (auto &pixel : view) { - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { pixel[channel_index] = dist(twister); before_diffusion[channel_index] += pixel[channel_index]; @@ -43,7 +43,7 @@ void heat_conservation_test(std::uint32_t seed) double after_diffusion[num_channels] = {0}; for (const auto &pixel : output_view) { - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { after_diffusion[channel_index] += pixel[channel_index]; } @@ -74,7 +74,7 @@ void convergence_to_mean_test(std::uint32_t seed) double before_diffusion[num_channels] = {0}; for (auto &pixel : view) { - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { pixel[channel_index] = dist(twister); mean_before_diffusion[channel_index] += pixel[channel_index]; @@ -88,7 +88,7 @@ void convergence_to_mean_test(std::uint32_t seed) for (auto &pixel : view) { - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { double difference = pixel[channel_index] - mean_before_diffusion[channel_index]; before_diffusion[channel_index] += difference * difference; @@ -110,7 +110,7 @@ void convergence_to_mean_test(std::uint32_t seed) for (auto &pixel : output_view) { - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { mean_after_diffusion[channel_index] += pixel[channel_index]; } @@ -123,7 +123,7 @@ void convergence_to_mean_test(std::uint32_t seed) for (auto &pixel : view) { - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) + for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { double difference = pixel[channel_index] - mean_after_diffusion[channel_index]; after_diffusion[channel_index] += difference * difference; From ce3bee4668fa000dec2153871639d42ceea2703b Mon Sep 17 00:00:00 2001 From: Olzhas Date: Thu, 25 Jun 2020 01:00:55 +0600 Subject: [PATCH 19/48] Use transform Rather than manually writing out each entry, use std::transform --- include/boost/gil/image_processing/diffusion.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index d2f2ddefe2..7db138e18a 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -104,8 +104,8 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns minus(result(x - 1, y), current), // west minus(result(x + 1, y), current) // east }; - std::array diffusivity{g(nabla[0]), g(nabla[1]), g(nabla[2]), - g(nabla[3])}; + std::array diffusivity; + std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), g); std::array product; std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), product.begin(), multiply); From 43bc160faef069888a597b7db79ba93e1af7754d Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Fri, 26 Jun 2020 00:23:20 +0600 Subject: [PATCH 20/48] Squash the damn accumulate bug The accumulate part was wrong, it multiplied by delta_t on every sum, which is wrong --- .../boost/gil/image_processing/diffusion.hpp | 30 +++++++++++-------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 7db138e18a..28c01c4d7d 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -52,14 +52,15 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns auto result = view(result_image); computation_image scratch_result_image(width + 2, height + 2, zero_pixel); auto scratch_result = view(scratch_result_image); - transform_pixels(input, subimage_view(result, 1, 1, width, height), [](const input_pixel_type& pixel) { - pixel_type converted; - for (std::size_t i = 0; i < num_channels{}; ++i) - { - converted[i] = pixel[i]; - } - return converted; - }); + transform_pixels(input, subimage_view(result, 1, 1, width, height), + [](const input_pixel_type& pixel) { + pixel_type converted; + for (std::size_t i = 0; i < num_channels{}; ++i) + { + converted[i] = pixel[i]; + } + return converted; + }); const auto minus = [](const pixel_type& lhs, const pixel_type& rhs) { pixel_type result_pixel; @@ -73,9 +74,8 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns }; const auto plus_delta_t = [delta_t](const pixel_type& lhs, const pixel_type& rhs) { pixel_type result_pixel; - static_transform(lhs, rhs, result_pixel, [delta_t](channel_type ilhs, channel_type irhs) { - return (ilhs + irhs) * delta_t; - }); + static_transform(lhs, rhs, result_pixel, + [delta_t](channel_type ilhs, channel_type irhs) { return (ilhs + irhs); }); return result_pixel; }; @@ -111,10 +111,14 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns multiply); auto sum = std::accumulate(product.begin(), product.end(), zero_pixel, plus_delta_t); - static_transform(result(x, y), sum, scratch_result(x, y), std::plus{}); + static_transform(sum, sum, + [delta_t](channel_type value) { return value * delta_t; }); + static_transform(result(x, y), sum, scratch_result(x, y), + std::plus{}); } } - std::swap(result, scratch_result); + using std::swap; + swap(result, scratch_result); } copy_pixels(subimage_view(result, 1, 1, width, height), output); From 149b07ddc25724a68e111db860fbe52cba152606 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Fri, 26 Jun 2020 02:52:31 +0600 Subject: [PATCH 21/48] Use 8 way nabla compute This is just different discretization of Laplace operator. https://en.wikipedia.org/wiki/Discrete_Laplace_operator --- .../boost/gil/image_processing/diffusion.hpp | 30 ++-- .../anisotropic_diffusion.cpp | 170 +++++++++--------- 2 files changed, 105 insertions(+), 95 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 28c01c4d7d..979350e722 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -72,7 +72,7 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns static_transform(lhs, rhs, result_pixel, std::multiplies{}); return result_pixel; }; - const auto plus_delta_t = [delta_t](const pixel_type& lhs, const pixel_type& rhs) { + const auto plus = [delta_t](const pixel_type& lhs, const pixel_type& rhs) { pixel_type result_pixel; static_transform(lhs, rhs, result_pixel, [delta_t](channel_type ilhs, channel_type irhs) { return (ilhs + irhs); }); @@ -98,19 +98,29 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns const auto x = relative_x + 1; const auto y = relative_y + 1; const auto current = result(x, y); - std::array nabla{ - minus(result(x, y - 1), current), // north - minus(result(x, y + 1), current), // south - minus(result(x - 1, y), current), // west - minus(result(x + 1, y), current) // east + constexpr std::size_t point_count = 8; + std::array nabla{ + minus(result(x, y - 1), current), // north + minus(result(x, y + 1), current), // south + minus(result(x - 1, y), current), // west + minus(result(x + 1, y), current), // east + minus(result(x - 1, y - 1), current), // north west + minus(result(x + 1, y - 1), current), // north east + minus(result(x - 1, y + 1), current), // south west + minus(result(x + 1, y + 1), current) // south east }; - std::array diffusivity; + std::array diffusivity; std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), g); - std::array product; + std::array product; std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), product.begin(), multiply); - auto sum = - std::accumulate(product.begin(), product.end(), zero_pixel, plus_delta_t); + + for (std::size_t i = 4; i < point_count; ++i) + { + static_transform(product[i], product[i], + [](channel_type value) { return value / 2; }); + } + auto sum = std::accumulate(product.begin(), product.end(), zero_pixel, plus); static_transform(sum, sum, [delta_t](channel_type value) { return value * delta_t; }); static_transform(result(x, y), sum, scratch_result(x, y), diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 06e558166b..907528bc58 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -28,7 +28,7 @@ void heat_conservation_test(std::uint32_t seed) auto view = gil::view(image); constexpr std::ptrdiff_t num_channels = gil::num_channels::value; double before_diffusion[num_channels] = {0}; - for (auto &pixel : view) + for (auto& pixel : view) { for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { @@ -39,9 +39,9 @@ void heat_conservation_test(std::uint32_t seed) gil::gray32f_image_t output(32, 32); auto output_view = gil::view(output); - gil::anisotropic_diffusion(view, output_view, 10, 1 / 7.0f, 30); + gil::anisotropic_diffusion(view, output_view, 10, 1 / 6.0f, 5); double after_diffusion[num_channels] = {0}; - for (const auto &pixel : output_view) + for (const auto& pixel : output_view) { for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { @@ -60,83 +60,83 @@ void heat_conservation_test(std::uint32_t seed) std::cout << '\n'; } -template -void convergence_to_mean_test(std::uint32_t seed) -{ - std::mt19937 twister(seed); - std::uniform_int_distribution dist; - - const std::size_t size = 32; - ImageType image(size, size); - auto view = gil::view(image); - constexpr std::ptrdiff_t num_channels = gil::num_channels::value; - double mean_before_diffusion[num_channels] = {0}; - double before_diffusion[num_channels] = {0}; - for (auto &pixel : view) - { - for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - pixel[channel_index] = dist(twister); - mean_before_diffusion[channel_index] += pixel[channel_index]; - } - } - - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - mean_before_diffusion[channel_index] /= size * size; - } - - for (auto &pixel : view) - { - for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - double difference = pixel[channel_index] - mean_before_diffusion[channel_index]; - before_diffusion[channel_index] += difference * difference; - } - } - - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - before_diffusion[channel_index] /= size * size; - before_diffusion[channel_index] = std::sqrt(before_diffusion[channel_index]); - } - - gil::gray32f_image_t output(32, 32); - auto output_view = gil::view(output); - gil::anisotropic_diffusion(view, output_view, 10, 1 / 7.0f, 30); - - double mean_after_diffusion[num_channels] = {0}; - double after_diffusion[num_channels] = {0}; - - for (auto &pixel : output_view) - { - for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - mean_after_diffusion[channel_index] += pixel[channel_index]; - } - } - - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - mean_after_diffusion[channel_index] /= size * size; - } - - for (auto &pixel : view) - { - for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - double difference = pixel[channel_index] - mean_after_diffusion[channel_index]; - after_diffusion[channel_index] += difference * difference; - } - } - - for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) - { - after_diffusion[channel_index] /= size * size; - after_diffusion[channel_index] = std::sqrt(after_diffusion[channel_index]); - BOOST_TEST(before_diffusion[channel_index] > after_diffusion[channel_index]); - } -} +// template +// void convergence_to_mean_test(std::uint32_t seed) +// { +// std::mt19937 twister(seed); +// std::uniform_int_distribution dist; + +// const std::size_t size = 32; +// ImageType image(size, size); +// auto view = gil::view(image); +// constexpr std::ptrdiff_t num_channels = gil::num_channels::value; +// double mean_before_diffusion[num_channels] = {0}; +// double before_diffusion[num_channels] = {0}; +// for (auto& pixel : view) +// { +// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// pixel[channel_index] = dist(twister); +// mean_before_diffusion[channel_index] += pixel[channel_index]; +// } +// } + +// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// mean_before_diffusion[channel_index] /= size * size; +// } + +// for (auto& pixel : view) +// { +// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// double difference = pixel[channel_index] - mean_before_diffusion[channel_index]; +// before_diffusion[channel_index] += difference * difference; +// } +// } + +// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// before_diffusion[channel_index] /= size * size; +// before_diffusion[channel_index] = std::sqrt(before_diffusion[channel_index]); +// } + +// gil::gray32f_image_t output(32, 32); +// auto output_view = gil::view(output); +// gil::anisotropic_diffusion(view, output_view, 10, 1 / 7.0f, 30); + +// double mean_after_diffusion[num_channels] = {0}; +// double after_diffusion[num_channels] = {0}; + +// for (auto& pixel : output_view) +// { +// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// mean_after_diffusion[channel_index] += pixel[channel_index]; +// } +// } + +// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// mean_after_diffusion[channel_index] /= size * size; +// } + +// for (auto& pixel : view) +// { +// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// double difference = pixel[channel_index] - mean_after_diffusion[channel_index]; +// after_diffusion[channel_index] += difference * difference; +// } +// } + +// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) +// { +// after_diffusion[channel_index] /= size * size; +// after_diffusion[channel_index] = std::sqrt(after_diffusion[channel_index]); +// BOOST_TEST(before_diffusion[channel_index] > after_diffusion[channel_index]); +// } +// } int main() { @@ -149,12 +149,12 @@ int main() std::cout << "rgb8 test:\n"; heat_conservation_test(seed); - std::cout << "convergence to mean test:\n" - << "gray8 test:\n"; + // std::cout << "convergence to mean test:\n" + // << "gray8 test:\n"; - convergence_to_mean_test(seed); - std::cout << "rgb8 test:\n"; - convergence_to_mean_test(seed); + // convergence_to_mean_test(seed); + // std::cout << "rgb8 test:\n"; + // convergence_to_mean_test(seed); } return boost::report_errors(); From 4fe247045874ce457eb802213c157f2b60be29a7 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 27 Jun 2020 10:44:49 +0600 Subject: [PATCH 22/48] Move to strategy --- example/anisotropic_diffusion.cpp | 4 +- .../boost/gil/image_processing/diffusion.hpp | 155 ++++++++++-------- .../anisotropic_diffusion.cpp | 4 +- 3 files changed, 95 insertions(+), 68 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index cfee1a05cf..1c9200e689 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -36,7 +36,7 @@ void gray_version(std::string const& input_path, std::string const& output_path, gil::gray32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); - gil::anisotropic_diffusion(gray_view, output_view, iteration_count, delta_t, kappa); + gil::anisotropic_diffusion(gray_view, output_view, iteration_count, {kappa, delta_t}); double sum_after = 0; gil::for_each_pixel(output_view, [&sum_after](gil::gray32f_pixel_t p) { sum_after += p[0]; }); @@ -73,7 +73,7 @@ void rgb_version(const std::string& input_path, const std::string& output_path, gil::rgb32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); - gil::anisotropic_diffusion(gray_view, output_view, iteration_count, delta_t, kappa); + gil::anisotropic_diffusion(gray_view, output_view, iteration_count, {kappa, delta_t}); double sum_after[3] = {}; gil::for_each_pixel(output_view, [&sum_after](gil::rgb32f_pixel_t p) { sum_after[0] += p[0]; diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 979350e722..015744eb33 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -18,6 +18,88 @@ #include namespace boost { namespace gil { + +template +struct diffusion_strategy +{ + using channel_type = typename channel_type::type; + channel_type kappa; + channel_type delta_t; + template + PixelType operator()(View padded_view, point_t origin) + { + // seems like lambda capture doesn't work for member variables in C++11 + const auto local_delta_t = delta_t; + const auto local_kappa = kappa; + + const auto minus = [](const PixelType& lhs, const PixelType& rhs) { + PixelType result_pixel; + static_transform(lhs, rhs, result_pixel, std::minus{}); + return result_pixel; + }; + const auto multiply = [](const PixelType& lhs, const PixelType& rhs) { + PixelType result_pixel; + static_transform(lhs, rhs, result_pixel, std::multiplies{}); + return result_pixel; + }; + const auto plus = [](const PixelType& lhs, const PixelType& rhs) { + PixelType result_pixel; + static_transform(lhs, rhs, result_pixel, + [](channel_type ilhs, channel_type irhs) { return (ilhs + irhs); }); + return result_pixel; + }; + + const auto g = [local_kappa](const PixelType& input_pixel) { + PixelType result_pixel(input_pixel); + static_transform(input_pixel, result_pixel, [local_kappa](channel_type value) { + value /= local_kappa; + return std::exp(-value * value); + }); + + return result_pixel; + }; + const auto x = origin.x; + const auto y = origin.y; + const auto current = padded_view(x, y); + constexpr std::size_t point_count = 8; + std::array nabla{ + minus(padded_view(x, y - 1), current), // north + minus(padded_view(x, y + 1), current), // south + minus(padded_view(x - 1, y), current), // west + minus(padded_view(x + 1, y), current), // east + minus(padded_view(x - 1, y - 1), current), // north west + minus(padded_view(x + 1, y - 1), current), // north east + minus(padded_view(x - 1, y + 1), current), // south west + minus(padded_view(x + 1, y + 1), current) // south east + }; + std::array diffusivity; + std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), g); + std::array product; + std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), product.begin(), multiply); + + for (std::size_t i = 4; i < point_count; ++i) + { + static_transform(product[i], product[i], [](channel_type value) { return value / 2; }); + } + const auto zero_pixel = []() { + PixelType pixel; + for (std::size_t i = 0; i < num_channels{}; ++i) + { + pixel[i] = 0; + } + + return pixel; + }(); + auto sum = std::accumulate(product.begin(), product.end(), zero_pixel, plus); + static_transform(sum, sum, + [local_delta_t](channel_type value) { return value * local_delta_t; }); + PixelType result_pixel; + static_transform(padded_view(x, y), sum, result_pixel, std::plus{}); + + return result_pixel; + } +}; + /// \brief Performs diffusion according to Perona-Malik equation /// /// WARNING: Output channel type must be floating point, @@ -27,15 +109,17 @@ namespace boost { namespace gil { /// edge boundaries and can work as an edge detector if suitable /// iteration count is set and grayscale image view is used /// as an input -template +template > void anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, - typename channel_type::type delta_t, - typename channel_type::type kappa) + ComputationStrategy strategy) { using input_pixel_type = typename InputView::value_type; using pixel_type = typename OutputView::value_type; using computation_image = image; - using channel_type = typename channel_type::type; + const auto width = input.width(); + const auto height = input.height(); + const point_t dims(width, height); const auto zero_pixel = []() { pixel_type pixel; for (std::size_t i = 0; i < num_channels{}; ++i) @@ -45,9 +129,6 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns return pixel; }(); - const auto width = input.width(); - const auto height = input.height(); - const point_t dims(width, height); computation_image result_image(width + 2, height + 2, zero_pixel); auto result = view(result_image); computation_image scratch_result_image(width + 2, height + 2, zero_pixel); @@ -62,69 +143,15 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns return converted; }); - const auto minus = [](const pixel_type& lhs, const pixel_type& rhs) { - pixel_type result_pixel; - static_transform(lhs, rhs, result_pixel, std::minus{}); - return result_pixel; - }; - const auto multiply = [](const pixel_type& lhs, const pixel_type& rhs) { - pixel_type result_pixel; - static_transform(lhs, rhs, result_pixel, std::multiplies{}); - return result_pixel; - }; - const auto plus = [delta_t](const pixel_type& lhs, const pixel_type& rhs) { - pixel_type result_pixel; - static_transform(lhs, rhs, result_pixel, - [delta_t](channel_type ilhs, channel_type irhs) { return (ilhs + irhs); }); - return result_pixel; - }; - - const auto g = [kappa](const pixel_type& input_pixel) { - pixel_type result_pixel(input_pixel); - static_transform(input_pixel, result_pixel, [kappa](channel_type value) { - value /= kappa; - return std::exp(-value * value); - }); - - return result_pixel; - }; - for (unsigned int iteration = 0; iteration < num_iter; ++iteration) { for (std::ptrdiff_t relative_y = 0; relative_y < height; ++relative_y) { for (std::ptrdiff_t relative_x = 0; relative_x < width; ++relative_x) { - const auto x = relative_x + 1; - const auto y = relative_y + 1; - const auto current = result(x, y); - constexpr std::size_t point_count = 8; - std::array nabla{ - minus(result(x, y - 1), current), // north - minus(result(x, y + 1), current), // south - minus(result(x - 1, y), current), // west - minus(result(x + 1, y), current), // east - minus(result(x - 1, y - 1), current), // north west - minus(result(x + 1, y - 1), current), // north east - minus(result(x - 1, y + 1), current), // south west - minus(result(x + 1, y + 1), current) // south east - }; - std::array diffusivity; - std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), g); - std::array product; - std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), product.begin(), - multiply); - - for (std::size_t i = 4; i < point_count; ++i) - { - static_transform(product[i], product[i], - [](channel_type value) { return value / 2; }); - } - auto sum = std::accumulate(product.begin(), product.end(), zero_pixel, plus); - static_transform(sum, sum, - [delta_t](channel_type value) { return value * delta_t; }); - static_transform(result(x, y), sum, scratch_result(x, y), - std::plus{}); + auto x = relative_x + 1; + auto y = relative_y + 1; + scratch_result(x, y) = strategy(result, {x, y}); } } using std::swap; diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 907528bc58..cb1c4861da 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -39,7 +39,7 @@ void heat_conservation_test(std::uint32_t seed) gil::gray32f_image_t output(32, 32); auto output_view = gil::view(output); - gil::anisotropic_diffusion(view, output_view, 10, 1 / 6.0f, 5); + gil::anisotropic_diffusion(view, output_view, 10, {5, 1 / 6.0f}); double after_diffusion[num_channels] = {0}; for (const auto& pixel : output_view) { @@ -55,7 +55,7 @@ void heat_conservation_test(std::uint32_t seed) std::abs(after_diffusion[channel_index] - before_diffusion[channel_index]) / after_diffusion[channel_index] * 100.0; std::cout << percentage << ' '; - BOOST_TEST(percentage < 10.0); + BOOST_TEST(percentage < 8.5); } std::cout << '\n'; } From 8628367a9496594edf1f4dbb2524489b224f3ab5 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 30 Jun 2020 01:34:53 +0600 Subject: [PATCH 23/48] opencv strategy iteration --- .../boost/gil/image_processing/diffusion.hpp | 29 +++++++++++-------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 015744eb33..6deb91e65b 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -35,6 +35,7 @@ struct diffusion_strategy const auto minus = [](const PixelType& lhs, const PixelType& rhs) { PixelType result_pixel; static_transform(lhs, rhs, result_pixel, std::minus{}); + return result_pixel; }; const auto multiply = [](const PixelType& lhs, const PixelType& rhs) { @@ -51,10 +52,16 @@ struct diffusion_strategy const auto g = [local_kappa](const PixelType& input_pixel) { PixelType result_pixel(input_pixel); - static_transform(input_pixel, result_pixel, [local_kappa](channel_type value) { - value /= local_kappa; - return std::exp(-value * value); - }); + channel_type sum = 0; + static_for_each(result_pixel, + [&sum](channel_type channel_value) { sum += std::abs(channel_value); }); + channel_type sigma = local_kappa * num_channels{} * channel_type(255); + channel_type exponent = (sum * sum) / (sigma * sigma); + static_fill(result_pixel, std::exp(exponent)); + // static_transform(scratch_pixel, result_pixel, [local_kappa](channel_type value) { + // value /= local_kappa * num_channels{} * channel_type(255); + // return std::exp(-value * value); + // }); return result_pixel; }; @@ -77,16 +84,14 @@ struct diffusion_strategy std::array product; std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), product.begin(), multiply); - for (std::size_t i = 4; i < point_count; ++i) - { - static_transform(product[i], product[i], [](channel_type value) { return value / 2; }); - } + // for (std::size_t i = 4; i < point_count; ++i) + // { + // static_transform(product[i], product[i], [](channel_type value) { return value / 2; + // }); + // } const auto zero_pixel = []() { PixelType pixel; - for (std::size_t i = 0; i < num_channels{}; ++i) - { - pixel[i] = 0; - } + static_fill(pixel, channel_type(0)); return pixel; }(); From db8a32613af7a032373d0cc94c9dc132a3957e89 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 6 Jul 2020 03:56:16 +0600 Subject: [PATCH 24/48] Segment into function objects This commit implements correct diffusion and separates the whole process into multiple pieces --- example/anisotropic_diffusion.cpp | 20 +- .../boost/gil/image_processing/diffusion.hpp | 305 +++++++++++++----- 2 files changed, 230 insertions(+), 95 deletions(-) diff --git a/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp index 1c9200e689..03defb7566 100644 --- a/example/anisotropic_diffusion.cpp +++ b/example/anisotropic_diffusion.cpp @@ -21,7 +21,7 @@ namespace gil = boost::gil; void gray_version(std::string const& input_path, std::string const& output_path, - unsigned int iteration_count, float kappa, float delta_t) + unsigned int iteration_count, float kappa) { gil::gray8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -36,7 +36,8 @@ void gray_version(std::string const& input_path, std::string const& output_path, gil::gray32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); - gil::anisotropic_diffusion(gray_view, output_view, iteration_count, {kappa, delta_t}); + // gil::anisotropic_diffusion(gray_view, output_view, iteration_count, {kappa, delta_t}); + gil::default_anisotropic_diffusion(gray_view, output_view, iteration_count, kappa); double sum_after = 0; gil::for_each_pixel(output_view, [&sum_after](gil::gray32f_pixel_t p) { sum_after += p[0]; }); @@ -52,7 +53,7 @@ void gray_version(std::string const& input_path, std::string const& output_path, } void rgb_version(const std::string& input_path, const std::string& output_path, - unsigned int iteration_count, float kappa, float delta_t) + unsigned int iteration_count, float kappa) { gil::rgb8_image_t input; gil::read_image(input_path, input, gil::png_tag{}); @@ -73,7 +74,8 @@ void rgb_version(const std::string& input_path, const std::string& output_path, gil::rgb32f_image_t output(gray.dimensions()); auto output_view = gil::view(output); - gil::anisotropic_diffusion(gray_view, output_view, iteration_count, {kappa, delta_t}); + // gil::anisotropic_diffusion(gray_view, output_view, iteration_count, {kappa, delta_t}); + gil::default_anisotropic_diffusion(gray_view, output_view, iteration_count, kappa); double sum_after[3] = {}; gil::for_each_pixel(output_view, [&sum_after](gil::rgb32f_pixel_t p) { sum_after[0] += p[0]; @@ -98,12 +100,11 @@ void rgb_version(const std::string& input_path, const std::string& output_path, int main(int argc, char* argv[]) { - if (argc != 7) + if (argc != 6) { std::cerr << "usage: " << argv[0] << " " - " " - "\n"; + " \n"; return -1; } std::string input_path = argv[1]; @@ -112,14 +113,13 @@ int main(int argc, char* argv[]) unsigned int iteration_count = static_cast(std::stoul(argv[4])); float kappa = std::stof(argv[5]); - float delta_t = std::stof(argv[6]); if (colorspace == "gray") { - gray_version(input_path, output_path, iteration_count, kappa, delta_t); + gray_version(input_path, output_path, iteration_count, kappa); } else if (colorspace == "rgb") { - rgb_version(input_path, output_path, iteration_count, kappa, delta_t); + rgb_version(input_path, output_path, iteration_count, kappa); } else { diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 6deb91e65b..2ff5b89a96 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -5,6 +5,7 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // +#include "boost/gil/detail/math.hpp" #include #include #include @@ -18,93 +19,214 @@ #include namespace boost { namespace gil { +namespace diffusion { +struct perona_malik_diffusivity +{ + double kappa; + template + Pixel operator()(Pixel input) + { + using channel_type = typename channel_type::type; + // C++11 doesn't seem to capture members + auto local_kappa = kappa; + static_transform(input, input, [local_kappa](channel_type value) { + value /= local_kappa; + return std::exp(-std::abs(value)); + }); + + return input; + } +}; + +struct gaussian_diffusivity +{ + double kappa; + template + Pixel operator()(Pixel input) + { + using channel_type = typename channel_type::type; + // C++11 doesn't seem to capture members + auto local_kappa = kappa; + static_transform(input, input, [local_kappa](channel_type value) { + value /= local_kappa; + return std::exp(-value * value); + }); + + return input; + } +}; + +struct wide_regions_diffusivity +{ + double kappa; + template + Pixel operator()(Pixel input) + { + using channel_type = typename channel_type::type; + // C++11 doesn't seem to capture members + auto local_kappa = kappa; + static_transform(input, input, [local_kappa](channel_type value) { + value /= local_kappa; + return 1.0 / (1.0 + value * value); + }); + + return input; + } +}; +struct more_wide_regions_diffusivity +{ + double kappa; + template + Pixel operator()(Pixel input) + { + using channel_type = typename channel_type::type; + // C++11 doesn't seem to capture members + auto local_kappa = kappa; + static_transform(input, input, [local_kappa](channel_type value) { + value /= local_kappa; + return 1.0 / std::sqrt((1.0 + value * value)); + }); + + return input; + } +}; +} // namespace diffusion + +namespace laplace_function { template -struct diffusion_strategy +using stencil_type = std::array; + +struct stencil_4points +{ + double delta_t = 0.25; + + template + stencil_type compute_laplace(SubImageView view, + point_t origin) + { + stencil_type stencil; + auto out = stencil.begin(); + auto current = view(origin); + using channel_type = typename channel_type::type; + std::array offsets{point_t{0, -1}, point_t{-1, 0}, point_t{+1, 0}, + point_t{0, 1}}; + for (auto offset : offsets) + { + ++out; // skip 45 degree values; + static_transform(view(origin.x + offset.x, origin.y + offset.y), current, *out++, + std::minus{}); + } + } + + template + Pixel reduce(const stencil_type& stencil) + { + auto first = stencil.begin(); + auto last = stencil.end(); + using channel_type = typename channel_type::value_type; + auto result = []() { + Pixel zero_pixel; + static_fill(zero_pixel, channel_type(0)); + return zero_pixel; + }(); + + for (; first != last; ++first) + { + static_transform(result, *first, result, std::plus{}); + } + Pixel delta_t_pixel; + static_fill(delta_t_pixel, delta_t); + static_transform(result, delta_t_pixel, result, std::multiplies{}); + + return result; + } +}; + +struct stencil_8points_standard { - using channel_type = typename channel_type::type; - channel_type kappa; - channel_type delta_t; - template - PixelType operator()(View padded_view, point_t origin) + double delta_t = 0.125; + + template + stencil_type compute_laplace(SubImageView view, + point_t origin) { - // seems like lambda capture doesn't work for member variables in C++11 - const auto local_delta_t = delta_t; - const auto local_kappa = kappa; - - const auto minus = [](const PixelType& lhs, const PixelType& rhs) { - PixelType result_pixel; - static_transform(lhs, rhs, result_pixel, std::minus{}); - - return result_pixel; - }; - const auto multiply = [](const PixelType& lhs, const PixelType& rhs) { - PixelType result_pixel; - static_transform(lhs, rhs, result_pixel, std::multiplies{}); - return result_pixel; - }; - const auto plus = [](const PixelType& lhs, const PixelType& rhs) { - PixelType result_pixel; - static_transform(lhs, rhs, result_pixel, - [](channel_type ilhs, channel_type irhs) { return (ilhs + irhs); }); - return result_pixel; - }; - - const auto g = [local_kappa](const PixelType& input_pixel) { - PixelType result_pixel(input_pixel); - channel_type sum = 0; - static_for_each(result_pixel, - [&sum](channel_type channel_value) { sum += std::abs(channel_value); }); - channel_type sigma = local_kappa * num_channels{} * channel_type(255); - channel_type exponent = (sum * sum) / (sigma * sigma); - static_fill(result_pixel, std::exp(exponent)); - // static_transform(scratch_pixel, result_pixel, [local_kappa](channel_type value) { - // value /= local_kappa * num_channels{} * channel_type(255); - // return std::exp(-value * value); - // }); - - return result_pixel; - }; - const auto x = origin.x; - const auto y = origin.y; - const auto current = padded_view(x, y); - constexpr std::size_t point_count = 8; - std::array nabla{ - minus(padded_view(x, y - 1), current), // north - minus(padded_view(x, y + 1), current), // south - minus(padded_view(x - 1, y), current), // west - minus(padded_view(x + 1, y), current), // east - minus(padded_view(x - 1, y - 1), current), // north west - minus(padded_view(x + 1, y - 1), current), // north east - minus(padded_view(x - 1, y + 1), current), // south west - minus(padded_view(x + 1, y + 1), current) // south east - }; - std::array diffusivity; - std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), g); - std::array product; - std::transform(nabla.begin(), nabla.end(), diffusivity.begin(), product.begin(), multiply); - - // for (std::size_t i = 4; i < point_count; ++i) - // { - // static_transform(product[i], product[i], [](channel_type value) { return value / 2; - // }); - // } - const auto zero_pixel = []() { - PixelType pixel; - static_fill(pixel, channel_type(0)); - - return pixel; + stencil_type stencil; + auto out = stencil.begin(); + auto current = view(origin); + using channel_type = typename channel_type::type; + std::array offsets{point_t{-1, -1}, point_t{0, -1}, point_t{+1, -1}, + point_t{-1, 0}, point_t{+1, 0}, point_t{-1, +1}, + point_t{0, +1}, point_t{+1, +1}}; + for (auto offset : offsets) + { + static_transform(view(origin.x + offset.x, origin.y + offset.y), current, *out++, + std::minus{}); + } + + return stencil; + } + + template + Pixel reduce(const stencil_type& stencil) + { + auto first = stencil.begin(); + auto last = stencil.end(); + using channel_type = typename channel_type::type; + auto result = []() { + Pixel zero_pixel; + static_fill(zero_pixel, channel_type(0)); + return zero_pixel; }(); - auto sum = std::accumulate(product.begin(), product.end(), zero_pixel, plus); - static_transform(sum, sum, - [local_delta_t](channel_type value) { return value * local_delta_t; }); - PixelType result_pixel; - static_transform(padded_view(x, y), sum, result_pixel, std::plus{}); + auto half = std::next(first, 4); + for (; first != half; ++first) + { + static_transform(result, *first, result, std::plus{}); + } - return result_pixel; + for (first = half; first != last; ++first) + { + Pixel half_pixel; + static_fill(half_pixel, channel_type(1 / 2.0)); + static_transform(*first, half_pixel, half_pixel, std::multiplies{}); + static_transform(result, half_pixel, result, std::plus{}); + } + + Pixel delta_t_pixel; + static_fill(delta_t_pixel, delta_t); + static_transform(result, delta_t_pixel, result, std::multiplies{}); + + return result; + } +}; +} // namespace laplace_function + +namespace brightness_function { +using laplace_function::stencil_type; +struct identity +{ + template + stencil_type operator()(const stencil_type& stencil) + { + return stencil; } }; +// TODO: Figure out how to implement color gradient brightness, as it +// seems to need dx and dy using sobel or scharr kernels + +// TODO: Implement luminance based brightness function + +} // namespace brightness_function + +template +void default_anisotropic_diffusion(const InputView& input, const OutputView& output, + unsigned int num_iter, double kappa) +{ + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_8points_standard{}, + brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); +} + /// \brief Performs diffusion according to Perona-Malik equation /// /// WARNING: Output channel type must be floating point, @@ -115,22 +237,23 @@ struct diffusion_strategy /// iteration count is set and grayscale image view is used /// as an input template > + typename LaplaceStrategy = laplace_function::stencil_8points_standard, + typename BrightnessFunction = brightness_function::identity, + typename DiffusivityFunction = diffusion::gaussian_diffusivity> void anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, - ComputationStrategy strategy) + LaplaceStrategy laplace, BrightnessFunction brightness, + DiffusivityFunction diffusivity) { using input_pixel_type = typename InputView::value_type; using pixel_type = typename OutputView::value_type; + using channel_type = typename channel_type::type; using computation_image = image; const auto width = input.width(); const auto height = input.height(); const point_t dims(width, height); const auto zero_pixel = []() { pixel_type pixel; - for (std::size_t i = 0; i < num_channels{}; ++i) - { - pixel[i] = 0; - } + static_fill(pixel, static_cast(0)); return pixel; }(); @@ -156,7 +279,19 @@ void anisotropic_diffusion(const InputView& input, const OutputView& output, uns { auto x = relative_x + 1; auto y = relative_y + 1; - scratch_result(x, y) = strategy(result, {x, y}); + auto stencil = laplace.compute_laplace(result, point_t(x, y)); + auto brightness_stencil = brightness(stencil); + laplace_function::stencil_type diffusivity_stencil; + std::transform(brightness_stencil.begin(), brightness_stencil.end(), + diffusivity_stencil.begin(), diffusivity); + laplace_function::stencil_type product_stencil; + std::transform(stencil.begin(), stencil.end(), diffusivity_stencil.begin(), + product_stencil.begin(), [](pixel_type lhs, pixel_type rhs) { + static_transform(lhs, rhs, lhs, std::multiplies{}); + return lhs; + }); + static_transform(result(x, y), laplace.reduce(product_stencil), + scratch_result(x, y), std::plus{}); } } using std::swap; From 4469ff2c3022a7d1d11b794d51e1552f58a6695b Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 6 Jul 2020 03:58:32 +0600 Subject: [PATCH 25/48] Fix test --- .../anisotropic_diffusion.cpp | 80 +------------------ 1 file changed, 1 insertion(+), 79 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index cb1c4861da..06392b70cd 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -39,7 +39,7 @@ void heat_conservation_test(std::uint32_t seed) gil::gray32f_image_t output(32, 32); auto output_view = gil::view(output); - gil::anisotropic_diffusion(view, output_view, 10, {5, 1 / 6.0f}); + gil::default_anisotropic_diffusion(view, output_view, 10, 5); double after_diffusion[num_channels] = {0}; for (const auto& pixel : output_view) { @@ -60,84 +60,6 @@ void heat_conservation_test(std::uint32_t seed) std::cout << '\n'; } -// template -// void convergence_to_mean_test(std::uint32_t seed) -// { -// std::mt19937 twister(seed); -// std::uniform_int_distribution dist; - -// const std::size_t size = 32; -// ImageType image(size, size); -// auto view = gil::view(image); -// constexpr std::ptrdiff_t num_channels = gil::num_channels::value; -// double mean_before_diffusion[num_channels] = {0}; -// double before_diffusion[num_channels] = {0}; -// for (auto& pixel : view) -// { -// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// pixel[channel_index] = dist(twister); -// mean_before_diffusion[channel_index] += pixel[channel_index]; -// } -// } - -// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// mean_before_diffusion[channel_index] /= size * size; -// } - -// for (auto& pixel : view) -// { -// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// double difference = pixel[channel_index] - mean_before_diffusion[channel_index]; -// before_diffusion[channel_index] += difference * difference; -// } -// } - -// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// before_diffusion[channel_index] /= size * size; -// before_diffusion[channel_index] = std::sqrt(before_diffusion[channel_index]); -// } - -// gil::gray32f_image_t output(32, 32); -// auto output_view = gil::view(output); -// gil::anisotropic_diffusion(view, output_view, 10, 1 / 7.0f, 30); - -// double mean_after_diffusion[num_channels] = {0}; -// double after_diffusion[num_channels] = {0}; - -// for (auto& pixel : output_view) -// { -// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// mean_after_diffusion[channel_index] += pixel[channel_index]; -// } -// } - -// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// mean_after_diffusion[channel_index] /= size * size; -// } - -// for (auto& pixel : view) -// { -// for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// double difference = pixel[channel_index] - mean_after_diffusion[channel_index]; -// after_diffusion[channel_index] += difference * difference; -// } -// } - -// for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) -// { -// after_diffusion[channel_index] /= size * size; -// after_diffusion[channel_index] = std::sqrt(after_diffusion[channel_index]); -// BOOST_TEST(before_diffusion[channel_index] > after_diffusion[channel_index]); -// } -// } - int main() { for (std::uint32_t seed = 0; seed < 1000; ++seed) From 946a073cf4827144e4fc9fd3d8b4668968c23904 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 6 Jul 2020 04:01:04 +0600 Subject: [PATCH 26/48] Fix distribution type bug --- test/core/image_processing/anisotropic_diffusion.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 06392b70cd..e1337a9329 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -22,7 +22,7 @@ template void heat_conservation_test(std::uint32_t seed) { std::mt19937 twister(seed); - std::uniform_int_distribution dist; + std::uniform_int_distribution dist; ImageType image(32, 32); auto view = gil::view(image); @@ -32,7 +32,7 @@ void heat_conservation_test(std::uint32_t seed) { for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { - pixel[channel_index] = dist(twister); + pixel[channel_index] = static_cast(dist(twister)); before_diffusion[channel_index] += pixel[channel_index]; } } From fcf118d6989bec80bbc2979d246ad842f2c6db19 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 6 Jul 2020 04:08:41 +0600 Subject: [PATCH 27/48] Fix distribution range --- test/core/image_processing/anisotropic_diffusion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index e1337a9329..fb203fc568 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -22,7 +22,7 @@ template void heat_conservation_test(std::uint32_t seed) { std::mt19937 twister(seed); - std::uniform_int_distribution dist; + std::uniform_int_distribution dist(0, std::numeric_limits::max()); ImageType image(32, 32); auto view = gil::view(image); From cf0ae18d189f87bed4212c2f5e3fab37589434ff Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Mon, 6 Jul 2020 04:25:07 +0600 Subject: [PATCH 28/48] Drastically decrease error threshold Found a bug inside tests, now it computes correctly, allowing to reduce error threshold drastically --- .../image_processing/anisotropic_diffusion.cpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index fb203fc568..2fe93a7c4b 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -18,7 +18,7 @@ namespace gil = boost::gil; -template +template void heat_conservation_test(std::uint32_t seed) { std::mt19937 twister(seed); @@ -37,7 +37,7 @@ void heat_conservation_test(std::uint32_t seed) } } - gil::gray32f_image_t output(32, 32); + OutputImageType output(32, 32); auto output_view = gil::view(output); gil::default_anisotropic_diffusion(view, output_view, 10, 5); double after_diffusion[num_channels] = {0}; @@ -52,10 +52,13 @@ void heat_conservation_test(std::uint32_t seed) for (std::ptrdiff_t channel_index = 0; channel_index < num_channels; ++channel_index) { const auto percentage = - std::abs(after_diffusion[channel_index] - before_diffusion[channel_index]) / - after_diffusion[channel_index] * 100.0; + before_diffusion[channel_index] != 0 + ? std::abs(after_diffusion[channel_index] - before_diffusion[channel_index]) / + after_diffusion[channel_index] * 100.0 + : std::abs(after_diffusion[channel_index] - before_diffusion[channel_index]) / + after_diffusion[channel_index] * 100.0; std::cout << percentage << ' '; - BOOST_TEST(percentage < 8.5); + BOOST_TEST(percentage < 1); } std::cout << '\n'; } @@ -67,9 +70,9 @@ int main() std::cout << "seed: " << seed << '\n'; std::cout << "conservation of heat test:\n" << "gray8 test:\n"; - heat_conservation_test(seed); + heat_conservation_test(seed); std::cout << "rgb8 test:\n"; - heat_conservation_test(seed); + heat_conservation_test(seed); // std::cout << "convergence to mean test:\n" // << "gray8 test:\n"; From 1ce4264be47db6021fca9b422d01320f76bac28b Mon Sep 17 00:00:00 2001 From: Olzhas Date: Mon, 6 Jul 2020 13:08:36 +0600 Subject: [PATCH 29/48] Surround cout with ifdef --- .../image_processing/anisotropic_diffusion.cpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 2fe93a7c4b..b2d5143c92 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -57,29 +57,30 @@ void heat_conservation_test(std::uint32_t seed) after_diffusion[channel_index] * 100.0 : std::abs(after_diffusion[channel_index] - before_diffusion[channel_index]) / after_diffusion[channel_index] * 100.0; +#ifdef BOOST_GIL_TEST_DEBUG std::cout << percentage << ' '; +#endif BOOST_TEST(percentage < 1); } +#ifdef BOOST_GIL_TEST_DEBUG std::cout << '\n'; +#endif } int main() { for (std::uint32_t seed = 0; seed < 1000; ++seed) { +#ifdef BOOST_GIL_TEST_DEBUG std::cout << "seed: " << seed << '\n'; std::cout << "conservation of heat test:\n" << "gray8 test:\n"; +#endif heat_conservation_test(seed); +#ifdef BOOST_GIL_TEST_DEBUG std::cout << "rgb8 test:\n"; +#endif heat_conservation_test(seed); - - // std::cout << "convergence to mean test:\n" - // << "gray8 test:\n"; - - // convergence_to_mean_test(seed); - // std::cout << "rgb8 test:\n"; - // convergence_to_mean_test(seed); } return boost::report_errors(); From e6a6548b38b21de7c4fc5192f154183fd5ebb160 Mon Sep 17 00:00:00 2001 From: Olzhas Date: Tue, 7 Jul 2020 01:56:56 +0600 Subject: [PATCH 30/48] Matlab version --- .../boost/gil/image_processing/diffusion.hpp | 35 +++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 2ff5b89a96..0af0c31213 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -219,6 +219,41 @@ struct identity } // namespace brightness_function +enum class matlab_connectivity { + minimal, + maximal +}; + +enum class matlab_conduction_method { + exponential, + quadratic +}; + +template +void matlab_anisotropic_diffusion(const InputView& input, const OutputView& output, + unsigned int num_iter, double kappa, matlab_connectivity connectivity, matlab_conduction_method conduction_method) +{ + if (connectivity == matlab_connectivity::minimal) { + if (conduction_method == matlab_conduction_method::exponential) { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); + } else if (conduction_method == matlab_conduction_method::quadratic) { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); + } else { + throw std::logic_error("unhandled conduction method found"); + } + } else if (connectivity == matlab_connectivity::maximal) { + if (conduction_method == matlab_conduction_method::exponential) { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); + } else if (conduction_method == matlab_conduction_method::quadratic) { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); + } else { + throw std::logic_error("unhandled conduction method found"); + } + } else { + throw std::logic_error("unhandled connectivity found"); + } +} + template void default_anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, double kappa) From c6534fac24debc3d5db92dbfe7d614114a00f8ad Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 7 Jul 2020 02:08:16 +0600 Subject: [PATCH 31/48] Add classic version --- .../boost/gil/image_processing/diffusion.hpp | 68 ++++++++++++++----- 1 file changed, 52 insertions(+), 16 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 0af0c31213..7c045a9438 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -219,37 +219,73 @@ struct identity } // namespace brightness_function -enum class matlab_connectivity { +enum class matlab_connectivity +{ minimal, maximal }; -enum class matlab_conduction_method { +enum class matlab_conduction_method +{ exponential, quadratic }; +template +void classic_anisotropic_diffusion(const InputView& input, const OutputView& output, + unsigned int num_iter, double kappa) +{ + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + brightness_function::identity{}, + diffusion::perona_malik_diffusivity{kappa}); +} + template void matlab_anisotropic_diffusion(const InputView& input, const OutputView& output, - unsigned int num_iter, double kappa, matlab_connectivity connectivity, matlab_conduction_method conduction_method) + unsigned int num_iter, double kappa, + matlab_connectivity connectivity, + matlab_conduction_method conduction_method) { - if (connectivity == matlab_connectivity::minimal) { - if (conduction_method == matlab_conduction_method::exponential) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); - } else if (conduction_method == matlab_conduction_method::quadratic) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); - } else { + if (connectivity == matlab_connectivity::minimal) + { + if (conduction_method == matlab_conduction_method::exponential) + { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + brightness_function::identity{}, + diffusion::gaussian_diffusivity{kappa}); + } + else if (conduction_method == matlab_conduction_method::quadratic) + { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + brightness_function::identity{}, + diffusion::gaussian_diffusivity{kappa}); + } + else + { throw std::logic_error("unhandled conduction method found"); } - } else if (connectivity == matlab_connectivity::maximal) { - if (conduction_method == matlab_conduction_method::exponential) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); - } else if (conduction_method == matlab_conduction_method::quadratic) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); - } else { + } + else if (connectivity == matlab_connectivity::maximal) + { + if (conduction_method == matlab_conduction_method::exponential) + { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + brightness_function::identity{}, + diffusion::gaussian_diffusivity{kappa}); + } + else if (conduction_method == matlab_conduction_method::quadratic) + { + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + brightness_function::identity{}, + diffusion::gaussian_diffusivity{kappa}); + } + else + { throw std::logic_error("unhandled conduction method found"); } - } else { + } + else + { throw std::logic_error("unhandled connectivity found"); } } From 4fadad0e2a2129b627f67d741069d2e9c043b203 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 7 Jul 2020 02:43:24 +0600 Subject: [PATCH 32/48] Add diffusivity tests --- .../boost/gil/image_processing/diffusion.hpp | 15 +++++++- .../anisotropic_diffusion.cpp | 35 +++++++++++++++++++ 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 7c045a9438..43d4051449 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -215,7 +215,20 @@ struct identity // TODO: Figure out how to implement color gradient brightness, as it // seems to need dx and dy using sobel or scharr kernels -// TODO: Implement luminance based brightness function +struct rgb_luminance +{ + using pixel_type = rgb32f_pixel_t; + stencil_type operator()(const stencil_type& stencil) + { + stencil_type output; + std::transform(stencil.begin(), stencil.end(), output.begin(), [](const pixel_type& pixel) { + float32_t luminance = 0.2126f * pixel[0] + 0.7152f * pixel[1] + 0.0722f * pixel[2]; + pixel_type result_pixel; + static_fill(result_pixel, luminance); + }); + return output; + } +}; } // namespace brightness_function diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index b2d5143c92..3b93bed3f3 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -5,6 +5,7 @@ // 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 @@ -18,6 +19,21 @@ namespace gil = boost::gil; +template +void diffusion_function_check(DiffusionFunction diffusion) +{ + using limits = std::numeric_limits; + using channel_type = typename gil::channel_type::type; + for (channel_type value = limits::min(); value <= limits::max(); ++value) + { + PixelType pixel; + gil::static_fill(pixel, value); + auto diffusivity_value = diffusion(pixel); + gil::static_for_each(diffusivity_value, + [](channel_type value) { BOOST_TEST(0 <= value && value <= 1.0); }); + } +} + template void heat_conservation_test(std::uint32_t seed) { @@ -83,5 +99,24 @@ int main() heat_conservation_test(seed); } + for (double kappa = 5; kappa <= 70; ++kappa) + { + diffusion_function_check( + gil::diffusion::perona_malik_diffusivity{kappa}); + diffusion_function_check(gil::diffusion::gaussian_diffusivity{kappa}); + diffusion_function_check( + gil::diffusion::wide_regions_diffusivity{kappa}); + diffusion_function_check( + gil::diffusion::more_wide_regions_diffusivity{kappa}); + + diffusion_function_check( + gil::diffusion::perona_malik_diffusivity{kappa}); + diffusion_function_check(gil::diffusion::gaussian_diffusivity{kappa}); + diffusion_function_check( + gil::diffusion::wide_regions_diffusivity{kappa}); + diffusion_function_check( + gil::diffusion::more_wide_regions_diffusivity{kappa}); + } + return boost::report_errors(); } From 370d1b37248c8e64e63c926f5e74c340096a1349 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 7 Jul 2020 03:05:23 +0600 Subject: [PATCH 33/48] Add identity brightness tests --- .../boost/gil/image_processing/diffusion.hpp | 1 + .../anisotropic_diffusion.cpp | 21 +++++++++++++++++++ 2 files changed, 22 insertions(+) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 43d4051449..571ec8076a 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -225,6 +225,7 @@ struct rgb_luminance float32_t luminance = 0.2126f * pixel[0] + 0.7152f * pixel[1] + 0.0722f * pixel[2]; pixel_type result_pixel; static_fill(result_pixel, luminance); + return result_pixel; }); return output; } diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 3b93bed3f3..23e66f246d 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -34,6 +34,25 @@ void diffusion_function_check(DiffusionFunction diffusion) } } +void brightness_function_test() +{ + std::vector rgb_pixels{ + gil::rgb32f_pixel_t(0, 11, 14), gil::rgb32f_pixel_t(2, 117, 200), + gil::rgb32f_pixel_t(223, 2, 180), gil::rgb32f_pixel_t(250, 254, 100), + gil::rgb32f_pixel_t(255, 255, 255), + }; + for (const auto& pixel : rgb_pixels) + { + boost::gil::laplace_function::stencil_type stencil; + std::fill(stencil.begin(), stencil.end(), pixel); + auto brightness_stencil = boost::gil::brightness_function::identity{}(stencil); + for (const auto& brightness_pixel : brightness_stencil) + { + BOOST_TEST(pixel == brightness_pixel); + } + } +} + template void heat_conservation_test(std::uint32_t seed) { @@ -118,5 +137,7 @@ int main() gil::diffusion::more_wide_regions_diffusivity{kappa}); } + brightness_function_test(); + return boost::report_errors(); } From acb240b173c373b01e51ad477696a45dc1323181 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 7 Jul 2020 03:25:21 +0600 Subject: [PATCH 34/48] Tested all brightness functions --- .../anisotropic_diffusion.cpp | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 23e66f246d..cbe014cf23 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -14,6 +14,7 @@ #include +#include #include #include @@ -51,6 +52,38 @@ void brightness_function_test() BOOST_TEST(pixel == brightness_pixel); } } + + std::vector corresponding_luminance_values{8.878f, 98.5436f, 61.8362f, 242.0308f, 255}; + std::size_t index = 0; + for (const auto& pixel : rgb_pixels) + { + boost::gil::laplace_function::stencil_type stencil; + std::fill(stencil.begin(), stencil.end(), pixel); + auto brightness_stencil = boost::gil::brightness_function::rgb_luminance{}(stencil); + for (const auto& brightness_pixel : brightness_stencil) + { + gil::static_for_each(brightness_pixel, [&corresponding_luminance_values, + index](boost::gil::float32_t value) { + BOOST_TEST(std::abs(value - corresponding_luminance_values[index]) < 1.0f); + }); + } + ++index; + } + + std::vector gray_pixels{ + gil::gray32f_pixel_t(11), gil::gray32f_pixel_t(0), gil::gray32f_pixel_t(255), + gil::gray32f_pixel_t(123), gil::gray32f_pixel_t(17), + }; + for (const auto& pixel : gray_pixels) + { + boost::gil::laplace_function::stencil_type stencil; + std::fill(stencil.begin(), stencil.end(), pixel); + auto brightness_stencil = boost::gil::brightness_function::identity{}(stencil); + for (const auto& brightness_pixel : brightness_stencil) + { + BOOST_TEST(pixel == brightness_pixel); + } + } } template From cd03a9a4c89178f43b29246a7a10fa0cbd4d5716 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Tue, 7 Jul 2020 03:47:19 +0600 Subject: [PATCH 35/48] Add Jamfile entries --- example/Jamfile | 1 + test/core/image_processing/Jamfile | 1 + 2 files changed, 2 insertions(+) diff --git a/example/Jamfile b/example/Jamfile index fc33ce319c..02bb33163c 100644 --- a/example/Jamfile +++ b/example/Jamfile @@ -21,6 +21,7 @@ project local sources = adaptive_threshold.cpp affine.cpp + anisotropic_diffusion.cpp convolution.cpp convolve2d.cpp dynamic_image.cpp diff --git a/test/core/image_processing/Jamfile b/test/core/image_processing/Jamfile index d9593f0793..f2f547f560 100644 --- a/test/core/image_processing/Jamfile +++ b/test/core/image_processing/Jamfile @@ -19,3 +19,4 @@ run hessian.cpp ; run sobel_scharr.cpp ; run box_filter.cpp ; run median_filter.cpp ; +run anisotropic_diffusion ; From 59396633ec9bd9a778be47eaad9fcd1630f04102 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Thu, 9 Jul 2020 15:41:23 +0600 Subject: [PATCH 36/48] Failing laplace function tests --- .../boost/gil/image_processing/diffusion.hpp | 2 +- .../anisotropic_diffusion.cpp | 142 ++++++++++++++++++ 2 files changed, 143 insertions(+), 1 deletion(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 571ec8076a..cf40ec7c88 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -124,7 +124,7 @@ struct stencil_4points { auto first = stencil.begin(); auto last = stencil.end(); - using channel_type = typename channel_type::value_type; + using channel_type = typename channel_type::type; auto result = []() { Pixel zero_pixel; static_fill(zero_pixel, channel_type(0)); diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index cbe014cf23..7e491edbf9 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -5,6 +5,7 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // +#include "boost/gil/algorithm.hpp" #include #include #include @@ -135,6 +136,145 @@ void heat_conservation_test(std::uint32_t seed) #endif } +template +void test_stencil_pixels(const gil::laplace_function::stencil_type& stencil, + const PixelType& reference) +{ + for (const auto& stencil_entry : stencil) + { + BOOST_TEST(stencil_entry == reference); + + gil::static_for_each(stencil_entry, + [](gil::float32_t value) { std::cout << value << ' '; }); + std::cout << '\n'; + } +} + +template +void test_stencil(const gil::laplace_function::stencil_type& stencil, + const std::array& expected_values) +{ + for (std::size_t i = 0; i < 8; ++i) + { + PixelType expected_pixel; + gil::static_fill(expected_pixel, expected_values[i]); + BOOST_TEST(stencil[i] == expected_pixel); + } +} + +void laplace_functions_test() +{ + // sanity checks + auto zero_rgb_pixel = []() { + gil::rgb32f_pixel_t zero_pixel; + gil::static_fill(zero_pixel, 0); + return zero_pixel; + }(); + auto zero_gray_pixel = []() { + gil::gray32f_pixel_t zero_pixel; + gil::static_fill(zero_pixel, 0); + return zero_pixel; + }(); + + auto image_size = gil::point_t(16, 16); + gil::rgb32f_image_t rgb_image(image_size, zero_rgb_pixel); + gil::gray32f_image_t gray_image(image_size, zero_gray_pixel); + auto rgb = gil::view(rgb_image); + auto gray = gil::view(gray_image); + + // gil::for_each_pixel(gil::view(zero_rgb_image), [](const gil::rgb32f_pixel_t& pixel) { + + // gil::static_for_each(pixel, [](gil::float32_t value) { + // BOOST_TEST(value == 0); + // }); + // }); + // gil::for_each_pixel(gil::view(zero_rgb_image), [](const gil::gray8_pixel_t& pixel) { + // gil::static_for_each(pixel, [](gil::float32_t value) { + // BOOST_TEST(value == 0); + // }); + // }); + auto _4way = gil::laplace_function::stencil_4points{}; + auto _8way = gil::laplace_function::stencil_8points_standard{}; + for (std::ptrdiff_t y = 1; y < image_size.y - 1; ++y) + { + for (std::ptrdiff_t x = 1; x < image_size.x - 1; ++x) + { + auto rgb_4way_stencil = _4way.compute_laplace(rgb, {x, y}); + auto gray_4way_stencil = _4way.compute_laplace(gray, {x, y}); + + auto rgb_8way_stencil = _8way.compute_laplace(rgb, {x, y}); + auto gray_8way_stencil = _8way.compute_laplace(gray, {x, y}); + + // test_stencil_pixels(rgb_4way_stencil, zero_rgb_pixel); + test_stencil_pixels(rgb_8way_stencil, zero_rgb_pixel); + // test_stencil_pixels(gray_4way_stencil, zero_gray_pixel); + test_stencil_pixels(gray_8way_stencil, zero_gray_pixel); + } + } + + // a predefined case + // 5 11 2 + // 17 25 58 + // 90 31 40 + + using rgb_pixel = gil::rgb32f_pixel_t; + rgb(1, 1) = rgb_pixel(5, 5, 5); + rgb(1, 2) = rgb_pixel(11, 11, 11); + rgb(1, 3) = rgb_pixel(2, 2, 2); + rgb(2, 1) = rgb_pixel(17, 17, 17); + rgb(2, 2) = rgb_pixel(25, 25, 25); + rgb(2, 3) = rgb_pixel(58, 58, 58); + rgb(3, 1) = rgb_pixel(90, 90, 90); + rgb(3, 2) = rgb_pixel(31, 31, 31); + rgb(3, 3) = rgb_pixel(40, 40, 40); + + using gray_pixel = gil::gray32f_pixel_t; + gray(1, 1) = gray_pixel(5); + gray(1, 2) = gray_pixel(11); + gray(1, 3) = gray_pixel(2); + gray(2, 1) = gray_pixel(17); + gray(2, 2) = gray_pixel(25); + gray(2, 3) = gray_pixel(58); + gray(3, 1) = gray_pixel(90); + gray(3, 2) = gray_pixel(31); + gray(3, 3) = gray_pixel(40); + + // 4 way stencil should be + // 0 -14 0 + // -8 33 + // 0 6 0 + + auto test_point = gil::point_t(2, 2); + std::array _4way_expected{0, -14, 0, -8, 33, 0, 6, 0}; + auto rgb_4way = _4way.compute_laplace(rgb, test_point); + // test_stencil(rgb_4way, _4way_expected); + auto gray_4way = _4way.compute_laplace(gray, test_point); + // test_stencil(gray_4way, _4way_expected); + + // 8 way stencil should be + // -25 -14 -23 + // -8 33 + // 65 6 15 + + std::array _8way_expected{-25, -14, -23, -8, 33, 65, 6, 15}; + auto rgb_8way = _8way.compute_laplace(rgb, test_point); + test_stencil(rgb_8way, _8way_expected); + auto gray_8way = _8way.compute_laplace(gray, test_point); + test_stencil(gray_8way, _8way_expected); + + // reduce result for 4 way should be -14 - 8 + 6 + 33 = 17 + auto rgb_reduced_4way = _4way.reduce(rgb_4way); + gil::static_for_each(rgb_reduced_4way, [](gil::float32_t value) { BOOST_TEST(value == 17); }); + auto gray_reduced_4way = _4way.reduce(gray_4way); + gil::static_for_each(gray_reduced_4way, [](gil::float32_t value) { BOOST_TEST(value == 17); }); + + // reduce result for 8 way should be (-25-23+65+15)*0.5+(-14-8+33+6) = 16+17=33 + auto rgb_reduced_8way = _8way.reduce(rgb_8way); + gil::static_for_each(rgb_reduced_8way, [](gil::float32_t value) { BOOST_TEST(value == 33); }); + auto gray_reduced_8way = _4way.reduce(gray_8way); + gil::static_for_each(gray_reduced_8way, [](gil::float32_t value) { BOOST_TEST(value == 33); }); +} + int main() { for (std::uint32_t seed = 0; seed < 1000; ++seed) @@ -172,5 +312,7 @@ int main() brightness_function_test(); + laplace_functions_test(); + return boost::report_errors(); } From 8c60cb40be30b54b4bf885cd51d009d9218a18e3 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 11 Jul 2020 02:02:05 +0600 Subject: [PATCH 37/48] Squash all of the bugs --- .../boost/gil/image_processing/diffusion.hpp | 27 +++--- .../anisotropic_diffusion.cpp | 85 ++++++++----------- 2 files changed, 48 insertions(+), 64 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index cf40ec7c88..5d1aa655dc 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -109,14 +109,17 @@ struct stencil_4points auto out = stencil.begin(); auto current = view(origin); using channel_type = typename channel_type::type; - std::array offsets{point_t{0, -1}, point_t{-1, 0}, point_t{+1, 0}, - point_t{0, 1}}; + std::array offsets{point_t{0, -1}, point_t{+1, 0}, point_t{0, +1}, + point_t{-1, 0}}; + typename SubImageView::value_type zero_pixel; + static_fill(zero_pixel, 0); for (auto offset : offsets) { - ++out; // skip 45 degree values; + *out++ = zero_pixel; static_transform(view(origin.x + offset.x, origin.y + offset.y), current, *out++, std::minus{}); } + return stencil; } template @@ -155,9 +158,9 @@ struct stencil_8points_standard auto out = stencil.begin(); auto current = view(origin); using channel_type = typename channel_type::type; - std::array offsets{point_t{-1, -1}, point_t{0, -1}, point_t{+1, -1}, - point_t{-1, 0}, point_t{+1, 0}, point_t{-1, +1}, - point_t{0, +1}, point_t{+1, +1}}; + std::array offsets{point_t{-1, -1}, point_t{0, -1}, point_t{+1, -1}, + point_t{+1, 0}, point_t{+1, +1}, point_t{0, +1}, + point_t{-1, +1}, point_t{-1, 0}}; for (auto offset : offsets) { static_transform(view(origin.x + offset.x, origin.y + offset.y), current, *out++, @@ -170,25 +173,23 @@ struct stencil_8points_standard template Pixel reduce(const stencil_type& stencil) { - auto first = stencil.begin(); - auto last = stencil.end(); using channel_type = typename channel_type::type; auto result = []() { Pixel zero_pixel; static_fill(zero_pixel, channel_type(0)); return zero_pixel; }(); - auto half = std::next(first, 4); - for (; first != half; ++first) + for (std::size_t index : {1u, 3u, 5u, 7u}) { - static_transform(result, *first, result, std::plus{}); + static_transform(result, stencil[index], result, std::plus{}); } - for (first = half; first != last; ++first) + for (std::size_t index : {0u, 2u, 4u, 6u}) { Pixel half_pixel; static_fill(half_pixel, channel_type(1 / 2.0)); - static_transform(*first, half_pixel, half_pixel, std::multiplies{}); + static_transform(stencil[index], half_pixel, half_pixel, + std::multiplies{}); static_transform(result, half_pixel, result, std::plus{}); } diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 7e491edbf9..1084e567a5 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -31,8 +31,9 @@ void diffusion_function_check(DiffusionFunction diffusion) PixelType pixel; gil::static_fill(pixel, value); auto diffusivity_value = diffusion(pixel); - gil::static_for_each(diffusivity_value, - [](channel_type value) { BOOST_TEST(0 <= value && value <= 1.0); }); + gil::static_for_each(diffusivity_value, [](channel_type channel_value) { + BOOST_TEST(0 <= channel_value && channel_value <= 1.0); + }); } } @@ -143,10 +144,6 @@ void test_stencil_pixels(const gil::laplace_function::stencil_type& s for (const auto& stencil_entry : stencil) { BOOST_TEST(stencil_entry == reference); - - gil::static_for_each(stencil_entry, - [](gil::float32_t value) { std::cout << value << ' '; }); - std::cout << '\n'; } } @@ -182,17 +179,6 @@ void laplace_functions_test() auto rgb = gil::view(rgb_image); auto gray = gil::view(gray_image); - // gil::for_each_pixel(gil::view(zero_rgb_image), [](const gil::rgb32f_pixel_t& pixel) { - - // gil::static_for_each(pixel, [](gil::float32_t value) { - // BOOST_TEST(value == 0); - // }); - // }); - // gil::for_each_pixel(gil::view(zero_rgb_image), [](const gil::gray8_pixel_t& pixel) { - // gil::static_for_each(pixel, [](gil::float32_t value) { - // BOOST_TEST(value == 0); - // }); - // }); auto _4way = gil::laplace_function::stencil_4points{}; auto _8way = gil::laplace_function::stencil_8points_standard{}; for (std::ptrdiff_t y = 1; y < image_size.y - 1; ++y) @@ -205,9 +191,9 @@ void laplace_functions_test() auto rgb_8way_stencil = _8way.compute_laplace(rgb, {x, y}); auto gray_8way_stencil = _8way.compute_laplace(gray, {x, y}); - // test_stencil_pixels(rgb_4way_stencil, zero_rgb_pixel); + test_stencil_pixels(rgb_4way_stencil, zero_rgb_pixel); test_stencil_pixels(rgb_8way_stencil, zero_rgb_pixel); - // test_stencil_pixels(gray_4way_stencil, zero_gray_pixel); + test_stencil_pixels(gray_4way_stencil, zero_gray_pixel); test_stencil_pixels(gray_8way_stencil, zero_gray_pixel); } } @@ -219,24 +205,24 @@ void laplace_functions_test() using rgb_pixel = gil::rgb32f_pixel_t; rgb(1, 1) = rgb_pixel(5, 5, 5); - rgb(1, 2) = rgb_pixel(11, 11, 11); - rgb(1, 3) = rgb_pixel(2, 2, 2); - rgb(2, 1) = rgb_pixel(17, 17, 17); + rgb(2, 1) = rgb_pixel(11, 11, 11); + rgb(3, 1) = rgb_pixel(2, 2, 2); + rgb(1, 2) = rgb_pixel(17, 17, 17); rgb(2, 2) = rgb_pixel(25, 25, 25); - rgb(2, 3) = rgb_pixel(58, 58, 58); - rgb(3, 1) = rgb_pixel(90, 90, 90); - rgb(3, 2) = rgb_pixel(31, 31, 31); + rgb(3, 2) = rgb_pixel(58, 58, 58); + rgb(1, 3) = rgb_pixel(90, 90, 90); + rgb(2, 3) = rgb_pixel(31, 31, 31); rgb(3, 3) = rgb_pixel(40, 40, 40); using gray_pixel = gil::gray32f_pixel_t; gray(1, 1) = gray_pixel(5); - gray(1, 2) = gray_pixel(11); - gray(1, 3) = gray_pixel(2); - gray(2, 1) = gray_pixel(17); + gray(2, 1) = gray_pixel(11); + gray(3, 1) = gray_pixel(2); + gray(1, 2) = gray_pixel(17); gray(2, 2) = gray_pixel(25); - gray(2, 3) = gray_pixel(58); - gray(3, 1) = gray_pixel(90); - gray(3, 2) = gray_pixel(31); + gray(3, 2) = gray_pixel(58); + gray(1, 3) = gray_pixel(90); + gray(2, 3) = gray_pixel(31); gray(3, 3) = gray_pixel(40); // 4 way stencil should be @@ -245,49 +231,46 @@ void laplace_functions_test() // 0 6 0 auto test_point = gil::point_t(2, 2); - std::array _4way_expected{0, -14, 0, -8, 33, 0, 6, 0}; + std::array _4way_expected{0, -14, 0, 33, 0, 6, 0, -8}; auto rgb_4way = _4way.compute_laplace(rgb, test_point); - // test_stencil(rgb_4way, _4way_expected); + test_stencil(rgb_4way, _4way_expected); auto gray_4way = _4way.compute_laplace(gray, test_point); - // test_stencil(gray_4way, _4way_expected); + test_stencil(gray_4way, _4way_expected); // 8 way stencil should be - // -25 -14 -23 + // -20 -14 -23 // -8 33 // 65 6 15 - std::array _8way_expected{-25, -14, -23, -8, 33, 65, 6, 15}; + std::array _8way_expected{-20, -14, -23, 33, 15, 6, 65, -8}; auto rgb_8way = _8way.compute_laplace(rgb, test_point); test_stencil(rgb_8way, _8way_expected); auto gray_8way = _8way.compute_laplace(gray, test_point); test_stencil(gray_8way, _8way_expected); - // reduce result for 4 way should be -14 - 8 + 6 + 33 = 17 + // reduce result for 4 way should be (-14 - 8 + 6 + 33) * 0.25 = 4.25 auto rgb_reduced_4way = _4way.reduce(rgb_4way); - gil::static_for_each(rgb_reduced_4way, [](gil::float32_t value) { BOOST_TEST(value == 17); }); + gil::static_for_each(rgb_reduced_4way, + [](gil::float32_t value) { BOOST_TEST(value == 4.25f); }); auto gray_reduced_4way = _4way.reduce(gray_4way); - gil::static_for_each(gray_reduced_4way, [](gil::float32_t value) { BOOST_TEST(value == 17); }); + gil::static_for_each(gray_reduced_4way, + [](gil::float32_t value) { BOOST_TEST(value == 4.25f); }); - // reduce result for 8 way should be (-25-23+65+15)*0.5+(-14-8+33+6) = 16+17=33 + // reduce result for 8 way should be ((-20-23+15+65)*0.5+(-14+33+6-8)) * 0.125 = (18.5+17) * + // 0.125 = 4.4375 auto rgb_reduced_8way = _8way.reduce(rgb_8way); - gil::static_for_each(rgb_reduced_8way, [](gil::float32_t value) { BOOST_TEST(value == 33); }); - auto gray_reduced_8way = _4way.reduce(gray_8way); - gil::static_for_each(gray_reduced_8way, [](gil::float32_t value) { BOOST_TEST(value == 33); }); + gil::static_for_each(rgb_reduced_8way, + [](gil::float32_t value) { BOOST_TEST(value == 4.4375); }); + auto gray_reduced_8way = _8way.reduce(gray_8way); + gil::static_for_each(gray_reduced_8way, + [](gil::float32_t value) { BOOST_TEST(value == 4.4375); }); } int main() { for (std::uint32_t seed = 0; seed < 1000; ++seed) { -#ifdef BOOST_GIL_TEST_DEBUG - std::cout << "seed: " << seed << '\n'; - std::cout << "conservation of heat test:\n" - << "gray8 test:\n"; -#endif heat_conservation_test(seed); -#ifdef BOOST_GIL_TEST_DEBUG - std::cout << "rgb8 test:\n"; -#endif heat_conservation_test(seed); } From abf7be739b5bfe0886d2605fc9b7708d6cb298ab Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 11 Jul 2020 02:04:46 +0600 Subject: [PATCH 38/48] Rename laplace stencils Laplace stencils are now the same as in mathematical notation --- .../boost/gil/image_processing/diffusion.hpp | 18 +++++++++--------- .../image_processing/anisotropic_diffusion.cpp | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 5d1aa655dc..c5c2e7d2c3 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -97,7 +97,7 @@ namespace laplace_function { template using stencil_type = std::array; -struct stencil_4points +struct stencil_5points { double delta_t = 0.25; @@ -146,7 +146,7 @@ struct stencil_4points } }; -struct stencil_8points_standard +struct stencil_9points_standard { double delta_t = 0.125; @@ -250,7 +250,7 @@ template void classic_anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, double kappa) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, diffusion::perona_malik_diffusivity{kappa}); } @@ -265,13 +265,13 @@ void matlab_anisotropic_diffusion(const InputView& input, const OutputView& outp { if (conduction_method == matlab_conduction_method::exponential) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); } else if (conduction_method == matlab_conduction_method::quadratic) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); } @@ -284,13 +284,13 @@ void matlab_anisotropic_diffusion(const InputView& input, const OutputView& outp { if (conduction_method == matlab_conduction_method::exponential) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); } else if (conduction_method == matlab_conduction_method::quadratic) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_4points{}, + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); } @@ -309,7 +309,7 @@ template void default_anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, double kappa) { - anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_8points_standard{}, + anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_9points_standard{}, brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); } @@ -323,7 +323,7 @@ void default_anisotropic_diffusion(const InputView& input, const OutputView& out /// iteration count is set and grayscale image view is used /// as an input template void anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 1084e567a5..9b61c5b45c 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -179,8 +179,8 @@ void laplace_functions_test() auto rgb = gil::view(rgb_image); auto gray = gil::view(gray_image); - auto _4way = gil::laplace_function::stencil_4points{}; - auto _8way = gil::laplace_function::stencil_8points_standard{}; + auto _4way = gil::laplace_function::stencil_5points{}; + auto _8way = gil::laplace_function::stencil_9points_standard{}; for (std::ptrdiff_t y = 1; y < image_size.y - 1; ++y) { for (std::ptrdiff_t x = 1; x < image_size.x - 1; ++x) From 6321b527948d03385fa1860778b6dd6bbd3d8dfe Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 11 Jul 2020 02:10:53 +0600 Subject: [PATCH 39/48] Add docs for laplace stencil indexing --- include/boost/gil/image_processing/diffusion.hpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index c5c2e7d2c3..84107a250b 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -94,6 +94,11 @@ struct more_wide_regions_diffusivity } // namespace diffusion namespace laplace_function { +// The functions assume clockwise enumeration of stencil points, as such +// NW North NE 0 1 2 (-1, -1) (0, -1) (+1, -1) +// West East ===> 7 3 ===> (-1, 0) (+1, 0) +// SW South SE 6 5 4 (-1, +1) (0, +1) (+1, +1) + template using stencil_type = std::array; From 8106bee9bc24c196150dcbc17977f1ce17d77773 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 11 Jul 2020 02:17:40 +0600 Subject: [PATCH 40/48] Use indexing and add offsets function The new function will provide a uniform way to generate stencils by making sure directions are indexed properly --- .../boost/gil/image_processing/diffusion.hpp | 30 ++++++++++++------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 84107a250b..492386d303 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -99,6 +99,12 @@ namespace laplace_function { // West East ===> 7 3 ===> (-1, 0) (+1, 0) // SW South SE 6 5 4 (-1, +1) (0, +1) (+1, +1) +std::array get_directed_offsets() +{ + return {point_t{-1, -1}, point_t{0, -1}, point_t{+1, -1}, point_t{+1, 0}, + point_t{+1, +1}, point_t{0, +1}, point_t{-1, +1}, point_t{-1, 0}}; +} + template using stencil_type = std::array; @@ -110,19 +116,23 @@ struct stencil_5points stencil_type compute_laplace(SubImageView view, point_t origin) { - stencil_type stencil; - auto out = stencil.begin(); auto current = view(origin); + stencil_type stencil; using channel_type = typename channel_type::type; - std::array offsets{point_t{0, -1}, point_t{+1, 0}, point_t{0, +1}, - point_t{-1, 0}}; + std::array offsets(get_directed_offsets()); typename SubImageView::value_type zero_pixel; static_fill(zero_pixel, 0); - for (auto offset : offsets) + for (std::size_t index = 0; index < offsets.size(); ++index) { - *out++ = zero_pixel; - static_transform(view(origin.x + offset.x, origin.y + offset.y), current, *out++, - std::minus{}); + if (index % 2 != 0) + { + static_transform(view(origin.x + offsets[index].x, origin.y + offsets[index].y), + current, stencil[index], std::minus{}); + } + else + { + stencil[index] = zero_pixel; + } } return stencil; } @@ -163,9 +173,7 @@ struct stencil_9points_standard auto out = stencil.begin(); auto current = view(origin); using channel_type = typename channel_type::type; - std::array offsets{point_t{-1, -1}, point_t{0, -1}, point_t{+1, -1}, - point_t{+1, 0}, point_t{+1, +1}, point_t{0, +1}, - point_t{-1, +1}, point_t{-1, 0}}; + std::array offsets(get_directed_offsets()); for (auto offset : offsets) { static_transform(view(origin.x + offset.x, origin.y + offset.y), current, *out++, From e82f2c2f09cb205b09d9febad5830934062d6bc9 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 11 Jul 2020 02:22:50 +0600 Subject: [PATCH 41/48] Add only required stencil points The 5 points Laplace stencil is now adding only required points and not assuming that others are zero --- include/boost/gil/image_processing/diffusion.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 492386d303..e12f6a1998 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -149,9 +149,9 @@ struct stencil_5points return zero_pixel; }(); - for (; first != last; ++first) + for (std::size_t index : {1u, 3u, 5u, 7u}) { - static_transform(result, *first, result, std::plus{}); + static_transform(result, stencil[index], result, std::plus{}); } Pixel delta_t_pixel; static_fill(delta_t_pixel, delta_t); From f8e382dc6f7803255288c2b402b7ccb0572903cb Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Sat, 11 Jul 2020 02:34:29 +0600 Subject: [PATCH 42/48] Add file prefix for Jamfile in tests --- test/core/image_processing/Jamfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/core/image_processing/Jamfile b/test/core/image_processing/Jamfile index f2f547f560..be4c755c02 100644 --- a/test/core/image_processing/Jamfile +++ b/test/core/image_processing/Jamfile @@ -19,4 +19,4 @@ run hessian.cpp ; run sobel_scharr.cpp ; run box_filter.cpp ; run median_filter.cpp ; -run anisotropic_diffusion ; +run anisotropic_diffusion.cpp ; From aaebbdb1386f367b91ec3b1c6826b6757fd8951c Mon Sep 17 00:00:00 2001 From: Olzhas Date: Thu, 6 Aug 2020 20:21:11 +0600 Subject: [PATCH 43/48] Rename diffusivity to conductivity --- .../boost/gil/image_processing/diffusion.hpp | 24 +++++++++---------- .../anisotropic_diffusion.cpp | 16 ++++++------- 2 files changed, 20 insertions(+), 20 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index e12f6a1998..99e28320d5 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -19,8 +19,8 @@ #include namespace boost { namespace gil { -namespace diffusion { -struct perona_malik_diffusivity +namespace conductivity { +struct perona_malik_conductivity { double kappa; template @@ -38,7 +38,7 @@ struct perona_malik_diffusivity } }; -struct gaussian_diffusivity +struct gaussian_conductivity { double kappa; template @@ -56,7 +56,7 @@ struct gaussian_diffusivity } }; -struct wide_regions_diffusivity +struct wide_regions_conductivity { double kappa; template @@ -74,7 +74,7 @@ struct wide_regions_diffusivity } }; -struct more_wide_regions_diffusivity +struct more_wide_regions_conductivity { double kappa; template @@ -265,7 +265,7 @@ void classic_anisotropic_diffusion(const InputView& input, const OutputView& out { anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, - diffusion::perona_malik_diffusivity{kappa}); + conductivity::perona_malik_conductivity{kappa}); } template @@ -280,13 +280,13 @@ void matlab_anisotropic_diffusion(const InputView& input, const OutputView& outp { anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, - diffusion::gaussian_diffusivity{kappa}); + conductivity::gaussian_conductivity{kappa}); } else if (conduction_method == matlab_conduction_method::quadratic) { anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, - diffusion::gaussian_diffusivity{kappa}); + conductivity::gaussian_conductivity{kappa}); } else { @@ -299,13 +299,13 @@ void matlab_anisotropic_diffusion(const InputView& input, const OutputView& outp { anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, - diffusion::gaussian_diffusivity{kappa}); + conductivity::gaussian_conductivity{kappa}); } else if (conduction_method == matlab_conduction_method::quadratic) { anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_5points{}, brightness_function::identity{}, - diffusion::gaussian_diffusivity{kappa}); + conductivity::gaussian_conductivity{kappa}); } else { @@ -323,7 +323,7 @@ void default_anisotropic_diffusion(const InputView& input, const OutputView& out unsigned int num_iter, double kappa) { anisotropic_diffusion(input, output, num_iter, laplace_function::stencil_9points_standard{}, - brightness_function::identity{}, diffusion::gaussian_diffusivity{kappa}); + brightness_function::identity{}, conductivity::gaussian_conductivity{kappa}); } /// \brief Performs diffusion according to Perona-Malik equation @@ -338,7 +338,7 @@ void default_anisotropic_diffusion(const InputView& input, const OutputView& out template + typename DiffusivityFunction = conductivity::gaussian_conductivity> void anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter, LaplaceStrategy laplace, BrightnessFunction brightness, DiffusivityFunction diffusivity) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index 9b61c5b45c..b329aba28d 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -277,20 +277,20 @@ int main() for (double kappa = 5; kappa <= 70; ++kappa) { diffusion_function_check( - gil::diffusion::perona_malik_diffusivity{kappa}); - diffusion_function_check(gil::diffusion::gaussian_diffusivity{kappa}); + gil::conductivity::perona_malik_conductivity{kappa}); + diffusion_function_check(gil::conductivity::gaussian_conductivity{kappa}); diffusion_function_check( - gil::diffusion::wide_regions_diffusivity{kappa}); + gil::conductivity::wide_regions_conductivity{kappa}); diffusion_function_check( - gil::diffusion::more_wide_regions_diffusivity{kappa}); + gil::conductivity::more_wide_regions_conductivity{kappa}); diffusion_function_check( - gil::diffusion::perona_malik_diffusivity{kappa}); - diffusion_function_check(gil::diffusion::gaussian_diffusivity{kappa}); + gil::conductivity::perona_malik_conductivity{kappa}); + diffusion_function_check(gil::conductivity::gaussian_conductivity{kappa}); diffusion_function_check( - gil::diffusion::wide_regions_diffusivity{kappa}); + gil::conductivity::wide_regions_conductivity{kappa}); diffusion_function_check( - gil::diffusion::more_wide_regions_diffusivity{kappa}); + gil::conductivity::more_wide_regions_conductivity{kappa}); } brightness_function_test(); From b5434da2efe48842860118ec428053360668ba05 Mon Sep 17 00:00:00 2001 From: Olzhas Date: Thu, 6 Aug 2020 20:26:49 +0600 Subject: [PATCH 44/48] Rename _4way=>s4way and _8way=>s8way --- .../anisotropic_diffusion.cpp | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index b329aba28d..ea4a31801b 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -179,17 +179,17 @@ void laplace_functions_test() auto rgb = gil::view(rgb_image); auto gray = gil::view(gray_image); - auto _4way = gil::laplace_function::stencil_5points{}; - auto _8way = gil::laplace_function::stencil_9points_standard{}; + auto s4way = gil::laplace_function::stencil_5points{}; + auto s8way = gil::laplace_function::stencil_9points_standard{}; for (std::ptrdiff_t y = 1; y < image_size.y - 1; ++y) { for (std::ptrdiff_t x = 1; x < image_size.x - 1; ++x) { - auto rgb_4way_stencil = _4way.compute_laplace(rgb, {x, y}); - auto gray_4way_stencil = _4way.compute_laplace(gray, {x, y}); + auto rgb_4way_stencil = s4way.compute_laplace(rgb, {x, y}); + auto gray_4way_stencil = s4way.compute_laplace(gray, {x, y}); - auto rgb_8way_stencil = _8way.compute_laplace(rgb, {x, y}); - auto gray_8way_stencil = _8way.compute_laplace(gray, {x, y}); + auto rgb_8way_stencil = s8way.compute_laplace(rgb, {x, y}); + auto gray_8way_stencil = s8way.compute_laplace(gray, {x, y}); test_stencil_pixels(rgb_4way_stencil, zero_rgb_pixel); test_stencil_pixels(rgb_8way_stencil, zero_rgb_pixel); @@ -232,9 +232,9 @@ void laplace_functions_test() auto test_point = gil::point_t(2, 2); std::array _4way_expected{0, -14, 0, 33, 0, 6, 0, -8}; - auto rgb_4way = _4way.compute_laplace(rgb, test_point); + auto rgb_4way = s4way.compute_laplace(rgb, test_point); test_stencil(rgb_4way, _4way_expected); - auto gray_4way = _4way.compute_laplace(gray, test_point); + auto gray_4way = s4way.compute_laplace(gray, test_point); test_stencil(gray_4way, _4way_expected); // 8 way stencil should be @@ -243,25 +243,25 @@ void laplace_functions_test() // 65 6 15 std::array _8way_expected{-20, -14, -23, 33, 15, 6, 65, -8}; - auto rgb_8way = _8way.compute_laplace(rgb, test_point); + auto rgb_8way = s8way.compute_laplace(rgb, test_point); test_stencil(rgb_8way, _8way_expected); - auto gray_8way = _8way.compute_laplace(gray, test_point); + auto gray_8way = s8way.compute_laplace(gray, test_point); test_stencil(gray_8way, _8way_expected); // reduce result for 4 way should be (-14 - 8 + 6 + 33) * 0.25 = 4.25 - auto rgb_reduced_4way = _4way.reduce(rgb_4way); + auto rgb_reduced_4way = s4way.reduce(rgb_4way); gil::static_for_each(rgb_reduced_4way, [](gil::float32_t value) { BOOST_TEST(value == 4.25f); }); - auto gray_reduced_4way = _4way.reduce(gray_4way); + auto gray_reduced_4way = s4way.reduce(gray_4way); gil::static_for_each(gray_reduced_4way, [](gil::float32_t value) { BOOST_TEST(value == 4.25f); }); // reduce result for 8 way should be ((-20-23+15+65)*0.5+(-14+33+6-8)) * 0.125 = (18.5+17) * // 0.125 = 4.4375 - auto rgb_reduced_8way = _8way.reduce(rgb_8way); + auto rgb_reduced_8way = s8way.reduce(rgb_8way); gil::static_for_each(rgb_reduced_8way, [](gil::float32_t value) { BOOST_TEST(value == 4.4375); }); - auto gray_reduced_8way = _8way.reduce(gray_8way); + auto gray_reduced_8way = s8way.reduce(gray_8way); gil::static_for_each(gray_reduced_8way, [](gil::float32_t value) { BOOST_TEST(value == 4.4375); }); } From c77ef2d1fc3ad274d443988592946c85aac1f169 Mon Sep 17 00:00:00 2001 From: Olzhas Date: Thu, 6 Aug 2020 20:30:51 +0600 Subject: [PATCH 45/48] Capture this by value instead of local Since this cannot be captured by ref, it is captured by value instead. The change is to not declare redundant variables --- .../boost/gil/image_processing/diffusion.hpp | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 99e28320d5..0a3f6369fe 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -28,9 +28,8 @@ struct perona_malik_conductivity { using channel_type = typename channel_type::type; // C++11 doesn't seem to capture members - auto local_kappa = kappa; - static_transform(input, input, [local_kappa](channel_type value) { - value /= local_kappa; + static_transform(input, input, [this](channel_type value) { + value /= kappa; return std::exp(-std::abs(value)); }); @@ -46,9 +45,8 @@ struct gaussian_conductivity { using channel_type = typename channel_type::type; // C++11 doesn't seem to capture members - auto local_kappa = kappa; - static_transform(input, input, [local_kappa](channel_type value) { - value /= local_kappa; + static_transform(input, input, [this](channel_type value) { + value /= kappa; return std::exp(-value * value); }); @@ -64,9 +62,8 @@ struct wide_regions_conductivity { using channel_type = typename channel_type::type; // C++11 doesn't seem to capture members - auto local_kappa = kappa; - static_transform(input, input, [local_kappa](channel_type value) { - value /= local_kappa; + static_transform(input, input, [this](channel_type value) { + value /= kappa; return 1.0 / (1.0 + value * value); }); @@ -82,9 +79,8 @@ struct more_wide_regions_conductivity { using channel_type = typename channel_type::type; // C++11 doesn't seem to capture members - auto local_kappa = kappa; - static_transform(input, input, [local_kappa](channel_type value) { - value /= local_kappa; + static_transform(input, input, [this](channel_type value) { + value /= kappa; return 1.0 / std::sqrt((1.0 + value * value)); }); From b0f81f03225779b931f0c33eb90f79a67714c349 Mon Sep 17 00:00:00 2001 From: Olzhas Date: Thu, 6 Aug 2020 21:27:54 +0600 Subject: [PATCH 46/48] Add docs for Laplacian stencils --- .../boost/gil/image_processing/diffusion.hpp | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp index 0a3f6369fe..72f55a2ae4 100644 --- a/include/boost/gil/image_processing/diffusion.hpp +++ b/include/boost/gil/image_processing/diffusion.hpp @@ -89,12 +89,23 @@ struct more_wide_regions_conductivity }; } // namespace diffusion +/** + \brief contains discrete approximations of 2D Laplacian operator +*/ namespace laplace_function { // The functions assume clockwise enumeration of stencil points, as such // NW North NE 0 1 2 (-1, -1) (0, -1) (+1, -1) // West East ===> 7 3 ===> (-1, 0) (+1, 0) // SW South SE 6 5 4 (-1, +1) (0, +1) (+1, +1) +/** + \brief This function makes sure all Laplace functions enumerate + values in the same order and direction. + + The first element is difference North West direction, second in North, + and so on in clockwise manner. Leave element as zero if it is not + to be computed. +*/ std::array get_directed_offsets() { return {point_t{-1, -1}, point_t{0, -1}, point_t{+1, -1}, point_t{+1, 0}, @@ -104,6 +115,11 @@ std::array get_directed_offsets() template using stencil_type = std::array; +/** + \brief 5 point stencil approximation of Laplacian + + Only main 4 directions are non-zero, the rest are zero +*/ struct stencil_5points { double delta_t = 0.25; @@ -157,6 +173,12 @@ struct stencil_5points } }; +/** + \brief 9 point stencil approximation of Laplacian + + This is full 8 way approximation, though diagonal + elements are halved during reduction. +*/ struct stencil_9points_standard { double delta_t = 0.125; From 6c9b272fbfcdcaceedefa19455b3af17ca106bb1 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Wed, 26 Aug 2020 01:24:03 +0600 Subject: [PATCH 47/48] Use random_value from test_fixture.hpp Also decrease seed range to make test duration manageable. --- .../image_processing/anisotropic_diffusion.cpp | 15 +++++++++------ test/core/test_fixture.hpp | 17 +++++++++-------- 2 files changed, 18 insertions(+), 14 deletions(-) diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp index ea4a31801b..58a5c9d6ce 100644 --- a/test/core/image_processing/anisotropic_diffusion.cpp +++ b/test/core/image_processing/anisotropic_diffusion.cpp @@ -5,6 +5,7 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // +#include "../test_fixture.hpp" #include "boost/gil/algorithm.hpp" #include #include @@ -91,8 +92,8 @@ void brightness_function_test() template void heat_conservation_test(std::uint32_t seed) { - std::mt19937 twister(seed); - std::uniform_int_distribution dist(0, std::numeric_limits::max()); + gil::test::fixture::random_value dist(seed, 0, + std::numeric_limits::max()); ImageType image(32, 32); auto view = gil::view(image); @@ -102,7 +103,7 @@ void heat_conservation_test(std::uint32_t seed) { for (std::size_t channel_index = 0; channel_index < num_channels; ++channel_index) { - pixel[channel_index] = static_cast(dist(twister)); + pixel[channel_index] = static_cast(dist()); before_diffusion[channel_index] += pixel[channel_index]; } } @@ -268,7 +269,7 @@ void laplace_functions_test() int main() { - for (std::uint32_t seed = 0; seed < 1000; ++seed) + for (std::uint32_t seed = 0; seed < 100; ++seed) { heat_conservation_test(seed); heat_conservation_test(seed); @@ -278,7 +279,8 @@ int main() { diffusion_function_check( gil::conductivity::perona_malik_conductivity{kappa}); - diffusion_function_check(gil::conductivity::gaussian_conductivity{kappa}); + diffusion_function_check( + gil::conductivity::gaussian_conductivity{kappa}); diffusion_function_check( gil::conductivity::wide_regions_conductivity{kappa}); diffusion_function_check( @@ -286,7 +288,8 @@ int main() diffusion_function_check( gil::conductivity::perona_malik_conductivity{kappa}); - diffusion_function_check(gil::conductivity::gaussian_conductivity{kappa}); + diffusion_function_check( + gil::conductivity::gaussian_conductivity{kappa}); diffusion_function_check( gil::conductivity::wide_regions_conductivity{kappa}); diffusion_function_check( diff --git a/test/core/test_fixture.hpp b/test/core/test_fixture.hpp index 819b2cac22..7d0606dc12 100644 --- a/test/core/test_fixture.hpp +++ b/test/core/test_fixture.hpp @@ -8,8 +8,8 @@ #ifndef BOOST_GIL_TEST_CORE_TEST_FIXTURE_HPP #define BOOST_GIL_TEST_CORE_TEST_FIXTURE_HPP -#include #include +#include #include #include @@ -60,19 +60,20 @@ template struct random_value { static_assert(std::is_integral::value, "T must be integral type"); - static constexpr auto range_min = std::numeric_limits::min(); - static constexpr auto range_max = std::numeric_limits::max(); - random_value() : rng_(rd_()), uid_(range_min, range_max) {} + random_value(T range_min = std::numeric_limits::min(), + T range_max = std::numeric_limits::max()) + : uid_(range_min, range_max) + {} + + random_value(std::uint32_t seed, T minimum, T maximum) : rng_(seed), uid_(minimum, maximum) + {} T operator()() { - auto value = uid_(rng_); - BOOST_ASSERT(range_min <= value && value <= range_max); - return static_cast(value); + return uid_(rng_); } - std::random_device rd_; std::mt19937 rng_; std::uniform_int_distribution::type> uid_; }; From fbad958e8342145b64f64d41c2f0c79494347fc2 Mon Sep 17 00:00:00 2001 From: Olzhas Zhumabek Date: Wed, 26 Aug 2020 01:26:39 +0600 Subject: [PATCH 48/48] Add range_min and max to random_value --- test/core/test_fixture.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/core/test_fixture.hpp b/test/core/test_fixture.hpp index 7d0606dc12..d7a35abe9a 100644 --- a/test/core/test_fixture.hpp +++ b/test/core/test_fixture.hpp @@ -74,6 +74,16 @@ struct random_value return uid_(rng_); } + T range_min() const noexcept + { + return uid_.a(); + } + + T range_max() const noexcept + { + return uid_.b(); + } + std::mt19937 rng_; std::uniform_int_distribution::type> uid_; };