Add rasterizer support for ellipse

example/Jamfile

@@ -32,6 +32,7 @@ local sources =
     mandelbrot.cpp
     morphology.cpp
     packed_pixel.cpp
+    rasterizer_ellipse.cpp
     resize.cpp
     sobel_scharr.cpp
     threshold.cpp

example/rasterizer_ellipse.cpp

@@ -0,0 +1,44 @@
+//
+// Copyright 2021 Prathamesh Tagore <prathameshtagore@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/extension/io/jpeg.hpp>
+#include <boost/gil.hpp>
+
+namespace gil = boost::gil;
+
+int main()
+{
+    // Syntax for usage :- 
+    // auto rasterizer = gil::midpoint_elliptical_rasterizer{};
+    // rasterizer(img_view, colour, center, semi-axes_length);
+    // Where
+    // img_view : gil view of the image on which ellipse is to be drawn.
+    // colour : Vector containing channel intensity values for img_view. Number of colours 
+    // provided must be equal to the number of channels present in img_view.
+    // center : Array containing positive integer x co-ordinate and y co-ordinate of the center
+    // respectively.
+    // semi-axes_length : Array containing positive integer lengths of horizontal semi-axis
+    // and vertical semi-axis respectively.
+
+    gil::gray8_image_t gray_buffer_image(256, 256);
+    auto gray_elliptical_rasterizer = gil::midpoint_elliptical_rasterizer{};
+    gray_elliptical_rasterizer(view(gray_buffer_image), {128}, {128, 128}, {100, 50});
+
+    gil::rgb8_image_t rgb_buffer_image(256, 256);
+    auto rgb_elliptical_rasterizer = gil::midpoint_elliptical_rasterizer{};
+    rgb_elliptical_rasterizer(view(rgb_buffer_image), {0, 0, 255}, {128, 128}, {50, 100});
+
+    gil::rgb8_image_t rgb_buffer_image_out_of_bound(256, 256);
+    auto rgb_elliptical_rasterizer_out_of_bound = gil::midpoint_elliptical_rasterizer{};
+    rgb_elliptical_rasterizer_out_of_bound(view(rgb_buffer_image_out_of_bound), {255, 0, 0},
+        {100, 100}, {160, 160});
+
+    gil::write_view("rasterized_ellipse_gray.jpg", view(gray_buffer_image), gil::jpeg_tag{});
+    gil::write_view("rasterized_ellipse_rgb.jpg", view(rgb_buffer_image), gil::jpeg_tag{});
+    gil::write_view("rasterized_ellipse_rgb_out_of_bound.jpg", view(rgb_buffer_image_out_of_bound),
+        gil::jpeg_tag{});
+}

include/boost/gil.hpp

@@ -41,6 +41,7 @@
 #include <boost/gil/premultiply.hpp>
 #include <boost/gil/promote_integral.hpp>
 #include <boost/gil/rasterization/circle.hpp>
+#include <boost/gil/rasterization/ellipse.hpp>
 #include <boost/gil/rasterization/line.hpp>
 #include <boost/gil/rgb.hpp>
 #include <boost/gil/rgba.hpp>

include/boost/gil/rasterization/ellipse.hpp

