How to define if a point is inside or outside a 3D pyramidal area?

Hi Dears,

In a 3D dimensions, I'm working on a camera placement and coverage project.

I'm looking for calculating or simply defining (function or formula) if a point is being covered or not by the camera Field Of View.

The camera Field Of View is supposed to be like a Pyramid form with a rectangular base such as we could have the horizontal perimeter greater than the vertical perimeter.

We can consider that the two angles (AFB adn AFD) that forms the vertical and horizontal base are supposed to be known.

The distance between the for extremities and the top of the pyramid (FA, FB, FC, FD) is also known.

F is supposed to be the center of the camera.

FG = to be defined.

Thanks in advance for the support and any suggestions.

If any other information are requested, I'll be happy to share it.

(This is an optimization problem)

Regards.

Sami

[Jan]

Easier to read:

What have you tried so far?

Obtain the height of G over the plane ABCD. If it is inside [0, Fz]: Construct A'B'C'D' in this height and check if G is inside the rectangle. This is cheaper than calculating the view angles.

[Sami Zemma]

Hi Jan,

Thank you for the suggestion and for image corrections :). Nice to find some helpers.

I think the answer @Karim is pretty much easier and exploitable :).

Now that I know the "pointlocation" function, can I use it in a defined area (like a room or imported image) with more thane one camera ? it means more than one Fiel Of View.

Answer 1

[Karim]

Hi, asuming you want to use as many standard matlab functions, the easiest method i can think is by constructing two terahedronds that form the pyramid. Then you can use the function "pointlocation" to know if a point is inside or not. See below for a demonstration.

% gather the points in a grid
MyGrid = [4 0 0; 4 4 0; 2 4 0; 2 0 0; 3 3 3];

% split the pyramid into two tetrahedrons, and create the triangulation
TR = triangulation([1 2 4 5; 2 3 4 5], MyGrid);

% plot the pyarmid and its triangulation
figure
subplot(1,2,1)
patch('Faces', [1 2 3 4; 1 2 5 nan; 2 3 5 nan; 3 4 5 nan; 4 1 5 nan],'Vertices',MyGrid,'FaceColor','r','FaceAlpha',0.5)
title('Pyramid view')
axis equal
grid on
view(3)

subplot(1,2,2)
tetramesh(TR)
title('Triangulation')
axis equal
grid on
view(3)

% now create a bunch of random points in 3D
MyPoints = [rand(50,1)*3+1 rand(50,1)*6-2 rand(50,1)*4-1];

% check if points are inside or not
ImInside = ~isnan(pointLocation(TR,MyPoints));

% plot patch, points that are indside in green and outside in red
figure
hold on
patch('Faces', [1 2 3 4; 1 2 5 nan; 2 3 5 nan; 3 4 5 nan; 4 1 5 nan],'Vertices',MyGrid,'FaceColor','y','FaceAlpha',0.25)
scatter3(MyPoints(ImInside,1),MyPoints(ImInside,2),MyPoints(ImInside,3),'g','filled')
scatter3(MyPoints(~ImInside,1),MyPoints(~ImInside,2),MyPoints(~ImInside,3),'r','filled')
hold off
axis equal
grid on
view(3)

[Sami Zemma]

Hi @Karim,

Many thanks for the answer. It really helped me.

Now that I know the "pointlocation" function, can I use it in a defined area (like a room or imported image) with more thane one camera ? it means more than one Fiel Of View.

Regards.

Sami

[Karim]

Yes, just include as many tetrahedronds in the triangulation as you want.

Note you can evaluate as many points as you want at the same time, the "pointlocation" function is extremly efficient. I have used it on very large datasets.

[Karim]

@Sami Zemma does one of the solutions fit your need? It's usefull to "accept" an answer so that it might be helpfull for other people looking in the answer section.

[Jan]

@Sami Zemma: Or wait a little bit with accepting, because the different solutions are interesting. I do learn from it.

[Sami Zemma]

Hello Karim,

Actually, I'm trying to adapt the solution by creating view cameras initially with a known position in a defined area (like a grid of 10m Width x 10m Length x 5m Height).

I also saw the "plotCoverageArea" function which uses a bird's eye sensor. It seems to be good too, because it has the ability to initialize the orientation, which is very important to know.

The link to Bird's Eye plot : https://fr.mathworks.com/help/driving/ref/birdseyeplot.plotcoveragearea.html

