boost/polygon/detail/polygon_simplify.hpp

   1 // Copyright 2011, Andrew Ross
   2 //
   3 // Use, modification and distribution are subject to the Boost Software License,
   4 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   5 // http://www.boost.org/LICENSE_1_0.txt).
   6 #ifndef BOOST_POLYGON_DETAIL_SIMPLIFY_HPP
   7 #define BOOST_POLYGON_DETAIL_SIMPLIFY_HPP
   8 #include <vector>
   9
  10 namespace boost { namespace polygon { namespace detail { namespace simplify_detail {
  11
  12   // Does a simplification/optimization pass on the polygon.  If a given
  13   // vertex lies within "len" of the line segment joining its neighbor
  14   // vertices, it is removed.
  15   template <typename T> //T is a model of point concept
  16   std::size_t simplify(std::vector<T>& dst, const std::vector<T>& src,
  17                        typename coordinate_traits<
  18                        typename point_traits<T>::coordinate_type
  19                        >::coordinate_distance len)
  20   {
  21     using namespace boost::polygon;
  22     typedef typename point_traits<T>::coordinate_type coordinate_type;
  23     typedef typename coordinate_traits<coordinate_type>::area_type ftype;
  24     typedef typename std::vector<T>::const_iterator iter;
  25
  26     std::vector<T> out;
  27     out.reserve(src.size());
  28     dst = src;
  29     std::size_t final_result = 0;
  30     std::size_t orig_size = src.size();
  31
  32     //I can't use == if T doesn't provide it, so use generic point concept compare
  33     bool closed = equivalence(src.front(), src.back());
  34
  35     //we need to keep smoothing until we don't find points to remove
  36     //because removing points in the first iteration through the
  37     //polygon may leave it in a state where more removal is possible
  38     bool not_done = true;
  39     while(not_done) {
  40       if(dst.size() < 3) {
  41         dst.clear();
  42         return orig_size;
  43       }
  44
  45       // Start with the second, test for the last point
  46       // explicitly, and exit after looping back around to the first.
  47       ftype len2 = ftype(len) * ftype(len);
  48       for(iter prev=dst.begin(), i=prev+1, next; /**/; i = next) {
  49         next = i+1;
  50         if(next == dst.end())
  51           next = dst.begin();
  52
  53         // points A, B, C
  54         ftype ax = x(*prev), ay = y(*prev);
  55         ftype bx = x(*i), by = y(*i);
  56         ftype cx = x(*next), cy = y(*next);
  57
  58         // vectors AB, BC and AC:
  59         ftype abx = bx-ax, aby = by-ay;
  60         ftype bcx = cx-bx, bcy = cy-by;
  61         ftype acx = cx-ax, acy = cy-ay;
  62
  63         // dot products
  64         ftype ab_ab = abx*abx + aby*aby;
  65         ftype bc_bc = bcx*bcx + bcy*bcy;
  66         ftype ac_ac = acx*acx + acy*acy;
  67         ftype ab_ac = abx*acx + aby*acy;
  68
  69         // projection of AB along AC
  70         ftype projf = ab_ac / ac_ac;
  71         ftype projx = acx * projf, projy = acy * projf;
  72
  73         // perpendicular vector from the line AC to point B (i.e. AB - proj)
  74         ftype perpx = abx - projx, perpy = aby - projy;
  75
  76         // Squared fractional distance of projection. FIXME: can
  77         // remove this division, the decisions below can be made with
  78         // just the sign of the quotient and a check to see if
  79         // abs(numerator) is greater than abs(divisor).
  80         ftype f2 = (projx*acx + projy*acx) / ac_ac;
  81
  82         // Square of the relevant distance from point B:
  83         ftype dist2;
  84         if     (f2 < 0) dist2 = ab_ab;
  85         else if(f2 > 1) dist2 = bc_bc;
  86         else            dist2 = perpx*perpx + perpy*perpy;
  87
  88         if(dist2 > len2) {
  89           prev = i; // bump prev, we didn't remove the segment
  90           out.push_back(*i);
  91         }
  92
  93         if(i == dst.begin())
  94           break;
  95       }
  96       std::size_t result = dst.size() - out.size();
  97       if(result == 0) {
  98         not_done = false;
  99       } else {
 100         final_result += result;
 101         dst = out;
 102         out.clear();
 103       }
 104     } //end of while loop
 105     if(closed) {
 106       //if the input was closed we want the output to be closed
 107       --final_result;
 108       dst.push_back(dst.front());
 109     }
 110     return final_result;
 111   }
 112
 113
 114 }}}}
 115
 116 #endif