Add implementation of Hough transforms (#512)

Support construction from step_size, step_count, and a function for angles

Implement angle and radious version of Hough line transform and adds a demo
with static line that goes over secondary diagonal.

Implement incremental line raster
Implement naive line raster
Implement Bresenham line raster
Leave only Bresenham line rasterization

Naive and incremental algorithms were removed because they are supposed
to produce the same results anyway.
The reason for diverging results is inaccuracy of floating point numbers

Add circle rendering through trigonometric functions, using
arctan(1 / (radius + 1)) as minimal angle step.

Trigonometric circle rasterizer does not follow circle equation, but still
produces very round shapes.
A new testing methodology needs to be devised for this rasterizer.

The new version accepts start and points inclusively and tries to use
canonic representation during computations.

Slope decided to be is (diff_y + 1) / (diff_x + 1).

example/hough_transform_circle.cpp

@@ -0,0 +1,55 @@
+// Boost.GIL (Generic Image Library) - tests
+//
+// 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.hpp>
+#include <boost/gil/extension/io/png.hpp>
+
+#include <iostream>
+#include <limits>
+#include <vector>
+
+namespace gil = boost::gil;
+
+int main()
+{
+    const std::size_t size = 128;
+    gil::gray8_image_t input_image(size, size);
+    auto input = gil::view(input_image);
+
+    const std::ptrdiff_t circle_radius = 16;
+    const gil::point_t circle_center = {64, 64};
+    const auto rasterizer = gil::midpoint_circle_rasterizer{};
+    std::vector<gil::point_t> circle_points(rasterizer.point_count(circle_radius));
+    rasterizer(circle_radius, circle_center, circle_points.begin());
+    for (const auto& point : circle_points)
+    {
+        input(point) = std::numeric_limits<gil::uint8_t>::max();
+    }
+
+    const auto radius_parameter =
+        gil::hough_parameter<std::ptrdiff_t>::from_step_count(circle_radius, 3, 3);
+    const auto x_parameter =
+        gil::hough_parameter<std::ptrdiff_t>::from_step_count(circle_center.x, 3, 3);
+    const auto y_parameter =
+        gil::hough_parameter<std::ptrdiff_t>::from_step_count(circle_center.x, 3, 3);
+
+    std::vector<gil::gray16_image_t> parameter_space_images(
+        radius_parameter.step_count,
+        gil::gray16_image_t(x_parameter.step_count, y_parameter.step_count));
+    std::vector<gil::gray16_view_t> parameter_space_views(parameter_space_images.size());
+    std::transform(parameter_space_images.begin(), parameter_space_images.end(),
+                   parameter_space_views.begin(),
+                   [](gil::gray16_image_t& img)
+                   {
+                       return gil::view(img);
+                   });
+
+    gil::hough_circle_transform_brute(input, radius_parameter, x_parameter, y_parameter,
+                                      parameter_space_views.begin(), rasterizer);
+    std::cout << parameter_space_views[3](3, 3) << '\n';
+}

example/hough_transform_line.cpp

@@ -0,0 +1,71 @@
+// Boost.GIL (Generic Image Library) - tests
+//
+// 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.hpp>
+#include <boost/gil/extension/io/png.hpp>
+#include <cmath>
+#include <cstddef>
+#include <iostream>
+
+namespace gil = boost::gil;
+
+int main()
+{
+    std::ptrdiff_t size = 32;
+    gil::gray16_image_t input_image(size, size);
+    auto input_view = gil::view(input_image);
+
+    // fill secondary diagonal with ones
+    // do note that origin is located at upper left,
+    // not bottom left as in usual plots
+    for (std::ptrdiff_t i = 0; i < size; ++i)
+    {
+        input_view(i, size - i - 1) = 1;
+    }
+
+    // print vertically flipped for better understanding of origin location
+    for (std::ptrdiff_t y = size - 1; y >= 0; --y)
+    {
+        for (std::ptrdiff_t x = 0; x < size; ++x)
+        {
+            std::cout << input_view(x, y)[0] << ' ';
+        }
+        std::cout << '\n';
+    }
+
+    double minimum_theta_step = std::atan(1.0 / size);
+    // this is the expected theta
+    double _45_degrees = gil::detail::pi / 4;
+    double _5_degrees = gil::detail::pi / 36;
+    std::size_t step_count = 5;
+    auto theta_parameter =
+        gil::make_theta_parameter(_45_degrees, _5_degrees, input_view.dimensions());
+    auto expected_radius = static_cast<std::ptrdiff_t>(std::round(std::cos(_45_degrees) * size));
+    auto radius_parameter =
+        gil::hough_parameter<std::ptrdiff_t>::from_step_size(expected_radius, 7, 1);
+    gil::gray32_image_t accumulator_array_image(theta_parameter.step_count,
+                                                radius_parameter.step_count);
+    auto accumulator_array = gil::view(accumulator_array_image);
+    gil::hough_line_transform(input_view, accumulator_array, theta_parameter, radius_parameter);
+    std::cout << "expecting maximum at theta=" << _45_degrees << " and radius=" << expected_radius
+              << '\n';
+    for (std::size_t theta_index = 0; theta_index < theta_parameter.step_count; ++theta_index)
+    {
+        for (std::size_t radius_index = 0; radius_index < radius_parameter.step_count;
+             ++radius_index)
+        {
+            double current_theta =
+                theta_parameter.start_point + theta_index * theta_parameter.step_size;
+            std::ptrdiff_t current_radius =
+                radius_parameter.start_point + radius_parameter.step_size * radius_index;
+            std::cout << "theta: " << current_theta << " radius: " << current_radius
+                      << " accumulated value: " << accumulator_array(theta_index, radius_index)[0]
+                      << '\n';
+        }
+    }
+}

example/rasterizer_circle.cpp

@@ -0,0 +1,33 @@
+// Boost.GIL (Generic Image Library) - tests
+//
+// 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.hpp>
+#include <boost/gil/extension/io/png.hpp>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+namespace gil = boost::gil;
+
+int main()
+{
+    const std::ptrdiff_t size = 256;
+    gil::gray8_image_t buffer_image(size, size);
+    auto buffer = gil::view(buffer_image);
+
+    const std::ptrdiff_t radius = 64;
+    const auto rasterizer = gil::trigonometric_circle_rasterizer{};
+    std::vector<gil::point_t> circle_points(rasterizer.point_count(radius));
+    rasterizer(radius, {128, 128}, circle_points.begin());
+    for (const auto& point : circle_points)
+    {
+        buffer(point) = std::numeric_limits<gil::uint8_t>::max();
+    }
+
+    gil::write_view("circle.png", buffer, gil::png_tag{});
+}

example/rasterizer_line.cpp

@@ -0,0 +1,42 @@
+// Boost.GIL (Generic Image Library) - tests
+//
+// 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.hpp>
+#include <boost/gil/extension/io/png.hpp>
+
+#include <limits>
+#include <vector>
+
+namespace gil = boost::gil;
+
+const std::ptrdiff_t size = 256;
+
+void line_bresenham(std::ptrdiff_t width, std::ptrdiff_t height, const std::string& output_name)
+{
+    const auto rasterizer = gil::bresenham_line_rasterizer{};
+    std::vector<gil::point_t> line_points(rasterizer.point_count(width, height));
+
+    gil::gray8_image_t image(size, size);
+    auto view = gil::view(image);
+
+    rasterizer({0, 0}, {width - 1, height - 1}, line_points.begin());
+    for (const auto& point : line_points)
+    {
+        view(point) = std::numeric_limits<gil::uint8_t>::max();
+    }
+
+    gil::write_view(output_name, view, gil::png_tag{});
+}
+
+int main()
+{
+    line_bresenham(256, 256, "line-bresenham-256-256.png");
+    line_bresenham(256, 128, "line-bresenham-256-128.png");
+    line_bresenham(256, 1, "line-bresenham-256-1.png");
+    line_bresenham(1, 256, "line-bresenham-1-256.png");
+}

include/boost/gil.hpp

@@ -24,6 +24,10 @@
 #include <boost/gil/gray.hpp>
 #include <boost/gil/histogram.hpp>
 #include <boost/gil/image.hpp>
+#include <boost/gil/image_processing/hough_parameter.hpp>
+#include <boost/gil/image_processing/hough_transform.hpp>
+#include <boost/gil/image_processing/scaling.hpp>
+#include <boost/gil/image_processing/threshold.hpp>
 #include <boost/gil/image_view.hpp>
 #include <boost/gil/image_view_factory.hpp>
 #include <boost/gil/iterator_from_2d.hpp>
@@ -38,13 +42,13 @@
 #include <boost/gil/point.hpp>
 #include <boost/gil/position_iterator.hpp>
 #include <boost/gil/premultiply.hpp>
+#include <boost/gil/rasterization/circle.hpp>
+#include <boost/gil/rasterization/line.hpp>
 #include <boost/gil/rgb.hpp>
 #include <boost/gil/rgba.hpp>
 #include <boost/gil/step_iterator.hpp>
 #include <boost/gil/typedefs.hpp>
 #include <boost/gil/utilities.hpp>
 #include <boost/gil/virtual_locator.hpp>
-#include <boost/gil/image_processing/scaling.hpp>
-#include <boost/gil/image_processing/threshold.hpp>
 
 #endif

include/boost/gil/image_processing/hough_parameter.hpp

@@ -0,0 +1,112 @@
+// Boost.GIL (Generic Image Library) - tests
+//
+// 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)
+//
+#ifndef BOOST_GIL_IMAGE_PROCESSING_HOUGH_PARAMETER_HPP
+#define BOOST_GIL_IMAGE_PROCESSING_HOUGH_PARAMETER_HPP
+
+#include <boost/gil/point.hpp>
+#include <cmath>
+#include <cstddef>
+
+namespace boost
+{
+namespace gil
+{
+/// \ingroup HoughTransform
+/// \brief A type to encapsulate Hough transform parameter range
+///
+/// This type provides a way to express value range for a parameter
+/// as well as some factory functions to simplify initialization
+template <typename T>
+struct hough_parameter
+{
+    T start_point;
+    T step_size;
+    std::size_t step_count;
+
+    /// \ingroup HoughTransform
+    /// \brief Create Hough parameter from value neighborhood and step count
+    ///
+    /// This function will take start_point as middle point, and in both
+    /// directions will try to walk half_step_count times until distance of
+    /// neighborhood is reached
+    static hough_parameter<T> from_step_count(T start_point, T neighborhood,
+                                              std::size_t half_step_count)
+    {
+        T step_size = neighborhood / half_step_count;
+        std::size_t step_count = half_step_count * 2 + 1;
+        // explicitly fill out members, as aggregate init will error out with narrowing
+        hough_parameter<T> parameter;
+        parameter.start_point = start_point - neighborhood;
+        parameter.step_size = step_size;
+        parameter.step_count = step_count;
+        return parameter;
+    }
+
+    /// \ingroup HoughTransform
+    /// \brief Create Hough parameter from value neighborhood and step size
+    ///
+    /// This function will take start_point as middle point, and in both
+    /// directions will try to walk step_size at a time until distance of
+    /// neighborhood is reached
+    static hough_parameter<T> from_step_size(T start_point, T neighborhood, T step_size)
+    {
+        std::size_t step_count =
+            2 * static_cast<std::size_t>(std::floor(neighborhood / step_size)) + 1;
+        // do not use step_size - neighborhood, as step_size might not allow
+        // landing exactly on that value when starting from start_point
+        // also use parentheses on step_count / 2 because flooring is exactly
+        // what we want
+
+        // explicitly fill out members, as aggregate init will error out with narrowing
+        hough_parameter<T> parameter;
+        parameter.start_point = start_point - step_size * (step_count / 2);
+        parameter.step_size = step_size;
+        parameter.step_count = step_count;
+        return parameter;
+    }
+};
+
+/// \ingroup HoughTransform
+/// \brief Calculate minimum angle which would be observable if walked on a circle
+///
+/// When drawing a circle or moving around a point in circular motion, it is
+/// important to not do too many steps, but also to not have disconnected
+/// trajectory. This function will calculate the minimum angle that is observable
+/// when walking on a circle or tilting a line.
+/// WARNING: do keep in mind IEEE 754 quirks, e.g. no-associativity,
+/// no-commutativity and precision. Do not expect expressions that are
+/// mathematically the same to produce the same values
+inline double minimum_angle_step(point_t dimensions)
+{
+    auto longer_dimension = dimensions.x > dimensions.y ? dimensions.x : dimensions.y;
+    return std::atan2(1, longer_dimension);
+}
+
+/// \ingroup HoughTransform
+/// \brief Create a Hough transform parameter with optimal angle step
+///
+/// Due to computational intensity and noise sensitivity of Hough transform,
+/// having any candidates missed or computed again is problematic. This function
+/// will properly encapsulate optimal value range around approx_angle with amplitude of
+/// neighborhood in each direction.
+/// WARNING: do keep in mind IEEE 754 quirks, e.g. no-associativity,
+/// no-commutativity and precision. Do not expect expressions that are
+/// mathematically the same to produce the same values
+inline hough_parameter<double> make_theta_parameter(double approx_angle, double neighborhood,
+                                                    point_t dimensions)
+{
+    auto angle_step = minimum_angle_step(dimensions);
+
+    // std::size_t step_count =
+    //     2 * static_cast<std::size_t>(std::floor(neighborhood / angle_step)) + 1;
+    // return {approx_angle - angle_step * (step_count / 2), angle_step, step_count};
+    return hough_parameter<double>::from_step_size(approx_angle, neighborhood, angle_step);
+}
+}} // namespace boost::gil
+#endif