extracting voxels within 8 points / vertices

here's the thing: I have 8 points / vertices in a 3D volume. I want to know if a point x is contained within the box whose vertices are given by the 8 points.

Moreover, I know that 4 points lie on a plane parallel (let's call it A) to the y axes, and the other 4 points also lie on another plane parallel (B) to the y axes.

So given A and B, which are the 4x3 matrices of the vertices.

I compute the lines joining the two pair of vertices with lie on the boundary hyperplane of the box

    LinesJoiningVertices = 
            [A(2,:) - A(1,:), B(1,:) - A(1,:);
             A(1,:) - A(1,:), A(4,:) - A(1,:);
             A(3,:) - A(4,:), B(4,:) - A(4,:);
             A(2,:) - B(2,:), B(3,:) - B(2,:);
             B(2,:) - B(1,:), B(4,:) - B(1,:);
             A(1,:) - B(1,:), B(4,:) - B(1,:)]

and I compute the normals of each boundary hyperplane

   for i = 1:6
       normals(i,:) = cross(LinesJoiningVertices(i,1:3), LinesJoiningVertices(i,4:6))
   end

Theoretically, a point x within the box should have a negative dot product with each normal (shouldn't it?). Practically, it doesn't work. I take a point x I know to be within the box, and the answer is greater that 0 for the first, the third and the fifth normal.

Is there an error in my code? Is there an error in my math?

Solution

FIRST OFF, let me tell you something.

With two parallel planes any three dimentional polygon containing two subsets of four points lying on one of each of these two planes can be identified with 6 points only, each 3 of them can represent the plane wich contains one of this polygone's faces.

for any 6 points, {A_0,A_2,A_3} on the plane parallel to the plane constitued by {B_0,B_1,B_2}, a condition is needed for the remaining two points that allows a polygon, or else, we can be faced to any inconsistently distributed coordinates which forms a random non-linear shape or a circle or anything. this illustration shows which two-by-two same-colored lines are bound to be parallel to give the remaining points as an intersection of these lines.

Sencondo, lets proceed to calculate given last informations as a background.

As a beginning, the points' coordinates can be given either generated by the system:

syms x y z
P = [x,y,z]
left_x=0;right_x=4;

for j=0:1, for k=0:1, A((j)*2+k+1,1:3)=[left_x;floor(10*rand());floor(10*rand())];end,end
for j=0:1, for k=0:1, B((j)*2+k+1,1:3)=[right_x;floor(10*rand());floor(10*rand())];end,end

or, typed manually

syms x y z
P = [x,y,z]


left_x=input('enter first x ');
for j=0:1, for k=0:1*~j,  A((j)*2+k+1,1:3)=[left_x;input(['print y' char((j)*2+k+49) ' ']);input(['print z' char((j)*2+k+49) ' '])];end,end,
right_x=input('enter second x ');
for j=0:1, for k=0:1*~j,  B((j)*2+k+1,1:3)=[right_x;input(['print y' char((j)*2+k+49) ' ']);input(['print z' char((j)*2+k+49) ' '])];end,end,

2- now we calculate which of these points are nearer to the y-axix (closer), which are nearer to x-axix (lower), to find equations of 4 virtual planes constituted by each 3 points.

close_right_point=B((B(:,3)==min(B(:,3))),:)

far_right_point=B((B(:,3)==max(B(:,3))),:)

if(numel(close_right_point(:,1))==1)
temp=B((B(:,3)==min(B(B(:,3)~=min(B(:,3)),3))) ,:);
middler_points(2,:)=temp(1,:);
else
middler_points(2,:)=close_right_point(2,:);
end




closer_point=middler_points(middler_points(:,3)==min(middler_points(1,3),middler_points(2,3)),:)
if(numel(closer_point(:,1))==1)
further_point=middler_points(middler_points(:,3)==max(middler_points(1,3),middler_points(2,3)),:)
else
further_point=closer_point(2,:)
end


close_normal = cross(closer_point(1,:)-close_right_point(1,:), closer_point(1,:)-close_left_point(1,:))

close_plane=dot(close_normal,P-closer_point)



far_normal = cross(further_point-far_right_point(1,:), further_point-far_left_point(1,:))

far_plane=dot(far_normal,P-further_point)









low_left_point=A((A(:,2)==min(A(:,2))),:)


if(numel(low_left_point(:,1))==1)
temp=A((A(:,2)==min(A(A(:,2)~=min(A(:,2)),2))),:);
average_points(1,:)=temp(1,:);
else
average_points(1,:)=low_left_point(2,:);
end

high_left_point=A((A(:,2)==max(A(:,2))),:)




low_right_point=B((B(:,2)==min(B(:,2))),:)

high_right_point=B((B(:,2)==max(B(:,2))),:)


if(numel(low_right_point(:,1))==1)
temp=B((B(:,2)==min(B(B(:,2)~=min(B(:,2)),2))),:);
average_points(2,:)=temp(1,:);
else
average_points(2,:)=low_right_point(2,:);
end



lower_point=average_points(average_points(:,2)==min(average_points(1,2),average_points(2,2)),:)
if(numel(lower_point(:,1))==1)
higher_point=average_points(average_points(:,2)==max(average_points(1,2),average_points(2,2)),:)
else
higher_point=lower_point(2,:);
end



low_normal = cross( lower_point(1,:)-low_right_point(1,:), lower_point(1,:)-low_left_point(1,:))

low_plane=dot(low_normal,P-lower_point)



high_normal = cross( higher_point-high_right_point(1,:), higher_point-high_left_point(1,:))

high_plane=dot(high_normal,P-higher_point)

note that, normals are perpendicular vectors to the plane they represent, the equation of any plane is extracted from the normal vector.

this is how the diagram is shown, see the parallel planes and 6 dynamic points that sufficiently constitue any 3D polygone.

3- having almost Finished with that, keeping on to next stage, where we receive coordinates of a desired point by user-input:

X=input('enter X ')
Y=input('enter Y ')
Z=input('enter Z ')

then, checking out if the point dont go beyond the extreme right, or left, relating to the two parallel planes with y-axis

if (X>right_x)
'out-bounds'
return;
end


if (X<left_x)
'out-bounds'
return;
end

then, after drawing a virtual line crossing this point parallelising y-axis, we make sure that this point stays always between the segment formed by this line's intersection with the higher and lower virual planes as this simulation shows

if Z> str2num(char(regexp(evalc(['x=' num2str(X) ';y=' num2str(Y) ';' char(solve(far_plane,z))]),'(\d+.\d*)','match')))
'out-bounds'
return;
end


if Z< str2num(char(regexp(evalc(['x=' num2str(X) ';y=' num2str(Y) ';' char(solve(close_plane,z))]),'(\d+.\d*)','match')))
'out-bounds'
return;
end


if Y> str2num(char(regexp(evalc(['x=' num2str(X) ';z=' num2str(Z) ';' char(solve(high_plane,y))]),'(\d+.\d*)','match')))
'out-of-bounds'
return;
end


if Y< str2num(char(regexp(evalc(['x=' num2str(X) ';z=' num2str(Z) ';' char(solve(low_plane,y))]),'(\d+.\d*)','match')))
'out-of-bounds'
return;
end

4- Finalising all this, you can check the validity of this solution using this data-input set

then floatting point

1
2
3

success!, you can see/modify/check results using this simulation by clicking on any point, then change definition text-field with any parameter, you can drag points too.