But, you should know that I try to define an area (example 10x*10y*5z), then initially place 2 or 3 cameras with a certain orientation and coverage. Then I calculate the percentage of area initially covered by these cameras and finally I have to implement "Metaheuristic Optimization Algorithms" to try to maximize the coverage.

So, after implementing the optimization algorithms, I will have to recalculate the camera poses that will have actually changed and recalculate the covered area that will have normally increased as well.

[Karim]

One method to do this, is by rotating and translating the grid that makes up the triangulation (or directly the points in Jan's solution) see the example below. This would give you 6 design variables per camera, which you can update to obtain the desired coverage.

MyGrid = [4 0 0; 4 4 0; 2 4 0; 2 0 0; 3 3 3];
TR = triangulation([1 2 4 5; 2 3 4 5], MyGrid);

% rotate the grid about the x axis
RotX = @(x) [1 0 0; 0 cos(x) -sin(x); 0 sin(x) cos(x)];
RotGrid = (RotX( pi/6 ) * MyGrid')';
TR_rot = triangulation([1 2 4 5; 2 3 4 5], RotGrid);

% translate the rotated grid
RotTransGrid = RotGrid + repmat( [ 2 -3 3],size(RotGrid,1),1);
TR_rot_trans = triangulation([1 2 4 5; 2 3 4 5], RotTransGrid);

figure
subplot(1,3,1)
tetramesh(TR)
title('Start')
axis equal
grid on
view(3)

subplot(1,3,2)
tetramesh(TR_rot)
title('Rotated')
axis equal
grid on
view(3)

subplot(1,3,3)
tetramesh(TR_rot_trans)
title('Rotated and translated')
axis equal
grid on
view(3)

Answer 2

[Jan]

If you have the 4 points of the pyramid given as 1x3 vectors, A,B,C,D,F (by the way, where is E?):

% [UNDER CONSTRUCTION! RESULT IS BUGGY!]

Tethra1 = [A;D;C;F];
Tethra2 = [C;B;A;F];
Point   = [1,3,2];

isInside = isPointInSimplex(Tetra1, P) | isPointInSimplex(Tetra2, P);

function in = isPointInSimplex(T, P)
% INPUT:  T:  [4 x 3] matrix, coordinates of tetrahedon, 4 points in 3D
%         P:  [N x 3] matrix, coordinates of N points to check in 3D
% OUTPUT: in: TRUE if point is inside the tetrahedon.
%
% (C) Jan, 2022, CC BY-SA 3.0

Ax = T(1, 1);  Ay = T(1, 2);  Az = T(1, 3);
Bx = T(2, 1);  By = T(2, 2);  Bz = T(2, 3);
Cx = T(3, 1);  Cy = T(3, 2);  Cz = T(3, 3);
Dx = T(4, 1);  Dy = T(4, 2);  Dz = T(4, 3);

Px = P(:, 1);  Py = P(:, 2);  Pz = P(:, 3);

b1 = Px*Bz*Cy - Px*By*Cz + Py*Bx*Cz - Py*Cx*Bz - Bx*Pz*Cy + ...
     Pz*By*Cx + Px*By*Dz - Px*Bz*Dy - Py*Bx*Dz + Py*Bz*Dx + ...
     Bx*Pz*Dy - Pz*By*Dx - Px*Cy*Dz + Px*Cz*Dy + Py*Cx*Dz - ...
     Py*Dx*Cz - Pz*Cx*Dy + Pz*Cy*Dx + Bx*Cy*Dz - Bx*Cz*Dy - ...
     By*Cx*Dz + By*Dx*Cz + Cx*Bz*Dy - Bz*Cy*Dx;

b2 = Ax*Pz*Cy - Ax*Py*Cz + Ay*Px*Cz - Ay*Cx*Pz - Px*Az*Cy + ...
     Az*Py*Cx + Ax*Py*Dz - Ax*Pz*Dy - Ay*Px*Dz + Ay*Pz*Dx + ...
     Px*Az*Dy - Az*Py*Dx - Ax*Cy*Dz + Ax*Cz*Dy + Ay*Cx*Dz - ...
     Ay*Dx*Cz - Az*Cx*Dy + Az*Cy*Dx + Px*Cy*Dz - Px*Cz*Dy - ...
     Py*Cx*Dz + Py*Dx*Cz + Cx*Pz*Dy - Pz*Cy*Dx;

b3 = Ax*Bz*Py - Ax*By*Pz + Ay*Bx*Pz - Ay*Px*Bz - Bx*Az*Py + ...
     Az*By*Px + Ax*By*Dz - Ax*Bz*Dy - Ay*Bx*Dz + Ay*Bz*Dx + ...
     Bx*Az*Dy - Az*By*Dx - Ax*Py*Dz + Ax*Pz*Dy + Ay*Px*Dz - ...
     Ay*Dx*Pz - Az*Px*Dy + Az*Py*Dx + Bx*Py*Dz - Bx*Pz*Dy - ...
     By*Px*Dz + By*Dx*Pz + Px*Bz*Dy - Bz*Py*Dx;

b4 = Ax*Bz*Cy - Ax*By*Cz + Ay*Bx*Cz - Ay*Cx*Bz - Bx*Az*Cy + ...
     Az*By*Cx + Ax*By*Pz - Ax*Bz*Py - Ay*Bx*Pz + Ay*Bz*Px + ...
     Bx*Az*Py - Az*By*Px - Ax*Cy*Pz + Ax*Cz*Py + Ay*Cx*Pz - ...
     Ay*Px*Cz - Az*Cx*Py + Az*Cy*Px + Bx*Cy*Pz - Bx*Cz*Py - ...
     By*Cx*Pz + By*Px*Cz + Cx*Bz*Py - Bz*Cy*Px;

% [TO BE EDITED: Depends on order of points in the simplex]
% To exclude points on the surface replace ">=" by ">" :
in = (b1 >= 0) & (b2 >= 0) & (b3 >= 0) & (b4 >= 0);
end

[Sami Zemma]

Hi Jan,

Can we set the initial position of the pyramid randomely ?

Thanks

Answer 3

[Matt J]

Just camera project the point(s) onto the plane of A,B,C,D and use inpolygon to see if it lies between them. If the plane of ABCD will always be the xy plane, like in your example, it would simply be,

ABCD = [4 0 0; 4 4 0; 2 4 0; 2 0 0];
F=[3,3,3];


% now create a bunch of random points in 3D
P = [rand(50,1)*3+1 rand(50,1)*6-2 rand(50,1)*4-1];

t=-F(3)./(P(:,3)-F(3));

Projections= F+t.*(P-F);

isInside=inpolygon(Projections(:,1),Projections(:,2),ABCD(:,1),ABCD(:,2)) & t>=0 &t<=1;

[Sami Zemma]

Hi Matt,

No, the plane can actually change after.

As I was saying to @Karim,

Actually, I'm trying to adapt the solution by creating view cameras initially with a known position in a defined area (like a grid of 10m Width x 10m Length x 5m Height).

I also saw the "plotCoverageArea" function which uses a bird's eye sensor. It seems to be good too, because it has the ability to initialize the orientation, which is very important to know.

The link to Bird's Eye plot : https://fr.mathworks.com/help/driving/ref/birdseyeplot.plotcoveragearea.html

But, you should know that I try to define an area (example 10x*10y*5z), then initially place 2 or 3 cameras with a certain orientation and coverage. Then I calculate the percentage of area initially covered by these cameras and finally I have to implement "Metaheuristic Optimization Algorithms" to try to maximize the coverage.

So, after implementing the optimization algorithms, I will have to recalculate the camera poses that will have actually changed and recalculate the covered area that will have normally increased as well.

Any suggestion will be very helpfull.

Thanks.

[Matt J]

Any suggestion will be very helpfull.

Your additional comments do not change my suggestions if your goal is still to test whether a point lies within the field of a camera. Even if one or more of the cameras do not have ABCD aligned with the xy-plane, you can still project the point into the coordinates of the ABCD plane and follow the same approach.

Answer 4

[Matt J]

Download this FEX submission,

https://www.mathworks.com/matlabcentral/fileexchange/30892-analyze-n-dimensional-convex-polyhedra

and use it to find the linear inequalities defining the pyramid. Then test whether the candidate points satisfy the inequalities.

ABCD = [4 0 0; 4 4 0; 2 4 0; 2 0 0];
F=[3,3,3];

[Aineq,bineq]=vert2lcon([ABCD;F]);

% now create a bunch of random points in 3D
P = [rand(50,1)*3+1 rand(50,1)*6-2 rand(50,1)*4-1];

isInside=all(Aineq*P'<=bineq);

[Matt J]

Then I calculate the percentage of area initially covered by these cameras and finally I have to implement "Metaheuristic Optimization Algorithms" to try to maximize the coverage.

The same FEX submission mentioned above contains intersectionHull(). You can use it to calculate the polyhedron formed by the intersection of the cameras. Then you can use convhull or convhulln to compute its volume.