centroid_bashein_detmer.hpp 8.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
  3. // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
  4. // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
  5. // This file was modified by Oracle on 2015.
  6. // Modifications copyright (c) 2015 Oracle and/or its affiliates.
  7. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  8. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  9. // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  10. // Use, modification and distribution is subject to the Boost Software License,
  11. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  12. // http://www.boost.org/LICENSE_1_0.txt)
  13. #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
  14. #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
  15. #include <cstddef>
  16. #include <boost/math/special_functions/fpclassify.hpp>
  17. #include <boost/mpl/if.hpp>
  18. #include <boost/numeric/conversion/cast.hpp>
  19. #include <boost/type_traits/is_void.hpp>
  20. #include <boost/geometry/arithmetic/determinant.hpp>
  21. #include <boost/geometry/core/coordinate_type.hpp>
  22. #include <boost/geometry/core/point_type.hpp>
  23. #include <boost/geometry/strategies/centroid.hpp>
  24. #include <boost/geometry/util/select_coordinate_type.hpp>
  25. namespace boost { namespace geometry
  26. {
  27. // Note: when calling the namespace "centroid", it sometimes,
  28. // somehow, in gcc, gives compilation problems (confusion with function centroid).
  29. namespace strategy { namespace centroid
  30. {
  31. /*!
  32. \brief Centroid calculation using algorithm Bashein / Detmer
  33. \ingroup strategies
  34. \details Calculates centroid using triangulation method published by
  35. Bashein / Detmer
  36. \tparam Point point type of centroid to calculate
  37. \tparam PointOfSegment point type of segments, defaults to Point
  38. \tparam CalculationType \tparam_calculation
  39. \author Adapted from "Centroid of a Polygon" by
  40. Gerard Bashein and Paul R. Detmer<em>,
  41. in "Graphics Gems IV", Academic Press, 1994</em>
  42. \qbk{
  43. [heading See also]
  44. [link geometry.reference.algorithms.centroid.centroid_3_with_strategy centroid (with strategy)]
  45. }
  46. */
  47. /*
  48. \par Research notes
  49. The algorithm gives the same results as Oracle and PostGIS but
  50. differs from MySQL
  51. (tried 5.0.21 / 5.0.45 / 5.0.51a / 5.1.23).
  52. Without holes:
  53. - this: POINT(4.06923363095238 1.65055803571429)
  54. - geolib: POINT(4.07254 1.66819)
  55. - MySQL: POINT(3.6636363636364 1.6272727272727)'
  56. - PostGIS: POINT(4.06923363095238 1.65055803571429)
  57. - Oracle: 4.06923363095238 1.65055803571429
  58. - SQL Server: POINT(4.06923362245959 1.65055804168294)
  59. Statements:
  60. - \b MySQL/PostGIS: select AsText(Centroid(GeomFromText(
  61. '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
  62. ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))')))
  63. - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null,
  64. sdo_elem_info_array(1, 1003, 1), sdo_ordinate_array(
  65. 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
  66. ,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3))
  67. , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
  68. ,mdsys.sdo_dim_element('y',0,10,.00000005)))
  69. from dual
  70. - \b SQL Server 2008: select geometry::STGeomFromText(
  71. '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
  72. ,5.4 1.2,4.9 0.8,2.9 0.7,2 1.3))',0)
  73. .STCentroid()
  74. .STAsText()
  75. With holes:
  76. - this: POINT(4.04663 1.6349)
  77. - geolib: POINT(4.04675 1.65735)
  78. - MySQL: POINT(3.6090580503834 1.607573932092)
  79. - PostGIS: POINT(4.0466265060241 1.63489959839357)
  80. - Oracle: 4.0466265060241 1.63489959839357
  81. - SQL Server: POINT(4.0466264962959677 1.6348996057331333)
  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
  85. ,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)
  86. ,(4 2,4.2 1.4,4.8 1.9,4.4 2.2,4 2))')));
  87. - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null
  88. , sdo_elem_info_array(1, 1003, 1, 25, 2003, 1)
  89. , sdo_ordinate_array(2,1.3,2.4,1.7,2.8,1.8,3.4,1.2,3.7,1.6,3.4,
  90. 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,
  91. 4.8,1.9, 4.4,2.2, 4,2))
  92. , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
  93. ,mdsys.sdo_dim_element('y',0,10,.00000005)))
  94. from dual
  95. */
  96. template
  97. <
  98. typename Point,
  99. typename PointOfSegment = Point,
  100. typename CalculationType = void
  101. >
  102. class bashein_detmer
  103. {
  104. private :
  105. // If user specified a calculation type, use that type,
  106. // whatever it is and whatever the point-type(s) are.
  107. // Else, use the most appropriate coordinate type
  108. // of the two points, but at least double
  109. typedef typename
  110. boost::mpl::if_c
  111. <
  112. boost::is_void<CalculationType>::type::value,
  113. typename select_most_precise
  114. <
  115. typename select_coordinate_type
  116. <
  117. Point,
  118. PointOfSegment
  119. >::type,
  120. double
  121. >::type,
  122. CalculationType
  123. >::type calculation_type;
  124. /*! subclass to keep state */
  125. class sums
  126. {
  127. friend class bashein_detmer;
  128. std::size_t count;
  129. calculation_type sum_a2;
  130. calculation_type sum_x;
  131. calculation_type sum_y;
  132. public :
  133. inline sums()
  134. : count(0)
  135. , sum_a2(calculation_type())
  136. , sum_x(calculation_type())
  137. , sum_y(calculation_type())
  138. {}
  139. };
  140. public :
  141. typedef sums state_type;
  142. static inline void apply(PointOfSegment const& p1,
  143. PointOfSegment const& p2, sums& state)
  144. {
  145. /* Algorithm:
  146. For each segment:
  147. begin
  148. ai = x1 * y2 - x2 * y1;
  149. sum_a2 += ai;
  150. sum_x += ai * (x1 + x2);
  151. sum_y += ai * (y1 + y2);
  152. end
  153. return POINT(sum_x / (3 * sum_a2), sum_y / (3 * sum_a2) )
  154. */
  155. // Get coordinates and promote them to calculation_type
  156. calculation_type const x1 = boost::numeric_cast<calculation_type>(get<0>(p1));
  157. calculation_type const y1 = boost::numeric_cast<calculation_type>(get<1>(p1));
  158. calculation_type const x2 = boost::numeric_cast<calculation_type>(get<0>(p2));
  159. calculation_type const y2 = boost::numeric_cast<calculation_type>(get<1>(p2));
  160. calculation_type const ai = geometry::detail::determinant<calculation_type>(p1, p2);
  161. state.count++;
  162. state.sum_a2 += ai;
  163. state.sum_x += ai * (x1 + x2);
  164. state.sum_y += ai * (y1 + y2);
  165. }
  166. static inline bool result(sums const& state, Point& centroid)
  167. {
  168. calculation_type const zero = calculation_type();
  169. if (state.count > 0 && ! math::equals(state.sum_a2, zero))
  170. {
  171. calculation_type const v3 = 3;
  172. calculation_type const a3 = v3 * state.sum_a2;
  173. typedef typename geometry::coordinate_type
  174. <
  175. Point
  176. >::type coordinate_type;
  177. // Prevent NaN centroid coordinates
  178. if (boost::math::isfinite(a3))
  179. {
  180. // NOTE: above calculation_type is checked, not the centroid coordinate_type
  181. // which means that the centroid can still be filled with INF
  182. // if e.g. calculation_type is double and centroid contains floats
  183. set<0>(centroid,
  184. boost::numeric_cast<coordinate_type>(state.sum_x / a3));
  185. set<1>(centroid,
  186. boost::numeric_cast<coordinate_type>(state.sum_y / a3));
  187. return true;
  188. }
  189. }
  190. return false;
  191. }
  192. };
  193. #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  194. namespace services
  195. {
  196. // Register this strategy for rings and (multi)polygons, in two dimensions
  197. template <typename Point, typename Geometry>
  198. struct default_strategy<cartesian_tag, areal_tag, 2, Point, Geometry>
  199. {
  200. typedef bashein_detmer
  201. <
  202. Point,
  203. typename point_type<Geometry>::type
  204. > type;
  205. };
  206. } // namespace services
  207. #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  208. }} // namespace strategy::centroid
  209. }} // namespace boost::geometry
  210. #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP