File Exchange

image thumbnail

Surface Intersection

version 1.0.0.0 (373 KB) by Jaroslaw Tuszynski
Intersection of two triangulated surfaces

36 Downloads

Updated 01 Dec 2014

View License

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 .

Comments and Ratings (29)

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

jdiva6t9

XL

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

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)

Paul Zhang

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

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.

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 ?

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

Thank you

MarieLe

Excellent work! super useful. Thanks a lot!

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

Faez Alkadi

Annie Tran

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

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

Ben

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

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

I have the same problem as Marta ... the surface was generated with surf2patch and is rather well-behaved - I'm testing self-intersections ...

Marta

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

Excellent work. Thanks

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