version 1.0.0.0 (373 KB) by
Jaroslaw Tuszynski

Intersection of two triangulated surfaces

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 (2021). Surface Intersection (https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection), MATLAB Central File Exchange. Retrieved .

Created with
R2014b

Compatible with any release

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

KrystianIt seems that it limits Intersection to x=-35

OctavianSame here, works well (~3 min) with separating amygdala from the rest of cortex, not hippocampus (slightly larger surface)

Error using '

Out of memory. Type HELP MEMORY for your options.

Error in SurfaceIntersection (line 188)

dv1 = dv(:,surface1.faces(:,1))';

Jose A. Fernandez-RoldanNot working properly with isosurfaces.

Out of memory. Type "help memory" for your options.

Error in SurfaceIntersection (line 155)

du1 = du(:,surface2.faces(:,1));

KonradFor all who struggel with the "Array indices must be positive integers or logical values" error when calling the look-up table:

Try to 'cleanup' the surfaces like this:

% remove duplicate vertices

tol = 1e-10; % tolerance

[~, vui, vi] = unique(round(vertices./tol),'rows');

vertices= vertices(vui,:);

% update faces:

faces = vi(faces);

% remove duplicate faces

faces = unique(sort(faces,2),'rows');

% remove faces that contain a vertex multiple times (dont do this with tetrahedral surfaces!)

faces = faces(all(diff(faces,1,2)>0,2),:);

Lican WangjonasClaudio TomasiI 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 ?

darovajdiva6t9Hassaan MaqsoodFrederik ThönnißenXLI 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)

XLI 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 GuptaPaul ZhangSeems 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 ZhangSufficient 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 THIBAUDExcellent 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 DouthwaiteDownloaded, 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 LambAhmed ZankoorThank you

Donghua ZhaoMarieLeExcellent work! super useful. Thanks a lot!

Jacob Lynch AugustTrying to get this to work with R2017b, but running into many issues.

Faez AlkadiAnnie TranHi, 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

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

Shih-Jung HsuBenWorks 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.

BenAndrew BoppLU ZHANHello 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 golkErnst SchwartzI have the same problem as Marta ... the surface was generated with surf2patch and is rather well-behaved - I'm testing self-intersections ...

MartaI 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 TranExcellent work. Thanks