[ceph.git] / ceph / src / boost / boost / geometry / strategies / cartesian / centroid_bashein_detmer.hpp

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

// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.

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

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

// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.

// 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_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP


#include <cstddef>

#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/mpl/if.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <boost/type_traits/is_void.hpp>

#include <boost/geometry/arithmetic/determinant.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/point_type.hpp>
#include <boost/geometry/strategies/centroid.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/select_coordinate_type.hpp>


namespace boost { namespace geometry
{

// Note: when calling the namespace "centroid", it sometimes,
// somehow, in gcc, gives compilation problems (confusion with function centroid).

namespace strategy { namespace centroid
{


/*!
\brief Centroid calculation using algorithm Bashein / Detmer
\ingroup strategies
\details Calculates centroid using triangulation method published by
    Bashein / Detmer
\tparam Point point type of centroid to calculate
\tparam PointOfSegment point type of segments, defaults to Point
\tparam CalculationType \tparam_calculation

\author Adapted from  "Centroid of a Polygon" by
    Gerard Bashein and Paul R. Detmer<em>,
in "Graphics Gems IV", Academic Press, 1994</em>


\qbk{
[heading See also]
[link geometry.reference.algorithms.centroid.centroid_3_with_strategy centroid (with strategy)]
}
*/

/*
\par Research notes
The algorithm gives the same results as Oracle and PostGIS but
    differs from MySQL
(tried 5.0.21 / 5.0.45 / 5.0.51a / 5.1.23).

Without holes:
- this:       POINT(4.06923363095238 1.65055803571429)
- geolib:     POINT(4.07254 1.66819)
- MySQL:      POINT(3.6636363636364  1.6272727272727)'
- PostGIS:    POINT(4.06923363095238 1.65055803571429)
- Oracle:           4.06923363095238 1.65055803571429
- SQL Server: POINT(4.06923362245959 1.65055804168294)

Statements:
- \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
    'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2,3.7 1.6,3.4 2,4.1 3,5.3 2.6
        ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))')))
- \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null,
        sdo_elem_info_array(1, 1003, 1), sdo_ordinate_array(
            2,1.3,2.4,1.7,2.8,1.8,3.4,1.2,3.7,1.6,3.4,2,4.1,3,5.3,2.6
            ,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3))
        , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
        ,mdsys.sdo_dim_element('y',0,10,.00000005)))
        from dual
- \b SQL Server 2008: select geometry::STGeomFromText(
    'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2,3.7 1.6,3.4 2,4.1 3,5.3 2.6
        ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))',0)
                .STCentroid()
                .STAsText()

With holes:
- this:       POINT(4.04663 1.6349)
- geolib:     POINT(4.04675 1.65735)
- MySQL:      POINT(3.6090580503834 1.607573932092)
- PostGIS:    POINT(4.0466265060241 1.63489959839357)
- Oracle:           4.0466265060241 1.63489959839357
- SQL Server: POINT(4.0466264962959677 1.6348996057331333)

Statements:
- \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
    'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2
        ,3.7 1.6,3.4 2,4.1 3,5.3 2.6,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3)
        ,(4 2,4.2 1.4,4.8 1.9,4.4 2.2,4 2))')));
- \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null
        , sdo_elem_info_array(1, 1003, 1, 25, 2003, 1)
        , sdo_ordinate_array(2,1.3,2.4,1.7,2.8,1.8,3.4,1.2,3.7,1.6,3.4,
        2,4.1,3,5.3,2.6,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3,4,2, 4.2,1.4,
        4.8,1.9, 4.4,2.2, 4,2))
        , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
        ,mdsys.sdo_dim_element('y',0,10,.00000005)))
        from dual
 */
template
<
    typename Point,
    typename PointOfSegment = Point,
    typename CalculationType = void
>
class bashein_detmer
{
private :
    // If user specified a calculation type, use that type,
    //   whatever it is and whatever the point-type(s) are.
    // Else, use the most appropriate coordinate type
    //    of the two points, but at least double
    typedef typename
        boost::mpl::if_c
        <
            boost::is_void<CalculationType>::type::value,
            typename select_most_precise
            <
                typename select_coordinate_type
                    <
                        Point,
                        PointOfSegment
                    >::type,
                double
            >::type,
            CalculationType
        >::type calculation_type;

