How to define if a point is inside or outside a 3D pyramidal area?
12 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
2 Commenti
Risposte (4)
Karim
il 7 Lug 2022
Modificato: Karim
il 7 Lug 2022
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)
6 Commenti
Karim
il 11 Lug 2022
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)
Jan
il 8 Lug 2022
Modificato: Jan
il 9 Lug 2022
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
Matt J
il 8 Lug 2022
Modificato: Matt J
il 8 Lug 2022
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;
2 Commenti
Matt J
il 11 Lug 2022
Modificato: Matt J
il 11 Lug 2022
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.
Matt J
il 8 Lug 2022
Modificato: Matt J
il 11 Lug 2022
Download this FEX submission,
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);
1 Commento
Matt J
il 11 Lug 2022
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.
Vedere anche
Categorie
Scopri di più su MATLAB Support Package for USB Webcams in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!