strategies/cartesian/centroid_weighted_length.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
0004 // Copyright (c) 2009-2015 Barend Gehrels, Amsterdam, the Netherlands.
0005
0006 // This file was modified by Oracle on 2015-2023.
0007 // Modifications copyright (c) 2015-2023, Oracle and/or its affiliates.
0008 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010
0011 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
0012 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
0013
0014 // Use, modification and distribution is subject to the Boost Software License,
0015 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0016 // http://www.boost.org/LICENSE_1_0.txt)
0017
0018 #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_WEIGHTED_LENGTH_HPP
0019 #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_WEIGHTED_LENGTH_HPP
0020
0021 #include <boost/math/special_functions/fpclassify.hpp>
0022 #include <boost/numeric/conversion/cast.hpp>
0023
0024 #include <boost/geometry/algorithms/assign.hpp>
0025
0026 #include <boost/geometry/arithmetic/arithmetic.hpp>
0027
0028 // Helper geometry
0029 #include <boost/geometry/geometries/point.hpp>
0030
0031 #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
0032 #include <boost/geometry/strategies/centroid.hpp>
0033
0034 #include <boost/geometry/util/algorithm.hpp>
0035 #include <boost/geometry/util/select_most_precise.hpp>
0036
0037
0038 namespace boost { namespace geometry
0039 {
0040
0041 namespace strategy { namespace centroid
0042 {
0043
0044 template
0045 <
0046     typename Ignored1 = void,
0047     typename Ignored2 = void,
0048     typename CalculationType = void
0049 >
0050 class weighted_length
0051 {
0052 private :
0053     typedef geometry::strategy::distance::pythagoras<CalculationType> pythagoras_strategy;
0054
0055     template <typename GeometryPoint, typename ResultPoint>
0056     struct calculation_type
0057     {
0058         // Below the distance between two GeometryPoints is calculated.
0059         // ResultPoint is taken into account by passing them together here.
0060         typedef typename pythagoras_strategy::template calculation_type
0061             <
0062                 GeometryPoint, ResultPoint
0063             >::type type;
0064     };
0065
0066     template <typename GeometryPoint, typename ResultPoint>
0067     class sums
0068     {
0069         friend class weighted_length;
0070         template <typename, typename> friend struct set_sum_div_length;
0071
0072         typedef typename calculation_type<GeometryPoint, ResultPoint>::type calc_type;
0073         typedef typename geometry::model::point
0074             <
0075                 calc_type,
0076                 geometry::dimension<ResultPoint>::value,
0077                 cs::cartesian
0078             > work_point;
0079
0080         calc_type length;
0081         work_point average_sum;
0082
0083     public:
0084         inline sums()
0085             : length(calc_type())
0086         {
0087             geometry::assign_zero(average_sum);
0088         }
0089     };
0090
0091 public :
0092     template <typename GeometryPoint, typename ResultPoint>
0093     struct state_type
0094     {
0095         typedef sums<GeometryPoint, ResultPoint> type;
0096     };
0097
0098     template <typename GeometryPoint, typename ResultPoint>
0099     static inline void apply(GeometryPoint const& p1, GeometryPoint const& p2,
0100                              sums<GeometryPoint, ResultPoint>& state)
0101     {
0102         typedef typename calculation_type<GeometryPoint, ResultPoint>::type distance_type;
0103
0104         distance_type const d = pythagoras_strategy::apply(p1, p2);
0105         state.length += d;
0106
0107         distance_type const d_half = d / distance_type(2);
0108         geometry::detail::for_each_dimension<ResultPoint>([&](auto dimension)
0109         {
0110             distance_type const coord1 = get<dimension>(p1);
0111             distance_type const coord2 = get<dimension>(p2);
0112             distance_type const wm = (coord1 + coord2) * d_half; // weighted median
0113             set<dimension>(state.average_sum, get<dimension>(state.average_sum) + wm);
0114         });
0115     }
0116
0117     template <typename GeometryPoint, typename ResultPoint>
0118     static inline bool result(sums<GeometryPoint, ResultPoint> const& state,
0119                               ResultPoint& centroid)
0120     {
0121         typedef typename calculation_type<GeometryPoint, ResultPoint>::type distance_type;
0122
0123         distance_type const zero = distance_type();
0124         if (! geometry::math::equals(state.length, zero)
0125             && boost::math::isfinite(state.length)) // Prevent NaN centroid coordinates
0126         {
0127             // NOTE: above distance_type is checked, not the centroid coordinate_type
0128             // which means that the centroid can still be filled with INF
0129             // if e.g. distance_type is double and centroid contains floats
0130             geometry::detail::for_each_dimension<ResultPoint>([&](auto dimension)
0131             {
0132                 typedef typename geometry::coordinate_type<ResultPoint>::type coordinate_type;
0133                 geometry::set<dimension>(
0134                     centroid,
0135                     boost::numeric_cast<coordinate_type>(
0136                         geometry::get<dimension>(state.average_sum) / state.length
0137                     )
0138                 );
0139             });
0140             return true;
0141         }
0142
0143         return false;
0144     }
0145 };
0146
0147 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0148
0149 namespace services
0150 {
0151
0152
0153 // Register this strategy for linear geometries, in all dimensions
0154
0155 template <std::size_t N, typename Point, typename Geometry>
0156 struct default_strategy
0157 <
0158     cartesian_tag,
0159     linear_tag,
0160     N,
0161     Point,
0162     Geometry
0163 >
0164 {
0165     typedef weighted_length
0166         <
0167             Point,
0168             typename point_type<Geometry>::type
0169         > type;
0170 };
0171
0172
0173 } // namespace services
0174
0175
0176 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0177
0178
0179 }} // namespace strategy::centroid
0180
0181
0182 }} // namespace boost::geometry
0183
0184
0185 #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_WEIGHTED_LENGTH_HPP