Generic Iterator Based Douglas Peucker

douglas_peuker.hpp

#ifndef DOUGLASS_PEUKER_H_
#define DOUGLASS_PEUKER_H_

#include <stack>
#include <tuple>
#include <vector>


/// @{
/// @tparam InputIterator  Any STL style container
/// @tparam OutputIterator An output iterator or inserter.
/// @tparam Function  Computes the distance of a point from the surface.
/// @tparam Epsilon  Should be the same type as the return of Function
///
/// Generalized Douglas-Peucker algorithm for path simplification.
///
/// @code
/// #include <cassert>
/// #include <cmath>
/// #include <ctime>
///
/// struct Point {
///   double x;
///   double y;
/// };
///
/// void
/// circle(std::vector<Point>& data)
/// {
///   size_t Size = 360;
///   double delta = 2 * M_PI / Size;
///
///   for (size_t i = 0; i < Size; ++i) {
///     Point p;
///     p.x = sin(delta * i);
///     p.y = cos(delta * i);
///     data.push_back(p);
///   }
/// }
///
/// int
/// main()
/// {
///   using namespace std;
///
///   vector<Point> dataI;
///   circle(dataI);
///
///   vector<Point> dataO;
///   size_t size =
///     DouglasPeucker(dataI.begin(),
///                    dataI.end(),
///                    .05,
///                    std::back_inserter(dataO),
///                    [](const Point& elem, const Point& head, const Point& tail) -> double {
///                      auto n = (tail.x - head.x) * (head.y - elem.y)
///                             - (head.x - elem.x) * (tail.y - head.y);
///                      auto d = (tail.x - head.x) * (tail.x - head.x)
///                             + (tail.y - head.x) * (tail.y - head.y);
///
///                      return std::abs(n) / std::sqrt(d);
///                    });
///   assert(size == dataO.size());
/// }
/// @endcode

template <
  typename InputIterator,
  typename Function,
  typename Epsilon=typename std::result_of<Function>::value
  >
size_t
DouglasPeucker(const InputIterator& begin,
               const InputIterator& end,
               Epsilon epsilon,
               std::vector<bool>& included,
               Function computeDistance)
{
  using IteratorPair = std::pair<InputIterator, InputIterator>;
  using Stack = std::stack<IteratorPair>;
  using BoolVector = std::vector<bool>;

  // Initialize
  included.resize(std::distance(begin, end), true);
  Stack stack;
  stack.push(std::make_pair(begin, end));

  size_t size = std::distance(begin, end);
  while (!stack.empty()) {
    // {
    auto  elem = stack.top();
    auto& head = elem.first;
    auto& tail = elem.second;
    stack.pop();
    // }

    // {
    Epsilon deltaMax(0);
    InputIterator winner;

    for (auto iter = head + 1; iter < tail; ++iter) {
        auto offset = std::distance(begin, iter);
        if (included[offset]) {
          Epsilon delta = computeDistance(*iter, *head, *(tail - 1));

          if (delta > deltaMax) {
            winner = iter;
            deltaMax = delta;
          }
        }
    }

    if (deltaMax > epsilon) {
      stack.push(std::make_pair(head, winner));
      stack.push(std::make_pair(winner, tail));
    } else {
      for (auto iter = (head + 1); iter < tail; ++iter) {
        auto offset = std::distance(begin, iter);
        included[offset] = false;
        --size;
      }
    }
  }

  return size;
}

/// Varitation 1:
template <
  typename InputIterator,
  typename OutputIterator,
  typename Function,
  typename Epsilon=typename std::result_of<Function>::value
  >
size_t
DouglasPeucker(const InputIterator& begin,
               const InputIterator& end,
               Epsilon epsilon,
               OutputIterator output,
               Function computeDistance)
{
  std::vector<bool> included;
  size_t size = DouglasPeucker(begin, end, epsilon, included, computeDistance);

  for (auto iter = begin; iter != end; ++iter) {
    auto offset = std::distance(begin, iter);
    if (included[offset]) {
      *(output++) = *iter;
    }
  }

  return size;
}

/// Varitation 2:
template <
  typename InputIterator,
  typename OutputContainer,
  typename Function,
  typename Epsilon=typename std::result_of<Function>::value
  >
size_t
DouglasPeuckerIndices(const InputIterator& begin,
                      const InputIterator& end,
                      Epsilon epsilon,
                      OutputContainer output,
                      Function computeDistance)
{
  using Container = typename OutputContainer::container_type;
  static_assert(
      std::is_integral<typename Container::value_type>::value,
      "Output container should be of an integral type"
  );

  std::vector<bool> included;
  size_t size = DouglasPeucker(begin, end, epsilon, included, computeDistance);

  for (auto iter = begin; iter != end; ++iter) {
    auto offset = std::distance(begin, iter);
    if (included[offset]) {
      *(output++) = offset;
    }
  }

  return size;
}
/// @}

#endif // DOUGLASS_PEUKER_H_