From daa9d3fcd071ef5b7b2bffdd1a62fd79cbffdaa7 Mon Sep 17 00:00:00 2001
From: Barend Gehrels <barend@xs4all.nl>
Date: Mon, 27 Oct 2014 23:27:31 +0100
Subject: [PATCH] [buffer][fix] fix round join strategy, it should not use acos
 for non-normalized dot product. Now we use atan2 to calculate the angles

---
 .../cartesian/buffer_join_round.hpp           | 41 ++++++++-----------
 1 file changed, 17 insertions(+), 24 deletions(-)

diff --git a/include/boost/geometry/strategies/cartesian/buffer_join_round.hpp b/include/boost/geometry/strategies/cartesian/buffer_join_round.hpp
index 9e467c85a..7e1ad838f 100644
--- a/include/boost/geometry/strategies/cartesian/buffer_join_round.hpp
+++ b/include/boost/geometry/strategies/cartesian/buffer_join_round.hpp
@@ -69,34 +69,27 @@ private :
                 DistanceType const& buffer_distance,
                 RangeOut& range_out) const
     {
-        PromotedType dx1 = get<0>(perp1) - get<0>(vertex);
-        PromotedType dy1 = get<1>(perp1) - get<1>(vertex);
-        PromotedType dx2 = get<0>(perp2) - get<0>(vertex);
-        PromotedType dy2 = get<1>(perp2) - get<1>(vertex);
+        PromotedType const dx1 = get<0>(perp1) - get<0>(vertex);
+        PromotedType const dy1 = get<1>(perp1) - get<1>(vertex);
+        PromotedType const dx2 = get<0>(perp2) - get<0>(vertex);
+        PromotedType const dy2 = get<1>(perp2) - get<1>(vertex);
 
-        BOOST_ASSERT(buffer_distance != 0);
-
-        dx1 /= buffer_distance;
-        dy1 /= buffer_distance;
-        dx2 /= buffer_distance;
-        dy2 /= buffer_distance;
-
-        PromotedType angle_diff = acos(dx1 * dx2 + dy1 * dy2);
-
-        PromotedType two = 2.0;
-        PromotedType steps = m_points_per_circle;
-        int n = boost::numeric_cast<int>(steps * angle_diff
-                    / (two * geometry::math::pi<PromotedType>()));
-
-        if (n <= 1)
-        {
-            return;
-        }
+        PromotedType const two = 2.0;
 
         PromotedType const angle1 = atan2(dy1, dx1);
-        PromotedType diff = angle_diff / PromotedType(n);
-        PromotedType a = angle1 - diff;
+        PromotedType angle2 = atan2(dy2, dx2);
+        while (angle2 > angle1)
+        {
+            angle2 -= 2 * M_PI;
+        }
 
+        // Divide the angle into an integer amount of steps to make it
+        // visually correct also for a low number of points / circle
+        int const n = static_cast<int>((m_points_per_circle * (angle1 - angle2)
+                            / (two * geometry::math::pi<PromotedType>())));
+
+        PromotedType const diff = (angle1 - angle2) / static_cast<PromotedType>(n);
+        PromotedType a = angle1 - diff;
         for (int i = 0; i < n - 1; i++, a -= diff)
         {
             Point p;