Asked by Musa Mustapha
on 26 Jun 2019

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 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:

Jon
on 27 Jun 2019

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.

Sign in to comment.

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

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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.