polygoncutting.hh

#ifndef DUNE_MMESH_GRID_POLYGONCUTTING_HH
#define DUNE_MMESH_GRID_POLYGONCUTTING_HH
 
#include <algorithm>
#include <iostream>
 
namespace Dune {
template <class Scalar, class Point>
class PolygonCutting {
 public:
  // implements the Sutherland-Hodgman algorithm to determine the polygon that
  // emerges from clipping a subject polygon with respect to the edges of a
  // clip polygon, subject polygon and clip polygon as well as the resulting
  // polygon are represented by lists of their corners points with the points
  // arranged in a consecutive counterclockwise order
  template <typename Polygon>
  static std::vector<Point> sutherlandHodgman(const Polygon& subjectPolygon,
                                              const Polygon& clipPolygonInput) {
    std::vector<Point> outputPolygon(subjectPolygon.begin(),
                                     subjectPolygon.end());
    std::vector<Point> clipPolygon(clipPolygonInput.begin(),
                                   clipPolygonInput.end());
 
    const int clipSize = clipPolygon.size();
 
    // iterate over the edges of the clipping polygon represented by
    // consecutive points (indices clipIdx and clipIdxNext) in clipPolygon
    for (int clipIdx = 0; clipIdx < clipSize; clipIdx++) {
      const int clipIdxNext = (clipIdx + 1) % clipSize;
 
      std::vector<Point> inputPolygon(outputPolygon);
      outputPolygon.clear();
      const int inputSize = inputPolygon.size();
 
      // iterate over the edges of the subject polygon represented by
      // consectutive points (indices subjIdx and subjIdxNext) and clip each
      // subject edge with respect to current clipping edge (indices clipIdx
      // and clipIdxNext)
      for (int inputIdx = 0; inputIdx < inputSize; inputIdx++) {
        const int inputIdxNext = (inputIdx + 1) % inputSize;
 
        // Case differentation: Determine the positions of the corner points
        // of the subject edge p1, p2 with respect to the current clipping
        // edge given as the line from q1 to q2.
        // This is done by evaluating the cross product
        // d := (q1 - q2) x (p_i - q_1) = [(q1 - q2) x p_i] + q2 x q1
        // for i = 1, 2.
        // d < 0 => p_i inside, d > 0 => p_i outside,
        // d = 0 => p_i on clipping edge
 
        Point deltaQ = clipPolygon[clipIdxNext] - clipPolygon[clipIdx];
 
        Scalar q2_cross_q1 =
            clipPolygon[clipIdx][1] * clipPolygon[clipIdxNext][0] -
            clipPolygon[clipIdx][0] * clipPolygon[clipIdxNext][1];
 
        Scalar firstPos = -deltaQ[0] * inputPolygon[inputIdx][1] +
                          deltaQ[1] * inputPolygon[inputIdx][0] + q2_cross_q1;
        Scalar secondPos = -deltaQ[0] * inputPolygon[inputIdxNext][1] +
                           deltaQ[1] * inputPolygon[inputIdxNext][0] +
                           q2_cross_q1;
 
        static constexpr double EPSILON = 1e-14;
        // case 1: first point outside, second point inside
        if (firstPos > EPSILON && secondPos < -EPSILON) {
          // add intersection point and second point to output polygon
          outputPolygon.push_back(lineIntersectionPoint(
              inputPolygon[inputIdx], inputPolygon[inputIdxNext],
              clipPolygon[clipIdx], clipPolygon[clipIdxNext]));
          outputPolygon.push_back(inputPolygon[inputIdxNext]);
        }
 
        // case 2: first point inside, second point outside
        else if (firstPos < -EPSILON && secondPos > EPSILON) {
          // add intersection point to output polygon
          outputPolygon.push_back(lineIntersectionPoint(
              inputPolygon[inputIdx], inputPolygon[inputIdxNext],
              clipPolygon[clipIdx], clipPolygon[clipIdxNext]));
        }
 
        // case 3: first point outside or on the edge, scond point outside
        else if (firstPos >= -EPSILON && secondPos > EPSILON) {
          // do nothing
          continue;
        }
 
        // case 4: first point inside, second point inside or on the edge
        else {
          // add second point to output polygon
          outputPolygon.push_back(inputPolygon[inputIdxNext]);
        }
      }
    }
 
    return outputPolygon;
  }
 
  // computes the point of intersection in 2d between the line given by the
  // points p1, p2 and the line defined by the points q1, q2
  static Point lineIntersectionPoint(const Point& p1, const Point& p2,
                                     const Point& q1, const Point& q2) {
    Point deltaP = p1 - p2;
    Point deltaQ = q1 - q2;
 
    Scalar scaleFactor = deltaP[0] * deltaQ[1] - deltaQ[0] * deltaP[1];
 
    deltaP *= q1[0] * q2[1] - q1[1] * q2[0];
    deltaQ *= p1[0] * p2[1] - p1[1] * p2[0];
 
    Point intersection = deltaQ - deltaP;
    intersection /= scaleFactor;
 
    return intersection;
  }
 
  // computes the area of a convex polygon by evaluating the shoelace formula,
  // the polygon is represented by a corner point list with the points arranged
  // consecutively in counterclockwise order (for clockwise order (-1)*area
  // will be returned)
  template <typename Points>
  static Scalar polygonArea(const Points& polygonPoints) {
    const int polygonSize = polygonPoints.size();
 
    if (polygonSize < 3) {
      return 0.0;
    }
 
    Scalar area(0.0);
 
    for (int i = 0; i < polygonSize; i++) {
      const int j = (i + 1) % polygonSize;
      area += polygonPoints[i][0] * polygonPoints[j][1];
      area -= polygonPoints[j][0] * polygonPoints[i][1];
    }
 
    return 0.5 * area;
  }
 
#if CGAL_INTERSECTION
  template <typename Triang>
  static Scalar intersectionVolumeCGAL(const Triang& subjectTriangle,
                                       const Triang& clipTriangle) {
    using Epeck = CGAL::Exact_predicates_exact_constructions_kernel;
    using PointCGAL = CGAL::Point_2<Epeck>;
    using Triangle = CGAL::Triangle_2<Epeck>;
 
    std::array<PointCGAL, 3> v;
    for (int i = 0; i < 3; ++i) {
      const auto& p = makePoint(subjectTriangle[i]);
      v[i] = PointCGAL(p.x(), p.y());
    }
    const Triangle triangle1(v[0], v[1], v[2]);
 
    for (int i = 0; i < 3; ++i) {
      const auto& p = makePoint(clipTriangle[i]);
      v[i] = PointCGAL(p.x(), p.y());
    }
    const Triangle triangle2(v[0], v[1], v[2]);
 
    const auto result = CGAL::intersection(triangle1, triangle2);
 
    using PointList = std::vector<PointCGAL>;
    using Polygon = CGAL::Polygon_2<Epeck>;
    if (result) {
      // intersection is triangle
      if (const Triangle* t = boost::get<Triangle>(&*result))
        return std::abs(CGAL::to_double(
            CGAL::area(t->vertex(0), t->vertex(1), t->vertex(2))));
 
      // intersection is polygon
      else if (const PointList* points = boost::get<PointList>(&*result)) {
        Polygon polygon(points->begin(), points->end());
        return std::abs(CGAL::to_double(polygon.area()));
      }
    }
 
    // there is no intersection
    return 0;
  }
#endif
};
}  // namespace Dune
 
#endif