[ceph.git] / ceph / src / boost / boost / geometry / algorithms / detail / direction_code.hpp

// Boost.Geometry (aka GGL, Generic Geometry Library)

// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.

// This file was modified by Oracle on 2015, 2017, 2019.
// Modifications copyright (c) 2015-2019 Oracle and/or its affiliates.

// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is 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_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP


#include <boost/geometry/core/access.hpp>
#include <boost/geometry/arithmetic/infinite_line_functions.hpp>
#include <boost/geometry/algorithms/detail/make/make.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_coordinate_type.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>

#include <boost/mpl/assert.hpp>


namespace boost { namespace geometry
{


#ifndef DOXYGEN_NO_DETAIL
namespace detail
{

template <typename CSTag>
struct direction_code_impl
{
    BOOST_MPL_ASSERT_MSG((false), NOT_IMPLEMENTED_FOR_THIS_CS, (CSTag));
};

template <>
struct direction_code_impl<cartesian_tag>
{
    template <typename Point1, typename Point2>
    static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
                            Point2 const& point)
    {
        typedef typename geometry::select_coordinate_type
            <
                Point1, Point2
            >::type calc_t;

        typedef model::infinite_line<calc_t> line_type;

        // point b is often equal to the specified point, check that first
        line_type const q = detail::make::make_infinite_line<calc_t>(segment_b, point);
        if (arithmetic::is_degenerate(q))
        {
            return 0;
        }

        line_type const p = detail::make::make_infinite_line<calc_t>(segment_a, segment_b);
        if (arithmetic::is_degenerate(p))
        {
            return 0;
        }

        // p extends a-b if direction is similar
        return arithmetic::similar_direction(p, q) ? 1 : -1;
    }
};

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,
                            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_MPL_ASSERT_MSG((boost::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();
        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);
        
        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));

        if ( is_b_pole && ((is_a_pole && math::sign(b1) == math::sign(a1))
                        || (is_p_pole && math::sign(b1) == math::sign(p1))) )
        {
            return 0;
        }

        // NOTE: as opposed to the implementation for cartesian CS
        // here point b is the origin

        calc_t const dlon1 = math::longitude_distance_signed<units_t, calc_t>(b0, a0);
        calc_t const dlon2 = math::longitude_distance_signed<units_t, calc_t>(b0, p0);

        bool is_antilon1 = false, is_antilon2 = false;
        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));

        int s1 = 0, s2 = 0;
        if (mx >= my)
        {
            s1 = dlon1 > 0 ? 1 : -1;
            s2 = dlon2 > 0 ? 1 : -1;
        }
        else
        {
            s1 = dlat1 > 0 ? 1 : -1;
            s2 = dlat2 > 0 ? 1 : -1;
        }

        return s1 == s2 ? -1 : 1;
    }

    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;
        static T const pi = constants::half_period();
        static T const c0 = 0;

        T res = lat2 - lat1;

        is_antilon = math::equals(math::abs(lon_ds), pi);
        if (is_antilon)
        {
            res = lat2 + lat1;
            if (res >= c0)
                res = pi - res;
            else
                res = -pi - res;
        }

        return res;
    }
};

template <>
struct direction_code_impl<spherical_polar_tag>
{
    template <typename Point1, typename Point2>
    static inline int apply(Point1 segment_a, Point1 segment_b,
                            Point2 p)
    {
        typedef 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<Point2>::type,
                typename cs_angular_units<Point2>::type
            > constants2;

        geometry::set<1>(segment_a,
            constants1::max_latitude() - geometry::get<1>(segment_a));
        geometry::set<1>(segment_b,
            constants1::max_latitude() - geometry::get<1>(segment_b));
        geometry::set<1>(p,
            constants2::max_latitude() - geometry::get<1>(p));

        return direction_code_impl
                <
                    spherical_equatorial_tag
                >::apply(segment_a, segment_b, p);
    }
};