    /*! subclass to keep state */
    class sums
    {
        friend class bashein_detmer;
        std::size_t count;
        calculation_type sum_a2;
        calculation_type sum_x;
        calculation_type sum_y;

    public :
        inline sums()
            : count(0)
            , sum_a2(calculation_type())
            , sum_x(calculation_type())
            , sum_y(calculation_type())
        {}
    };

public :
    typedef sums state_type;

    static inline void apply(PointOfSegment const& p1,
            PointOfSegment const& p2, sums& state)
    {
        /* Algorithm:
        For each segment:
        begin
            ai = x1 * y2 - x2 * y1;
            sum_a2 += ai;
            sum_x += ai * (x1 + x2);
            sum_y += ai * (y1 + y2);
        end
        return POINT(sum_x / (3 * sum_a2), sum_y / (3 * sum_a2) )
        */

        // Get coordinates and promote them to calculation_type
        calculation_type const x1 = boost::numeric_cast<calculation_type>(get<0>(p1));
        calculation_type const y1 = boost::numeric_cast<calculation_type>(get<1>(p1));
        calculation_type const x2 = boost::numeric_cast<calculation_type>(get<0>(p2));
        calculation_type const y2 = boost::numeric_cast<calculation_type>(get<1>(p2));
        calculation_type const ai = geometry::detail::determinant<calculation_type>(p1, p2);
        state.count++;
        state.sum_a2 += ai;
        state.sum_x += ai * (x1 + x2);
        state.sum_y += ai * (y1 + y2);
    }

    static inline bool result(sums const& state, Point& centroid)
    {
        calculation_type const zero = calculation_type();
        if (state.count > 0 && ! math::equals(state.sum_a2, zero))
        {
            calculation_type const v3 = 3;
            calculation_type const a3 = v3 * state.sum_a2;

            typedef typename geometry::coordinate_type
                <
                    Point
                >::type coordinate_type;

            // Prevent NaN centroid coordinates
            if (boost::math::isfinite(a3))
            {
                // NOTE: above calculation_type is checked, not the centroid coordinate_type
                // which means that the centroid can still be filled with INF
                // if e.g. calculation_type is double and centroid contains floats
                set<0>(centroid,
                    boost::numeric_cast<coordinate_type>(state.sum_x / a3));
                set<1>(centroid,
                    boost::numeric_cast<coordinate_type>(state.sum_y / a3));
                return true;
            }
        }

        return false;
    }

};

#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS

namespace services
{

// Register this strategy for rings and (multi)polygons, in two dimensions
template <typename Point, typename Geometry>
struct default_strategy<cartesian_tag, areal_tag, 2, Point, Geometry>
{
    typedef bashein_detmer
        <
            Point,
            typename point_type<Geometry>::type
        > type;
};


} // namespace services


#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS


}} // namespace strategy::centroid


}} // namespace boost::geometry


#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
Commit	Line	Data
7c673cae FG	1	// Boost.Geometry (aka GGL, Generic Geometry Library)
	2
	3	// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
	4	// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
	5	// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
	6
92f5a8d4 TL	7	// This file was modified by Oracle on 2015, 2018.
92f5a8d4 TL	8	// Modifications copyright (c) 2015, 2018, Oracle and/or its affiliates.
7c673cae FG	9
	10	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
	11
	12	// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
	13	// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
	14
	15	// Use, modification and distribution is subject to the Boost Software License,
	16	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	17	// http://www.boost.org/LICENSE_1_0.txt)
	18
	19	#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
	20	#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
	21
	22
	23	#include <cstddef>
	24
	25	#include <boost/math/special_functions/fpclassify.hpp>
	26	#include <boost/mpl/if.hpp>
	27	#include <boost/numeric/conversion/cast.hpp>
	28	#include <boost/type_traits/is_void.hpp>
	29
	30	#include <boost/geometry/arithmetic/determinant.hpp>
	31	#include <boost/geometry/core/coordinate_type.hpp>
	32	#include <boost/geometry/core/point_type.hpp>
	33	#include <boost/geometry/strategies/centroid.hpp>
