How to determine whether point lies within pyramid volume

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.

Answer 1

[Jan]

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:

https://www.mathworks.com/matlabcentral/answers/327990-generate-random-coordinates-inside-a-convex-polytope

[Jon]

I'm wondering what the goal is here and why you prefer (unaccepted earlier answer using alphaShape) this answer. While I agree Jan's approach is quite nice, it is not clear to me that there is an advantage to writing, and testing a new, special method that just works for regular pyramids when there is already a two line solution (as outlined in my earlier post) available using alphaShape. Unless this is for a homework assignment where you must develop your own algorithm, (in which case there may be some ethical questions about using Jan's solution) ,or a performance issue, I would generally recommend leveraging MATLAB's built in, presumably thoroughly tested solutions whenever possible. Especially, when as in this case, they are so simple.

[Jan]

@Jonathan: With alphaShape you must check explicitely, if the automagically selected Alpha value produce a convex shape. I admit, for this simple example it is hard to create a not convex hull, but this is not guaranteed, as far as I understand the documentation.

Musa has worked on this problem for 15 days now and mentioned, that this is a part of a Monte-Carlo-simulation. Therefore I think posting a working solution is fair. I've suggested John D'Erricio's constructive approach in another thread of Musa, but I admit that its implementation is not trivial. See here for an implementation: https://www.researchgate.net/publication/275348534_Picking_a_Uniformly_Random_Point_from_an_Arbitrary_Simplex ("Picking a Uniformly Random Point from an Arbitrary Simplex", Christian Grimme). For productive work, I'd prefer a constructive method, but maybe the costs choosing the points can be neglected.

The readers of this thread will find 3 votes for your alphaShape approach currently (one is from me). This is a clear statement about "what is te advantage of writing new code?"

[Jon]

Hi Jan,

Thanks for your explanation for why it is worth writing your own code for this. Based on your remarks regarding the possible non-convexity of the alpha shape, I looked further at the details of the alpha shape documentation. Seeing the possible pitfalls, and the considerable conceptual complexity of alpha shapes I think they may infact be overkill for this simple geometry. So I agree, it may be better to write a little, simple, easy to understand, and failsafe code like yours, when all of the underlying complexity of the two line built in solution is not needed. So, I go with your solution too.

As you point out, depending on the application if you really want to generate points inside of the pyramid, then a constructive approach would be very nice.

[Jun Liang]

Dear Jan,

Thanks for your code.

Can your code calculate a non-regular pyramid, like this?

IIf not, how do I change it?

Thanks you very much.

[Jan]

I do not see, how this is a pyramid.

[Jun Liang]

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.

Answer 2

[Jon]

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

[Jon]

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.

Some approaches and (non MATLAB) implementations are given in https://stackoverflow.com/questions/25179693/how-to-check-whether-the-point-is-in-the-tetrahedron-or-not

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.

Answer 3

[Torsten]

Another program for your issue:

https://de.mathworks.com/matlabcentral/fileexchange/10226-inhull