// if spherical_tag is passed then pick cs_tag based on Point1 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,
                            Point2 p)
    {
        return direction_code_impl
            <
                typename boost::mpl::if_c
                    <
                        boost::is_same
                            <
                                typename geometry::cs_tag<Point1>::type,
                                spherical_polar_tag
                            >::value,
                        spherical_polar_tag,
                        spherical_equatorial_tag
                    >::type
            >::apply(segment_a, segment_b, p);
    }
};

template <>
struct direction_code_impl<geographic_tag>
    : direction_code_impl<spherical_equatorial_tag>
{};

// Gives sense of direction for point p, collinear w.r.t. segment (a,b)
// Returns -1 if p goes backward w.r.t (a,b), so goes from b in direction of a
// 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,
                          Point2 const& p)
{
    return direction_code_impl<CSTag>::apply(segment_a, segment_b, p);
}


} // namespace detail
#endif //DOXYGEN_NO_DETAIL


}} // namespace boost::geometry

#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
Commit	Line	Data
b32b8144 FG	1	// Boost.Geometry (aka GGL, Generic Geometry Library)
	2
	3	// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
	4
92f5a8d4 TL	5	// This file was modified by Oracle on 2015, 2017, 2019.
92f5a8d4 TL	6	// Modifications copyright (c) 2015-2019 Oracle and/or its affiliates.
b32b8144 FG	7
	8	// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
	9	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
	10
	11	// Use, modification and distribution is subject to the Boost Software License,
	12	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	13	// http://www.boost.org/LICENSE_1_0.txt)
	14
92f5a8d4 TL	15	#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
92f5a8d4 TL	16	#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP
b32b8144 FG	17
	18
	19	#include <boost/geometry/core/access.hpp>
