polygon/detail/polygon_simplify.hpp

0001 // Copyright 2011, Andrew Ross
0002 //
0003 // Use, modification and distribution are subject to the Boost Software License,
0004 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0005 // http://www.boost.org/LICENSE_1_0.txt).
0006 #ifndef BOOST_POLYGON_DETAIL_SIMPLIFY_HPP
0007 #define BOOST_POLYGON_DETAIL_SIMPLIFY_HPP
0008 #include <vector>
0009
0010 namespace boost { namespace polygon { namespace detail { namespace simplify_detail {
0011
0012   // Does a simplification/optimization pass on the polygon.  If a given
0013   // vertex lies within "len" of the line segment joining its neighbor
0014   // vertices, it is removed.
0015   template <typename T> //T is a model of point concept
0016   std::size_t simplify(std::vector<T>& dst, const std::vector<T>& src,
0017                        typename coordinate_traits<
0018                        typename point_traits<T>::coordinate_type
0019                        >::coordinate_distance len)
0020   {
0021     using namespace boost::polygon;
0022     typedef typename point_traits<T>::coordinate_type coordinate_type;
0023     typedef typename coordinate_traits<coordinate_type>::area_type ftype;
0024     typedef typename std::vector<T>::const_iterator iter;
0025
0026     std::vector<T> out;
0027     out.reserve(src.size());
0028     dst = src;
0029     std::size_t final_result = 0;
0030     std::size_t orig_size = src.size();
0031
0032     //I can't use == if T doesn't provide it, so use generic point concept compare
0033     bool closed = equivalence(src.front(), src.back());
0034
0035     //we need to keep smoothing until we don't find points to remove
0036     //because removing points in the first iteration through the
0037     //polygon may leave it in a state where more removal is possible
0038     bool not_done = true;
0039     while(not_done) {
0040       if(dst.size() < 3) {
0041         dst.clear();
0042         return orig_size;
0043       }
0044
0045       // Start with the second, test for the last point
0046       // explicitly, and exit after looping back around to the first.
0047       ftype len2 = ftype(len) * ftype(len);
0048       for(iter prev=dst.begin(), i=prev+1, next; /**/; i = next) {
0049         next = i+1;
0050         if(next == dst.end())
0051           next = dst.begin();
0052
0053         // points A, B, C
0054         ftype ax = x(*prev), ay = y(*prev);
0055         ftype bx = x(*i), by = y(*i);
0056         ftype cx = x(*next), cy = y(*next);
0057
0058         // vectors AB, BC and AC:
0059         ftype abx = bx-ax, aby = by-ay;
0060         ftype bcx = cx-bx, bcy = cy-by;
0061         ftype acx = cx-ax, acy = cy-ay;
0062
0063         // dot products
0064         ftype ab_ab = abx*abx + aby*aby;
0065         ftype bc_bc = bcx*bcx + bcy*bcy;
0066         ftype ac_ac = acx*acx + acy*acy;
0067         ftype ab_ac = abx*acx + aby*acy;
0068
0069         // projection of AB along AC
0070         ftype projf = ab_ac / ac_ac;
0071         ftype projx = acx * projf, projy = acy * projf;
0072
0073         // perpendicular vector from the line AC to point B (i.e. AB - proj)
0074         ftype perpx = abx - projx, perpy = aby - projy;
0075
0076         // Squared fractional distance of projection. FIXME: can
0077         // remove this division, the decisions below can be made with
0078         // just the sign of the quotient and a check to see if
0079         // abs(numerator) is greater than abs(divisor).
0080         ftype f2 = (projx*acx + projy*acx) / ac_ac;
0081
0082         // Square of the relevant distance from point B:
0083         ftype dist2;
0084         if     (f2 < 0) dist2 = ab_ab;
0085         else if(f2 > 1) dist2 = bc_bc;
0086         else            dist2 = perpx*perpx + perpy*perpy;
0087
0088         if(dist2 > len2) {
0089           prev = i; // bump prev, we didn't remove the segment
0090           out.push_back(*i);
0091         }
0092
0093         if(i == dst.begin())
0094           break;
0095       }
0096       std::size_t result = dst.size() - out.size();
0097       if(result == 0) {
0098         not_done = false;
0099       } else {
0100         final_result += result;
0101         dst = out;
0102         out.clear();
0103       }
0104     } //end of while loop
0105     if(closed) {
0106       //if the input was closed we want the output to be closed
0107       --final_result;
0108       dst.push_back(dst.front());
0109     }
0110     return final_result;
0111   }
0112
0113
0114 }}}}
0115
0116 #endif