← Back to branch summary

~ubuntu-branches/ubuntu/saucy/merkaartor/saucy

« back to all changes in this revision

Viewing changes to include/builtin-ggl/ggl/strategies/cartesian/cart_centroid.hpp

Tags: upstream-0.15.3+svn20934
Import upstream version 0.15.3+svn20934

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/builtin-ggl/ggl/strategies/cartesian/cart_centroid.hpp
 
1
// Generic Geometry Library
 
2
//
 
3
// Copyright Barend Gehrels 1995-2009, Geodan Holding B.V. Amsterdam, the Netherlands.
 
4
// Copyright Bruno Lalande 2008, 2009
 
5
// Use, modification and distribution is subject to the Boost Software License,
 
6
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
 
7
// http://www.boost.org/LICENSE_1_0.txt)
 
8
 
 
9
#ifndef GGL_STRATEGY_CARTESIAN_CENTROID_HPP
 
10
#define GGL_STRATEGY_CARTESIAN_CENTROID_HPP
 
11
 
 
12
 
 
13
#include <boost/mpl/if.hpp>
 
14
#include <boost/numeric/conversion/cast.hpp>
 
15
#include <boost/type_traits.hpp>
 
16
 
 
17
#include <ggl/geometries/point_xy.hpp>
 
18
#include <ggl/geometries/segment.hpp>
 
19
 
 
20
#include <ggl/algorithms/assign.hpp>
 
21
 
 
22
#include <ggl/util/select_coordinate_type.hpp>
 
23
#include <ggl/util/copy.hpp>
 
24
 
 
25
 
 
26
namespace ggl
 
