How to cut a 3D object using a plane and deleting the part above the plane, but creating new faces/vertices based on the points to ensure that the geometry below the cut is maintained?

I am working with an stl file (faces, vertices, normals, etc) and have been able to establish a plane, find the coordinates of the points where the plane intersects the edges of the 3D objects, but am not able to modify the the stl file in such a way that I keep all the geometry intact below the plane cut but remove the geometry above the plane cut.

More clarification using an analogy: Imagine you had a donut in real life that you sliced using a knife to create 2 crescent shaped halves. Now I am trying to solve the same problem using matlab and an stl file. I have the donut, I have the knife (the plane), I know all the new vertices that should be formed when the knife cuts, I am just not able to cut it in such a way that the solid part of the donut is still present and the hole (well a semicircular hole after cutting) is also present.

Is there someone who has expertise in this area and could help me out please :) Thank you!

[Matt J]

Could you attach an example set of face/vertex data and the cutting plane equation?

[Anshuman Agrawal]

The plane is just defined as a row vector of a point, and 2 vectors. After that I found the edges and found the points on these edges/lines that intersected the plane. I have attached the data. Thank you!

[John D'Errico]

Note that you don't need TWO vectors to define a plane. In fact, that makes it more complex.

A plane is defined by a point (P) in the plane, and the normal vector (N) to the plane. Thus any point (X) in the plane satisfies the relation:

dot(X - P,N) = 0

The sign of that dot product, if the point is not in the plane is positive or negative. That tells you which side of the curtting plane the point lies on. Note that a dot product is most simply done using the * operator. This allows you to dot two vectors together, or a vector with an entire array of points.

For example:

XYZ = rand(100,3);
P = [.6 .6 .6];
N = [1 2 3];
N = N/norm(N); % normalized to unit length

So given a set of points in XYZ, we can tell if a point is in the plane, or off the plane, and which side of the plane the points lie on, simply by checking the sign of those dot products.

u = (XYZ - P)*N';

Note that the above uses a capability of MATLAB R2016b or later.

[Anshuman Agrawal]

Thank you for helping me! Greatly appreciate the help!

[Ismaeel Zaman]

Hi, would you be able to share how you obtained the intersection points from an STL file and a plane?

Answer 1

[John D'Errico]

Pretty simple really. Just the execution that takes some effort.

1. Identify any simplex that has all vertices above or below the cutting plane. Keep those uncut that are in the region of interest. Discard those that are all outside. The test to see if a vertex is above or below the cutting plane is as simple as a dot product, so just a matrix*vector multiply for all vertices at once.
2. All simplexes that remain have at least one vertex above, and one below the cutting plane. Operate on these simplexes one at a time.
3. Retain any vertices that are entirely in the region of interest, or on the plane itself. For any edges of the simplex, locate all edges that cross the cutting plane. Using simple linear interpolation, find the point of crossing.
4. Combine the set of vertices that are below the cutting plane with those that lie in the cutting plane from this simplex.
5. Tessellate that set, append the set of simplexes created to those that lie entirely below the cutting plane.

This operation is not that difficult. As I said, it just takes some code.

[Matt J]

3. Retain any vertices that are entirely in the region of interest, or on the plane itself. For any edges of the simplex, locate all edges that cross the cutting plane. Using simple linear interpolation, find the point of crossing.

My intersectionHull utility (Download), may simplify this step. Given the vertices, V, of any particular simplex from step 2, just do

I=intersectionHull('vert',V,'lcon',A,b)

where A*x<=b is the half-space, bounded by the cutting plane, that you want to keep. The vertices of the intersection of the half-space and the simplex will be given in I.vert which you can tesselate directly.

[John D'Errico]

Actually, the simplicial complex toolbox I wrote many years ago does it all. It has a planar truncation utility that slices through a simplicial complex in any number of dimensions, returning the truncated complex that lives on the designated side of a slicing plane. I never posted it on the FEX, because people don't seem to be able to use it without more coaching from me than I want to give.

xyz = rand(1000,3);
k = sqrt(sum((xyz - 0.5).^2,2)) < 0.5;
xyz = xyz(k,:);
sc = delaunays(xyz)
sc
sc = 
  struct with fields:

          domain: [553×3 double]
    tessellation: [3288×4 double]
           range: [553×0 double]
           
plotsc(sc)
axis equal
box on
grid on

So a quasi-spherical blob. I'm feeling too lazy to create a toroidal thing using the tube utility I have written, or perhaps by tessellating a torus in a cylindrical coordinate system, then converting to cartesian coordinates.

sc2 = planartruncate(sc,[.6 .6 .6],[1 2 3],'domain','below');
plotsc(sc2)
axis equal
box on
grid on

I've rotated it to show the cutting plane.

A different rotation shows the triangulation as it happens in the cutting plane.

As I said, all very doable. But until you start to understand how tessellations work, and how to work with them, it seems hard.

[Anshuman Agrawal]

Thank you so much for your time and help! This is pretty much exactly what I was trying to do. I assume the functions delaunays and planartruncate are part of your toolbox, hence was wondering if you could send me the toolbox since you have shown me what to do with it.

Thank you once again!

[hagege ruben]

I have a similar problem and would like to have this toolbox too

[Anu_M]

Hello,

Is your toolbox still available?

[Ismaeel Zaman]

Hi I would also like to use this tool box. Is there a download anywhere?