intersectPlaneSurf

Intersection points of an arbitrary surface with an arbitrary plane.

Description

out=intersectPlaneSurf(p0, v, exx, eyy, ezz)
Intersection of a surface with an arbitrary plane. The plane must be specified with "p0" which is a point that the plane includes and a normal vector "v" of that plane.
"exx", "eyy" and "ezz" is the surface coordinates such that can be used with surf command. the uotput "out" is a nx3 cell that contains n number of discreete intersection regions in 3d coordinaetes.

You can use provided example file to get the result that is used in screenshot.

Note: You have to download geom3d toolbox by David Legland from FEX to visualize results correctly, and "myfrenet" function which i provide that can use example file

• MATLAB release: MATLAB 7.6 (R2008a)

[Ben]

Hi Mehmet,

Thanks a lot for your code! But I got error message at the line

if minDist1<thr || minDist2<thr

The error message is

??? Operands to the || and && operators must be convertible to
logical scalar values.

This is due to that the minDist1 and minDist2 are empty. Did you meet that?

Best,
Ben

[Mehmet OZTURK]

Sorry, you can find necessary function from

http://people.sc.fsu.edu/~jburkardt/m_src/geometry/geometry.html

[Mehmet OZTURK]

Hi Ben, I think the following code solves your problem...

clear,clc

load intSurfPlane;

x = InstMesh.Vertices(1,:)';
y = InstMesh.Vertices(2,:)';
z = InstMesh.Vertices(3,:)';

p0=[0 0 0]; v=[1 0 0];

% plane=createPlane(p0,v); % createPlane from geom3d toolbox
[ A, B, C, D ] = plane_normal2imp_3d ( p0, v ); % create implicit plane function
segment_start=nan(3,round(size(InstMesh.Faces,2)/2));
segment_finish=segment_start;
count=1;
for s=1:size(InstMesh.Faces,2)
t(:,1)=InstMesh.Vertices(:,InstMesh.Faces(1,s));
t(:,2)=InstMesh.Vertices(:,InstMesh.Faces(2,s));
t(:,3)=InstMesh.Vertices(:,InstMesh.Faces(3,s));
[ num_int, pi ] = plane_imp_triangle_int_3d ( A, B, C, D, t );
if num_int==2
segment_start(:,count)=pi(:,1);
segment_finish(:,count)=pi(:,2);
count=count+1;
% hold on,plot3(pi(1,:),pi(2,:),pi(3,:))
end
end
segment_start(:,all(isnan(segment_start),1))=[]; % remove unused poritons
segment_finish(:,all(isnan(segment_finish),1))=[]; % remove unused poritons

thr=2; nol=1;
while ~isempty(segment_start)
lin{nol}=[segment_start(:,1) segment_finish(:,1)];
segment_start(:,1)=[];
segment_finish(:,1)=[];
while 1
testDistStart1=sum((lin{nol}(1,end)-segment_start(1,:)).^2 + ...
(lin{nol}(2,end)-segment_start(2,:)).^2 + ...
(lin{nol}(3,end)-segment_start(3,:)).^2,1);
testDistStart2=sum((lin{nol}(1,end)-segment_finish(1,:)).^2 + ...
(lin{nol}(2,end)-segment_finish(2,:)).^2 + ...
(lin{nol}(3,end)-segment_finish(3,:)).^2,1);
[minDist1, best_ind1]=min(testDistStart1);
[minDist2, best_ind2]=min(testDistStart2);
if minDist1<thr || minDist2<thr
if minDist1<minDist2
lin{nol}=[lin{nol} segment_finish(:,best_ind1)];
segment_start(:,best_ind1)=[];
segment_finish(:,best_ind1)=[];
else
lin{nol}=[lin{nol} segment_start(:,best_ind2)];
segment_start(:,best_ind2)=[];
segment_finish(:,best_ind2)=[];
end
else
nol=nol+1;
break
end
end
end

s = trisurf(InstMesh.Faces', x, y, z, 'EdgeColor', 'none', 'FaceColor', 'interp');
daspect([1 1 1]);
for m=1:size(lin)
hold on,plot3(lin{m}(1,:),lin{m}(2,:),lin{m}(3,:))
end
hold off

[Ben]

Hi Mehmet,

What if I have a triangulated surface mesh, where "X", "Y" and "Z" are the coordinates of the vertices? In this case, is there a way to use your code?

Thank you!
Ben

[Barry]

Excellent and very useful...

[Paula Saameno]

Very useful and easy to understand. Thanks!!

[Mehmet OZTURK]

"myfrenet.m".
http://www.mathworks.com/matlabcentral/fileexchange/22376-myfrenet

[Mehdi mzn]

I can not find "myfrenet" function there is no function with this name in geom3d toolbox or any where else.
Where can I find it?

[Erik Benkler]

There is a small error: on line 16 it should be sqrt instead of fsqrt. The result is exactly what I was looking for, but documentation and comments in the code could be better.