How to determine whether point lies within pyramid volume

Asked by Musa Mustapha on 26 Jun 2019

Latest activity Edited by Jon on 27 Jun 2019

54 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.

0 Comments

3 Answers

Link

Answer by Jan on 27 Jun 2019

Edited by Jan on 27 Jun 2019

Accepted Answer

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

3 Comments

Jon on 27 Jun 2019

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 on 27 Jun 2019

@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 on 27 Jun 2019

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.

Link

Answer by Jon on 26 Jun 2019

Edited by 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 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.

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.

Link

Answer by Torsten on 27 Jun 2019

Another program for your issue:

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

You are now following this question

Musa Mustapha (view profile)