annotate DEPENDENCIES/generic/include/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp @ 133:4acb5d8d80b6 tip

Don't fail environmental check if README.md exists (but .txt and no-suffix don't)
author Chris Cannam
date Tue, 30 Jul 2019 12:25:44 +0100
parents c530137014c0
children
rev   line source
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