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

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.

[Eric Bernard Dilger]

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

Answer 1

[Star Strider]

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.

.

[Eric Bernard Dilger]

Thank you! It helped a lot!

[Star Strider]

As always, my pleasure!

This was an interesting challenge!

[Eric Bernard Dilger]

Hello buddy! It's me again. Now, I'd like to plot the contourf of this 3D surface with a non-horizontal plane. I could plot the intersection between them, just defining the points of Z_matrix to be different than z_plane*ones(2).

Z_matrix = [0.3 0.3; 1.7 1.7];

surf(xlim, ylim, plane_z*Z_matrix, 'FaceColor','b', 'FaceAlpha',0.5) % Plot Plane At 'plane_z'

Unfortunatelly, I can't plot the contourf using the command below:

[c,h] = contourf(plotx, ploty, P, [plotx, ploty, plane_z*Z_matrix]);

Matlab is returning the following message:

"Error using horzcat. Dimensions of arrays being concatenated are not consistent."

Do you have any suggestion?

[Star Strider]

‘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.

Answer 2

[burak ergocmen]

With this code:

https://www.mathworks.com/matlabcentral/fileexchange/48613-surface-intersection?s_tid=srchtitle