@@ -0,0 +1,192 @@
+//
+// Copyright 2021 Prathamesh Tagore <prathameshtagore@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_RASTERIZATION_ELLIPSE_HPP
+#define BOOST_GIL_RASTERIZATION_ELLIPSE_HPP
+
+#include <array>
+#include <vector>
+#include <iostream>
+
+namespace boost { namespace gil {
+
+/// \defgroup EllipseRasterization
+/// \ingroup Rasterization
+/// \brief Ellipse rasterization algorithms.
+
+/// \ingroup EllipseRasterization
+/// \brief Performs ellipse rasterization using midpoint algorithm. Initially, program considers
+/// origin as center of ellipse and obtains first quadrant trajectory points. After that,
+/// it shifts origin to provided co-ordinates of center and then draws the curve.
+struct midpoint_elliptical_rasterizer
+{
+    /// \brief Returns a vector containing co-ordinates of first quadrant points which lie on
+    /// rasterizer trajectory of the ellipse.
+    /// \param semi_axes - Array containing half of lengths of horizontal and vertical axis
+    /// respectively.
+    auto obtain_trajectory(std::array<unsigned int, 2> const semi_axes)
+        -> std::vector<std::array<std::ptrdiff_t, 2>>
+    {
+        // Citation : J. Van Aken, "An Efficient Ellipse-Drawing Algorithm" in IEEE Computer
+        // Graphics and Applications, vol. 4, no. 09, pp. 24-35, 1984.
+        // doi: 10.1109/MCG.1984.275994
+        // keywords: {null}
+        // url: https://doi.ieeecomputersociety.org/10.1109/MCG.1984.275994
+        std::vector<std::array<std::ptrdiff_t, 2>> trajectory_points;
+        std::ptrdiff_t x = semi_axes[0], y = 0;
+
+        // Variables declared on following lines are temporary variables used for improving 
+        // performance since they help in converting all multiplicative operations inside the while
+        // loop into additive/subtractive operations.
+        long long int const t1 = semi_axes[0] * semi_axes[0];
+        long long int const t4 = semi_axes[1] * semi_axes[1];
+        long long int t2, t3, t5, t6, t8, t9;
+        t2 = 2 * t1, t3 = 2 * t2;
+        t5 = 2 * t4, t6 = 2 * t5;
+        long long int const t7 = semi_axes[0] * t5;
+        t8 = 2 * t7, t9 = 0;
+
+        // Following variables serve as decision parameters and help in choosing the right point
+        // to be included in rasterizer trajectory.
+        long long int d1, d2;
+        d1 = t2 - t7 + t4 / 2, d2 = t1 / 2 - t8 + t5;
+
+        while (d2 < 0)
+        {
+            trajectory_points.push_back({x, y});
+            y += 1;
+            t9 += t3;
+            if (d1 < 0)
+            {
+                d1 += t9 + t2;
+                d2 += t9;
+            }
+            else 
+            {
+                x -= 1;
+                t8 -= t6;
+                d1 += t9 + t2 - t8;
+                d2 += t5 + t9 - t8;
+            }
+        }
+        while (x >= 0)
+        {
+            trajectory_points.push_back({x, y});
+            x -= 1;
+            t8 -= t6;
+            if (d2 < 0)
+            {
+                y += 1;
+                t9 += t3;
+                d2 += t5 + t9 - t8;
+            }
+            else 
+            {
+                d2 += t5 - t8;
+            }
+        }
+        return trajectory_points;
+    }
+
+    /// \brief Fills pixels returned by function 'obtain_trajectory' as well as pixels
+    /// obtained from their reflection along major axis, minor axis and line passing through
+    /// center with slope -1 using colours provided by user.
+    /// \param view - Gil view of image on which the elliptical curve is to be drawn.
+    /// \param colour - Constant vector specifying colour intensity values for all channels present
+    ///                 in 'view'.
+    /// \param center - Constant array specifying co-ordinates of center of ellipse to be drawn.
+    /// \param trajectory_points - Constant vector specifying pixel co-ordinates of points lying 
+    ///                            on rasterizer trajectory.
+    /// \tparam View - Type of input image view.
+    template<typename View>
+    void draw_curve(View view, std::vector<unsigned int> const colour,
+        std::array<unsigned int, 2> const center,
+        std::vector<std::array<std::ptrdiff_t, 2>> const trajectory_points)
+    {
+        for (int i = 0, colour_index = 0; i < static_cast<int>(view.num_channels()); 
+            ++i, ++colour_index)
+        {
+            for (std::array<std::ptrdiff_t, 2> point : trajectory_points)
+            {
+                std::array<std::ptrdiff_t, 4> co_ords = {center[0] + point[0],
+                center[0] - point[0], center[1] + point[1], center[1] - point[1]
+                };
+                bool validity[4] = {0};
+                if (co_ords[0] < view.width())
+                {
+                    validity[0] = 1;
+                }
+                if (co_ords[1] >= 0 && co_ords[1] < view.width())
+                {
+                    validity[1] = 1;
+                }
+                if (co_ords[2] < view.height())
+                {
+                    validity[2] = 1;
+                }
+                if (co_ords[3] >= 0 && co_ords[3] < view.height())
+                {
+                    validity[3] = 1;
+                }
+                if (validity[0] && validity[2])
+                {
+                    nth_channel_view(view, i)(co_ords[0], co_ords[2])[0] = colour[colour_index];
+                }
+                if (validity[1] && validity[2])
+                {
+                    nth_channel_view(view, i)(co_ords[1], co_ords[2])[0] = colour[colour_index];
+                }
+                if (validity[1] && validity[3])
+                {
+                    nth_channel_view(view, i)(co_ords[1], co_ords[3])[0] = colour[colour_index];
+                }
+                if (validity[0] && validity[3])
+                {
+                    nth_channel_view(view, i)(co_ords[0], co_ords[3])[0] = colour[colour_index];
+                }
+            }
+        }
+    }
+
+    /// \brief Calls the function 'obtain_trajectory' and then passes obtained trajectory points
+    ///        in the function 'draw_curve' for drawing the desired ellipse.
+    /// \param view - Gil view of image on which the elliptical curve is to be drawn.
+    /// \param colour - Constant vector specifying colour intensity values for all channels present
+    ///                 in 'view'.
+    /// \param center - Array containing positive integer x co-ordinate and y co-ordinate of the
+    /// center respectively.
+    /// \param semi_axes - Array containing positive integer lengths of horizontal semi-axis
+    /// and vertical semi-axis respectively.
+    /// \tparam View - Type of input image view.
+    template<typename View>
+    void operator()(View view, std::vector<unsigned int> const colour,
+        std::array<unsigned int, 2> center, std::array<unsigned int, 2> const semi_axes)
+    {
+        --center[0], --center[1]; // For converting center co-ordinate values to zero based indexing.
+        if (colour.size() != view.num_channels())
+        {
+            throw std::length_error("Number of channels in given image is not equal to the "
+                "number of colours provided.");
+        }
+        if (center[0] + semi_axes[0] >= view.width() || center[1] + semi_axes[1] >= view.height()
+            || static_cast<int>(center[0] - semi_axes[0]) < 0 
+            || static_cast<int>(center[0] - semi_axes[0]) >= view.width()
+            || static_cast<int>(center[1] - semi_axes[1]) < 0
+            || static_cast<int>(center[1] - semi_axes[1]) >= view.height())
+        {
+            std::cout << "Image can't contain whole curve.\n"
+                "However, it will contain those parts of curve which can fit inside it.\n"
+                "Note : Image width = " << view.width() << " and Image height = " << 
+                view.height() << "\n";
+        }
+        std::vector<std::array<std::ptrdiff_t, 2>> trajectory_points = 
+            obtain_trajectory(semi_axes);
+        draw_curve(view, colour, center, trajectory_points);
+    }
+}; // midpoint elliptical rasterizer
+}} // namespace boost::gil
+#endif