27
{
 
28
 
 
29
namespace strategy
 
30
{
 
31
    namespace centroid
 
32
    {
 
33
        /*!
 
34
            \brief Centroid calculation
 
35
            \details Geolib original version,
 
36
            \par Template parameters and concepts: see bashein_detmer
 
37
            \author Barend and Maarten, 1995/1996
 
38
            \author Revised for templatized library, Barend Gehrels, 2007
 
39
            \note The results are slightly different from Bashein/Detmer, so probably slightly wrong.
 
40
        */
 
41
 
 
42
        template<typename PC, typename PS = PC>
 
43
        class geolib1995
 
44
        {
 
45
            private  :
 
46
                typedef typename select_most_precise
 
47
                    <
 
48
                        typename select_coordinate_type<PC, PS>::type,
 
49
                        double
 
50
                    >::type T;
 
51
 
 
52
                /*! subclass to keep state */
 
53
                struct sums
 
54
                {
 
55
                     PC sum_ms, sum_m;
 
56
 
 
57
                    sums()
 
58
                    {
 
59
                        assign_point(sum_m, T());
 
60
                        assign_point(sum_ms, T());
 
61
                    }
 
62
                    void centroid(PC& point)
 
63
                    {
 
64
                        point = sum_ms;
 
65
 
 
66
                        if (get<0>(sum_m) != 0 && get<1>(sum_m) != 0)
 
67
                        {
 
68
                            multiply_value(point, 0.5);
 
69
                            divide_point(point, sum_m);
 
70
                        }
 
71
                        else
 
72
                        {
 
73
                            // exception?
 
74
                        }
 
75
                    }
 
76
                };
 
77
 
 
78
            public :
 
79
                typedef sums state_type;
 
80
                inline bool operator()(const segment<const PS>& s, state_type& state) const
 
81
                {
 
82
                    /* Algorithm:
 
83
                    For each segment:
 
84
                    begin
 
85
                        dx = x2 - x1;
 
86
                        dy = y2 - y1;
 
87
                        sx = x2 + x1;
 
88
                        sy = y2 + y1;
 
89
                        mx = dx * sy;
 
90
                        my = sx * dy;
 
91
 
 
92
                        sum_mx += mx;
 
93
                        sum_my += my;
 
94
                        sum_msx += mx * sx;
 
95
                        sum_msy += my * sy;
 
96
                    end;
 
97
                    return POINT(0.5 * sum_msx / sum_mx, 0.5 * sum_msy / sum_my);
 
98
                    */
 
99
 
 
100
                    PS diff = s.second, sum = s.second;
 
101
                    subtract_point(diff, s.first);
 
102
                    add_point(sum, s.first);
 
103
 
 
104
                    // We might create an arithmatic operation for this.
 
105
                    PS m;
 
106
                    get<0>(m) = get<0>(diff) * get<1>(sum);
 
107
                    get<1>(m) = get<0>(sum) * get<1>(diff);
 
108
 
 
109
                    add_point(state.sum_m, m);
 
110
                    multiply_point(m, sum);
 
111
                    add_point(state.sum_ms, m);
 
112
 
 
113
                    return true;
 
114
                }
 
115
 
 
116
        };
 
117
 
 
118
 
 
119
        /*!
 
120
            \brief Centroid calculation using algorith Bashein / Detmer
 
121
            \details Calculates centroid using triangulation method published by Bashein / Detmer
 
122
            \tparam PC point type of centroid to calculate
 
123
            \tparam PS point type of segments, defaults to PC
 
124
            \par Concepts for PC and PS:
 
125
            - specialized point_traits class
 
126
            \author Adapted from  "Centroid of a Polygon" by Gerard Bashein and Paul R. Detmer<em>,
 
127
            in "Graphics Gems IV", Academic Press, 1994</em>
 
128
            \par Research notes
 
129
            The algorithm gives the same results as Oracle and PostGIS but differs from MySQL
 
130
            (tried 5.0.21 / 5.0.45 / 5.0.51a / 5.1.23).
 
131
 
 
132
            Without holes:
 
133
            - this:       POINT(4.06923363095238 1.65055803571429)
 
134
            - geolib:     POINT(4.07254 1.66819)
 
135
            - MySQL:      POINT(3.6636363636364  1.6272727272727)'
 
136
            - PostGIS:    POINT(4.06923363095238 1.65055803571429)
 
137
            - Oracle:           4.06923363095238 1.65055803571429
 
138
            - SQL Server: POINT(4.06923362245959 1.65055804168294)
 
139
 
 
140
            Statements:
 
141
            - \b MySQL/PostGIS: select AsText(Centroid(GeomFromText('POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2
 
142
                            ,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))')))
 
143
            - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null,
 
144
                    sdo_elem_info_array(1, 1003, 1), sdo_ordinate_array(2,1.3,2.4,1.7,2.8,1.8
 
145
                    ,3.4,1.2,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))
 
146
                    , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
 
147
                    ,mdsys.sdo_dim_element('y',0,10,.00000005)))
 
148
                    from dual
 
149
            - \b SQL Server 2008: select geometry::STGeomFromText('POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2
 
150
                            ,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))',0)
 
151
                            .STCentroid()
 
152
                            .STAsText()
 
153
 
 
154
            With holes:
 
155
            - this:       POINT(4.04663 1.6349)
 
156
            - geolib:     POINT(4.04675 1.65735)
 
157
            - MySQL:      POINT(3.6090580503834 1.607573932092)
 
158
            - PostGIS:    POINT(4.0466265060241 1.63489959839357)
 
159
            - Oracle:           4.0466265060241 1.63489959839357
 
160
            - SQL Server: POINT(4.0466264962959677 1.6348996057331333)
 
161
 
 
162
            Statements:
 
163
            - \b MySQL/PostGIS: select AsText(Centroid(GeomFromText('POLYGON((2 1.3,2.4 1.7,2.8 1.8,3.4 1.2
 
164
                    ,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)
 
165
                    ,(4 2,4.2 1.4,4.8 1.9,4.4 2.2,4 2))')));
 
166
            - \b Oracle: select sdo_geom.sdo_centroid(sdo_geometry(2003, null, null
 
167
                    , sdo_elem_info_array(1, 1003, 1, 25, 2003, 1)
 
168
                    , sdo_ordinate_array(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
 
169
                    ,2.6,5.4,1.2,4.9,0.8,2.9,0.7,2,1.3,4,2, 4.2,1.4, 4.8,1.9, 4.4,2.2, 4,2))
 
170
                    , mdsys.sdo_dim_array(mdsys.sdo_dim_element('x',0,10,.00000005)
 
171
                    ,mdsys.sdo_dim_element('y',0,10,.00000005)))
 
172
                    from dual
 
173
         */
 
174
        template
 
175
        <
 
176
            typename CentroidPointType,
 
177
            typename SegmentPointType = CentroidPointType,
 
178
            typename CalculationType = void
 
179
        >
 
180
        class bashein_detmer
 
181
        {
 
182
            private :
 
183
                // If user specified a calculation type, use that type,
 
184
                //   whatever it is and whatever the point-type(s) are.
 
185
                // Else, use the most appropriate coordinate type
 
186
                //    of the two points, but at least double
 
187
                typedef typename
 
188
                    boost::mpl::if_c
 
189
                    <
 
190
                        boost::is_void<CalculationType>::type::value,
 
191
                        typename select_most_precise
 
192
                        <
 
193
                            typename select_coordinate_type
 
194
                                <
 
195
                                    CentroidPointType,
 
196
                                    SegmentPointType
 
197
                                >::type,
 
198
                            double
 
199
                        >::type,
 
200
                        CalculationType
 
201
                    >::type calc_type;
 
202
 
 
203
                /*! subclass to keep state */
 
204
                struct sums
 
205
                {
 
206
                    calc_type sum_a2;
 
207
                    calc_type sum_x;
 
208
                    calc_type sum_y;
 
209
                    inline sums()
 
210
                        : sum_a2(calc_type())
 
211
                        , sum_x(calc_type())
 
212
                        , sum_y(calc_type())
 
213
                    {
 
214
                        typedef calc_type ct;
 
215
                        //std::cout << "-> calctype: " << typeid(ct).name()
 
216
                        //    << " size: " << sizeof(ct)
 
217
                        //    << " init: " << sum_a2
 
218
                        //    << std::endl;
 
219
                    }
 
220
 
 
221
                    inline void centroid(CentroidPointType& point)
 
222
                    {
 
223
                        if (sum_a2 != 0)
 
224
                        {
 
225
                            calc_type const v3 = 3;
 
226
                            calc_type const a3 = v3 * sum_a2;
 
227
 
 
228
                            typedef typename ggl::coordinate_type
 
229
                                <
 
230
                                    CentroidPointType
 
231
                                >::type coordinate_type;
 
232
 
 
233
                            set<0>(point,
 
234
                                boost::numeric_cast<coordinate_type>(sum_x / a3));
 
235
                            set<1>(point,
 
236
                                boost::numeric_cast<coordinate_type>(sum_y / a3));
 
237
                        }
 
238
                        else
 
239
                        {
 
240
                            // exception?
 
241
                        }
 
242
                    }
 
243
 
 
244
                };
 
245
 
 
246
            public :
 
247
                typedef sums state_type;
 
248
 
 
249
                inline bool operator()(segment<const SegmentPointType> const& s, state_type& state) const
 
250
                {
 
251
                    /* Algorithm:
 
252
                    For each segment:
 
253
                    begin
 
254
                        ai = x1 * y2 - x2 * y1;
 
255
                        sum_a2 += ai;
 
256
                        sum_x += ai * (x1 + x2);
 
257
                        sum_y += ai * (y1 + y2);
 
258
                    end
 
259
                    return POINT(sum_x / (3 * sum_a2), sum_y / (3 * sum_a2) )
 
260
                    */
 
261
 
 
262
                    // Get coordinates and promote them to calc_type
 
263
                    calc_type const x1 = boost::numeric_cast<calc_type>(get<0, 0>(s));
 
264
                    calc_type const y1 = boost::numeric_cast<calc_type>(get<0, 1>(s));
 
265
                    calc_type const x2 = boost::numeric_cast<calc_type>(get<1, 0>(s));
 
266
                    calc_type const y2 = boost::numeric_cast<calc_type>(get<1, 1>(s));
 
267
                    calc_type const ai = x1 * y2 - x2 * y1;
 
268
                    state.sum_a2 += ai;
 
269
                    state.sum_x += ai * (x1 + x2);
 
270
                    state.sum_y += ai * (y1 + y2);
 
271
 
 
272
                    return true;
 
273
                }
 
274
 
 
275
        };
 
276
    } // namespace centroid
 
277
 
 
278
} // namespace strategy
 
279
 
 
280
 
 
281
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 
282
template <typename P, typename PS>
 
283
struct strategy_centroid<cartesian_tag, P, PS>
 
284
{
 
285
    typedef strategy::centroid::bashein_detmer<P, PS> type;
 
286
};
 
287
#endif
 
288
 
 
289
} // namespace ggl
 
290
 
 
291
 
 
292
#endif // GGL_STRATEGY_CARTESIAN_CENTROID_HPP