Creating a projection of a surface on a plane for area calculation.

Hi!

Recently I've been working on a program that takes an object (a tetrahedron in my example) and rotates it to face a stratified sample of random points on a semi-sphere (which were generated with this wonderful person's code) one by one. I now want to calculate the area of the object's orthogonal projection on the x-z plane, and I'm not sure how to do it. In simpler terms, I basically want the shadow area of the object in the x-z plane. I've been trying to use this response to a similar question to guide me but it doesn't work for me, unless I'm doing it wrong. I posted my code below, I'm sure that it is not the most efficient way to tackle the problem, but it makes sense to me so far. Could someone please give me guidance on how to create a projection of the tetrahedron on the x-z plane? Thank you in advance!

Plot setup

view(0,0)
grid on
xlim([-1.02 1.02]);ylim([-1.02 1.02]);zlim([0.02 2.02])
xlabel('x-axis');ylabel('y-axis');zlabel('z-axis')

Create the pyramid

[x,y] = pol2cart(deg2rad(0:120:360),1);
x1 = [x*0; x; x*0];
y1 = [y*0; y; y*0];
z1 = [x*0+1; x*0; x*0];
s = surface(x1,y1,z1+1); % the pyramid was raised by one unit in the z-axis.

Get a 1000 random stratified sample of points on a sphere. (https://uk.mathworks.com/matlabcentral/fileexchange/37004-suite-of-functions-to-perform-uniform-sampling-of-a-sphere)

V=RandSampleSphere(1000,'stratified');

Run the pyramid rotation using the previously acquired points on a sphere.

for k = 1:length(V)/2
    delete(s) %reset the pyramid
    x1 = [x*0; x; x*0];
    y1 = [y*0; y; y*0];
    z1 = [x*0+1; x*0; x*0];
    s = surface(x1,y1,z1+1); %reestablish the pyramid (I did not know how to reset only the pyramid without resetting everything)
    alpha 0.5
    
    a = [V(k,1), V(k,2), 0];
    b = [0 0 2];
    g = cross(a,b); %find the cross product between a and b
    if norm(g)==0; continue; end %if the magnitude of the cross product is zero, then continue to the next iteration
    
    beta = acos(V(k,3));        
    rotate(s,g,-rad2deg(beta), [0 0 1])
    drawnow     
end

Answer 1

[Karim]

There is no need for rotations or even fine sampling, since it's a pyramid it has straight edges. If you project the points that make up these edges onto the plane you can figure out the projected surface and evaluate the area. See below for a demonstration.

% define the plane we want to project onto
po = [0 0 0]; % point of our plane, here simply the origin
no = [0 1 0]; % normal vector of the plane i.e. pointing along the y-axis

% set up the pyramid according to the OP's parameters and inputs
[x,y] = pol2cart(deg2rad(0:120:360),1);
x1 = [x*0; x; x*0];
y1 = [y*0; y; y*0];
z1 = [x*0+1; x*0; x*0] + 1; % the pyramid was raised by one unit in the z-axis.

% concatenate the pyarmiad grid
grid = [x1(:) y1(:) z1(:)];

% vectors from the pyramid grid to P0:
v = grid - repmat(po,numel(x1(:)),1);

% dot vector with the normal vector of the plane of intrest
% this represents the distance from the grid points to our plane along the normal
d = v(:,1)*no(1) + v(:,2)*no(2) + v(:,3)*no(3); 

% get the projected points
projected_grid = grid - d .* repmat(no,numel(x1(:)),1);

% make a plot of what we have so far
figure
subplot(1,2,1)
s = surface(x1,y1,z1); 
grid on
view(3)
title("Pyramid in 3D")
xlabel('x-axis');ylabel('y-axis');zlabel('z-axis')
subplot(1,2,2)
scatter3(projected_grid(:,1),projected_grid(:,2),projected_grid(:,3),'r','filled')
grid on
view(3)
title("Projected points in red")
xlabel('x-axis');ylabel('y-axis');zlabel('z-axis')

% as you can see we have all the projected points in the plane
% there are some redunant points due to the way the pyramid is setup, first
% extract the unique points
projected_grid = uniquetol(projected_grid,'ByRows',true);

% we only need the outer points of the region, any internal points won't be
% visible in the shadow, we can use the convex hull function to determine them
[U,S] = svd( bsxfun( @minus, projected_grid, mean(projected_grid)),0);  
[idx,area_convhull] = convhull(U*S(:,1:2));

% note that the convex hull function also returns the area, however just
% for fun we can take it a few steps further and also plot the shadow
area_convhull

area_convhull = 0.7500

% since the shadow is a 2d surface we can use polyshape to easily plot the
% shadow
shadow = polyshape(projected_grid(idx,1),projected_grid(idx,3));

figure
plot(shadow)
grid on
title("Pyramid shadow in the xz plane")
xlabel('x-axis'); ylabel('z-axis')

% as a double check we can also request the area from the polyshape:
area = shadow.area

area = 0.7500

[Denis Kiryukhin]

Thank you, this works perfectly! :)

[Denis Kiryukhin]

Hi again! Sorry to unaccept the answer. The projection works in the initial position of the pyramid, but when I rotate it, the resultant shadow doesn't look right. I attached some screenshots. The left is what the projection looks like in the figure viewer and the right is the shadow. It seems like I need to only plot the points with maximum/minimum values for both x and y?

Here's the updated code:

Setup

tic
po = [0 0 0]; % point of our plane, here simply the origin
no = [0 1 0]; % normal vector of the plane i.e. pointing along the y-axis
num = 100000;
areas = zeros(1000,1,'double');

Create the pyramid

[x,y] = pol2cart(deg2rad(0:120:360),1);
x1 = [x*0; x; x*0];
y1 = [y*0; y; y*0];
z1 = [x*0+sqrt(2); x*0; x*0];
s = surface(x1,y1,z1+1); % the pyramid was raised by one unit in the z-axis.

Get a random uniform sample of num points on a sphere. (https://uk.mathworks.com/matlabcentral/fileexchange/37004-suite-of-functions-to-perform-uniform-sampling-of-a-sphere)

V=RandSampleSphere(num,'uniform');

Run the pyramid rotation using the previously acquired points on a sphere.

for k = 1:length(V)
    drawnow

Get x-z projection of the shape

    %get data as a result of rotation
    XX = get(s,'xdata');
    YY = get(s,'ydata');
    ZZ = get(s,'zdata');
    % concat the pyramid grid
    gr = [XX(:) YY(:) ZZ(:)];
    v = gr - repmat(po,numel(XX(:)),1);
    % dot vector with the normal vector of the plane of intrest
    % this represents the distance from the grid points to our plane along the normal
    d = v(:,1)*no(1) + v(:,2)*no(2) + v(:,3)*no(3); 
    % get the projected points
    projected_grid = gr - d .* repmat(no,numel(XX(:)),1);
    projected_grid = unique(projected_grid,'rows','stable');
    % since the shadow is a 2d surface we can use polyshape
    shadow = polyshape(projected_grid(:,1),projected_grid(:,3)); 
    warning('off','last')
    plot(shadow)
    
    view(3);
    xlim([-2.02 2.02]);ylim([-2.02 2.02]);zlim([-1.02 3.02])
    xlabel('x-axis'); ylabel('y-axis'); zlabel('z-axis')
    areas(k) = shadow.area;

Rotate the pyramid

    delete(s) %reset the pyramid
    s = surface(x1,y1,z1+1); %reestablish the pyramid (I did not know how to reset only the pyramid without resetting everything)
    alpha 0.5

    a = [V(k,1), V(k,2), 0];
    b = [0 0 2];
    g = cross(a,b); %find the cross product between a and b
    if norm(g)==0; continue; end %if the magnitude of the cross product is zero, then continue to the next iteration
    
    beta = acos(V(k,3));        
    rotate(s,g,-rad2deg(beta), [0 0 1]);
end
mean(areas)
toc

[Denis Kiryukhin]

Thanks a lot, works flawlessly now!

Answer 2

[Matt J]

Could someone please give me guidance on how to create a projection of the tetrahedron on the x-z plane?

Just throw away the y-coordinates of all the vertices. Then you can use convhull to calculate its area.