test/core/rasterization/CMakeLists.txt

@@ -8,7 +8,8 @@
 #
 foreach(_name
     line
-    circle)
+    circle
+    ellipse)
     set(_test t_core_rasterization_${_name})
     set(_target test_core_rasterization_${_name})
 

test/core/rasterization/Jamfile

@@ -10,3 +10,4 @@ import testing ;
 
 run line.cpp ;
 run circle.cpp ;
+run ellipse.cpp ;

test/core/rasterization/ellipse.cpp

@@ -0,0 +1,83 @@
+//
+// Copyright 2021 Prathamesh Tagore <prathameshtagore@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/core/lightweight_test.hpp>
+#include <boost/gil.hpp>
+#include <array>
+#include <cmath>
+#include <vector>
+
+namespace gil = boost::gil;
+
+// This function utilizes the fact that sum of distances of a point on an ellipse from its foci 
+// is equal to the length of major axis of the ellipse.
+// Parameters b and a represent half of lengths of vertical and horizontal axis respectively.
+void test_rasterizer_follows_equation(
+    std::vector<std::array<std::ptrdiff_t, 2>> trajectory_points, float a, float b)
+{
+    float focus_x, focus_y;
+    if (a > b) // For horizontal ellipse
+    {
+        focus_x = a * std::sqrt(1 - b * b / (a * a));
+        focus_y = 0;
+    }
+    else // For vertical ellipse
+    {
+        focus_x = 0;
+        focus_y = b * std::sqrt(1 - a * a / (b * b));
+    }
+
+    for (auto trajectory_point : trajectory_points)
+    {
+        // To suppress conversion warnings from compiler
+        std::array<float, 2> point {
+            static_cast<float>(trajectory_point[0]), static_cast<float>(trajectory_point[1])};
+
+        double dist_sum = std::sqrt(std::pow(focus_x - point[0], 2) + 
+            std::pow(focus_y - point[1], 2)) + std::sqrt(std::pow( - focus_x - point[0], 2) + 
+            std::pow( - focus_y - point[1], 2));
+        if (a > b)
+        {
+            BOOST_TEST(std::abs(dist_sum - 2 * a) < 1);
+        }
+        else 
+        {
+            BOOST_TEST(std::abs(dist_sum - 2 * b) < 1);
+        }
+    }
+}
+
+// This function verifies that the difference between x co-ordinates and y co-ordinates for two
+// successive trajectory points is less than or equal to 1. This ensures that the curve is connected.
+void test_connectivity(std::vector<std::array<std::ptrdiff_t, 2>> points)
+{
+    for (std::size_t i = 1; i < points.size(); ++i)
+    {
+        std::ptrdiff_t diff_x = points[i][0] - points[i - 1][0];
+        std::ptrdiff_t diff_y = points[i][1] - points[i - 1][1];
+        BOOST_TEST_LE(diff_x, 1);
+        BOOST_TEST_LE(diff_y, 1);
+    }
+}
+
+// We verify all test cases for the portion of ellipse in first quadrant, since all other portions
+// can be constructed with simple reflection, they tend to be correct if first quadrant is verified.
+int main()
+{
+    for (float a = 1; a < 101; ++a)
+    {
+        for (float b = 1; b < 101; ++b)
+        {
+            auto rasterizer = gil::midpoint_elliptical_rasterizer{};
+            std::vector<std::array<std::ptrdiff_t, 2>> points = rasterizer.obtain_trajectory(
+                {static_cast<unsigned int>(a), static_cast<unsigned int>(b)});
+            test_rasterizer_follows_equation(points, a, b);
+            test_connectivity(points);
+        }
+    }
+    return boost::report_errors();
+}