/*
 * Software License Agreement (BSD License)
 *
 *  Copyright (c) 2010, Willow Garage, Inc.
 *  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions
 *  are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *   * Redistributions in binary form must reproduce the above
 *     copyright notice, this list of conditions and the following
 *     disclaimer in the documentation and/or other materials provided
 *     with the distribution.
 *   * Neither the name of Willow Garage, Inc. nor the names of its
 *     contributors may be used to endorse or promote products derived
 *     from this software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 *  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 *  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 *  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 *  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 *  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 *  POSSIBILITY OF SUCH DAMAGE.
 *
 */

#ifndef PCL_GEOMETRY_POLYGON_OPERATIONS_HPP_
#define PCL_GEOMETRY_POLYGON_OPERATIONS_HPP_

#include <pcl/geometry/polygon_operations.h>

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
void
pcl::approximatePolygon(const PlanarPolygon<PointT>& polygon,
                        PlanarPolygon<PointT>& approx_polygon,
                        float threshold,
                        bool refine,
                        bool closed)
{
  const Eigen::Vector4f& coefficients = polygon.getCoefficients();
  const typename pcl::PointCloud<PointT>::VectorType& contour = polygon.getContour();

  Eigen::Vector3f rotation_axis(coefficients[1], -coefficients[0], 0.0f);
  rotation_axis.normalize();

  float rotation_angle = acosf(coefficients[2]);
  Eigen::Affine3f transformation = Eigen::Translation3f(0, 0, coefficients[3]) *
                                   Eigen::AngleAxisf(rotation_angle, rotation_axis);

  typename pcl::PointCloud<PointT>::VectorType polygon2D(contour.size());
  for (std::size_t pIdx = 0; pIdx < polygon2D.size(); ++pIdx)
    polygon2D[pIdx].getVector3fMap() = transformation * contour[pIdx].getVector3fMap();

  typename pcl::PointCloud<PointT>::VectorType approx_polygon2D;
  approximatePolygon2D<PointT>(polygon2D, approx_polygon2D, threshold, refine, closed);

  typename pcl::PointCloud<PointT>::VectorType& approx_contour =
      approx_polygon.getContour();
  approx_contour.resize(approx_polygon2D.size());

  Eigen::Affine3f inv_transformation = transformation.inverse();
  for (std::size_t pIdx = 0; pIdx < approx_polygon2D.size(); ++pIdx)
    approx_contour[pIdx].getVector3fMap() =
        inv_transformation * approx_polygon2D[pIdx].getVector3fMap();
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
template <typename PointT>
void
pcl::approximatePolygon2D(const typename pcl::PointCloud<PointT>::VectorType& polygon,
                          typename pcl::PointCloud<PointT>::VectorType& approx_polygon,
                          float threshold,
                          bool refine,
                          bool closed)
{
  approx_polygon.clear();
  if (polygon.size() < 3)
    return;

  std::vector<std::pair<unsigned, unsigned>> intervals;
  std::pair<unsigned, unsigned> interval(0, 0);

  if (closed) {
    float max_distance = .0f;
    for (std::size_t idx = 1; idx < polygon.size(); ++idx) {
      float distance =
          (polygon[0].x - polygon[idx].x) * (polygon[0].x - polygon[idx].x) +
          (polygon[0].y - polygon[idx].y) * (polygon[0].y - polygon[idx].y);

      if (distance > max_distance) {
        max_distance = distance;
        interval.second = idx;
      }
    }

    for (std::size_t idx = 1; idx < polygon.size(); ++idx) {
      float distance = (polygon[interval.second].x - polygon[idx].x) *
                           (polygon[interval.second].x - polygon[idx].x) +
                       (polygon[interval.second].y - polygon[idx].y) *
                           (polygon[interval.second].y - polygon[idx].y);

      if (distance > max_distance) {
        max_distance = distance;
        interval.first = idx;
      }
    }

    if (max_distance < threshold * threshold)
      return;

    intervals.push_back(interval);
    std::swap(interval.first, interval.second);
    intervals.push_back(interval);
  }
  else {
    interval.first = 0;
    interval.second = static_cast<unsigned int>(polygon.size()) - 1;
    intervals.push_back(interval);
  }

  std::vector<unsigned> result;
  // recursively refine
  while (!intervals.empty()) {
    std::pair<unsigned, unsigned>& currentInterval = intervals.back();
    float line_x = polygon[currentInterval.first].y - polygon[currentInterval.second].y;
    float line_y = polygon[currentInterval.second].x - polygon[currentInterval.first].x;
    float line_d =
        polygon[currentInterval.first].x * polygon[currentInterval.second].y -
        polygon[currentInterval.first].y * polygon[currentInterval.second].x;

    float linelen = 1.0f / std::sqrt(line_x * line_x + line_y * line_y);

    line_x *= linelen;
    line_y *= linelen;
    line_d *= linelen;

    float max_distance = 0.0;
    unsigned first_index = currentInterval.first + 1;
    unsigned max_index = 0;

    // => 0-crossing
    if (currentInterval.first > currentInterval.second) {
      for (std::size_t idx = first_index; idx < polygon.size(); idx++) {
        float distance =
            std::abs(line_x * polygon[idx].x + line_y * polygon[idx].y + line_d);
        if (distance > max_distance) {
          max_distance = distance;
          max_index = idx;
        }
      }
      first_index = 0;
    }

    for (unsigned int idx = first_index; idx < currentInterval.second; idx++) {
      float distance =
          std::abs(line_x * polygon[idx].x + line_y * polygon[idx].y + line_d);
      if (distance > max_distance) {
        max_distance = distance;
        max_index = idx;
      }
    }

    if (max_distance > threshold) {
      std::pair<unsigned, unsigned> interval(max_index, currentInterval.second);
      currentInterval.second = max_index;
      intervals.push_back(interval);
    }
    else {
      result.push_back(currentInterval.second);
      intervals.pop_back();
    }
  }

  approx_polygon.reserve(result.size());
  if (refine) {
    std::vector<Eigen::Vector3f, Eigen::aligned_allocator<Eigen::Vector3f>> lines(
        result.size());
    std::reverse(result.begin(), result.end());
    for (std::size_t rIdx = 0; rIdx < result.size(); ++rIdx) {
      std::size_t nIdx = rIdx + 1;
      if (nIdx == result.size())
        nIdx = 0;

      Eigen::Vector2f centroid = Eigen::Vector2f::Zero();
      Eigen::Matrix2f covariance = Eigen::Matrix2f::Zero();
      std::size_t pIdx = result[rIdx];
      unsigned num_points = 0;
      if (pIdx > result[nIdx]) {
        num_points = static_cast<unsigned>(polygon.size()) - pIdx;
        for (; pIdx < polygon.size(); ++pIdx) {
          covariance.coeffRef(0) += polygon[pIdx].x * polygon[pIdx].x;
          covariance.coeffRef(1) += polygon[pIdx].x * polygon[pIdx].y;
          covariance.coeffRef(3) += polygon[pIdx].y * polygon[pIdx].y;
          centroid[0] += polygon[pIdx].x;
          centroid[1] += polygon[pIdx].y;
        }
        pIdx = 0;
      }

      num_points += result[nIdx] - pIdx;
      for (; pIdx < result[nIdx]; ++pIdx) {
        covariance.coeffRef(0) += polygon[pIdx].x * polygon[pIdx].x;
        covariance.coeffRef(1) += polygon[pIdx].x * polygon[pIdx].y;
        covariance.coeffRef(3) += polygon[pIdx].y * polygon[pIdx].y;
        centroid[0] += polygon[pIdx].x;
        centroid[1] += polygon[pIdx].y;
      }

      covariance.coeffRef(2) = covariance.coeff(1);

      float norm = 1.0f / float(num_points);
      centroid *= norm;
      covariance *= norm;
      covariance.coeffRef(0) -= centroid[0] * centroid[0];
      covariance.coeffRef(1) -= centroid[0] * centroid[1];
      covariance.coeffRef(3) -= centroid[1] * centroid[1];

      float eval;
      Eigen::Vector2f normal;
      eigen22(covariance, eval, normal);

      // select the one which is more "parallel" to the original line
      Eigen::Vector2f direction;
      direction[0] = polygon[result[nIdx]].x - polygon[result[rIdx]].x;
      direction[1] = polygon[result[nIdx]].y - polygon[result[rIdx]].y;
      direction.normalize();

      if (std::abs(direction.dot(normal)) > float(M_SQRT1_2)) {
        std::swap(normal[0], normal[1]);
        normal[0] = -normal[0];
      }

      // needs to be on the left side of the edge
      if (direction[0] * normal[1] < direction[1] * normal[0])
        normal *= -1.0;

      lines[rIdx].head<2>().matrix() = normal;
      lines[rIdx][2] = -normal.dot(centroid);
    }

    float threshold2 = threshold * threshold;
    for (std::size_t rIdx = 0; rIdx < lines.size(); ++rIdx) {
      std::size_t nIdx = rIdx + 1;
      if (nIdx == result.size())
        nIdx = 0;

      Eigen::Vector3f vertex = lines[rIdx].cross(lines[nIdx]);
      vertex /= vertex[2];
      vertex[2] = 0.0;

      PointT point;
      // test whether we need another edge since the intersection point is too far away
      // from the original vertex
      Eigen::Vector3f pq = polygon[result[nIdx]].getVector3fMap() - vertex;
      pq[2] = 0.0;

      float distance = pq.squaredNorm();
      if (distance > threshold2) {
        // test whether the old point is inside the new polygon or outside
        if ((pq[0] * lines[rIdx][0] + pq[1] * lines[rIdx][1] < 0.0) &&
            (pq[0] * lines[nIdx][0] + pq[1] * lines[nIdx][1] < 0.0)) {
          float distance1 = lines[rIdx][0] * polygon[result[nIdx]].x +
                            lines[rIdx][1] * polygon[result[nIdx]].y + lines[rIdx][2];
          float distance2 = lines[nIdx][0] * polygon[result[nIdx]].x +
                            lines[nIdx][1] * polygon[result[nIdx]].y + lines[nIdx][2];

          point.x = polygon[result[nIdx]].x - distance1 * lines[rIdx][0];
          point.y = polygon[result[nIdx]].y - distance1 * lines[rIdx][1];

          approx_polygon.push_back(point);

          vertex[0] = polygon[result[nIdx]].x - distance2 * lines[nIdx][0];
          vertex[1] = polygon[result[nIdx]].y - distance2 * lines[nIdx][1];
        }
      }
      point.getVector3fMap() = vertex;
      approx_polygon.push_back(point);
    }
  }
  else {
    // we have a new polygon in results, but inverted (clockwise <-> counter-clockwise)
    for (std::vector<unsigned>::reverse_iterator it = result.rbegin();
         it != result.rend();
         ++it)
      approx_polygon.push_back(polygon[*it]);
  }
}

#endif // PCL_GEOMETRY_POLYGON_OPERATIONS_HPP_