src/Utilities/pointInPolygon.cc
//------------------------------------------------------------------------------
// Test a polygon in 3-D.
// This version allows you to pass in an indirect addressing definition of the
// polygon vertices.
// Note -- does not currently work for polygons with more than one loop!
//------------------------------------------------------------------------------
Test a polygon in 3-D.
This version allows you to pass in an indirect addressing definition of the
polygon vertices.
Note -- does not currently work for polygons with more than one loop!
Test a polygon in 3 - D . This version allows you to pass in an indirect addressing definition of the polygon vertices . Note -- does not currently work for polygons with more than one loop!
bool pointInPolygon(const Dim<3>::Vector& p,
                    const vector<Dim<3>::Vector>& vertices,
                    const vector<unsigned>& ipoints,
                    const Dim<3>::Vector& normal,
                    const bool countBoundary,
                    const double tol) {
  const auto npts = ipoints.size();
  unsigned i, j, ik, jk;
  bool result = false;
  const double px = p.x();
  const double py = p.y();
  const double pz = p.z();
  double fxmin = std::numeric_limits<double>::max(),    fymin = std::numeric_limits<double>::max(),   fzmin =  std::numeric_limits<double>::max();
  double fxmax = std::numeric_limits<double>::lowest(), fymax = std::numeric_limits<double>::lowest(), fzmax = std::numeric_limits<double>::lowest();
  for (const auto i: ipoints) {
    fxmin = std::min(fxmin, vertices[i].x());
    fymin = std::min(fymin, vertices[i].y());
    fzmin = std::min(fzmin, vertices[i].z());
    fxmax = std::max(fxmax, vertices[i].x());
    fymax = std::max(fymax, vertices[i].y());
    fzmax = std::max(fzmax, vertices[i].z());
  }
  const auto fuzz = std::max(1.0e-5, tol)*std::max(fxmax - fxmin, std::max(fymax - fymin, fzmax - fzmin));
  fxmin -= fuzz;
  fymin -= fuzz;
  fzmin -= fuzz;
  fxmax += fuzz;
  fymax += fuzz;
  fzmax += fuzz;
  if (px >= fxmin and px <= fxmax and
      py >= fymin and py <= fymax and
      pz >= fzmin and pz <= fzmax) {
    if (pointOnPolygon(p, vertices, ipoints, tol)) return countBoundary;
    const double nmax = normal.maxAbsElement();
    if (std::abs(normal.x()) > 0.9*nmax) {
      for (ik = 0, jk = npts - 1; ik < npts; jk = ik++) {
        i = ipoints[ik];
        j = ipoints[jk];
        if ( ((vertices[i].z() > pz) != (vertices[j].z() > pz)) &&
             (py < (vertices[j].y() - vertices[i].y()) * (pz - vertices[i].z()) / (vertices[j].z() - vertices[i].z()) + vertices[i].y()) )
          result = not result;
      }
    } else if (std::abs(normal.y()) > 0.9*nmax) {
      for (ik = 0, jk = npts - 1; ik < npts; jk = ik++) {
        i = ipoints[ik];
        j = ipoints[jk];
        if ( ((vertices[i].x() > px) != (vertices[j].x() > px)) &&
             (pz < (vertices[j].z() - vertices[i].z()) * (px - vertices[i].x()) / (vertices[j].x() - vertices[i].x()) + vertices[i].z()) )
          result = not result;
      }
    } else {
      for (ik = 0, jk = npts - 1; ik < npts; jk = ik++) {
        i = ipoints[ik];
        j = ipoints[jk];
        if ( ((vertices[i].y() > py) != (vertices[j].y() > py)) &&
             (px < (vertices[j].x() - vertices[i].x()) * (py - vertices[i].y()) / (vertices[j].y() - vertices[i].y()) + vertices[i].x()) )
          result = not result;
      }
    }
  }
  return result;
}
bool pointInPolygon ( const Dim < 3 > :: Vector & p , const vector < Dim < 3 > :: Vector > & vertices , const vector < unsigned > & ipoints , const Dim < 3 > :: Vector & normal , const bool countBoundary , const double tol ) { const auto npts = ipoints . size ( ) ; unsigned i , j , ik , jk ; bool result = false ; const double px = p . x ( ) ; const double py = p . y ( ) ; const double pz = p . z ( ) ; double fxmin = std :: numeric_limits < double > :: max ( ) , fymin = std :: numeric_limits < double > :: max ( ) , fzmin = std :: numeric_limits < double > :: max ( ) ; double fxmax = std :: numeric_limits < double > :: lowest ( ) , fymax = std :: numeric_limits < double > :: lowest ( ) , fzmax = std :: numeric_limits < double > :: lowest ( ) ; for ( const auto i : ipoints ) { fxmin = std :: min ( fxmin , vertices [ i ] . x ( ) ) ; fymin = std :: min ( fymin , vertices [ i ] . y ( ) ) ; fzmin = std :: min ( fzmin , vertices [ i ] . z ( ) ) ; fxmax = std :: max ( fxmax , vertices [ i ] . x ( ) ) ; fymax = std :: max ( fymax , vertices [ i ] . y ( ) ) ; fzmax = std :: max ( fzmax , vertices [ i ] . z ( ) ) ; } const auto fuzz = std :: max ( 1.0e-5 , tol ) * std :: max ( fxmax - fxmin , std :: max ( fymax - fymin , fzmax - fzmin ) ) ; fxmin -= fuzz ; fymin -= fuzz ; fzmin -= fuzz ; fxmax += fuzz ; fymax += fuzz ; fzmax += fuzz ; if ( px >= fxmin and px <= fxmax and py >= fymin and py <= fymax and pz >= fzmin and pz <= fzmax ) { if ( pointOnPolygon ( p , vertices , ipoints , tol ) ) return countBoundary ; const double nmax = normal . maxAbsElement ( ) ; if ( std :: abs ( normal . x ( ) ) > 0.9 * nmax ) { for ( ik = 0 , jk = npts - 1 ; ik < npts ; jk = ik ++ ) { i = ipoints [ ik ] ; j = ipoints [ jk ] ; if ( ( ( vertices [ i ] . z ( ) > pz ) != ( vertices [ j ] . z ( ) > pz ) ) && ( py < ( vertices [ j ] . y ( ) - vertices [ i ] . y ( ) ) * ( pz - vertices [ i ] . z ( ) ) / ( vertices [ j ] . z ( ) - vertices [ i ] . z ( ) ) + vertices [ i ] . y ( ) ) ) result = not result ; } } else if ( std :: abs ( normal . y ( ) ) > 0.9 * nmax ) { for ( ik = 0 , jk = npts - 1 ; ik < npts ; jk = ik ++ ) { i = ipoints [ ik ] ; j = ipoints [ jk ] ; if ( ( ( vertices [ i ] . x ( ) > px ) != ( vertices [ j ] . x ( ) > px ) ) && ( pz < ( vertices [ j ] . z ( ) - vertices [ i ] . z ( ) ) * ( px - vertices [ i ] . x ( ) ) / ( vertices [ j ] . x ( ) - vertices [ i ] . x ( ) ) + vertices [ i ] . z ( ) ) ) result = not result ; } } else { for ( ik = 0 , jk = npts - 1 ; ik < npts ; jk = ik ++ ) { i = ipoints [ ik ] ; j = ipoints [ jk ] ; if ( ( ( vertices [ i ] . y ( ) > py ) != ( vertices [ j ] . y ( ) > py ) ) && ( px < ( vertices [ j ] . x ( ) - vertices [ i ] . x ( ) ) * ( py - vertices [ i ] . y ( ) ) / ( vertices [ j ] . y ( ) - vertices [ i ] . y ( ) ) + vertices [ i ] . x ( ) ) ) result = not result ; } } } return result ; }
Test a polygon in 3-D.
Test a polygon in 3 - D .
// typedef boost::tuple<double, double> BGPoint; // typedef boost::geometry::model::polygon<BGPoint> BGPolygon; // Prerequisites. // BEGIN_CONTRACT_SCOPE // { //   REQUIRE(ipoints.size() > 2); //   const double normmag = normal.magnitude(); //   Vector normi; //   for (i = 1; i != npts; ++i) { //     j = (i + 1) % npts; //     k = (i + 2) % npts; //     normi = (vertices[ipoints[j]] - vertices[ipoints[i]]).cross(vertices[ipoints[k]] - vertices[ipoints[i]]); //     REQUIRE2(fuzzyEqual(std::abs(normi.dot(normal)), normmag*normi.magnitude(), 1.0e-5), normi << " " << normal << " " << normi.dot(normal) << " " << normmag*normi.magnitude()); //   } //   // REQUIRE2(fuzzyEqual(pointPlaneDistance(p, vertices[ipoints[0]], normal.unitVector()), 0.0, 1.0e-3), pointPlaneDistance(p, vertices[ipoints[0]], normal.unitVector())); // } // END_CONTRACT_SCOPE // Find the bounding box for the polygon. // Check if the point is on the boundary (within tolerance). // Figure out which plane we're going to project to. // x plane -- use (y,z) coordinates. // vector<BGPoint> points; // for (i = 0; i != npts; ++i) points.push_back(BGPoint(vertices[ipoints[i]].y(), vertices[ipoints[i]].z())); // points.push_back(points[0]); // BGPolygon poly; // boost::geometry::append(poly, points); // result = boost::geometry::within(BGPoint(py, pz), poly); // y plane -- use (z,x) coordinates. // vector<BGPoint> points; // for (i = 0; i != npts; ++i) points.push_back(BGPoint(vertices[ipoints[i]].z(), vertices[ipoints[i]].x())); // points.push_back(points[0]); // BGPolygon poly; // boost::geometry::append(poly, points); // result = boost::geometry::within(BGPoint(pz, px), poly); // z plane -- use (x,y) coordinate. // vector<BGPoint> points; // for (i = 0; i != npts; ++i) points.push_back(BGPoint(vertices[ipoints[i]].x(), vertices[ipoints[i]].y())); // points.push_back(points[0]); // BGPolygon poly; // boost::geometry::append(poly, points); // result = boost::geometry::within(BGPoint(px, py), poly); // That's it.