extension/rasterization/circle.hpp

0001 //
0002 // Copyright 2020 Olzhas Zhumabek <anonymous.from.applecity@gmail.com>
0003 //
0004 // Use, modification and distribution are subject to the Boost Software License,
0005 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0006 // http://www.boost.org/LICENSE_1_0.txt)
0007 //
0008 #ifndef BOOST_GIL_EXTENSION_RASTERIZATION_CIRCLE_HPP
0009 #define BOOST_GIL_EXTENSION_RASTERIZATION_CIRCLE_HPP
0010
0011 #include <boost/gil/detail/math.hpp>
0012 #include <boost/gil/extension/rasterization/apply_rasterizer.hpp>
0013 #include <boost/gil/point.hpp>
0014
0015 #include <cmath>
0016 #include <cstddef>
0017 #include <vector>
0018
0019 namespace boost { namespace gil {
0020
0021 struct circle_rasterizer_t{};
0022
0023 /// \defgroup CircleRasterization
0024 /// \ingroup Rasterization
0025 /// \brief Circle rasterization algorithms
0026 ///
0027 /// The main problems are connectivity and equation following. Circle can be easily moved
0028 /// to new offset, and rotation has no effect on it (not recommended to do rotation).
0029
0030 /// \ingroup CircleRasterization
0031 /// \brief Rasterize trigonometric circle according to radius by sine and radius by cosine
0032 ///
0033 /// This rasterizer is the one used that is used in standard Hough circle transform in
0034 /// the books. It is also quite expensive to compute.
0035 /// WARNING: the product of this rasterizer does not follow circle equation, even though it
0036 /// produces quite round like shapes.
0037 struct trigonometric_circle_rasterizer
0038 {
0039     using type = circle_rasterizer_t;
0040
0041     /// \brief Creates a trigonometric circle rasterizer
0042     /// \param center_point - Point containing positive integer x co-ordinate and y co-ordinate of the
0043     /// center respectively.
0044     /// \param circle_radius - Radius of the circle
0045     trigonometric_circle_rasterizer(point_t center_point, std::ptrdiff_t circle_radius)
0046         : center(center_point), radius(circle_radius)
0047     {}
0048
0049     /// \brief Calculates minimum angle step that is distinguishable when walking on circle
0050     ///
0051     /// It is important to not have disconnected circle and to not compute unnecessarily,
0052     /// thus the result of this function is used when rendering.
0053     double minimum_angle_step() const noexcept
0054     {
0055         const auto diameter = radius * 2 - 1;
0056         return std::atan2(1.0, diameter);
0057     }
0058
0059     /// \brief Calculate the amount of points that rasterizer will output
0060     std::ptrdiff_t point_count() const noexcept
0061     {
0062         return 8 * static_cast<std::ptrdiff_t>(
0063                        std::round(detail::pi / 4 / minimum_angle_step()) + 1);
0064     }
0065
0066     /// \brief perform rasterization and output into d_first
0067     template <typename OutputIterator>
0068     void operator()(OutputIterator d_first) const
0069     {
0070         const double minimum_angle_step = std::atan2(1.0, radius);
0071         auto translate_mirror_points = [this, &d_first](point_t p) {
0072             *d_first++ = point_t{center.x + p.x, center.y + p.y};
0073             *d_first++ = point_t{center.x + p.x, center.y - p.y};
0074             *d_first++ = point_t{center.x - p.x, center.y + p.y};
0075             *d_first++ = point_t{center.x - p.x, center.y - p.y};
0076             *d_first++ = point_t{center.x + p.y, center.y + p.x};
0077             *d_first++ = point_t{center.x + p.y, center.y - p.x};
0078             *d_first++ = point_t{center.x - p.y, center.y + p.x};
0079             *d_first++ = point_t{center.x - p.y, center.y - p.x};
0080         };
0081         const std::ptrdiff_t iteration_count = point_count() / 8;
0082         double angle = 0;
0083         // do note that + 1 was done inside count estimation, thus <= is not needed, only <
0084         for (std::ptrdiff_t i = 0; i < iteration_count; ++i, angle += minimum_angle_step)
0085         {
0086             std::ptrdiff_t x = static_cast<std::ptrdiff_t>(std::round(radius * std::cos(angle)));
0087             std::ptrdiff_t y = static_cast<std::ptrdiff_t>(std::round(radius * std::sin(angle)));
0088             translate_mirror_points({x, y});
0089         }
0090     }
0091
0092     point_t center;
0093     std::ptrdiff_t radius;
0094 };
0095
0096 /// \ingroup CircleRasterization
0097 /// \brief Perform circle rasterization according to Midpoint algorithm
0098 ///
0099 /// This algorithm givess reasonable output and is cheap to compute.
0100 /// reference:
0101 /// https://en.wikipedia.org/wiki/Midpoint_circle_algorithm
0102 struct midpoint_circle_rasterizer
0103 {
0104     using type = circle_rasterizer_t;
0105
0106     /// \brief Creates a midpoint circle rasterizer
0107     /// \param center_point - Point containing positive integer x co-ordinate and y co-ordinate of the
0108     /// center respectively.
0109     /// \param circle_radius - Radius of the circle
0110     midpoint_circle_rasterizer(point_t center_point, std::ptrdiff_t circle_radius)
0111         : center(center_point), radius(circle_radius)
0112     {}
0113
0114     /// \brief Calculate the amount of points that rasterizer will output
0115     std::ptrdiff_t point_count() const noexcept
0116     {
0117         // the reason for pulling 8 out is so that when the expression radius * cos(45 degrees)
0118         // is used, it would yield the same result as here
0119         // + 1 at the end is because the point at radius itself is computed as well
0120         return 8 * static_cast<std::ptrdiff_t>(
0121                        std::round(radius * std::cos(boost::gil::detail::pi / 4)) + 1);
0122     }
0123
0124     /// \brief perform rasterization and output into d_first
0125     template <typename OutputIterator>
0126     void operator()(OutputIterator d_first) const
0127     {
0128         auto translate_mirror_points = [this, &d_first](point_t p) {
0129             *d_first++ = point_t{center.x + p.x, center.y + p.y};
0130             *d_first++ = point_t{center.x + p.x, center.y - p.y};
0131             *d_first++ = point_t{center.x - p.x, center.y + p.y};
0132             *d_first++ = point_t{center.x - p.x, center.y - p.y};
0133             *d_first++ = point_t{center.x + p.y, center.y + p.x};
0134             *d_first++ = point_t{center.x + p.y, center.y - p.x};
0135             *d_first++ = point_t{center.x - p.y, center.y + p.x};
0136             *d_first++ = point_t{center.x - p.y, center.y - p.x};
0137         };
0138         std::ptrdiff_t iteration_distance = point_count() / 8;
0139         std::ptrdiff_t y_current = radius;
0140         std::ptrdiff_t r_squared = radius * radius;
0141         translate_mirror_points({0, y_current});
0142         for (std::ptrdiff_t x = 1; x < iteration_distance; ++x)
0143         {
0144             std::ptrdiff_t midpoint = x * x + y_current * y_current - y_current - r_squared;
0145             if (midpoint > 0)
0146             {
0147                 --y_current;
0148             }
0149             translate_mirror_points({x, y_current});
0150         }
0151     }
0152
0153     point_t center;
0154     std::ptrdiff_t radius;
0155 };
0156
0157 namespace detail {
0158
0159 template <typename View, typename Rasterizer, typename Pixel>
0160 struct apply_rasterizer_op<View, Rasterizer, Pixel, circle_rasterizer_t>
0161 {
0162     void operator()(
0163         View const& view, Rasterizer const& rasterizer, Pixel const& pixel)
0164     {
0165         std::vector<point_t> trajectory(rasterizer.point_count());
0166         rasterizer(std::begin(trajectory));
0167
0168         for (auto const& point : trajectory)
0169         {
0170             view(point) = pixel;
0171         }
0172     }
0173 };
0174
0175 } //namespace detail
0176
0177 }} // namespace boost::gil
0178
0179 #endif