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/example/anisotropic_diffusion.cpp b/example/anisotropic_diffusion.cpp
new file mode 100644
index 0000000000..03defb7566
--- /dev/null
+++ b/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";
+    }
+}
diff --git a/include/boost/gil/image_processing/diffusion.hpp b/include/boost/gil/image_processing/diffusion.hpp
new file mode 100644
index 0000000000..72f55a2ae4
--- /dev/null
+++ b/include/boost/gil/image_processing/diffusion.hpp
@@ -0,0 +1,421 @@
+//
+// 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/detail/math.hpp"
+#include <boost/gil/algorithm.hpp>
+#include <boost/gil/color_base_algorithm.hpp>
+#include <boost/gil/image.hpp>
+#include <boost/gil/image_view.hpp>
+#include <boost/gil/image_view_factory.hpp>
+#include <boost/gil/pixel.hpp>
+#include <boost/gil/point.hpp>
+#include <boost/gil/typedefs.hpp>
+#include <functional>
+#include <numeric>
+#include <vector>
+
+namespace boost { namespace gil {
+namespace conductivity {
+struct perona_malik_conductivity
+{
+    double kappa;
+    template <typename Pixel>
+    Pixel operator()(Pixel input)
+    {
+        using channel_type = typename channel_type<Pixel>::type;
+        // C++11 doesn't seem to capture members
+        static_transform(input, input, [this](channel_type value) {
+            value /= kappa;
+            return std::exp(-std::abs(value));
+        });
+
+        return input;
+    }
+};
+
+struct gaussian_conductivity
+{
+    double kappa;
+    template <typename Pixel>
+    Pixel operator()(Pixel input)
+    {
+        using channel_type = typename channel_type<Pixel>::type;
+        // C++11 doesn't seem to capture members
+        static_transform(input, input, [this](channel_type value) {
+            value /= kappa;
+            return std::exp(-value * value);
+        });
+
+        return input;
+    }
+};
+
+struct wide_regions_conductivity
+{
+    double kappa;
+    template <typename Pixel>
+    Pixel operator()(Pixel input)
+    {
+        using channel_type = typename channel_type<Pixel>::type;
+        // C++11 doesn't seem to capture members
+        static_transform(input, input, [this](channel_type value) {
+            value /= kappa;
+            return 1.0 / (1.0 + value * value);
+        });
+
+        return input;
+    }
+};
+
+struct more_wide_regions_conductivity
+{
+    double kappa;
+    template <typename Pixel>
+    Pixel operator()(Pixel input)
+    {
+        using channel_type = typename channel_type<Pixel>::type;
+        // C++11 doesn't seem to capture members
+        static_transform(input, input, [this](channel_type value) {
+            value /= kappa;
+            return 1.0 / std::sqrt((1.0 + value * value));
+        });
+
+        return input;
+    }
+};
+} // 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<gil::point_t, 8> 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 <typename PixelType>
+using stencil_type = std::array<PixelType, 8>;
+
+/**
+    \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;
+
+    template <typename SubImageView>
+    stencil_type<typename SubImageView::value_type> compute_laplace(SubImageView view,
+                                                                    point_t origin)
+    {
+        auto current = view(origin);
+        stencil_type<typename SubImageView::value_type> stencil;
+        using channel_type = typename channel_type<typename SubImageView::value_type>::type;
+        std::array<gil::point_t, 8> offsets(get_directed_offsets());
+        typename SubImageView::value_type zero_pixel;
+        static_fill(zero_pixel, 0);
+        for (std::size_t index = 0; index < offsets.size(); ++index)
+        {
+            if (index % 2 != 0)
+            {
+                static_transform(view(origin.x + offsets[index].x, origin.y + offsets[index].y),
+                                 current, stencil[index], std::minus<channel_type>{});
+            }
+            else
+            {
+                stencil[index] = zero_pixel;
+            }
+        }
+        return stencil;
+    }
+
+    template <typename Pixel>
+    Pixel reduce(const stencil_type<Pixel>& stencil)
+    {
+        auto first = stencil.begin();
+        auto last = stencil.end();
+        using channel_type = typename channel_type<Pixel>::type;
+        auto result = []() {
+            Pixel zero_pixel;
+            static_fill(zero_pixel, channel_type(0));
+            return zero_pixel;
+        }();
+
+        for (std::size_t index : {1u, 3u, 5u, 7u})
+        {
+            static_transform(result, stencil[index], result, std::plus<channel_type>{});
+        }
+        Pixel delta_t_pixel;
+        static_fill(delta_t_pixel, delta_t);
+        static_transform(result, delta_t_pixel, result, std::multiplies<channel_type>{});
+
+        return result;
+    }
+};
+
+/**
+    \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;
+
+    template <typename SubImageView>
+    stencil_type<typename SubImageView::value_type> compute_laplace(SubImageView view,
+                                                                    point_t origin)
+    {
+        stencil_type<typename SubImageView::value_type> stencil;
+        auto out = stencil.begin();
+        auto current = view(origin);
+        using channel_type = typename channel_type<typename SubImageView::value_type>::type;
+        std::array<gil::point_t, 8> offsets(get_directed_offsets());
+        for (auto offset : offsets)
+        {
+            static_transform(view(origin.x + offset.x, origin.y + offset.y), current, *out++,
+                             std::minus<channel_type>{});
+        }
+
+        return stencil;
+    }
+
+    template <typename Pixel>
+    Pixel reduce(const stencil_type<Pixel>& stencil)
+    {
+        using channel_type = typename channel_type<Pixel>::type;
+        auto result = []() {
+            Pixel zero_pixel;
+            static_fill(zero_pixel, channel_type(0));
+            return zero_pixel;
+        }();
+        for (std::size_t index : {1u, 3u, 5u, 7u})
+        {
+            static_transform(result, stencil[index], result, std::plus<channel_type>{});
+        }
+
+        for (std::size_t index : {0u, 2u, 4u, 6u})
+        {
+            Pixel half_pixel;
+            static_fill(half_pixel, channel_type(1 / 2.0));
+            static_transform(stencil[index], half_pixel, half_pixel,
+                             std::multiplies<channel_type>{});
+            static_transform(result, half_pixel, result, std::plus<channel_type>{});
+        }
+
+        Pixel delta_t_pixel;
+        static_fill(delta_t_pixel, delta_t);
+        static_transform(result, delta_t_pixel, result, std::multiplies<channel_type>{});
+
+        return result;
+    }
+};
+} // namespace laplace_function
+
+namespace brightness_function {
+using laplace_function::stencil_type;
+struct identity
+{
+    template <typename Pixel>
+    stencil_type<Pixel> operator()(const stencil_type<Pixel>& 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
+
+struct rgb_luminance
+{
+    using pixel_type = rgb32f_pixel_t;
+    stencil_type<pixel_type> operator()(const stencil_type<pixel_type>& stencil)
+    {
+        stencil_type<pixel_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 result_pixel;
+        });
+        return output;
+    }
+};
+
+} // namespace brightness_function
+
+enum class matlab_connectivity
+{
+    minimal,
+    maximal
+};
+
+enum class matlab_conduction_method
+{
+    exponential,
+    quadratic
+};
+
+template <typename InputView, typename OutputView>
+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_5points{},
+                          brightness_function::identity{},
+                          conductivity::perona_malik_conductivity{kappa});
+}
+
+template <typename InputView, typename OutputView>
+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_5points{},
+                                  brightness_function::identity{},
+                                  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{},
+                                  conductivity::gaussian_conductivity{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_5points{},
+                                  brightness_function::identity{},
+                                  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{},
+                                  conductivity::gaussian_conductivity{kappa});
+        }
+        else
+        {
+            throw std::logic_error("unhandled conduction method found");
+        }
+    }
+    else
+    {
+        throw std::logic_error("unhandled connectivity found");
+    }
+}
+
+template <typename InputView, typename OutputView>
+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_9points_standard{},
+                          brightness_function::identity{}, conductivity::gaussian_conductivity{kappa});
+}
+
+/// \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 <typename InputView, typename OutputView,
+          typename LaplaceStrategy = laplace_function::stencil_9points_standard,
+          typename BrightnessFunction = brightness_function::identity,
+          typename DiffusivityFunction = conductivity::gaussian_conductivity>
+void anisotropic_diffusion(const InputView& input, const OutputView& output, unsigned int num_iter,
+                           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<pixel_type>::type;
+    using computation_image = image<pixel_type>;
+    const auto width = input.width();
+    const auto height = input.height();
+    const point_t dims(width, height);
+    const auto zero_pixel = []() {
+        pixel_type pixel;
+        static_fill(pixel, static_cast<channel_type>(0));
+
+        return pixel;
+    }();
+    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);
+    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<pixel_type>{}; ++i)
+                         {
+                             converted[i] = pixel[i];
+                         }
+                         return converted;
+                     });
+
+    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)
+            {
+                auto x = relative_x + 1;
+                auto y = relative_y + 1;
+                auto stencil = laplace.compute_laplace(result, point_t(x, y));
+                auto brightness_stencil = brightness(stencil);
+                laplace_function::stencil_type<pixel_type> diffusivity_stencil;
+                std::transform(brightness_stencil.begin(), brightness_stencil.end(),
+                               diffusivity_stencil.begin(), diffusivity);
+                laplace_function::stencil_type<pixel_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<channel_type>{});
+                                   return lhs;
+                               });
+                static_transform(result(x, y), laplace.reduce(product_stencil),
+                                 scratch_result(x, y), std::plus<channel_type>{});
+            }
+        }
+        using std::swap;
+        swap(result, scratch_result);
+    }
+
+    copy_pixels(subimage_view(result, 1, 1, width, height), output);
+}
+
+}} // namespace boost::gil
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/Jamfile b/test/core/image_processing/Jamfile
index d9593f0793..be4c755c02 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.cpp ;
diff --git a/test/core/image_processing/anisotropic_diffusion.cpp b/test/core/image_processing/anisotropic_diffusion.cpp
new file mode 100644
index 0000000000..58a5c9d6ce
--- /dev/null
+++ b/test/core/image_processing/anisotropic_diffusion.cpp
@@ -0,0 +1,304 @@
+//
+// 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 "../test_fixture.hpp"
+#include "boost/gil/algorithm.hpp"
+#include <boost/gil/color_base_algorithm.hpp>
+#include <boost/gil/image.hpp>
+#include <boost/gil/image_processing/diffusion.hpp>
+#include <boost/gil/image_view.hpp>
+#include <boost/gil/pixel.hpp>
+#include <boost/gil/typedefs.hpp>
+
+#include <boost/core/lightweight_test.hpp>
+
+#include <cmath>
+#include <limits>
+#include <random>
+
+namespace gil = boost::gil;
+
+template <typename PixelType, typename DiffusionFunction>
+void diffusion_function_check(DiffusionFunction diffusion)
+{
+    using limits = std::numeric_limits<gil::uint8_t>;
+    using channel_type = typename gil::channel_type<PixelType>::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 channel_value) {
+            BOOST_TEST(0 <= channel_value && channel_value <= 1.0);
+        });
+    }
+}
+
+void brightness_function_test()
+{
+    std::vector<gil::rgb32f_pixel_t> 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<gil::rgb32f_pixel_t> 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);
+        }
+    }
+
+    std::vector<float> 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<gil::rgb32f_pixel_t> 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<gil::gray32f_pixel_t> 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<gil::gray32f_pixel_t> 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 <typename ImageType, typename OutputImageType>
+void heat_conservation_test(std::uint32_t seed)
+{
+    gil::test::fixture::random_value<std::uint32_t> dist(seed, 0,
+                                                         std::numeric_limits<gil::uint8_t>::max());
+
+    ImageType image(32, 32);
+    auto view = gil::view(image);
+    constexpr std::ptrdiff_t num_channels = gil::num_channels<decltype(view)>::value;
+    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] = static_cast<gil::uint8_t>(dist());
+            before_diffusion[channel_index] += pixel[channel_index];
+        }
+    }
+
+    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};
+    for (const auto& pixel : output_view)
+    {
+        for (std::size_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 =
+            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;
+#ifdef BOOST_GIL_TEST_DEBUG
+        std::cout << percentage << ' ';
+#endif
+        BOOST_TEST(percentage < 1);
+    }
+#ifdef BOOST_GIL_TEST_DEBUG
+    std::cout << '\n';
+#endif
+}
+
+template <typename PixelType>
+void test_stencil_pixels(const gil::laplace_function::stencil_type<PixelType>& stencil,
+                         const PixelType& reference)
+{
+    for (const auto& stencil_entry : stencil)
+    {
+        BOOST_TEST(stencil_entry == reference);
+    }
+}
+
+template <typename PixelType>
+void test_stencil(const gil::laplace_function::stencil_type<PixelType>& stencil,
+                  const std::array<float, 8>& 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);
+
+    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 = s4way.compute_laplace(rgb, {x, y});
+            auto gray_4way_stencil = s4way.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);
+            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(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(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(2, 1) = gray_pixel(11);
+    gray(3, 1) = gray_pixel(2);
+    gray(1, 2) = gray_pixel(17);
+    gray(2, 2) = gray_pixel(25);
+    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
+    // 0 -14 0
+    // -8  <not computed for center> 33
+    // 0 6 0
+
+    auto test_point = gil::point_t(2, 2);
+    std::array<float, 8> _4way_expected{0, -14, 0, 33, 0, 6, 0, -8};
+    auto rgb_4way = s4way.compute_laplace(rgb, test_point);
+    test_stencil(rgb_4way, _4way_expected);
+    auto gray_4way = s4way.compute_laplace(gray, test_point);
+    test_stencil(gray_4way, _4way_expected);
+
+    // 8 way stencil should be
+    // -20 -14 -23
+    // -8 <not computed for center> 33
+    // 65 6 15
+
+    std::array<float, 8> _8way_expected{-20, -14, -23, 33, 15, 6, 65, -8};
+    auto rgb_8way = s8way.compute_laplace(rgb, test_point);
+    test_stencil(rgb_8way, _8way_expected);
+    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 = s4way.reduce(rgb_4way);
+    gil::static_for_each(rgb_reduced_4way,
+                         [](gil::float32_t value) { BOOST_TEST(value == 4.25f); });
+    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 = s8way.reduce(rgb_8way);
+    gil::static_for_each(rgb_reduced_8way,
+                         [](gil::float32_t value) { BOOST_TEST(value == 4.4375); });
+    auto gray_reduced_8way = s8way.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 < 100; ++seed)
+    {
+        heat_conservation_test<gil::gray8_image_t, gil::gray32f_image_t>(seed);
+        heat_conservation_test<gil::rgb8_image_t, gil::rgb32f_image_t>(seed);
+    }
+
+    for (double kappa = 5; kappa <= 70; ++kappa)
+    {
+        diffusion_function_check<gil::rgb32f_pixel_t>(
+            gil::conductivity::perona_malik_conductivity{kappa});
+        diffusion_function_check<gil::rgb32f_pixel_t>(
+            gil::conductivity::gaussian_conductivity{kappa});
+        diffusion_function_check<gil::rgb32f_pixel_t>(
+            gil::conductivity::wide_regions_conductivity{kappa});
+        diffusion_function_check<gil::rgb32f_pixel_t>(
+            gil::conductivity::more_wide_regions_conductivity{kappa});
+
+        diffusion_function_check<gil::gray32f_pixel_t>(
+            gil::conductivity::perona_malik_conductivity{kappa});
+        diffusion_function_check<gil::gray32f_pixel_t>(
+            gil::conductivity::gaussian_conductivity{kappa});
+        diffusion_function_check<gil::gray32f_pixel_t>(
+            gil::conductivity::wide_regions_conductivity{kappa});
+        diffusion_function_check<gil::gray32f_pixel_t>(
+            gil::conductivity::more_wide_regions_conductivity{kappa});
+    }
+
+    brightness_function_test();
+
+    laplace_functions_test();
+
+    return boost::report_errors();
+}
diff --git a/test/core/test_fixture.hpp b/test/core/test_fixture.hpp
index 819b2cac22..d7a35abe9a 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 <boost/gil/promote_integral.hpp>
 #include <boost/assert.hpp>
+#include <boost/gil/promote_integral.hpp>
 
 #include <cstdint>
 #include <initializer_list>
@@ -60,19 +60,30 @@ template <typename T>
 struct random_value
 {
     static_assert(std::is_integral<T>::value, "T must be integral type");
-    static constexpr auto range_min = std::numeric_limits<T>::min();
-    static constexpr auto range_max = std::numeric_limits<T>::max();
 
-    random_value() : rng_(rd_()), uid_(range_min, range_max) {}
+    random_value(T range_min = std::numeric_limits<T>::min(),
+                 T range_max = std::numeric_limits<T>::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<T>(value);
+        return uid_(rng_);
+    }
+
+    T range_min() const noexcept
+    {
+        return uid_.a();
+    }
+
+    T range_max() const noexcept
+    {
+        return uid_.b();
     }
 
-    std::random_device rd_;
     std::mt19937 rng_;
     std::uniform_int_distribution<typename gil::promote_integral<T>::type> uid_;
 };