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.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 : sums()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 apply(PointOfSegment const & p1,PointOfSegment const & p2,sums & state)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 result(sums const & state,Point & centroid)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 259