How to determine whether point lies within pyramid volume

15 views (last 30 days)
A square pyramid is described by five points E0, E1, E2, E3, E4. eg E0= [0,0,1];E1= [1,1,0]; E2= [-1,1,0]; E3= [-1,-1,0]; E4= [1,-1,0], how to determine if giving point p[x,y,z] is located within the pramid volume or not. i need the algorithms or the mathematical equation.

Accepted Answer

Jan
Jan on 27 Jun 2019
Edited: Jan on 27 Jun 2019
Because the pyramid is regular, the algorithm is easy:
E0 = [ 0, 0, 1];
E1 = [ 1, 1, 0];
E2 = [-1, 1, 0];
E3 = [-1, -1, 0];
E4 = [ 1, -1, 0];
% Define the parameters of the pyramid:
Pyramid.Z = E0(3);
Pyramid.Base = [min(E(:, 1)), max(E(:, 1)), min(E(:, 2)), max(E(:,2))]
% Choose random points in the cuboid surrounding the pyramid:
P = rand(10000, 3) .* [2, 2, 1] - [1, 1, 0];
% Select points inside the pyramid:
T = inPyramid(P, Pyramid);
P = P(T, :);
% Show the result:
figure;
plot3(P(:,1), P(:,2), P(:,3), '.')
hold on
plot3(E0(1), E0(2), E0(3), 'ro');
plot3(E1(1), E1(2), E1(3), 'ro');
plot3(E2(1), E2(2), E2(3), 'ro');
plot3(E3(1), E3(2), E3(3), 'ro');
plot3(E4(1), E4(2), E4(3), 'ro');
xlabel('x')
ylabel('y')
zlabel('z')
For a point with at Z==0, you have to check, if it is inside the square. This square is shrinking linearily with the height until Z==1.
function T = inPyramid(P, Pyramid)
% P: 3D position as [N x 3] matrix
% Pyramid: Struct, .Z: Height of pyramid
% .Base: minimum and maximum of X and Y positions at the
% base as [1 x 4] vector
S = Pyramid.Base .* (Pyramid.Z - P(:, 3)) / Pyramid.Z;
T = (P(:, 3) >= 0 & P(:, 3) <= Pyramid.Z) & ... % Inside Z range
(P(:, 1) >= S(:, 1) & P(:, 1) <= S(:, 2)) & ... % Inside X range
(P(:, 2) >= S(:, 3) & P(:, 2) <= S(:, 4)); % Inside Y range
end
This is elementary maths only.
A constructive approach instead of this rejection method:
  6 Comments
Jun Liang
Jun Liang on 29 Apr 2020
Dear Jan,
Thanks for your answer.
Select the same vertex of the cube (ie. origin of three coordinates -POS,NEU,NEG )as the vertex of the three pyramids.
The three faces that are not adjacent to the vertex are the bottom faces of the three pyramids.
Thanks you very much.

Sign in to comment.

More Answers (2)

Jon
Jon on 26 Jun 2019
Edited: Jon on 26 Jun 2019
Hi Musa,
You can first define your pyramid as an alphaShape, for example:
P = alphaShape([0 1 -1 1]',[0 1 -1 -1]',[1 0 0 0]')
Note the arguments to alphaShape are the x coordinates of each vertex, the y coordinates of each vertex and the z coordinates of each vertex. Note they have to be column vectors, that is why I applied the transpose ' operator.
You can then use alphaShape's inShape method to test if a point is inside the pyramid
for example the point [0.5, 0.5 0.5] is in the interior, you can test this using
isInside = inShape(P,[0.5,0.5,0.5])
which returns true
and the point [0.5, 0.5, 3] is not in the interior
so you can test this using
isInside = inShape(P,[0.5,0.5,3])
which returns false
  1 Comment
Jon
Jon on 27 Jun 2019
Edited: Jon on 27 Jun 2019
I took a quick look, but the Mathworks does not seem to provide any readily available documentation of the algorithm used by the inShape method of their alphaShape class.
For the special case of a tetrahedron (pyramid) I'm sure there are many approaches to determine whether a point is in the interior.
In case it is helpful, I will try to provide a further explanation for the idea behind some of the approaches which were more succintly presented in the link above.
First some background. Given a point, p, and plane P, an algebraic sign can be assigned to the outward normal from the plane to the point. One way to do this is using the vector cross product. So given a point p and plane P the point is either on one side of the plane or the other and it can be assigned a value of 1 or -1.
Now to determine if a point p is in the interior, we start at one of the four faces (planes) and determine if the point p is on the same side of the face as the remaining vertex, (the vertex that is not in the face) We repeat this for all four faces. If the point is on the same side as the remaining vertex for all four faces then it on the inside, otherwise it is not on the inside.

Sign in to comment.


Torsten
Torsten on 27 Jun 2019
Another program for your issue:
https://de.mathworks.com/matlabcentral/fileexchange/10226-inhull

Categories

Find more on Bounding Regions in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!