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