intersect_clipper.h

Go to the documentation of this file.

/*
This file is part of the Ristra portage project.
Please see the license file at the root of this repository, or at:
    https://github.com/laristra/portage/blob/master/LICENSE
*/

#ifndef PORTAGE_INTERSECT_INTERSECTCLIPPER_H_
#define PORTAGE_INTERSECT_INTERSECTCLIPPER_H_

#include <iostream>
#include <cmath>
#include <cfloat>
#include <algorithm>

// portage includes
#include "portage/intersect/clipper.hpp"
#include "portage/support/portage.h"
#include "wonton/support/Point.h"

namespace Portage {

inline
std::vector<double> areaAndMomentPolygon(std::vector<Wonton::Point<2>> const& poly){
  double area = 0;
  double cx = 0;
  double cy = 0;
  std::vector <double> ret;

  int lastIndex = poly.size()-1;

  for(int i=0; i < lastIndex; i++){
    double a = (poly[i][0]*poly[i+1][1]-poly[i+1][0]*poly[i][1]);
    area+=a;
    cx+= (poly[i][0]+poly[i+1][0])*a;
    cy+= (poly[i][1]+poly[i+1][1])*a;
  }

  //close the polygon
  double a = poly[lastIndex][0]*poly[0][1] - poly[0][0]*poly[lastIndex][1];
  area+= a;
  cx+= (poly[lastIndex][0]+poly[0][0])*a;
  cy+= (poly[lastIndex][1]+poly[0][1])*a;
  ret.emplace_back(.5*area);
  ret.emplace_back(cx/6.);
  ret.emplace_back(cy/6.);
  return ret;
}

template <typename SourceMeshType, typename TargetMeshType=SourceMeshType> class IntersectClipper
{

public:

typedef std::vector<Wonton::Point<2>> Poly;
typedef std::pair<double, Wonton::Point<2>> Moment;

IntersectClipper(const SourceMeshType &s, const TargetMeshType &t):sourceMeshWrapper(s), targetMeshWrapper(t){}

std::vector<std::vector<double> > operator() (const int cellA, const int cellB) const {
  Poly polyA, polyB;
  sourceMeshWrapper.cell_get_coordinates(cellA, &polyA);
  targetMeshWrapper.cell_get_coordinates(cellB, &polyB);
  double max_size_poly = 0;
  max_size_poly = IntersectClipper::updateMaxSize(polyA, max_size_poly);
  max_size_poly = IntersectClipper::updateMaxSize(polyB, max_size_poly);

  const ClipperLib::Path intPolyA = IntersectClipper::convertPoly2int(polyA, max_size_poly);
  const ClipperLib::Path intPolyB = IntersectClipper::convertPoly2int(polyB, max_size_poly);

  //Make clipper aware of polyA and polyB
  ClipperLib::Clipper clpr;
  clpr.AddPath(intPolyA, ClipperLib::ptSubject, true);
  clpr.AddPath(intPolyB, ClipperLib::ptClip, true);

  //Compute the intersection of all paths which have been added to clipper (in this case
  //polyA and polyB; use the Even-Odd winding rule in case of self-intersecting polygons
  //For more information on the winding rules see:
  //http://www.angusj.com/delphi/clipper/documentation/Docs/Units/ClipperLib/Types/PolyFillType.htm
  ClipperLib::Paths solution;
  clpr.Execute(ClipperLib::ctIntersection, solution,
               ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);

  std::vector<std::vector<Wonton::Point<2>>> intersectionList;
  for(auto const &i: solution){
    std::vector<Wonton::Point<2>> poly;
    for(auto const &j: i){
      //Build list of intersections and convert back to doubles
      poly.emplace_back(integer2real(j.X, max_size_poly), integer2real(j.Y, max_size_poly));
    }
    intersectionList.emplace_back(poly);
  }
  std::vector<std::vector<double>> moments;
  for(auto const &i: intersectionList){
    moments.emplace_back(areaAndMomentPolygon(i));
  }
  return moments;
}

IntersectClipper() = default;

IntersectClipper(const IntersectClipper &) = delete;

IntersectClipper & operator = (const IntersectClipper &) = delete;


private:

//We must use the same max size for all the polygons, so the number we are looking for is the maximum value in the set--all the X and Y values will be converted using this value
static double updateMaxSize( const std::vector<Wonton::Point<2>> poly, double max_size_poly){
  for(auto const &i: poly){
    double m = std::max(std::abs(i[0]), std::abs(i[1]));
    if (m > max_size_poly) max_size_poly = m;
  }
  return max_size_poly;
}

static long real2integer(double a, const double max_size){
  int exp;

  // We want to move the decimal point of the floating point number so that the last digit given is before the decimal point.
  // The alogrithm is to find the exponent (in the floating point representation) of the largest number in the polygons (either X or Y coordinate).
  //Multiply by the power of two which is number of digits in mantissa (in this case, 53) - exponent
  //lrint, then rounds to an integer

  //Compute the exponent
  frexp(max_size, &exp);
  //Size of long must be >= size double
  return lrint(ldexp(a, DBL_MANT_DIG-exp));
}

static double integer2real(long a, const double max_size){
  int exp;
  frexp(max_size, &exp);
  return ldexp(a, exp-DBL_MANT_DIG);
}

//Convert an entire polygon (specifiied as a std::vector<Wonton::Point>) to a std::vector<IntPoint>
static std::vector<ClipperLib::IntPoint> convertPoly2int(std::vector<Wonton::Point<2>> poly, double max_size_poly){
  std::vector<ClipperLib::IntPoint> intpoly(poly.size());
  std::transform(poly.begin(), poly.end(), intpoly.begin(), [max_size_poly](Wonton::Point<2> point){
      return ClipperLib::IntPoint(real2integer(point[0], max_size_poly), real2integer(point[1], max_size_poly));});
  return intpoly;
}

private:
const SourceMeshType &sourceMeshWrapper;
const TargetMeshType &targetMeshWrapper;

};  // class IntersectClipper

}  // namespace Portage

#endif  // PORTAGE_INTERSECT_INTERSECTCLIPPER_H_