92f5a8d4	34	#include <boost/geometry/util/math.hpp>
7c673cae FG	35	#include <boost/geometry/util/select_coordinate_type.hpp>
	36
	37
	38	namespace boost { namespace geometry
	39	{
	40
	41	// Note: when calling the namespace "centroid", it sometimes,
	42	// somehow, in gcc, gives compilation problems (confusion with function centroid).
	43
	44	namespace strategy { namespace centroid
	45	{
	46
	47
	48
	49	/*!
	50	\brief Centroid calculation using algorithm Bashein / Detmer
	51	\ingroup strategies
	52	\details Calculates centroid using triangulation method published by
	53	Bashein / Detmer
	54	\tparam Point point type of centroid to calculate
	55	\tparam PointOfSegment point type of segments, defaults to Point
	56	\tparam CalculationType \tparam_calculation
	57
	58	\author Adapted from "Centroid of a Polygon" by
	59	Gerard Bashein and Paul R. Detmer<em>,
	60	in "Graphics Gems IV", Academic Press, 1994</em>
	61
	62
	63	\qbk{
	64	[heading See also]
	65	[link geometry.reference.algorithms.centroid.centroid_3_with_strategy centroid (with strategy)]
	66	}
	67	*/
	68
	69	/*
	70	\par Research notes
	71	The algorithm gives the same results as Oracle and PostGIS but
	72	differs from MySQL
	73	(tried 5.0.21 / 5.0.45 / 5.0.51a / 5.1.23).
	74
	75	Without holes:
	76	- this: POINT(4.06923363095238 1.65055803571429)
	77	- geolib: POINT(4.07254 1.66819)
	78	- MySQL: POINT(3.6636363636364 1.6272727272727)'
	79	- PostGIS: POINT(4.06923363095238 1.65055803571429)
	80	- Oracle: 4.06923363095238 1.65055803571429
	81	- SQL Server: POINT(4.06923362245959 1.65055804168294)
	82
	83	Statements:
	84	- \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
	85	'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2,3.7 1.6,3.4 2,4.1 3,5.3 2.6
	86	,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))')))
	87	- \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null,
	88	sdo_elem_info_array(1, 1003, 1), sdo_ordinate_array(
	89	2,1.3,2.4,1.7,2.8,1.8,3.4,1.2,3.7,1.6,3.4,2,4.1,3,5.3,2.6
	90	,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3))
	91	, mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
	92	,mdsys.sdo_dim_element('y',0,10,.00000005)))
	93	from dual
	94	- \b SQL Server 2008: select geometry::STGeomFromText(
	95	'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2,3.7 1.6,3.4 2,4.1 3,5.3 2.6
	96	,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))',0)
	97	.STCentroid()
	98	.STAsText()
99
100	With holes:
101	- this: POINT(4.04663 1.6349)
102	- geolib: POINT(4.04675 1.65735)
103	- MySQL: POINT(3.6090580503834 1.607573932092)
104	- PostGIS: POINT(4.0466265060241 1.63489959839357)
105	- Oracle: 4.0466265060241 1.63489959839357
106	- SQL Server: POINT(4.0466264962959677 1.6348996057331333)
107
108	Statements:
109	- \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
110	'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2
111	,3.7 1.6,3.4 2,4.1 3,5.3 2.6,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3)
112	,(4 2,4.2 1.4,4.8 1.9,4.4 2.2,4 2))')));
113	- \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null
114	, sdo_elem_info_array(1, 1003, 1, 25, 2003, 1)
115	, sdo_ordinate_array(2,1.3,2.4,1.7,2.8,1.8,3.4,1.2,3.7,1.6,3.4,
116	2,4.1,3,5.3,2.6,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3,4,2, 4.2,1.4,
117	4.8,1.9, 4.4,2.2, 4,2))
118	, mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
119	,mdsys.sdo_dim_element('y',0,10,.00000005)))
120	from dual
121	*/
122	template
123	<
124	typename Point,
125	typename PointOfSegment = Point,
126	typename CalculationType = void
127	>
128	class bashein_detmer
129	{
130	private :
131	// If user specified a calculation type, use that type,
132	// whatever it is and whatever the point-type(s) are.
133	// Else, use the most appropriate coordinate type
134	// of the two points, but at least double
135	typedef typename
136	boost::mpl::if_c
137	<
138	boost::is_void<CalculationType>::type::value,
139	typename select_most_precise
140	<
141	typename select_coordinate_type
142	<
143	Point,
144	PointOfSegment
145	>::type,
146	double
147	>::type,
148	CalculationType
149	>::type calculation_type;
150
151	/! subclass to keep state /
152	class sums
153	{
154	friend class bashein_detmer;
155	std::size_t count;
156	calculation_type sum_a2;
157	calculation_type sum_x;
158	calculation_type sum_y;
159
160	public :
161	inline sums()
162	: count(0)
163	, sum_a2(calculation_type())
164	, sum_x(calculation_type())
165	, sum_y(calculation_type())
166	{}
167	};
168
169	public :
170	typedef sums state_type;
171
172	static inline void apply(PointOfSegment const& p1,
173	PointOfSegment const& p2, sums& state)
174	{
175	/* Algorithm:
176	For each segment:
177	begin
178	ai = x1 * y2 - x2 * y1;
179	sum_a2 += ai;
180	sum_x += ai * (x1 + x2);
181	sum_y += ai * (y1 + y2);
182	end
183	return POINT(sum_x / (3 * sum_a2), sum_y / (3 * sum_a2) )
184	*/
185
186	// Get coordinates and promote them to calculation_type
187	calculation_type const x1 = boost::numeric_cast<calculation_type>(get<0>(p1));
188	calculation_type const y1 = boost::numeric_cast<calculation_type>(get<1>(p1));
189	calculation_type const x2 = boost::numeric_cast<calculation_type>(get<0>(p2));
190	calculation_type const y2 = boost::numeric_cast<calculation_type>(get<1>(p2));
191	calculation_type const ai = geometry::detail::determinant<calculation_type>(p1, p2);
192	state.count++;
193	state.sum_a2 += ai;
194	state.sum_x += ai * (x1 + x2);
195	state.sum_y += ai * (y1 + y2);
196	}
197
198	static inline bool result(sums const& state, Point& centroid)
199	{
200	calculation_type const zero = calculation_type();
201	if (state.count > 0 && ! math::equals(state.sum_a2, zero))
202	{
203	calculation_type const v3 = 3;
204	calculation_type const a3 = v3 * state.sum_a2;
205
206	typedef typename geometry::coordinate_type
207	<
208	Point
209	>::type coordinate_type;
210
211	// Prevent NaN centroid coordinates
212	if (boost::math::isfinite(a3))
213	{
214	// NOTE: above calculation_type is checked, not the centroid coordinate_type
215	// which means that the centroid can still be filled with INF
216	// if e.g. calculation_type is double and centroid contains floats
217	set<0>(centroid,
218	boost::numeric_cast<coordinate_type>(state.sum_x / a3));
219	set<1>(centroid,
220	boost::numeric_cast<coordinate_type>(state.sum_y / a3));
221	return true;
222	}
223	}
224
225	return false;
226	}
227
228	};
229
230	#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
231
232	namespace services
233	{
234
235	// Register this strategy for rings and (multi)polygons, in two dimensions
236	template <typename Point, typename Geometry>
237	struct default_strategy<cartesian_tag, areal_tag, 2, Point, Geometry>
238	{
239	typedef bashein_detmer
240	<
241	Point,
242	typename point_type<Geometry>::type
243	> type;
244	};
245
246
247	} // namespace services
248
249
250	#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
251
252
253	}} // namespace strategy::centroid
254
255
256	}} // namespace boost::geometry
257
258
259	#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP