← Back to branch summary

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

« back to all changes in this revision

Viewing changes to include/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/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