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