ceph/src/boost/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp

   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
   7 // This file was modified by Oracle on 2015, 2018.
   8 // Modifications copyright (c) 2015, 2018, Oracle and/or its affiliates.
   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>
  34 #include <boost/geometry/util/math.hpp>
  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