Creating a projection of a surface on a plane for area calculation.
조회 수: 91 (최근 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
댓글 수: 0
채택된 답변
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
% 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
댓글 수: 3
참고 항목
카테고리
Help Center 및 File Exchange에서 Surface and Mesh Plots에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!