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