How to extract intersection between two 3D shapes

I have a two 3D shapes, one defined by a patch object, the other by a simple analytical function (in my case a sphere):

[verts, faces, cindex] = teapotGeometry;
sphere_origin = [1 1 1];
sphere_radius = 1;
figure;
p = patch('Faces',faces,'Vertices',verts,'FaceVertexCData',cindex,'FaceColor','interp');
sph = drawSphere([sphere_origin sphere_radius]); % using matGeom (https://github.com/mattools/matGeom)

I'd like a way to slice the patch into 2 or more sub-objects that contains the part of the teapot outside the sphere, and the one that is inside.

My main problems are that I don't know:

1. How to find the exact intersection between the sphere and the 3D shape (not just the closest vertice, but the extrapolated one). (DONE: using surfaceintersection as @darova suggested)
2. How to add the faces to the extrapolated vertices in a way that conserve the integrity of the teapot if i were to plot only the sub-objects found.

I suppose the right way to start is to find the faces that are intersecting the sphere like this:

out_pnt = find(vecnorm(verts-sphere_origin,2,2)>=sphere_radius);
out_faces = any(ismember(faces,out_pnt),2);

But going from there to finding the exact set of points that define the intersection is really difficult for me right now.

Any part of the answer could really help me, thanks in advance.

[darova]

Please show the data. Maybe it's possible to find intersection using contour

For more complicated cases use surfaceIntersection

[hagege ruben]

Thanks for you'r quick response!

The code should run as is given that you installed matGeom

[darova]

Can you please just attach the data? I don't want to install anything

[hagege ruben]

Ok, so instead of using drawSphere you can do the following:

hold on;
[X,Y,Z] = sphere(32);
sph = surf(gca,X.*sphere_radius+sphere_origin(1),Y.*sphere_radius+sphere_origin(2),Z.*sphere_radius+sphere_origin(3));
view([135 0])
axis equal

[darova]

You shuold use surfaceIntersection

[hagege ruben]

For anyone who is interested in doing something similar, I tried the following:

tea.faces = faces;
tea.vertices = verts;
DT = delaunayTriangulation(Z(:).*sphere_radius+sphere_origin(3), Y(:).*sphere_radius+sphere_origin(2), X(:).*sphere_radius+sphere_origin(1));
[sphere.faces, sphere.vertices] = freeBoundary(DT);
[intersect, Surf] = SurfaceIntersection(tea,sphere);
patch('Faces',Surf.faces,'Vertices',Surf.vertices,'LineWidth',2)

It solved the first part of my problem, thank you very much!

Now do you know how I can "slice" the teapot so that I have the part outside of the sphere with a hole in the shape of the sphere and the part inside the sphere? (the two parts should touch each other)

[darova]

Im using older MATLAB version and don't have teapot. Can you show it? What are you trying to achieve?

[hagege ruben]

This is a matlab built-in, the file has been created in 2014 so you must have a very old version

I've attached the file, if you need help reproducing another step becasuse of the version of you'r matlab I'd be happy to help.

My goal is to divide the "teapot" (or any object really) into parts that are divided by the intersection of another object (sphere in this case).

[hagege ruben]

This seems like the function I need, too bad it's not available...

[Joseph Olson]

Seems like the inShape function is what you want!

https://www.mathworks.com/matlabcentral/answers/496500-intersection-volume-of-two-3d-alphashapes

Answer 1

[DGM]

Normally, this would be my answer:

https://www.mathworks.com/matlabcentral/answers/2169008-are-there-matlab-functions-or-scripts-to-perform-boolean-operations-on-triangulated-solids

I was going to make another demo for this specific example, but there are a few problems. The given teapot model is represented using quadrilaterals instead of triangles, but that's fairly easy to convert. The bigger problem is that it's not a closed surface, and it contains at least one non-manifold vertex. That means we can't use CGAL on it unless we do a bunch of work to fix the model.

Using FEX #48613 can be done, but all it will give you is a list of edges in 3D. It won't return solid geometry. That doesn't seem like the intended goal.

EDIT:

I dug up a few alternate teapot models. A lot of them aren't valid solid geometry, and a lot of those that are have some other issues (extremely high-aspect triangles, near-coincident vertices) which result in problems due to the relatively limited precision of exported STLs. There are cleaner models than this one, but this one is smaller.

% given parameters
sphere_origin = [1.75 1 1]*3;
sphere_radius = 2.5;

% say you have some teapot geometry
T = stlread('utahteapot.stl'); % a different version (may still be problematic for CGAL)
Ft = T.ConnectivityList;
Vt = T.Points;

% and you have some sphere geometry
[X,Y,Z] = sphere(32); % using sphere to get gridded XYZ data
[Fs,Vs] = mesh2tri(X,Y,Z,'f'); % convert to triangulated FV data (FEX #28327)
Vs = Vs*sphere_radius + sphere_origin; % scale and offset

% preview display
patch('faces',Ft,'vertices',Vt,'facecolor','w','edgecolor','none');
hold on;
patch('faces',Fs,'vertices',Vs,'facecolor','w','edgecolor','none');
view(3); view(144,44); camlight; 
axis equal; grid on
xlabel('X'); ylabel('Y'); zlabel('Z')

% do the CSG
[Fout Vout] = scad_op(Ft,Vt,Fs,Vs,'difference');

% plot the result
clf
patch('faces',Fout,'vertices',Vout,'facecolor','w','edgecolor','none');
view(3); view(144,44); camlight; 
axis equal; grid on
xlabel('X'); ylabel('Y'); zlabel('Z')

If instead, you want the volume enclosed by their intersection, you can do that by specifying 'intersection' instead of 'difference'.

Again, all of this relies on external tools. The model file and scad_op() are attached. You'd still need OpenSCAD for it to work.