92f5a8d4 TL	20	#include <boost/geometry/arithmetic/infinite_line_functions.hpp>
92f5a8d4 TL	21	#include <boost/geometry/algorithms/detail/make/make.hpp>
b32b8144 FG	22	#include <boost/geometry/util/math.hpp>
	23	#include <boost/geometry/util/select_coordinate_type.hpp>
	24	#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
	25
	26	#include <boost/mpl/assert.hpp>
	27
	28
	29	namespace boost { namespace geometry
	30	{
	31
	32
	33	#ifndef DOXYGEN_NO_DETAIL
	34	namespace detail
	35	{
	36
92f5a8d4	37	template <typename CSTag>
b32b8144 FG	38	struct direction_code_impl
	39	{
	40	BOOST_MPL_ASSERT_MSG((false), NOT_IMPLEMENTED_FOR_THIS_CS, (CSTag));
	41	};
	42
92f5a8d4 TL	43	template <>
92f5a8d4 TL	44	struct direction_code_impl<cartesian_tag>
b32b8144 FG	45	{
	46	template <typename Point1, typename Point2>
	47	static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
92f5a8d4	48	Point2 const& point)
b32b8144 FG	49	{
	50	typedef typename geometry::select_coordinate_type
	51	<
	52	Point1, Point2
	53	>::type calc_t;
	54
92f5a8d4 TL	55	typedef model::infinite_line<calc_t> line_type;
	56
	57	// point b is often equal to the specified point, check that first
	58	line_type const q = detail::make::make_infinite_line<calc_t>(segment_b, point);
	59	if (arithmetic::is_degenerate(q))
b32b8144 FG	60	{
	61	return 0;
	62	}
	63
92f5a8d4 TL	64	line_type const p = detail::make::make_infinite_line<calc_t>(segment_a, segment_b);
92f5a8d4 TL	65	if (arithmetic::is_degenerate(p))
b32b8144	66	{
92f5a8d4	67	return 0;
b32b8144 FG	68	}
b32b8144 FG	69
92f5a8d4 TL	70	// p extends a-b if direction is similar
92f5a8d4 TL	71	return arithmetic::similar_direction(p, q) ? 1 : -1;
b32b8144 FG	72	}
	73	};
	74
92f5a8d4 TL	75	template <>
92f5a8d4 TL	76	struct direction_code_impl<spherical_equatorial_tag>
b32b8144 FG	77	{
	78	template <typename Point1, typename Point2>
	79	static inline int apply(Point1 const& segment_a, Point1 const& segment_b,
	80	Point2 const& p)
	81	{
	82	typedef typename coordinate_type<Point1>::type coord1_t;
	83	typedef typename coordinate_type<Point2>::type coord2_t;
92f5a8d4 TL	84	typedef typename cs_angular_units<Point1>::type units_t;
92f5a8d4 TL	85	typedef typename cs_angular_units<Point2>::type units2_t;
b32b8144 FG	86	BOOST_MPL_ASSERT_MSG((boost::is_same<units_t, units2_t>::value),
	87	NOT_IMPLEMENTED_FOR_DIFFERENT_UNITS,
	88	(units_t, units2_t));
	89
	90	typedef typename geometry::select_coordinate_type <Point1, Point2>::type calc_t;
	91	typedef math::detail::constants_on_spheroid<coord1_t, units_t> constants1;
	92	typedef math::detail::constants_on_spheroid<coord2_t, units_t> constants2;
	93	static coord1_t const pi_half1 = constants1::max_latitude();
	94	static coord2_t const pi_half2 = constants2::max_latitude();
	95	static calc_t const c0 = 0;
	96
	97	coord1_t const a0 = geometry::get<0>(segment_a);
	98	coord1_t const a1 = geometry::get<1>(segment_a);
	99	coord1_t const b0 = geometry::get<0>(segment_b);
	100	coord1_t const b1 = geometry::get<1>(segment_b);
	101	coord2_t const p0 = geometry::get<0>(p);
	102	coord2_t const p1 = geometry::get<1>(p);
	103
	104	if ( (math::equals(b0, a0) && math::equals(b1, a1))
	105	\|\| (math::equals(b0, p0) && math::equals(b1, p1)) )
	106	{
	107	return 0;
	108	}
	109
	110	bool const is_a_pole = math::equals(pi_half1, math::abs(a1));
	111	bool const is_b_pole = math::equals(pi_half1, math::abs(b1));
	112	bool const is_p_pole = math::equals(pi_half2, math::abs(p1));
	113
	114	if ( is_b_pole && ((is_a_pole && math::sign(b1) == math::sign(a1))
	115	\|\| (is_p_pole && math::sign(b1) == math::sign(p1))) )
	116	{
	117	return 0;
	118	}
	119
	120	// NOTE: as opposed to the implementation for cartesian CS
	121	// here point b is the origin
	122
	123	calc_t const dlon1 = math::longitude_distance_signed<units_t, calc_t>(b0, a0);
	124	calc_t const dlon2 = math::longitude_distance_signed<units_t, calc_t>(b0, p0);
	125
	126	bool is_antilon1 = false, is_antilon2 = false;
	127	calc_t const dlat1 = latitude_distance_signed<units_t, calc_t>(b1, a1, dlon1, is_antilon1);
	128	calc_t const dlat2 = latitude_distance_signed<units_t, calc_t>(b1, p1, dlon2, is_antilon2);
	129
	130	calc_t mx = is_a_pole \|\| is_b_pole \|\| is_p_pole ?
	131	c0 :
	132	(std::min)(is_antilon1 ? c0 : math::abs(dlon1),
	133	is_antilon2 ? c0 : math::abs(dlon2));
	134	calc_t my = (std::min)(math::abs(dlat1),
	135	math::abs(dlat2));
	136
	137	int s1 = 0, s2 = 0;
	138	if (mx >= my)
	139	{
	140	s1 = dlon1 > 0 ? 1 : -1;
	141	s2 = dlon2 > 0 ? 1 : -1;
	142	}
	143	else
	144	{
	145	s1 = dlat1 > 0 ? 1 : -1;
	146	s2 = dlat2 > 0 ? 1 : -1;
	147	}
	148
	149	return s1 == s2 ? -1 : 1;
150	}
151
152	template <typename Units, typename T>
153	static inline T latitude_distance_signed(T const& lat1, T const& lat2, T const& lon_ds, bool & is_antilon)
154	{
155	typedef math::detail::constants_on_spheroid<T, Units> constants;
156	static T const pi = constants::half_period();
157	static T const c0 = 0;
158
159	T res = lat2 - lat1;
160
161	is_antilon = math::equals(math::abs(lon_ds), pi);
162	if (is_antilon)
163	{
164	res = lat2 + lat1;
165	if (res >= c0)
166	res = pi - res;
167	else
168	res = -pi - res;
169	}
170
171	return res;
172	}
173	};
174
92f5a8d4 TL	175	template <>
92f5a8d4 TL	176	struct direction_code_impl<spherical_polar_tag>
b32b8144 FG	177	{
	178	template <typename Point1, typename Point2>
	179	static inline int apply(Point1 segment_a, Point1 segment_b,
	180	Point2 p)
	181	{
	182	typedef math::detail::constants_on_spheroid
	183	<
	184	typename coordinate_type<Point1>::type,
92f5a8d4	185	typename cs_angular_units<Point1>::type
b32b8144 FG	186	> constants1;
	187	typedef math::detail::constants_on_spheroid
	188	<
	189	typename coordinate_type<Point2>::type,
92f5a8d4	190	typename cs_angular_units<Point2>::type
b32b8144 FG	191	> constants2;
	192
	193	geometry::set<1>(segment_a,
	194	constants1::max_latitude() - geometry::get<1>(segment_a));
	195	geometry::set<1>(segment_b,
	196	constants1::max_latitude() - geometry::get<1>(segment_b));
	197	geometry::set<1>(p,
	198	constants2::max_latitude() - geometry::get<1>(p));
	199
	200	return direction_code_impl
	201	<
92f5a8d4	202	spherical_equatorial_tag
b32b8144 FG	203	>::apply(segment_a, segment_b, p);
	204	}
	205	};
	206
92f5a8d4 TL	207	// if spherical_tag is passed then pick cs_tag based on Point1 type
	208	// with spherical_equatorial_tag as the default
	209	template <>
	210	struct direction_code_impl<spherical_tag>
	211	{
	212	template <typename Point1, typename Point2>
	213	static inline int apply(Point1 segment_a, Point1 segment_b,
	214	Point2 p)
	215	{
	216	return direction_code_impl
	217	<
	218	typename boost::mpl::if_c
	219	<
	220	boost::is_same
	221	<
	222	typename geometry::cs_tag<Point1>::type,
	223	spherical_polar_tag
	224	>::value,
	225	spherical_polar_tag,
	226	spherical_equatorial_tag
	227	>::type
	228	>::apply(segment_a, segment_b, p);
	229	}
	230	};
	231
	232	template <>
	233	struct direction_code_impl<geographic_tag>
	234	: direction_code_impl<spherical_equatorial_tag>
b32b8144 FG	235	{};
	236
	237	// Gives sense of direction for point p, collinear w.r.t. segment (a,b)
	238	// Returns -1 if p goes backward w.r.t (a,b), so goes from b in direction of a
	239	// Returns 1 if p goes forward, so extends (a,b)
	240	// Returns 0 if p is equal with b, or if (a,b) is degenerate
	241	// Note that it does not do any collinearity test, that should be done before
92f5a8d4	242	template <typename CSTag, typename Point1, typename Point2>
b32b8144 FG	243	inline int direction_code(Point1 const& segment_a, Point1 const& segment_b,
	244	Point2 const& p)
	245	{
92f5a8d4	246	return direction_code_impl<CSTag>::apply(segment_a, segment_b, p);
b32b8144 FG	247	}
	248
	249
	250	} // namespace detail
	251	#endif //DOXYGEN_NO_DETAIL
	252
	253
	254	}} // namespace boost::geometry
	255
92f5a8d4	256	#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECTION_CODE_HPP