Function calculates intersection of any two triangulated surfaces using triangle/triangle intersection algorithm proposed by Tomas Möller (1997) and implemented as highly vectorized MATLAB code. The algorithm was expanded to include calculation of the intersection surface, in addition to boolean matrix cataloging which triangle from one surface intersects with which triangle in the other surface. Function can be used for contour line calculations and can handle surfaces residing on the same plane.
Jaroslaw Tuszynski (2020). Surface Intersection (https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection), MATLAB Central File Exchange. Retrieved .
Create scripts with code, output, and formatted text in a single executable document.
Claudio Tomasi (view profile)
I would like to use it for 2D meshes. For now I just putted a 0 z-coordinate in each point and I get the results, but when the meshes differ too much I get errors like "The condition input argument must be a scalar logical.
Error in SurfaceIntersection>TriangleIntersection2D (line 531)
assert(max(faces(:))<=size(vertices,1), 'Bad surface definition')"
Can it be used for 2D meshes or is only for 3D ?
darova (view profile)
jdiva6t9 (view profile)
Hassaan Maqsood (view profile)
Frederik Thönnißen (view profile)
XL (view profile)
I used Surface intersection function to get the intersection points between the cylinder (I created cylinder mesh from my input data) and a horizontal planes; however, I am getting more points than expected. I used the code below, any suggestion on how to solve this?
% Below are expanded to 2D version (you have 13 rows and 18 columns for your design);
[x,z] = meshgrid(0:1:5, 0:1:2);
% generate mesh for 2D version,
DT = delaunay(x,z);
%Below are 3D points I would like to generate mesh
x1 = [9.75 9.73 9.70 -25.70 -25.78 -25.83 9.72 9.69 9.62 40.69 40.62 40.59 40.70 40.65 40.59 9.75 9.73 9.70];
z1 = [36.82 49.24 61.73 36.62 49.06 61.54 36.80 49.21 61.69 36.95 49.38 61.87 36.96 49.39 61.88 36.82 49.24 61.73];
y1 = [-111.48 -111.47 -111.47 -82.26 -82.33 -82.40 -39.80 -39.69 -39.66 -57.66 -57.66 -57.60 -93.47 -93.42 -93.41 -111.48 -111.47 -111.47];
trimesh(DT,x1,y1,z1);
DT = delaunayTriangulation(x1(:), y1(:), z1(:));
hold on;
% trimesh(tri,x,y, z);
[Surface1.faces, Surface1.vertices] = freeBoundary(DT);
%%%SURFACE 2%%%%% The Horizontal Cutting Surface
Surface2=[];
z=45;
Surface2.vertices((1:3),:) = [-40, 0, z; -40, -140, z; 160, -80, z];
Surface2.faces(1,:) = (1:3);
[~, Surf12] = SurfaceIntersection(Surface1, Surface2);
IntersectXYZ=Surf12.vertices; % Coordinates of Intersection pints
Surf12.vertices(:,3) = -1; % project the contour lines on a single plane
% plot the results
S=Surface2; trisurf(S.faces, S.vertices(:,1),S.vertices(:,2),S.vertices(:,3),'FaceAlpha', 0.25, 'FaceColor', 'b')
S=Surf12;
line([S.vertices(S.edges(:,1),1), S.vertices(S.edges(:,2),1)]',...
[S.vertices(S.edges(:,1),2), S.vertices(S.edges(:,2),2)]',...
[S.vertices(S.edges(:,1),3), S.vertices(S.edges(:,2),3)]','Color', 'r');
axis equal
plot3(IntersectXYZ(:,1), IntersectXYZ(:,2), IntersectXYZ(:,3), 'o');
view(170, 20)
XL (view profile)
I used Surface intersection function to get the intersection points between the cylinder (I created cylinder mesh from my input data) and a horizontal planes; however, I am getting more points than expected. I used the code below, any suggestion on how to solve this?
% Below are expanded to 2D version (you have 13 rows and 18 columns for your design);
[x,z] = meshgrid(0:1:18, 0:1:12);
% generate mesh for 2D version,
tri = delaunay(x,z);
%Below are 3D points I would like to generate mesh
x1= [ x1 values];
y1= [ y1 values];
z1= [ z1 values];
Plot the results for 3D version
trimesh(tri,x1,y1,z1);
DT = delaunayTriangulation(x1(:), y1(:), z1(:));
trimesh(tri,x1,y1,z1)
hold on;
% trimesh(tri,x,y, z);
[Surface1.faces, Surface1.vertices] = freeBoundary(DT);
%%%SURFACE 2%%%%% The Horizontal Cutting Surface
Surface2=[];
z=180;
Surface2.vertices((1:3),:) = [-40, 0, z; -40, -140, z; 160, -80, z];
Surface2.faces(1,:) = (1:3);
[~, Surf12] = SurfaceIntersection(Surface1, Surface2);
IntersectXYZ=Surf12.vertices; % Coordinates of Intersection pints
Surf12.vertices(:,3) = -1; % project the contour lines on a single plane
% plot the results
% figure(1); clf; hold on
% S=Surface1; trisurf(S.faces, S.vertices(:,1),S.vertices(:,2),S.vertices(:,3),'FaceAlpha', 0.5, 'FaceColor', 'r')
S=Surface2; trisurf(S.faces, S.vertices(:,1),S.vertices(:,2),S.vertices(:,3),'FaceAlpha', 0.25, 'FaceColor', 'b')
S=Surf12;
line([S.vertices(S.edges(:,1),1), S.vertices(S.edges(:,2),1)]',...
[S.vertices(S.edges(:,1),2), S.vertices(S.edges(:,2),2)]',...
[S.vertices(S.edges(:,1),3), S.vertices(S.edges(:,2),3)]','Color', 'r');
axis equal
view(170, 20)
Chirag Gupta (view profile)
Paul Zhang (view profile)
Seems to not handle surfaces with 0 area well. But you can make an approximation by adding 1e-6 perturbation to some vertices. For example if intSurface1 represents a 1d curve and triangle faces all have index [i j j], you can do the following to make intSurfacet an approximation of the 1d curve that's has basically the same intersection. surface5 is whatever surface you want to intersect with a curve. This hack may or may not be appropriate for your use case.
eps = 1e-6;
intSurfacet = rmfield(intSurface1,'edges');
nV = size(intSurfacet.vertices,1);
intSurfacet.vertices = [intSurfacet.vertices; intSurfacet.vertices(intSurfacet.faces(:,3),:) + (rand(size(intSurfacet.vertices(intSurfacet.faces(:,3),:)))-.5)*eps];
intSurfacet.faces(:,3) = nV+1:size(intSurfacet.vertices,1);
[intMatrix_1_1, intSurface_1_1] = SurfaceIntersection(intSurfacet, surface5);
Paul Zhang (view profile)
Sufficient for surface intersection. When I take the intersection surface (a curve) and try to intersect it with another surface I get the error James Douthwaite mentions below. Would be nice to have a fix.
Sébastien THIBAUD (view profile)
Excellent but i've tried on a non closed surface and i didn't obtained the complete intersection. Do you have an idea to solve this problem ?
James Douthwaite (view profile)
Downloaded, other than the example, I recieve the following error :
Array indices must be positive integers or logical values.
Error in SurfaceIntersection>TriangleIntersection3D_Moller (line 353)
a1 = lut(sign(dv)*[9; 3; 1] + 14); % calculate the key
and call the look-up table
Its clear if the geometries need more conditioning than the .vertices and .face structure.
Nikolas Lamb (view profile)
Ahmed Zankoor (view profile)
Thank you
Donghua Zhao (view profile)
MarieLe (view profile)
Excellent work! super useful. Thanks a lot!
Jacob Lynch August (view profile)
Trying to get this to work with R2017b, but running into many issues.
Faez Alkadi (view profile)
Annie Tran (view profile)
Hi, I have problem with self-intersection test too.
I defined two surfaces (each include three points). it works when two planes are exactly the same, i.e as following.
verticesA = [-0.268473,-0.515378,0;
-0.268473,-0.515378,1;
2.731527,-0.515378,0];
verticesB = [-0.268473,-0.515378,0;
-0.268473,-0.515378,1;
2.731527,-0.515378,0];
When I change the Z -coordinate of point 2 in surface B
verticesB = [-0.268473,-0.515378,0;
-0.268473,-0.515378,2.5;
2.731527,-0.515378,0];
=> it throws error: 'Bad surface definition'
Any idea? Thanks
DF (view profile)
Fantastic functions, works great in MATLAB R2014a and R2014b. Thank you!
Shih-Jung Hsu (view profile)
Ben (view profile)
Works well.
A little heavy on memory if you input larger meshes, and seems to return a few duplicate edges in the output, but otherwise runs smoothly and reliably.
Ben (view profile)
Andrew Bopp (view profile)
LU ZHAN (view profile)
Hello guys, i have the same issue said:
Error in SurfaceIntersection>TriangleIntersection3D_Moller (line 357)
But then i try to get rid of all the NaN values in my surface data. It works!
Ehsan golk (view profile)
Ernst Schwartz (view profile)
I have the same problem as Marta ... the surface was generated with surf2patch and is rather well-behaved - I'm testing self-intersections ...
Marta (view profile)
I can't get this to work...
"Subscript indices must either be real positive integers or logicals.
Error in SurfaceIntersection>TriangleIntersection3D_Moller (line 357)
a1 = lut(sign(dv)*[9; 3; 1] + 14); % calculate the key and call the look-up table
Error in SurfaceIntersection (line 218)
[intMatrix(tMsk), intSurface] = TriangleIntersection3D_Moller(..."
The 'rapid' algorithm also crashes. Can you help?
Thien Tran (view profile)
Excellent work. Thanks