[fix] support different geometry types

This fixes #1139 , and more required code changes.

Including unit test.

Author: barendgehrels

doc/release_notes.qbk

@@ -30,6 +30,7 @@
 [*Solved issues]
 
 * [@https://github.com/boostorg/geometry/issues/1100 1100] Fix for union
+* [@https://github.com/boostorg/geometry/issues/1139 1139] Fix for different geometry types
 
 [*Breaking changes]
 

include/boost/geometry/algorithms/detail/direction_code.hpp

@@ -46,16 +46,16 @@ struct direction_code_impl
 template <>
 struct direction_code_impl<cartesian_tag>
 {
-    template <typename Point1, typename Point2>
-    static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
+    template <typename PointSegmentA, typename PointSegmentB, typename Point2>
+    static inline int apply(PointSegmentA const& segment_a, PointSegmentB const& segment_b,
                             Point2 const& point)
     {
-        typedef typename geometry::select_coordinate_type
+        using calc_t = typename geometry::select_coordinate_type
             <
-                Point1, Point2
-            >::type calc_t;
+                PointSegmentA, PointSegmentB, Point2
+            >::type;
 
-        typedef model::infinite_line<calc_t> line_type;
+        using line_type = model::infinite_line<calc_t>;
 
         // Situation and construction of perpendicular line
         //
@@ -87,42 +87,59 @@ struct direction_code_impl<cartesian_tag>
 template <>
 struct direction_code_impl<spherical_equatorial_tag>
 {
-    template <typename Point1, typename Point2>
-    static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
+    template <typename PointSegmentA, typename PointSegmentB, typename Point2>
+    static inline int apply(PointSegmentA const& segment_a, PointSegmentB const& segment_b,
                             Point2 const& p)
     {
-        typedef typename coordinate_type<Point1>::type coord1_t;
-        typedef typename coordinate_type<Point2>::type coord2_t;
-        typedef typename cs_angular_units<Point1>::type units_t;
-        typedef typename cs_angular_units<Point2>::type units2_t;
-        BOOST_GEOMETRY_STATIC_ASSERT(
-            (std::is_same<units_t, units2_t>::value),
-            "Not implemented for different units.",
-            units_t, units2_t);
-
-        typedef typename geometry::select_coordinate_type <Point1, Point2>::type calc_t;
-        typedef math::detail::constants_on_spheroid<coord1_t, units_t> constants1;
-        typedef math::detail::constants_on_spheroid<coord2_t, units_t> constants2;
-        static coord1_t const pi_half1 = constants1::max_latitude();
-        static coord2_t const pi_half2 = constants2::max_latitude();
+        {
+            using units_sa_t =  typename cs_angular_units<PointSegmentA>::type;
+            using units_sb_t =  typename cs_angular_units<PointSegmentB>::type;
+            using units_p_t = typename cs_angular_units<Point2>::type;
+            BOOST_GEOMETRY_STATIC_ASSERT(
+                (std::is_same<units_sa_t, units_sb_t>::value),
+                "Not implemented for different units.",
+                units_sa_t, units_sb_t);
+            BOOST_GEOMETRY_STATIC_ASSERT(
+                (std::is_same<units_sa_t, units_p_t>::value),
+                "Not implemented for different units.",
+                units_sa_t, units_p_t);
+        }
+
+        using coor_sa_t = typename coordinate_type<PointSegmentA>::type;
+        using coor_sb_t = typename coordinate_type<PointSegmentB>::type;
+        using coor_p_t = typename coordinate_type<Point2>::type;
+
+        // Declare unit type (equal for all types) and calc type (coerced to most precise)
+        using units_t = typename cs_angular_units<Point2>::type;
+        using calc_t = typename geometry::select_coordinate_type
+            <
+                PointSegmentA, PointSegmentB, Point2
+            >::type;
+        using constants_sa_t = math::detail::constants_on_spheroid<coor_sa_t, units_t>;
+        using constants_sb_t = math::detail::constants_on_spheroid<coor_sb_t, units_t>;
+        using constants_p_t = math::detail::constants_on_spheroid<coor_p_t, units_t>;
+
+        static coor_sa_t const pi_half_sa = constants_sa_t::max_latitude();
+        static coor_sb_t const pi_half_sb = constants_sb_t::max_latitude();
+        static coor_p_t const pi_half_p = constants_p_t::max_latitude();
         static calc_t const c0 = 0;
 
-        coord1_t const a0 = geometry::get<0>(segment_a);
-        coord1_t const a1 = geometry::get<1>(segment_a);
-        coord1_t const b0 = geometry::get<0>(segment_b);
-        coord1_t const b1 = geometry::get<1>(segment_b);
-        coord2_t const p0 = geometry::get<0>(p);
-        coord2_t const p1 = geometry::get<1>(p);
-        
+        coor_sa_t const a0 = geometry::get<0>(segment_a);
+        coor_sa_t const a1 = geometry::get<1>(segment_a);
+        coor_sb_t const b0 = geometry::get<0>(segment_b);
+        coor_sb_t const b1 = geometry::get<1>(segment_b);
+        coor_p_t const p0 = geometry::get<0>(p);
+        coor_p_t const p1 = geometry::get<1>(p);
+
         if ( (math::equals(b0, a0) && math::equals(b1, a1))
           || (math::equals(b0, p0) && math::equals(b1, p1)) )
         {
             return 0;
         }
 
-        bool const is_a_pole = math::equals(pi_half1, math::abs(a1));
-        bool const is_b_pole = math::equals(pi_half1, math::abs(b1));
-        bool const is_p_pole = math::equals(pi_half2, math::abs(p1));
+        bool const is_a_pole = math::equals(pi_half_sa, math::abs(a1));
+        bool const is_b_pole = math::equals(pi_half_sb, math::abs(b1));
+        bool const is_p_pole = math::equals(pi_half_p, math::abs(p1));
 
         if ( is_b_pole && ((is_a_pole && math::sign(b1) == math::sign(a1))
                         || (is_p_pole && math::sign(b1) == math::sign(p1))) )
@@ -140,12 +157,12 @@ struct direction_code_impl<spherical_equatorial_tag>
         calc_t const dlat1 = latitude_distance_signed<units_t, calc_t>(b1, a1, dlon1, is_antilon1);
         calc_t const dlat2 = latitude_distance_signed<units_t, calc_t>(b1, p1, dlon2, is_antilon2);
 
-        calc_t mx = is_a_pole || is_b_pole || is_p_pole ?
-                    c0 :
-                    (std::min)(is_antilon1 ? c0 : math::abs(dlon1),
-                               is_antilon2 ? c0 : math::abs(dlon2));
-        calc_t my = (std::min)(math::abs(dlat1),
-                               math::abs(dlat2));
+        calc_t const mx = is_a_pole || is_b_pole || is_p_pole
+            ? c0
+            : (std::min)(is_antilon1 ? c0 : math::abs(dlon1),
+                         is_antilon2 ? c0 : math::abs(dlon2));
+        calc_t const my = (std::min)(math::abs(dlat1),
+                                     math::abs(dlat2));
 
         int s1 = 0, s2 = 0;
         if (mx >= my)
@@ -165,7 +182,7 @@ struct direction_code_impl<spherical_equatorial_tag>
     template <typename Units, typename T>
     static inline T latitude_distance_signed(T const& lat1, T const& lat2, T const& lon_ds, bool & is_antilon)
     {
-        typedef math::detail::constants_on_spheroid<T, Units> constants;
+        using constants = math::detail::constants_on_spheroid<T, Units>;
         static T const pi = constants::half_period();
         static T const c0 = 0;
 
@@ -188,27 +205,27 @@ struct direction_code_impl<spherical_equatorial_tag>
 template <>
 struct direction_code_impl<spherical_polar_tag>
 {
-    template <typename Point1, typename Point2>
-    static inline int apply(Point1 segment_a, Point1 segment_b,
+    template <typename PointSegmentA, typename PointSegmentB, typename Point2>
+    static inline int apply(PointSegmentA segment_a, PointSegmentB segment_b,
                             Point2 p)
     {
-        typedef math::detail::constants_on_spheroid
+        using constants_sa_t = math::detail::constants_on_spheroid
             <
-                typename coordinate_type<Point1>::type,
-                typename cs_angular_units<Point1>::type
-            > constants1;
-        typedef math::detail::constants_on_spheroid
+                typename coordinate_type<PointSegmentA>::type,
+                typename cs_angular_units<PointSegmentA>::type
+            >;
+        using constants_p_t = math::detail::constants_on_spheroid
             <
                 typename coordinate_type<Point2>::type,
                 typename cs_angular_units<Point2>::type
-            > constants2;
+            >;
 
         geometry::set<1>(segment_a,
-            constants1::max_latitude() - geometry::get<1>(segment_a));
+            constants_sa_t::max_latitude() - geometry::get<1>(segment_a));
         geometry::set<1>(segment_b,
-            constants1::max_latitude() - geometry::get<1>(segment_b));
+            constants_sa_t::max_latitude() - geometry::get<1>(segment_b));
         geometry::set<1>(p,
-            constants2::max_latitude() - geometry::get<1>(p));
+            constants_p_t::max_latitude() - geometry::get<1>(p));
 
         return direction_code_impl
                 <
@@ -217,13 +234,13 @@ struct direction_code_impl<spherical_polar_tag>
     }
 };
 
-// if spherical_tag is passed then pick cs_tag based on Point1 type
+// if spherical_tag is passed then pick cs_tag based on PointSegmentA type
 // with spherical_equatorial_tag as the default
 template <>
 struct direction_code_impl<spherical_tag>
 {
-    template <typename Point1, typename Point2>
-    static inline int apply(Point1 segment_a, Point1 segment_b,
+    template <typename PointSegmentA, typename PointSegmentB, typename Point2>
+    static inline int apply(PointSegmentA segment_a, PointSegmentB segment_b,
                             Point2 p)
     {
         return direction_code_impl
@@ -232,7 +249,7 @@ struct direction_code_impl<spherical_tag>
                     <
                         std::is_same
                             <
-                                typename geometry::cs_tag<Point1>::type,
+                                typename geometry::cs_tag<PointSegmentA>::type,
                                 spherical_polar_tag
                             >::value,
                         spherical_polar_tag,
@@ -252,8 +269,10 @@ struct direction_code_impl<geographic_tag>
 // Returns 1 if p goes forward, so extends (a,b)
 // Returns 0 if p is equal with b, or if (a,b) is degenerate
 // Note that it does not do any collinearity test, that should be done before
-template <typename CSTag, typename Point1, typename Point2>
-inline int direction_code(Point1 const& segment_a, Point1 const& segment_b,
+// In some cases the "segment" consists of different source points, and therefore
+// their types might differ.
+template <typename CSTag, typename PointSegmentA, typename PointSegmentB, typename Point2>
+inline int direction_code(PointSegmentA const& segment_a, PointSegmentB const& segment_b,
                           Point2 const& p)
 {
     return direction_code_impl<CSTag>::apply(segment_a, segment_b, p);

include/boost/geometry/algorithms/detail/make/make.hpp

@@ -19,10 +19,10 @@ namespace boost { namespace geometry
 namespace detail { namespace make
 {
 
-template <typename Type, typename Coordinate>
+template <typename Type, typename Coordinate1, typename Coordinate2>
 inline
-model::infinite_line<Type> make_infinite_line(Coordinate const& x1,
-    Coordinate const& y1, Coordinate const& x2, Coordinate const& y2)
+model::infinite_line<Type> make_infinite_line(Coordinate1 const& x1,
+    Coordinate1 const& y1, Coordinate2 const& x2, Coordinate2 const& y2)
 {
     model::infinite_line<Type> result;
     result.a = y1 - y2;
@@ -31,9 +31,9 @@ model::infinite_line<Type> make_infinite_line(Coordinate const& x1,
     return result;
 }
 
-template <typename Type, typename Point>
+template <typename Type, typename PointA, typename PointB>
 inline
-model::infinite_line<Type> make_infinite_line(Point const& a, Point const& b)
+model::infinite_line<Type> make_infinite_line(PointA const& a, PointB const& b)
 {
     return make_infinite_line<Type>(geometry::get<0>(a), geometry::get<1>(a),
         geometry::get<0>(b), geometry::get<1>(b));
@@ -49,9 +49,9 @@ model::infinite_line<Type> make_infinite_line(Segment const& segment)
         geometry::get<1, 1>(segment));
 }
 
-template <typename Type, typename Point>
+template <typename Type, typename PointA, typename PointB, typename PointC>
 inline
-model::infinite_line<Type> make_perpendicular_line(Point const& a, Point const& b, Point const& c)
+model::infinite_line<Type> make_perpendicular_line(PointA const& a, PointB const& b, PointC const& c)
 {
     // https://www.math-only-math.com/equation-of-a-line-perpendicular-to-a-line.html
     model::infinite_line<Type> const line = make_infinite_line<Type>(a, b);

include/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp

@@ -23,6 +23,7 @@
 #include <boost/static_assert.hpp>
 
 #include <boost/geometry/algorithms/append.hpp>
+#include <boost/geometry/algorithms/detail/convert_point_to_point.hpp>
 #include <boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp>
 #include <boost/geometry/algorithms/detail/equals/point_point.hpp>
 
@@ -99,7 +100,16 @@ inline void append_no_dups_or_spikes(Range& range, Point const& point,
         return;
     }
 
-    traits::push_back<Range>::apply(range, point);
+    auto append = [](auto& r, const auto& p)
+    {
+        using point_t = typename boost::range_value<Range>::type;
+        point_t rp;
+        geometry::detail::conversion::convert_point_to_point(p, rp);
+        traits::push_back<Range>::apply(r, std::move(rp));
+    };
+
+    append(range, point);
+
 
     // If a point is equal, or forming a spike, remove the pen-ultimate point
     // because this one caused the spike.
@@ -115,7 +125,7 @@ inline void append_no_dups_or_spikes(Range& range, Point const& point,
     {
         // Use the Concept/traits, so resize and append again
         traits::resize<Range>::apply(range, boost::size(range) - 2);
-        traits::push_back<Range>::apply(range, point);
+        append(range, point);
     }
 }
 

include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp

@@ -192,13 +192,16 @@ template<typename VerifyPolicy>
 struct turn_info_verification_functions
 {
     template <typename Point1, typename Point2>
-    static inline typename geometry::coordinate_type<Point1>::type
-            distance_measure(Point1 const& a, Point2 const& b)
+    static inline
+    typename select_coordinate_type<Point1, Point2>::type
+    distance_measure(Point1 const& a, Point2 const& b)
     {
         // TODO: revise this using comparable distance for various
         // coordinate systems
-        auto const dx = get<0>(a) - get<0>(b);
-        auto const dy = get<1>(a) - get<1>(b);
+        using coor_t = typename select_coordinate_type<Point1, Point2>::type;
+
+        coor_t const dx = get<0>(a) - get<0>(b);
+        coor_t const dy = get<1>(a) - get<1>(b);
         return dx * dx + dy * dy;
     }
 
@@ -1032,8 +1035,9 @@ struct collinear : public base_turn_handler
             return false;
         }
 
-        auto const dm = fun::distance_measure(info.intersections[1],
-                arrival_p == 1 ? range_q.at(1) : range_p.at(1));
+       auto const dm = arrival_p == 1
+              ? fun::distance_measure(info.intersections[1], range_q.at(1))
+              : fun::distance_measure(info.intersections[1], range_p.at(1));
         decltype(dm) const zero = 0;
         return math::equals(dm, zero);
     }

include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp

@@ -127,10 +127,10 @@ struct less_false
     }
 };
 
-template <typename Point, typename SideStrategy, typename LessOnSame, typename Compare>
+template <typename PointOrigin, typename PointTurn, typename SideStrategy, typename LessOnSame, typename Compare>
 struct less_by_side
 {
-    less_by_side(const Point& p1, const Point& p2, SideStrategy const& strategy)
+    less_by_side(const PointOrigin& p1, const PointTurn& p2, SideStrategy const& strategy)
         : m_origin(p1)
         , m_turn_point(p2)
         , m_strategy(strategy)
@@ -209,8 +209,8 @@ struct less_by_side
     }
 
 private :
-    Point const& m_origin;
-    Point const& m_turn_point;
+    PointOrigin const& m_origin;
+    PointTurn const& m_turn_point;
     SideStrategy const& m_strategy;
 };
 
@@ -379,7 +379,8 @@ public :
         }
     }
 
-    void apply(Point const& turn_point)
+    template <typename PointTurn>
+    void apply(PointTurn const& turn_point)
     {
         // We need three compare functors:
         // 1) to order clockwise (union) or counter clockwise (intersection)
@@ -388,8 +389,8 @@ public :
         //    to give colinear points
 
         // Sort by side and assign rank
-        less_by_side<Point, SideStrategy, less_by_index, Compare> less_unique(m_origin, turn_point, m_strategy);
-        less_by_side<Point, SideStrategy, less_false, Compare> less_non_unique(m_origin, turn_point, m_strategy);
+        less_by_side<Point, PointTurn, SideStrategy, less_by_index, Compare> less_unique(m_origin, turn_point, m_strategy);
+        less_by_side<Point, PointTurn, SideStrategy, less_false, Compare> less_non_unique(m_origin, turn_point, m_strategy);
 
         std::sort(m_ranked_points.begin(), m_ranked_points.end(), less_unique);