algorithms/detail/calculate_point_order.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2019-2021, Oracle and/or its affiliates.
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Licensed under the Boost Software License version 1.0.
0009 // http://www.boost.org/users/license.html
0010
0011 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_CALCULATE_POINT_ORDER_HPP
0012 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_CALCULATE_POINT_ORDER_HPP
0013
0014
0015 #include <vector>
0016
0017 #include <boost/range/size.hpp>
0018
0019 #include <boost/geometry/algorithms/area.hpp>
0020 #include <boost/geometry/core/point_order.hpp>
0021 #include <boost/geometry/core/radian_access.hpp>
0022 #include <boost/geometry/core/static_assert.hpp>
0023 #include <boost/geometry/strategies/geographic/point_order.hpp>
0024 #include <boost/geometry/util/math.hpp>
0025 #include <boost/geometry/util/range.hpp>
0026
0027
0028 namespace boost { namespace geometry
0029 {
0030
0031 namespace detail
0032 {
0033
0034 template <typename Iter, typename CalcT>
0035 struct clean_point
0036 {
0037     explicit clean_point(Iter const& iter)
0038         : m_iter(iter), m_azi(0), m_razi(0), m_azi_diff(0)
0039         , m_is_azi_valid(false), m_is_azi_diff_valid(false)
0040     {}
0041
0042     decltype(auto) ref() const
0043     {
0044         return *m_iter;
0045     }
0046
0047     CalcT const& azimuth() const
0048     {
0049         return m_azi;
0050     }
0051
0052     CalcT const& reverse_azimuth() const
0053     {
0054         return m_razi;
0055     }
0056
0057     CalcT const& azimuth_difference() const
0058     {
0059         return m_azi_diff;
0060     }
0061
0062     void set_azimuths(CalcT const& azi, CalcT const& razi)
0063     {
0064         m_azi = azi;
0065         m_razi = razi;
0066         m_is_azi_valid = true;
0067     }
0068
0069     void set_azimuth_invalid()
0070     {
0071         m_is_azi_valid = false;
0072     }
0073
0074     bool is_azimuth_valid() const
0075     {
0076         return m_is_azi_valid;
0077     }
0078
0079     void set_azimuth_difference(CalcT const& diff)
0080     {
0081         m_azi_diff = diff;
0082         m_is_azi_diff_valid = true;
0083     }
0084
0085     void set_azimuth_difference_invalid()
0086     {
0087         m_is_azi_diff_valid = false;
0088     }
0089
0090     bool is_azimuth_difference_valid() const
0091     {
0092         return m_is_azi_diff_valid;
0093     }
0094
0095 private:
0096     Iter m_iter;
0097     CalcT m_azi;
0098     CalcT m_razi;
0099     CalcT m_azi_diff;
0100     // NOTE: these flags could be removed and replaced with some magic number
0101     //       assigned to the above variables, e.g. CalcT(1000).
0102     bool m_is_azi_valid;
0103     bool m_is_azi_diff_valid;
0104 };
0105
0106 struct calculate_point_order_by_azimuth
0107 {
0108     template <typename Ring, typename Strategy>
0109     static geometry::order_selector apply(Ring const& ring, Strategy const& strategy)
0110     {
0111         typedef typename boost::range_iterator<Ring const>::type iter_t;
0112         typedef typename Strategy::template result_type<Ring>::type calc_t;
0113         typedef clean_point<iter_t, calc_t> clean_point_t;
0114
0115         calc_t const zero = 0;
0116         calc_t const pi = math::pi<calc_t>();
0117
0118         std::size_t const count = boost::size(ring);
0119         if (count < 3)
0120         {
0121             return geometry::order_undetermined;
0122         }
0123
0124         // non-duplicated, non-spike points
0125         std::vector<clean_point_t> cleaned;
0126         cleaned.reserve(count);
0127
0128         for (iter_t it = boost::begin(ring); it != boost::end(ring); ++it)
0129         {
0130             // Add point
0131             cleaned.push_back(clean_point_t(it));
0132
0133             while (cleaned.size() >= 3)
0134             {
0135                 auto it0 = cleaned.end() - 3;
0136                 auto it1 = cleaned.end() - 2;
0137                 auto it2 = cleaned.end() - 1;
0138
0139                 calc_t diff;
0140                 if (get_or_calculate_azimuths_difference(*it0, *it1, *it2, diff, strategy)
0141                     && ! math::equals(math::abs(diff), pi))
0142                 {
0143                     // neither duplicate nor a spike - difference already stored
0144                     break;
0145                 }
0146                 else
0147                 {
0148                     // spike detected
0149                     // TODO: angles have to be invalidated only if spike is detected
0150                     // for duplicates it'd be ok to leave them
0151                     it0->set_azimuth_invalid();
0152                     it0->set_azimuth_difference_invalid();
0153                     it2->set_azimuth_difference_invalid();
0154                     cleaned.erase(it1);
0155                 }
0156             }
0157         }
0158
0159         // filter-out duplicates and spikes at the front and back of cleaned
0160         auto cleaned_b = cleaned.begin();
0161         auto cleaned_e = cleaned.end();
0162         std::size_t cleaned_count = cleaned.size();
0163         bool found = false;
0164         do
0165         {
0166             found = false;
0167             while(cleaned_count >= 3)
0168             {
0169                 auto it0 = cleaned_e - 2;
0170                 auto it1 = cleaned_e - 1;
0171                 auto it2 = cleaned_b;
0172                 auto it3 = cleaned_b + 1;
0173
0174                 calc_t diff = 0;
0175                 if (! get_or_calculate_azimuths_difference(*it0, *it1, *it2, diff, strategy)
0176                     || math::equals(math::abs(diff), pi))
0177                 {
0178                     // spike at the back
0179                     // TODO: angles have to be invalidated only if spike is detected
0180                     // for duplicates it'd be ok to leave them
0181                     it0->set_azimuth_invalid();
0182                     it0->set_azimuth_difference_invalid();
0183                     it2->set_azimuth_difference_invalid();
0184                     --cleaned_e;
0185                     --cleaned_count;
0186                     found = true;
0187                 }
0188                 else if (! get_or_calculate_azimuths_difference(*it1, *it2, *it3, diff, strategy)
0189                          || math::equals(math::abs(diff), pi))
0190                 {
0191                     // spike at the front
0192                     // TODO: angles have to be invalidated only if spike is detected
0193                     // for duplicates it'd be ok to leave them
0194                     it1->set_azimuth_invalid();
0195                     it1->set_azimuth_difference_invalid();
0196                     it3->set_azimuth_difference_invalid();
0197                     ++cleaned_b;
0198                     --cleaned_count;
0199                     found = true;
0200                 }
0201                 else
0202                 {
0203                     break;
0204                 }
0205             }
0206         }
0207         while (found);
0208
0209         if (cleaned_count < 3)
0210         {
0211             return geometry::order_undetermined;
0212         }
0213
0214         // calculate the sum of external angles
0215         calc_t angles_sum = zero;
0216         for (auto it = cleaned_b; it != cleaned_e; ++it)
0217         {
0218             auto it0 = (it == cleaned_b ? cleaned_e - 1 : it - 1);
0219             auto it2 = (it == cleaned_e - 1 ? cleaned_b : it + 1);
0220
0221             calc_t diff = 0;
0222             get_or_calculate_azimuths_difference(*it0, *it, *it2, diff, strategy);
0223
0224             angles_sum += diff;
0225         }
0226
0227 #ifdef BOOST_GEOMETRY_DEBUG_POINT_ORDER
0228         std::cout << angles_sum  << " for " << geometry::wkt(ring) << std::endl;
0229 #endif
0230
0231         return angles_sum == zero ? geometry::order_undetermined
0232              : angles_sum > zero  ? geometry::clockwise
0233                                   : geometry::counterclockwise;
0234     }
0235
0236 private:
0237     template <typename Iter, typename T, typename Strategy>
0238     static bool get_or_calculate_azimuths_difference(clean_point<Iter, T> & p0,
0239                                                      clean_point<Iter, T> & p1,
0240                                                      clean_point<Iter, T> const& p2,
0241                                                      T & diff,
0242                                                      Strategy const& strategy)
0243     {
0244         if (p1.is_azimuth_difference_valid())
0245         {
0246             diff = p1.azimuth_difference();
0247             return true;
0248         }
0249
0250         T azi1, razi1, azi2, razi2;
0251         if (get_or_calculate_azimuths(p0, p1, azi1, razi1, strategy)
0252             && get_or_calculate_azimuths(p1, p2, azi2, razi2, strategy))
0253         {
0254             diff = strategy.apply(p0.ref(), p1.ref(), p2.ref(), razi1, azi2);
0255             p1.set_azimuth_difference(diff);
0256             return true;
0257         }
0258         return false;
0259     }
0260
0261     template <typename Iter, typename T, typename Strategy>
0262     static bool get_or_calculate_azimuths(clean_point<Iter, T> & p0,
0263                                           clean_point<Iter, T> const& p1,
0264                                           T & azi, T & razi,
0265                                           Strategy const& strategy)
0266     {
0267         if (p0.is_azimuth_valid())
0268         {
0269             azi = p0.azimuth();
0270             razi = p0.reverse_azimuth();
0271             return true;
0272         }
0273
0274         if (strategy.apply(p0.ref(), p1.ref(), azi, razi))
0275         {
0276             p0.set_azimuths(azi, razi);
0277             return true;
0278         }
0279
0280         return false;
0281     }
0282 };
0283
0284 struct calculate_point_order_by_area
0285 {
0286     template <typename Ring, typename Strategy>
0287     static geometry::order_selector apply(Ring const& ring, Strategy const& strategy)
0288     {
0289         auto const result = detail::area::ring_area::apply(
0290                                 ring,
0291                                 // TEMP - in the future (umbrella) strategy will be passed
0292                                 geometry::strategies::area::services::strategy_converter
0293                                     <
0294                                         decltype(strategy.get_area_strategy())
0295                                     >::get(strategy.get_area_strategy()));
0296
0297         decltype(result) const zero = 0;
0298         return result == zero ? geometry::order_undetermined
0299              : result > zero  ? geometry::clockwise
0300                               : geometry::counterclockwise;
0301     }
0302 };
0303
0304 } // namespace detail
0305
0306 namespace dispatch
0307 {
0308
0309 template
0310 <
0311     typename Strategy,
0312     typename VersionTag = typename Strategy::version_tag
0313 >
0314 struct calculate_point_order
0315 {
0316     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
0317         "Not implemented for this VersionTag.",
0318         VersionTag);
0319 };
0320
0321 template <typename Strategy>
0322 struct calculate_point_order<Strategy, strategy::point_order::area_tag>
0323     : geometry::detail::calculate_point_order_by_area
0324 {};
0325
0326 template <typename Strategy>
0327 struct calculate_point_order<Strategy, strategy::point_order::azimuth_tag>
0328     : geometry::detail::calculate_point_order_by_azimuth
0329 {};
0330
0331
0332 } // namespace dispatch
0333
0334 namespace detail
0335 {
0336
0337 template <typename Ring, typename Strategy>
0338 inline geometry::order_selector calculate_point_order(Ring const& ring, Strategy const& strategy)
0339 {
0340     concepts::check<Ring const>();
0341
0342     return dispatch::calculate_point_order<Strategy>::apply(ring, strategy);
0343 }
0344
0345 template <typename Ring>
0346 inline geometry::order_selector calculate_point_order(Ring const& ring)
0347 {
0348     typedef typename strategy::point_order::services::default_strategy
0349         <
0350             typename geometry::cs_tag<Ring>::type
0351         >::type strategy_type;
0352
0353     concepts::check<Ring const>();
0354
0355     return dispatch::calculate_point_order<strategy_type>::apply(ring, strategy_type());
0356 }
0357
0358
0359 } // namespace detail
0360
0361 }} // namespace boost::geometry
0362
0363
0364 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_CALCULATE_POINT_ORDER_HPP