src/Utilities/pointInPolyhedron.cc
//------------------------------------------------------------------------------
// The method in question -- test if a point is on the interior of a polyhedron.
//------------------------------------------------------------------------------
The method in question -- test if a point is on the interior of a polyhedron.
The method in question -- test if a point is on the interior of a polyhedron .
bool pointInPolyhedron(const Dim<3>::Vector& p,
                       const Dim<3>::FacetedVolume& polyhedron,
                       const bool countBoundary,
                       const double tol) {
  typedef Dim<3>::Vector Vector;
  typedef Dim<2>::Vector Vector2d;
  const Vector rayhat(0, 0, 1);    
  const auto pend = p + Vector(0, 0, 1e60);
  const auto& vertices = polyhedron.vertices();
  const auto& facets = polyhedron.facets();
  const auto  nfacets = facets.size();
  const Vector2d pxy(p.x(), p.y()),
                 pzx(p.z(), p.x()),
                 pyz(p.y(), p.z());
  if (not testPointInBox(p, polyhedron.xmin(), polyhedron.xmax(), tol)) return false;
  if (countBoundary and polyhedron.distance(p) < tol) return true;
  int numxy = 0;
#pragma omp parallel
  {
    int numxy_thread = 0;
#pragma omp for
    for (auto ifacet = 0u; ifacet < nfacets; ++ifacet) {
      const auto& facet = facets[ifacet];
      const auto& ipts = facet.ipoints();
      const auto  npts = ipts.size();
      const auto& fhat = facet.normal();
      const auto  fdotr = fhat.dot(rayhat);
      if (std::abs(fdotr) <= tol) {
        if (std::abs((vertices[ipts[0]] - p).dot(fhat)) <= tol) {
          size_t i = 0;
          while (i < npts and
                 segmentSegmentDistance(vertices[ipts[i]], vertices[ipts[(i + 1) % npts]], p, pend) > tol) ++i;
          if (i < npts) ++numxy_thread;
        }
      } else {
        const auto si = fhat.dot(vertices[ipts[0]] - p)*safeInvVar(fdotr);
        CHECK2(std::abs((p + si*rayhat - vertices[ipts[0]]).dot(fhat)) < tol*std::max(1.0, 1e-3*(p - vertices[ipts[0]]).magnitude()),
               "Point in plane?  " << p << " " << vertices[ipts[0]] << " " << rayhat << " " << si << " : " << std::abs((p + si*rayhat - vertices[ipts[0]]).dot(fhat)) << " !< " << tol);
        if (si > 0.0) {
          vector<Vector2d> vxy;
          for (size_t i = 0; i < npts; ++i) {
            vxy.push_back(Vector2d(vertices[ipts[i]].x(), vertices[ipts[i]].y()));
          }
          vxy.push_back(vxy[0]);
          if (pointInPolygon(pxy, vxy)) {
            ++numxy_thread;
          }
        }
      }
    }
#pragma omp critical
    {
      numxy += numxy_thread;
    }
  }
  return ((numxy % 2) == 1);
}
bool pointInPolyhedron ( const Dim < 3 > :: Vector & p , const Dim < 3 > :: FacetedVolume & polyhedron , const bool countBoundary , const double tol ) { typedef Dim < 3 > :: Vector Vector ; typedef Dim < 2 > :: Vector Vector2d ; const Vector rayhat ( 0 , 0 , 1 ) ; const auto pend = p + Vector ( 0 , 0 , 1e60 ) ; const auto & vertices = polyhedron . vertices ( ) ; const auto & facets = polyhedron . facets ( ) ; const auto nfacets = facets . size ( ) ; const Vector2d pxy ( p . x ( ) , p . y ( ) ) , pzx ( p . z ( ) , p . x ( ) ) , pyz ( p . y ( ) , p . z ( ) ) ; if ( not testPointInBox ( p , polyhedron . xmin ( ) , polyhedron . xmax ( ) , tol ) ) return false ; if ( countBoundary and polyhedron . distance ( p ) < tol ) return true ; int numxy = 0 ; #pragma  omp parallel { int numxy_thread = 0 ; #pragma  omp for for ( auto ifacet = 0u ; ifacet < nfacets ; ++ ifacet ) { const auto & facet = facets [ ifacet ] ; const auto & ipts = facet . ipoints ( ) ; const auto npts = ipts . size ( ) ; const auto & fhat = facet . normal ( ) ; const auto fdotr = fhat . dot ( rayhat ) ; if ( std :: abs ( fdotr ) <= tol ) { if ( std :: abs ( ( vertices [ ipts [ 0 ] ] - p ) . dot ( fhat ) ) <= tol ) { size_t i = 0 ; while ( i < npts and segmentSegmentDistance ( vertices [ ipts [ i ] ] , vertices [ ipts [ ( i + 1 ) % npts ] ] , p , pend ) > tol ) ++ i ; if ( i < npts ) ++ numxy_thread ; } } else { const auto si = fhat . dot ( vertices [ ipts [ 0 ] ] - p ) * safeInvVar ( fdotr ) ; CHECK2 ( std :: abs ( ( p + si * rayhat - vertices [ ipts [ 0 ] ] ) . dot ( fhat ) ) < tol * std :: max ( 1.0 , 1e-3 * ( p - vertices [ ipts [ 0 ] ] ) . magnitude ( ) ) , " " << p << " " << vertices [ ipts [ 0 ] ] << " " << rayhat << " " << si << " " << std :: abs ( ( p + si * rayhat - vertices [ ipts [ 0 ] ] ) . dot ( fhat ) ) << " " << tol ) ; if ( si > 0.0 ) { vector < Vector2d > vxy ; for ( size_t i = 0 ; i < npts ; ++ i ) { vxy . push_back ( Vector2d ( vertices [ ipts [ i ] ] . x ( ) , vertices [ ipts [ i ] ] . y ( ) ) ) ; } vxy . push_back ( vxy [ 0 ] ) ; if ( pointInPolygon ( pxy , vxy ) ) { ++ numxy_thread ; } } } } #pragma  omp critical { numxy += numxy_thread ; } } return ( ( numxy % 2 ) == 1 ) ; }
The method in question -- test if a point is on the interior of a polyhedron.
The method in question -- test if a point is on the interior of a polyhedron .
// Project straight up in z. // First reject any points outside the bounding box. // Are we checking the surface? // if (countBoundary and pointOnPolyhedron(p, polyhedron, tol)) return true; // Keep count of how many facets we intersect. // Walk the facets. // Check if the facet is parallel to ray. // If they are, we only care if the ray is in the plane of the facet. // The ray and plane are parallel, and the ray is in the plane. // Check if the ray crosses any of the facet lines -- if so, counts as an intersection. // Find the intersection point of the line corresponding to our ray and the plane of the facet. // What we're solving for here is parametric s where the line is represented as //   pline = p + s*rayhat // We only need to check this facet if si is positive, indicating the the intersection // point is in the positive ray direction, i.e., above p in z. // Copy the vertices as XY pairs. // Does a ray straight up the +z direction intersect this projected polygon? // cerr << "Intersect facet: " << si << " " << pxy << endl; // for (size_t i = 0; i < npts; ++i) cerr << " " << vxy[i]; // cerr << endl; // If the ray passed through an even number of intersections the point is external to the  // polyhedron.