Add Perona–Malik anisotropic diffusion algorithm (#500)

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

The accumulate part was wrong, it multiplied by delta_t on every sum, which is wrong

Use 8 way nabla compute.
This is just different discretization of Laplace operator
https://en.wikipedia.org/wiki/Discrete_Laplace_operator

Laplace stencils are now the same as in mathematical notation

The new function will provide a uniform way to generate stencils
by making sure directions are indexed properly

Add only required stencil points:
The 5 points Laplace stencil is now adding only required points and not assuming that others are zero

example/Jamfile

@@ -21,6 +21,7 @@ project
 local sources =
     adaptive_threshold.cpp
     affine.cpp
+    anisotropic_diffusion.cpp
     convolution.cpp
     convolve2d.cpp
     dynamic_image.cpp

example/anisotropic_diffusion.cpp

@@ -0,0 +1,128 @@
+//
+// Copyright 2020 Olzhas Zhumabek <anonymous.from.applecity@gmail.com>
+//
+// 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 <boost/gil/algorithm.hpp>
+#include <boost/gil/extension/io/png.hpp>
+#include <boost/gil/image.hpp>
+#include <boost/gil/image_processing/diffusion.hpp>
+#include <boost/gil/image_view.hpp>
+#include <boost/gil/io/typedefs.hpp>
+#include <boost/gil/typedefs.hpp>
+
+#include <cmath>
+#include <iostream>
+#include <string>
+#include <type_traits>
+
+namespace gil = boost::gil;
+
+void gray_version(std::string const& input_path, std::string const& output_path,
+                  unsigned int iteration_count, float kappa)
+{
+    gil::gray8_image_t input;
+    gil::read_image(input_path, input, gil::png_tag{});
+    auto input_view = gil::view(input);
+
+    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]; });
+    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, 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]; });
+
+    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<gil::uint8_t>(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, float 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]);
+    });
+    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, 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];
+        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<gil::uint8_t>(p[0]), static_cast<gil::uint8_t>(p[1]),
+                                 static_cast<gil::uint8_t>(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[])
+{
+    if (argc != 6)
+    {
+        std::cerr << "usage: " << argv[0]
+                  << " <input.png> <output.png>"
+                     " <colorspace: gray|rgb> <positive iteration count> <positive kappa~30>\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<unsigned int>(std::stoul(argv[4]));
+    float kappa = std::stof(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";
+    }
+}