File Exchange

## Surface Intersection

version 1.0.0.0 (373 KB) by Jaroslaw Tuszynski

### Jaroslaw Tuszynski (view profile)

Intersection of two triangulated surfaces

Updated 01 Dec 2014

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.

### Cite As

Jaroslaw Tuszynski (2020). Surface Intersection (https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection), MATLAB Central File Exchange. Retrieved .

Claudio Tomasi

### 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)

Can it be used for 2D meshes or is only for 3D ?

darova

jdiva6t9

Hassaan Maqsood

### Hassaan Maqsood (view profile)

Frederik Thönnißen

XL

### 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

### 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

Paul Zhang

### 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

### 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

### 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

### James Douthwaite (view profile)

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

Ahmed Zankoor

Thank you

Donghua Zhao

MarieLe

### MarieLe (view profile)

Excellent work! super useful. Thanks a lot!

Jacob Lynch August

### Jacob Lynch August (view profile)

Trying to get this to work with R2017b, but running into many issues.

Annie Tran

### 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

### DF (view profile)

Fantastic functions, works great in MATLAB R2014a and R2014b. Thank you!

Shih-Jung Hsu

Ben

### 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

Andrew Bopp

LU ZHAN

### 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

Ernst Schwartz

### 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

### 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

### Thien Tran (view profile)

Excellent work. Thanks

##### MATLAB Release Compatibility
Created with R2014b
Compatible with any release
##### Platform Compatibility
Windows macOS Linux