strategies/cartesian/centroid_bashein_detmer.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004 // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
0005 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
0006
0007 // This file was modified by Oracle on 2015-2021.
0008 // Modifications copyright (c) 2015-2021, Oracle and/or its affiliates.
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010
0011 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
0012 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
0013
0014 // Use, modification and distribution is subject to the Boost Software License,
0015 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0016 // http://www.boost.org/LICENSE_1_0.txt)
0017
0018 #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
0019 #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
0020
0021
0022 #include <cstddef>
0023
0024 #include <boost/math/special_functions/fpclassify.hpp>
0025 #include <boost/numeric/conversion/cast.hpp>
0026
0027 #include <boost/geometry/arithmetic/determinant.hpp>
0028 #include <boost/geometry/core/coordinate_type.hpp>
0029 #include <boost/geometry/core/point_type.hpp>
0030 #include <boost/geometry/strategies/centroid.hpp>
0031 #include <boost/geometry/util/math.hpp>
0032 #include <boost/geometry/util/select_most_precise.hpp>
0033
0034
0035 namespace boost { namespace geometry
0036 {
0037
0038 // Note: when calling the namespace "centroid", it sometimes,
0039 // somehow, in gcc, gives compilation problems (confusion with function centroid).
0040
0041 namespace strategy { namespace centroid
0042 {
0043
0044
0045
0046 /*!
0047 \brief Centroid calculation using algorithm Bashein / Detmer
0048 \ingroup strategies
0049 \details Calculates centroid using triangulation method published by
0050     Bashein / Detmer
0051 \tparam Point point type of centroid to calculate
0052 \tparam PointOfSegment point type of segments, defaults to Point
0053 \tparam CalculationType \tparam_calculation
0054
0055 \author Adapted from  "Centroid of a Polygon" by
0056     Gerard Bashein and Paul R. Detmer<em>,
0057 in "Graphics Gems IV", Academic Press, 1994</em>
0058
0059
0060 \qbk{
0061 [heading See also]
0062 [link geometry.reference.algorithms.centroid.centroid_3_with_strategy centroid (with strategy)]
0063 }
0064 */
0065
0066 /*
0067 \par Research notes
0068 The algorithm gives the same results as Oracle and PostGIS but
0069     differs from MySQL
0070 (tried 5.0.21 / 5.0.45 / 5.0.51a / 5.1.23).
0071
0072 Without holes:
0073 - this:       POINT(4.06923363095238 1.65055803571429)
0074 - geolib:     POINT(4.07254 1.66819)
0075 - MySQL:      POINT(3.6636363636364  1.6272727272727)'
0076 - PostGIS:    POINT(4.06923363095238 1.65055803571429)
0077 - Oracle:           4.06923363095238 1.65055803571429
0078 - SQL Server: POINT(4.06923362245959 1.65055804168294)
0079
0080 Statements:
0081 - \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
0082     '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
0083         ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))')))
0084 - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null,
0085         sdo_elem_info_array(1, 1003, 1), sdo_ordinate_array(
0086             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
0087             ,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3))
0088         , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
0089         ,mdsys.sdo_dim_element('y',0,10,.00000005)))
0090         from dual
0091 - \b SQL Server 2008: select geometry::STGeomFromText(
0092     '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
0093         ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))',0)
0094                 .STCentroid()
0095                 .STAsText()
0096
0097 With holes:
0098 - this:       POINT(4.04663 1.6349)
0099 - geolib:     POINT(4.04675 1.65735)
0100 - MySQL:      POINT(3.6090580503834 1.607573932092)
0101 - PostGIS:    POINT(4.0466265060241 1.63489959839357)
0102 - Oracle:           4.0466265060241 1.63489959839357
0103 - SQL Server: POINT(4.0466264962959677 1.6348996057331333)
0104
0105 Statements:
0106 - \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
0107     'POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2
0108         ,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)
0109         ,(4 2,4.2 1.4,4.8 1.9,4.4 2.2,4 2))')));
0110 - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null
0111         , sdo_elem_info_array(1, 1003, 1, 25, 2003, 1)
0112         , sdo_ordinate_array(2,1.3,2.4,1.7,2.8,1.8,3.4,1.2,3.7,1.6,3.4,
0113         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,
0114         4.8,1.9, 4.4,2.2, 4,2))
0115         , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
0116         ,mdsys.sdo_dim_element('y',0,10,.00000005)))
0117         from dual
0118  */
0119 template
0120 <
0121     typename Ignored1 = void,
0122     typename Ignored2 = void,
0123     typename CalculationType = void
0124 >
0125 class bashein_detmer
0126 {
0127 private :
0128     // If user specified a calculation type, use that type,
0129     //   whatever it is and whatever the point-type(s) are.
0130     // Else, use the most appropriate coordinate type
0131     //    of the two points, but at least double
0132     template <typename GeometryPoint, typename ResultPoint>
0133     struct calculation_type
0134         : std::conditional
0135             <
0136                 std::is_void<CalculationType>::value,
0137                 typename select_most_precise
0138                     <
0139                         typename coordinate_type<GeometryPoint>::type,
0140                         typename coordinate_type<ResultPoint>::type,
0141                         double
0142                     >::type,
0143                 CalculationType
0144             >
0145     {};
0146
0147     /*! subclass to keep state */
0148     template <typename GeometryPoint, typename ResultPoint>
0149     class sums
0150     {
0151         typedef typename calculation_type<GeometryPoint, ResultPoint>::type calc_type;
0152
0153         friend class bashein_detmer;
0154         std::size_t count;
0155         calc_type sum_a2;
0156         calc_type sum_x;
0157         calc_type sum_y;
0158
0159     public :
0160         inline sums()
0161             : count(0)
0162             , sum_a2(calc_type())
0163             , sum_x(calc_type())
0164             , sum_y(calc_type())
0165         {}
0166     };
0167
0168 public :
0169     template <typename GeometryPoint, typename ResultPoint>
0170     struct state_type
0171     {
0172         typedef sums<GeometryPoint, ResultPoint> type;
0173     };
0174
0175     template <typename GeometryPoint, typename ResultPoint>
0176     static inline void apply(GeometryPoint const& p1, GeometryPoint const& p2,
0177                              sums<GeometryPoint, ResultPoint>& state)
0178     {
0179         /* Algorithm:
0180         For each segment:
0181         begin
0182             ai = x1 * y2 - x2 * y1;
0183             sum_a2 += ai;
0184             sum_x += ai * (x1 + x2);
0185             sum_y += ai * (y1 + y2);
0186         end
0187         return POINT(sum_x / (3 * sum_a2), sum_y / (3 * sum_a2) )
0188         */
0189
0190         typedef typename calculation_type<GeometryPoint, ResultPoint>::type calc_type;
0191
0192         // Get coordinates and promote them to calculation_type
0193         calc_type const x1 = boost::numeric_cast<calc_type>(get<0>(p1));
0194         calc_type const y1 = boost::numeric_cast<calc_type>(get<1>(p1));
0195         calc_type const x2 = boost::numeric_cast<calc_type>(get<0>(p2));
0196         calc_type const y2 = boost::numeric_cast<calc_type>(get<1>(p2));
0197         calc_type const ai = geometry::detail::determinant<calc_type>(p1, p2);
0198         state.count++;
0199         state.sum_a2 += ai;
0200         state.sum_x += ai * (x1 + x2);
0201         state.sum_y += ai * (y1 + y2);
0202     }
0203
0204     template <typename GeometryPoint, typename ResultPoint>
0205     static inline bool result(sums<GeometryPoint, ResultPoint> const& state,
0206                               ResultPoint& centroid)
0207     {
0208         typedef typename calculation_type<GeometryPoint, ResultPoint>::type calc_type;
0209
0210         calc_type const zero = calc_type();
0211         if (state.count > 0 && ! math::equals(state.sum_a2, zero))
0212         {
0213             calc_type const v3 = 3;
0214             calc_type const a3 = v3 * state.sum_a2;
0215
0216             typedef typename geometry::coordinate_type
0217                 <
0218                     ResultPoint
0219                 >::type coordinate_type;
0220
0221             // Prevent NaN centroid coordinates
0222             if (boost::math::isfinite(a3))
0223             {
0224                 // NOTE: above calculation_type is checked, not the centroid coordinate_type
0225                 // which means that the centroid can still be filled with INF
0226                 // if e.g. calculation_type is double and centroid contains floats
0227                 set<0>(centroid,
0228                     boost::numeric_cast<coordinate_type>(state.sum_x / a3));
0229                 set<1>(centroid,
0230                     boost::numeric_cast<coordinate_type>(state.sum_y / a3));
0231                 return true;
0232             }
0233         }
0234
0235         return false;
0236     }
0237
0238 };
0239
0240 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0241
0242 namespace services
0243 {
0244
0245 // Register this strategy for rings and (multi)polygons, in two dimensions
0246 template <typename Point, typename Geometry>
0247 struct default_strategy<cartesian_tag, areal_tag, 2, Point, Geometry>
0248 {
0249     typedef bashein_detmer
0250         <
0251             Point,
0252             typename point_type<Geometry>::type
0253         > type;
0254 };
0255
0256
0257 } // namespace services
0258
0259
0260 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0261
0262
0263 }} // namespace strategy::centroid
0264
0265
0266 }} // namespace boost::geometry
0267
0268
0269 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP