# How do I plot an intersection of a surface 3D and a plane?

조회 수: 66(최근 30일)
Eric Bernard Dilger 2021년 4월 25일
댓글: Star Strider 2021년 5월 31일
I'm interested to plot the intersection of a 3D surface with a 2D plane (x,y) for differents values of the height z. I need to plot the upper view of this intersection and represent its points as a dark pixel of value1 if they're below the intersection, and as a white pixel of value 0 if they're above. I also need to store this data points in a matricial form so that it represents the (x,y) cutting plane. Any suggestions?
The code below it's a midpoint displacement algorithm that I'm using to plot 3D surfaces.
##### 댓글 수: 1표시숨기기 없음
Eric Bernard Dilger 2021년 4월 25일
H=0.5;
nmax=6;
length=2^nmax+1;
step=2^nmax;
halfstep=step/2;
x=linspace(0,1,length);
y=x;
[x,y]=meshgrid(x,y);
z=zeros(size(x));
for n=1 : nmax
range=2^(-2*n*H);
for i = 1 : step : length - step
for j = 1 : step : length - step
z(i+halfstep,j)=0.5*(z(i,j)+z(i+step,j))+(1-2*rand)*range; % (1+0,5 ; j)
z(i,j+halfstep)=0.5*(z(i,j)+z(i,j+step))+(1-2*rand)*range; % (i ; j+0,5)
z(i+step,j+halfstep)=0.5*(z(i+step,j)+z(i+step,j+step))+(1-2*rand)*range; % (i+1 ; j+0,5)
z(i+halfstep,j+step)=0.5*(z(i,j+step)+z(i+step,j+step))+(1-2*rand)*range; % (i+0.5 ; j+1)
z(i+halfstep,j+halfstep)=0.25*(z(i,j)+z(i+step,j)+z(i,j+step)+z(i+step,j+step))+(1-2*rand)*range;
end
end
step=step/2;
halfstep=halfstep/2;
p=1;
for i=1 : 2^n +1
q=1;
for j = 1 : 2^n + 1
plotx(i,j) = x(p,q);
ploty(i,j) = y(p,q);
plotz(i,j) = z(p,q);
q = q + step;
end
p = p + step;
end
surf(plotx,ploty,plotz,'facecolor','w')
axis([0 1 0 1 -0.5 0.5])
axis off
shg
pause(0.1)
end

댓글을 달려면 로그인하십시오.

### 채택된 답변

Star Strider 2021년 4월 25일
Use the contourf function to get the intersection of the plane and the surface, and plot the result —
H=0.5;
nmax=6;
length=2^nmax+1;
step=2^nmax;
halfstep=step/2;
x=linspace(0,1,length);
y=x;
[x,y]=meshgrid(x,y);
z=zeros(size(x));
for n=1 : nmax
range=2^(-2*n*H);
for i = 1 : step : length - step
for j = 1 : step : length - step
z(i+halfstep,j)=0.5*(z(i,j)+z(i+step,j))+(1-2*rand)*range; % (1+0,5 ; j)
z(i,j+halfstep)=0.5*(z(i,j)+z(i,j+step))+(1-2*rand)*range; % (i ; j+0,5)
z(i+step,j+halfstep)=0.5*(z(i+step,j)+z(i+step,j+step))+(1-2*rand)*range; % (i+1 ; j+0,5)
z(i+halfstep,j+step)=0.5*(z(i,j+step)+z(i+step,j+step))+(1-2*rand)*range; % (i+0.5 ; j+1)
z(i+halfstep,j+halfstep)=0.25*(z(i,j)+z(i+step,j)+z(i,j+step)+z(i+step,j+step))+(1-2*rand)*range;
end
end
step=step/2;
halfstep=halfstep/2;
p=1;
for i=1 : 2^n +1
q=1;
for j = 1 : 2^n + 1
plotx(i,j) = x(p,q);
ploty(i,j) = y(p,q);
plotz(i,j) = z(p,q);
q = q + step;
end
p = p + step;
end
surf(plotx,ploty,plotz,'facecolor','w')
axis([0 1 0 1 -0.5 0.5])
% axis off
shg
pause(0.1)
end
hold on
plane_z = +0.15;
surf(xlim, ylim, plane_z*ones(2), 'FaceColor','b', 'FaceAlpha',0.5) % Plot Plane At 'plane_z'
hold off
xlabel('X \rightarrow', 'Rotation',20)
ylabel('\leftarrow Y', 'Rotation',-25)
zlabel('Z \rightarrow')
ZL = zlim;
figure
[c,h] = contourf(plotx, ploty, plotz, [ZL(1) plane_z]); % Contour Plot At 'plane_z'
colormap(bone(2));
colorbar
axis('equal')
grid
xlabel('X')
ylabel('Y')
Levels = plane_z;
idx = find(c(1,:) == Levels);
Len = c(2,idx);
for k = 1:numel(idx) % Retrieve Contour (x,y) Data
xc{k} = c(1,idx+1:Len(k));
yc{k} = c(2,idx+1:Len(k));
end
Contour_Information = sprintf('Number of different contours at %.2f = %d\n', plane_z, numel(idx))
Contour_Information =
'Number of different contours at 0.15 = 7 '
I altered the code slightly to show the plane in the original surf plot in order to demonstrate that the contour plot creates it correctly.
See the contourf documentation on Contour Matrix M to understand how to get the coordinates of the contours themselves. (This is not trivial, however it is not difficult. I included code that does that and appears to work here.)
I will let you experiment with getting the image and saving it. See the documentation on imwrite for one approach.
.
##### 댓글 수: 4표시숨기기 이전 댓글 수: 3
Star Strider 2021년 5월 31일
Do you have any suggestion?
Not immediately. I would need to have the rest of the code to see if I can understand what the probllem is.

댓글을 달려면 로그인하십시오.

### Community Treasure Hunt

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

Start Hunting!

Translated by