MATLAB Answers

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

24 ビュー (過去 30 日間)
Eric Bernard Dilger
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
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
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 件のコメント
Star Strider
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.

サインインしてコメントする。

その他の回答 (0 件)

Community Treasure Hunt

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

Start Hunting!

Translated by