Determine if 3d point belongs to polyhedron

I want to determine whether a 3d point belongs to a polyhedron. This polyhedron is a rectangular pyramid which extends out from the origin. I was thinking of doing by defining the normal to each side of the polyhedron and then taking the distance between each side and the point. But I'm kind of stuck what to do after that. I was thinking of this plus a combination of inequalities would lead me to an answer.

[Dyuman Joshi]

You can make an alphaShape of the polyhedron and use inShape to check whether the point is inside or not.

[John D'Errico]

This is a viable idea, but remember to supply an alpha of inf. An alpha shape when alpha is infinite is the same as a convex hull, so the prediction of whether a point lies inside the domain will be correct.

Given that, you can also just use other tools. For example, you can build a delaunayTriangulation of the pyramid. The boundary surface of a delaunay triangulation is simply the convex hull. Then you could use a tool like isInterior to learn if a point lies inside or not. This does not rely on a value of alpha for an alpha shape.

Finally, yes, you can convert the polyhedron into a set of inequalities. @Matt J has a tool on the FEX called vert2lcon. It will provide a set of inequalities you can use to perform the test. If ALL of the inequalities are satisfied, then the point lies inside the polyhedron.

[Dyuman Joshi]

@John D'Errico, the documentation of isInterior mentions only working with 2D triangulations.

How would that work here?

[John D'Errico]

Sorry. pointLocation is the tool, not isInterior.

[Dyuman Joshi]

Ah, neat. I didn't know about pointLocation.

Answer 1

[John D'Errico]

Since I'm going to show how to do it using a triangulation anyway, I might as well post it as an answer instead of a comment.

xyz = [0 0 0;1 0 0;0 1 0;1 1 0;.5 .5 1];
tri = delaunayTriangulation(xyz)

tri = 
  delaunayTriangulation with properties:

              Points: [5×3 double]
    ConnectivityList: [2×4 double]
         Constraints: []

help pointLocation

--- help for triangulation/pointLocation ---

  pointLocation  Triangle or tetrahedron containing specified point
  TI = pointLocation(T, QP)  returns the index of the triangle/tetrahedron
       enclosing the query point QP. The matrix QP contains the coordinates
       of the query points. QP is a mpts-by-ndim matrix where mpts is the 
       number of query points and 2 <= ndim <= 3.
       TI is a column vector of triangle or tetrahedron IDs corresponding to the 
       row numbers of the triangulation connectivity matrix T.ConnectivityList. 
       The triangle/tetrahedron enclosing the point QP(k,:) is TI(k). 
       Returns NaN for points not located in a triangle or tetrahedron of T.
  
       TI = pointLocation(T, QX,QY) and TI = pointLocation (T, QX,QY,QZ)
       allow the query points to be specified in alternative column vector 
       format when working in 2D and 3D.
  
       [TI, BC] = pointLocation(T,...) returns, in addition, the Barycentric
       coordinates BC. BC is a mpts-by-ndim matrix, each row BC(i,:) represents
       the Barycentric coordinates of QP(i,:) with respect to the enclosing TI(i).
 
     Example 1:
         % Point Location in 2D
         T = triangulation([1 2 4; 1 4 3; 2 4 5],[0 0; 2 0; 0 1; 1 1; 2 1]);
         % Find the triangle that contains the following query point
         QP = [1, 0.5];
         triplot(T,'-b*'), hold on
         plot(QP(:,1),QP(:,2),'ro')
         % The query point QP is located in a triangle with index TI = 1
         TI = pointLocation(T, QP)
 
     Example 2:
         % Point Location in 3D plus barycentric coordinate evaluation
         % for a Delaunay triangulation
         x = rand(10,1); y = rand(10,1); z = rand(10,1);
         DT = delaunayTriangulation(x,y,z);
         % Find the tetrahedra that contain the following query points
         QP = [0.25 0.25 0.25; 0.5 0.5 0.5]
         [TI, BC] = pointLocation(DT, QP)
 
     See also triangulation, triangulation.nearestNeighbor, delaunayTriangulation.

    Documentation for triangulation/pointLocation
       doc triangulation/pointLocation

    Other uses of pointLocation

       DelaunayTri/pointLocation

Now we can use that tool to identify a point as inside or out.

pointLocation(tri,[.5 .5 .5;1 2 3])

ans = 2×1

2 NaN

If it returns a number, then the point lies inside the triangulation. And it even tells you in which simplex the point falls. Note that in the case of a point on the surface, it should return a true result, but sometimes floating point trash can be an issue. If the point lies outside the triangulation, then you will get a NaN. So the test can be most simply written as

~isnan(pointLocation(tri,[.5 .5 .5;1 2 3]))

ans = 2×1 logical array
   1
   0

True is inside, false is outside.

Could you use an alpha shape? Again, yes. But you need to tell it the value of alpha.

S = alphaShape(xyz,inf)

S = 
  alphaShape with properties:

             Points: [5×3 double]
              Alpha: Inf
      HoleThreshold: 0
    RegionThreshold: 0

Alpha shapes have many nice capabilities.

methods(S)

Methods for class alphaShape:

alphaShape          alphaTriangulation  boundaryFacets      inShape             numRegions          plot                volume              
alphaSpectrum       area                criticalAlpha       nearestNeighbor     perimeter           surfaceArea         

One of them is plot.

plot(S)

But now you can also use the inShape tool to identify a point inside.

inShape(S,[.5 .5 .5;1 2 3])

ans = 2×1 logical array
   1
   0

Both approaches will suffice. An alpha shape does have one virtue, in that if the polyhedron is not convex, then it will still work, as long as you are able to choose an alpha that works. However, that can sometimes be problematic.