Are there MATLAB functions or scripts to perform boolean operations on triangulated solids?

It is a very specific question. There is no need for more details.

Answer 1

[DGM]

If your solids can each be expressed as a manifold STL, then I don't see why you can't use external tools to do the job. I normally use OpenSCAD to do CSG, but I'm pretty sure there are command-line tools that can also do basic operations. OpenSCAD can be used from the command line as well, but that's not how I use it. If you want the operation to be part of a MATLAB script, just use system() to deal with the external tools.

For example, this is the difference of two unit cubes, each composed of a minimal 12 triangles.

% a triangulated unit cube
Q.vertices = [0 0 0; 1 0 0; 1 1 0; 0 1 0; ...
              0 0 1; 1 0 1; 1 1 1; 0 1 1]; 
Q.faces = [3 2 1; 1 4 3; 5 6 7; 7 8 5; 1 2 6; 6 5 1; ...
           3 4 8; 8 7 3; 4 1 5; 5 8 4; 2 3 6; 3 7 6];
TR = triangulation(Q.faces,Q.vertices);
stlwrite(TR,'cube1.stl')

% the same thing, but translated by [1 -1 1]*0.5
Q.vertices = Q.vertices + [1 -1 1]*0.5;
TR = triangulation(Q.faces,Q.vertices);
stlwrite(TR,'cube2.stl')

% create a simple SCAD model 
% this could also just be written externally and reused
% writing scripts via fprintf() is ridiculously cumbersome
% though it demonstrates that names, etc. can be programmatically incorporated
fid = fopen('mymodel.scad','w');
fprintf(fid,'difference(){\n\timport("cube1.stl");\n\timport("cube2.stl");\n}');
fclose(fid);

% compute the composite geometry, output as a new STL
[status,result] = system('openscad -o output.stl mymodel.scad');

% read the new STL as a triangulation object and display it
TRnew = stlread('output.stl');
hs = trisurf(TRnew,'facecolor',[1 1 1]*0.7)
camlight
axis equal
view(-22,35)

If you don't want to use external tools, supposedly GPT has CSG tools (e.g. mesh_boolean()), though I've never used them. You will need to be able to compile the mex tools in order to use CGAL stuff.

https://www.mathworks.com/matlabcentral/fileexchange/49692-gptoolbox

That would arguably be a more direct way of doing it, but in the end, you're using the same external library (CGAL) either way.

See also:

https://www.mathworks.com/matlabcentral/answers/462934-make-valid-stl-file-from-overlapping-spheres

https://www.mathworks.com/matlabcentral/answers/561428-how-to-extract-intersection-between-two-3d-shapes

https://www.mathworks.com/matlabcentral/answers/349939-need-help-with-mesh-boolean-in-matlab

[Jahir Pabon]

Thank you. This is very informative.

[Paul]

@DGM

I've never seen this type of thing before. Can you explain what the "difference of two unit cubes" means?

Here are the cubes in question

Q.vertices = [0 0 0; 1 0 0; 1 1 0; 0 1 0; ...
              0 0 1; 1 0 1; 1 1 1; 0 1 1]; 
Q.faces = [3 2 1; 1 4 3; 5 6 7; 7 8 5; 1 2 6; 6 5 1; ...
           3 4 8; 8 7 3; 4 1 5; 5 8 4; 2 3 6; 3 7 6];
TR1 = triangulation(Q.faces,Q.vertices);
% the same thing, but translated by [1 -1 1]*0.5
Q.vertices = Q.vertices + [1 -1 1]*0.5;
TR2 = triangulation(Q.faces,Q.vertices);
figure
hold on
trisurf(TR1)
trisurf(TR2);
view(3)

Based on your result, it looks like the difference, D = Cube1 - Cube2, removes the partial volume of Cube1 (green and dark blue) that is overlapped by the partial volume of Cube2 (yellow and light blue), so that that portion of the resulting difference (D) has 0 volume. But what about the partial volume of Cube2 that doesn't overlap any part of Cube1? That portion of Cube2 would become negative volume, for lack of a better term, in the resulting difference. Does that negative volume just get tossed away in this type of math?

[DGM]

In OpenSCAD, difference() takes the first child object and removes its intersection with each of the remaining child objects in the scope. In other words, it's the set difference -- the volume in A, and not in B, C, etc. Forgive the wrong syntax highlighting, but:

difference(){
    // take this cube
    cube([1,1,1]);
    
    // and remove the part which intersects with this cube
    translate([1,-1,1]*0.5)
        cube([1,1,1]);
    
    // and the part which intersects with this cylinder
    translate([0,1,0.75])
        cylinder(r=0.5,h=1,$fn=100);
}

https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/CSG_Modelling#difference

If you wanted to do the XOR (symmetric difference) of two cubes, you could, but it's a bit of extra work (and showing it is another step).

difference(){
    // the cavity created between these objects isn't visible
    objxor(){
        // take both of these cubes and remove their intersection
        cube([1,1,1]);
        
        translate([1,-1,1]*0.5)
            cube([1,1,1]);
    }
    
    // so cut a section through both of them using another rotated cube
    translate([0,1,-1])
        rotate([0,0,-135])
        cube([1,1,1]*3);
}

// there isn't a built-in object-level xor(),
// but it's common to make one.
module objxor(){
    difference(){
        union(){
            children(0);
            children(1);
        };
        intersection(){
            children(0);
            children(1);
        }
    }
}

The complication with using XOR is usually that it is going to create non-manifold edges where the remaining volumes touch. I'm not sure what the proper term for that form is. I guess it's a type of T-connected face defect. OpenSCAD will normally throw warnings during render if you're creating non-manifold geometry, but not in command-line mode.

EDIT: I was wrong. It does produce the non-manifold warnings in command-line mode, but the way I'm using system() doesn't display them.

Answer 2

[Matt J]

If you mean you want to take intersections and unions of meshes, then no, there are no MathWorks-authored commands to do that. However, there are FEX files that let you convert a mesh to a binary volume and back, e.g.,

https://www.mathworks.com/matlabcentral/fileexchange/27390-mesh-voxelisation?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/fileexchange/27733-converting-a-3d-logical-array-into-an-stl-surface-mesh?s_tid=srchtitle

so you might be able to do it by transfering to the binary image domain if you can tolerate some discretization error.

[Jahir Pabon]

Thank you Matt. What you propose, though a bit cumbersome, sounds like a possible alternative. I would like to preserve the simplicity of the original triangulations. But if there is nothing else in the Matlab world I will give that a try.

Answer 3

[John D'Errico]

Are there such codes in MATLAB? Yes. I wrote them. intersect, union. In 2-d and in 3-d. (Not sure if I wrote a setdiff variant. Probably not, as the problem I was working on would not need them.) Are they publically available? No. Sorry. The codes belong to my client. I'm not sure they are still in use, but I don't own them.

You want them to retain the original "simplicity" of the triangulations? Sorry, but no. The result tends to get quite complicated near the boundaries.

Could you write something? They were not that terrible to implement. Kind of slow if the tessellations were large, but that always happens. You can start with these operations on a pair of tetrahedra.

[Jahir Pabon]

Thank you John. Glad to hear hyou have done it. Not so glad to hear somebody else owns the work now.

By "preserving the simplicity" I meant, not ending up with a large set of triangles, which is what the Voxelization will return. I understand that some of the original triangles will be split into smaller ones in the process of computing intersections. That is expected.