Chris@16: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@16: 
Chris@16: // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
Chris@16: // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
Chris@16: // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
Chris@16: 
Chris@101: // This file was modified by Oracle on 2015.
Chris@101: // Modifications copyright (c) 2015 Oracle and/or its affiliates.
Chris@101: 
Chris@101: // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
Chris@101: 
Chris@16: // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
Chris@16: // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
Chris@16: 
Chris@16: // Use, modification and distribution is subject to the Boost Software License,
Chris@16: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
Chris@16: #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
Chris@16: 
Chris@16: 
Chris@101: #include <cstddef>
Chris@101: 
Chris@16: #include <boost/mpl/if.hpp>
Chris@16: #include <boost/numeric/conversion/cast.hpp>
Chris@16: #include <boost/type_traits.hpp>
Chris@16: 
Chris@16: #include <boost/geometry/arithmetic/determinant.hpp>
Chris@16: #include <boost/geometry/core/coordinate_type.hpp>
Chris@16: #include <boost/geometry/core/point_type.hpp>
Chris@16: #include <boost/geometry/strategies/centroid.hpp>
Chris@16: #include <boost/geometry/util/select_coordinate_type.hpp>
Chris@16: 
Chris@16: 
Chris@16: namespace boost { namespace geometry
Chris@16: {
Chris@16: 
Chris@16: // Note: when calling the namespace "centroid", it sometimes,
Chris@16: // somehow, in gcc, gives compilation problems (confusion with function centroid).
Chris@16: 
Chris@16: namespace strategy { namespace centroid
Chris@16: {
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: /*!
Chris@16: \brief Centroid calculation using algorithm Bashein / Detmer
Chris@16: \ingroup strategies
Chris@16: \details Calculates centroid using triangulation method published by
Chris@16:     Bashein / Detmer
Chris@16: \tparam Point point type of centroid to calculate
Chris@16: \tparam PointOfSegment point type of segments, defaults to Point
Chris@16: \tparam CalculationType \tparam_calculation
Chris@16: 
Chris@16: \author Adapted from  "Centroid of a Polygon" by
Chris@16:     Gerard Bashein and Paul R. Detmer<em>,
Chris@16: in "Graphics Gems IV", Academic Press, 1994</em>
Chris@16: 
Chris@16: 
Chris@16: \qbk{
Chris@16: [heading See also]
Chris@16: [link geometry.reference.algorithms.centroid.centroid_3_with_strategy centroid (with strategy)]
Chris@16: }
Chris@16: */
Chris@16: 
Chris@16: /*
Chris@16: \par Research notes
Chris@16: The algorithm gives the same results as Oracle and PostGIS but
Chris@16:     differs from MySQL
Chris@16: (tried 5.0.21 / 5.0.45 / 5.0.51a / 5.1.23).
Chris@16: 
Chris@16: Without holes:
Chris@16: - this:       POINT(4.06923363095238 1.65055803571429)
Chris@16: - geolib:     POINT(4.07254 1.66819)
Chris@16: - MySQL:      POINT(3.6636363636364  1.6272727272727)'
Chris@16: - PostGIS:    POINT(4.06923363095238 1.65055803571429)
Chris@16: - Oracle:           4.06923363095238 1.65055803571429
Chris@16: - SQL Server: POINT(4.06923362245959 1.65055804168294)
Chris@16: 
Chris@16: Statements:
Chris@16: - \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
Chris@16:     '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
Chris@16:         ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))')))
Chris@16: - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null,
Chris@16:         sdo_elem_info_array(1, 1003, 1), sdo_ordinate_array(
Chris@16:             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
Chris@16:             ,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3))
Chris@16:         , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
Chris@16:         ,mdsys.sdo_dim_element('y',0,10,.00000005)))
Chris@16:         from dual
Chris@16: - \b SQL Server 2008: select geometry::STGeomFromText(
Chris@16:     '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
Chris@16:         ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))',0)
Chris@16:                 .STCentroid()
Chris@16:                 .STAsText()
Chris@16: 
Chris@16: With holes:
Chris@16: - this:       POINT(4.04663 1.6349)
Chris@16: - geolib:     POINT(4.04675 1.65735)
Chris@16: - MySQL:      POINT(3.6090580503834 1.607573932092)
Chris@16: - PostGIS:    POINT(4.0466265060241 1.63489959839357)
Chris@16: - Oracle:           4.0466265060241 1.63489959839357
Chris@16: - SQL Server: POINT(4.0466264962959677 1.6348996057331333)
Chris@16: 
Chris@16: Statements:
Chris@16: - \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
Chris@16:     'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2
Chris@16:         ,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)
Chris@16:         ,(4 2,4.2 1.4,4.8 1.9,4.4 2.2,4 2))')));
Chris@16: - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null
Chris@16:         , sdo_elem_info_array(1, 1003, 1, 25, 2003, 1)
Chris@16:         , sdo_ordinate_array(2,1.3,2.4,1.7,2.8,1.8,3.4,1.2,3.7,1.6,3.4,
Chris@16:         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,
Chris@16:         4.8,1.9, 4.4,2.2, 4,2))
Chris@16:         , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
Chris@16:         ,mdsys.sdo_dim_element('y',0,10,.00000005)))
Chris@16:         from dual
Chris@16:  */
Chris@16: template
Chris@16: <
Chris@16:     typename Point,
Chris@16:     typename PointOfSegment = Point,
Chris@16:     typename CalculationType = void
Chris@16: >
Chris@16: class bashein_detmer
Chris@16: {
Chris@16: private :
Chris@16:     // If user specified a calculation type, use that type,
Chris@16:     //   whatever it is and whatever the point-type(s) are.
Chris@16:     // Else, use the most appropriate coordinate type
Chris@16:     //    of the two points, but at least double
Chris@16:     typedef typename
Chris@16:         boost::mpl::if_c
Chris@16:         <
Chris@16:             boost::is_void<CalculationType>::type::value,
Chris@16:             typename select_most_precise
Chris@16:             <
Chris@16:                 typename select_coordinate_type
Chris@16:                     <
Chris@16:                         Point,
Chris@16:                         PointOfSegment
Chris@16:                     >::type,
Chris@16:                 double
Chris@16:             >::type,
Chris@16:             CalculationType
Chris@16:         >::type calculation_type;
Chris@16: 
Chris@16:     /*! subclass to keep state */
Chris@16:     class sums
Chris@16:     {
Chris@16:         friend class bashein_detmer;
Chris@101:         std::size_t count;
Chris@16:         calculation_type sum_a2;
Chris@16:         calculation_type sum_x;
Chris@16:         calculation_type sum_y;
Chris@16: 
Chris@16:     public :
Chris@16:         inline sums()
Chris@16:             : count(0)
Chris@16:             , sum_a2(calculation_type())
Chris@16:             , sum_x(calculation_type())
Chris@16:             , sum_y(calculation_type())
Chris@16:         {}
Chris@16:     };
Chris@16: 
Chris@16: public :
Chris@16:     typedef sums state_type;
Chris@16: 
Chris@16:     static inline void apply(PointOfSegment const& p1,
Chris@16:             PointOfSegment const& p2, sums& state)
Chris@16:     {
Chris@16:         /* Algorithm:
Chris@16:         For each segment:
Chris@16:         begin
Chris@16:             ai = x1 * y2 - x2 * y1;
Chris@16:             sum_a2 += ai;
Chris@16:             sum_x += ai * (x1 + x2);
Chris@16:             sum_y += ai * (y1 + y2);
Chris@16:         end
Chris@16:         return POINT(sum_x / (3 * sum_a2), sum_y / (3 * sum_a2) )
Chris@16:         */
Chris@16: 
Chris@16:         // Get coordinates and promote them to calculation_type
Chris@16:         calculation_type const x1 = boost::numeric_cast<calculation_type>(get<0>(p1));
Chris@16:         calculation_type const y1 = boost::numeric_cast<calculation_type>(get<1>(p1));
Chris@16:         calculation_type const x2 = boost::numeric_cast<calculation_type>(get<0>(p2));
Chris@16:         calculation_type const y2 = boost::numeric_cast<calculation_type>(get<1>(p2));
Chris@16:         calculation_type const ai = geometry::detail::determinant<calculation_type>(p1, p2);
Chris@16:         state.count++;
Chris@16:         state.sum_a2 += ai;
Chris@16:         state.sum_x += ai * (x1 + x2);
Chris@16:         state.sum_y += ai * (y1 + y2);
Chris@16:     }
Chris@16: 
Chris@16:     static inline bool result(sums const& state, Point& centroid)
Chris@16:     {
Chris@16:         calculation_type const zero = calculation_type();
Chris@16:         if (state.count > 0 && ! math::equals(state.sum_a2, zero))
Chris@16:         {
Chris@16:             calculation_type const v3 = 3;
Chris@16:             calculation_type const a3 = v3 * state.sum_a2;
Chris@16: 
Chris@16:             typedef typename geometry::coordinate_type
Chris@16:                 <
Chris@16:                     Point
Chris@16:                 >::type coordinate_type;
Chris@16: 
Chris@16:             set<0>(centroid,
Chris@16:                 boost::numeric_cast<coordinate_type>(state.sum_x / a3));
Chris@16:             set<1>(centroid,
Chris@16:                 boost::numeric_cast<coordinate_type>(state.sum_y / a3));
Chris@16:             return true;
Chris@16:         }
Chris@16: 
Chris@16:         return false;
Chris@16:     }
Chris@16: 
Chris@16: };
Chris@16: 
Chris@16: #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@16: 
Chris@16: namespace services
Chris@16: {
Chris@16: 
Chris@16: // Register this strategy for rings and (multi)polygons, in two dimensions
Chris@16: template <typename Point, typename Geometry>
Chris@16: struct default_strategy<cartesian_tag, areal_tag, 2, Point, Geometry>
Chris@16: {
Chris@16:     typedef bashein_detmer
Chris@16:         <
Chris@16:             Point,
Chris@16:             typename point_type<Geometry>::type
Chris@16:         > type;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: } // namespace services
Chris@16: 
Chris@16: 
Chris@16: #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
Chris@16: 
Chris@16: 
Chris@16: }} // namespace strategy::centroid
Chris@16: 
Chris@16: 
Chris@16: }} // namespace boost::geometry
Chris@16: 
Chris@16: 
Chris@16: #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP