필터 지우기
필터 지우기

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

조회 수: 91 (최근 30일)
Denis Kiryukhin
Denis Kiryukhin 2022년 12월 30일
댓글: Denis Kiryukhin 2022년 12월 30일
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

채택된 답변

Karim
Karim 2022년 12월 30일
편집: Karim 2022년 12월 30일
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
  댓글 수: 3
Denis Kiryukhin
Denis Kiryukhin 2022년 12월 30일
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
Denis Kiryukhin 2022년 12월 30일
Thanks a lot, works flawlessly now!

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

추가 답변 (1개)

Matt J
Matt J 2022년 12월 30일
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.

카테고리

Help CenterFile Exchange에서 Surface and Mesh Plots에 대해 자세히 알아보기

제품


릴리스

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by