Mapping a 3D point inside a hexahedron to a unit cube

Question

Let A_4x8 be a hexahedron defined by 8 points in 3D homogeneous coordinates.

M_8x8 is a map that transforms the hexahedron to a unit cube, B_4x8.

If A = B . M, then

M = B^-1 . A

How can I map a 3D point, P_4x1 in A to B using M?

---

Example:

A = \begin{bmatrix} -0.909198 & -0.717041 & -0.745330 & -0.938879 & -1.435481 & -1.196454 & -1.231644 & -1.472402 \\ 0.867548 & 0.893142 & 0.633169 & 0.615845 & 1.071125 & 1.102962 & 0.779577 & 0.758028 \\ 2.259478 & 2.391340 & 2.435814 & 2.300889 & 2.942843 & 3.106868 & 3.162191 & 2.994355 \\ 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 \end{bmatrix}

B = \begin{bmatrix} 0.0 & 1.0 & 1.0 & 0.0 & 0.0 & 1.0 & 1.0 & 0.0\\ 1.0 & 1.0 & 0.0 & 0.0 & 1.0 & 1.0 & 0.0 & 0.0\\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 1.0 & 1.0 & 1.0\\ 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 & 1.0 \end{bmatrix}

P = \begin{bmatrix} -1.052645\\ 0.833213\\ 2.661413\\ 1.0 \end{bmatrix}

Answer 1

Let $A_i (x_i,y_i,z_i)^T (0 \leq i \leq 7) $ the 8 vertices of the hexahedron with their 3D coordinates (we don't need the 4th projective coordinate).

I don't need to assume that the hexahedron is a regular pyramid frutrum.

If the unit cube is defined by coordinates $\{0,1\}^3$, a convenient mapping is as follows :

$$\tag{1}(u,v,w) \in [0,1]^3 \ \ \mapsto \ \ \left\{\begin{array}{rl}(1-u)(1-v)(1-w) & \hspace{-0.5em}A_0 \ \ +\\ (1-u)(1-v)w & \hspace{-0.5em} A_1 \ \ +\\ (1-u)v(1-w) &\hspace{-0.5em}A_2 \ \ +\\ (1-u)vw & \hspace{-0.5em}A_3 \ \ +\\ u(1-v)(1-w)& \hspace{-0.5em}A_4 \ \ +\\ u(1-v)w& \hspace{-0.5em}A_5 \ \ +\\ uv(1-w)& \hspace{-0.5em}A_6 \ \ +\\ uvw&\hspace{-0.5em}A_7\end{array}\right.$$

Short explanation: 2 "ingredients":

1) Boolean logic (even if $u,v,w$ are reals): for example, vertex $A_6$ (and only this vertex) is "triggered" iff $(u,v,w)=(1,1,0)$ which is the binary representation of number 6.

2) Convex hull building: if the hexahedron is convex, the fact that $u,v,w \in [0,1]$ implies that each coefficient, from $(1-u)(1-v)(1-w)$ to $uvw$ is in $[0,1]$. Morever, one can prove that (1) is a bijection, because each point can be described in this way.

Remark: Of course, an important simplication occurs in the case of a pyramid frustrum.

---

Edit: How to build a function in the other way : Frustrum $\mapsto$ Cube.

As it is likely that you know projective transforms, here is a solution that uses them.

Let $B_1,B_2,B_3,B_4$ the basis of your frustrum. Let $S$ be the summit (apex) of the pyramid.

The projective matrix $P$ that maps Frustrum $\mapsto$ Cube is such that:

$$\underbrace{\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 & 1 & 1 & 0\\ 0 & 0 & 1 & 1\\ 1 & 1 & 1 & 0 \end{bmatrix}}_C=P\underbrace{\begin{bmatrix} x_1 & x_2 & x_3 & x_S\\ y_1 & y_2 & y_3 & y_S\\ z_1 & z_2 & z_3 & z_S\\ 1 & 1 & 1 & 1 \end{bmatrix}}_F$$

In particular $P$ maps $\begin{bmatrix} x_S\\ y_S\\ z_S\\ 1 \end{bmatrix}$ onto $\begin{bmatrix} 0\\ 0\\ 1\\ 0 \end{bmatrix}$ (apex is sent to the point at infinity - last coordinate $0$ - in the right direction ; this is what could be called "the opening of the cube")

Thus, one has to compute:

$$P=CF^{-1}$$

Let $P_{ij}$ be the entries of $P$ ; then the transformation equations are:

$$\begin{cases}X=\frac{P_{11}x+P_{12}y+P_{13}z+P_{14}}{P_{41}x+P_{42}y+P_{43}z+P_{44}}\\Y=\frac{P_{21}x+P_{22}y+P_{23}z+P_{24}}{P_{41}x+P_{42}y+P_{43}z+P_{44}}\\Z=\frac{P_{31}x+P_{32}y+P_{33}z+P_{34}}{P_{41}x+P_{42}y+P_{43}z+P_{44}}\\\end{cases}$$

Remark: you will probably need to take $\begin{bmatrix} kx_S\\ ky_S\\ kz_S\\ k \end{bmatrix}$ with a certain coefficient $k$ instead of $\begin{bmatrix} x_S\\ y_S\\ z_S\\ 1 \end{bmatrix}$ in order to take into account the heigth where the cut is done in the pyramid.

Here is a Matlab program that can help you:

clear all;

C = [0     1     1     0
     0     0     1     0
     0     0     0     1
     1     1     1     0];

H = 5;    %height of pyramid
k = 1.5;  % coefficient k has been obtained by trial and error
          % in order that(as an example) the frustrum's height is exactly 3 
x1 = 3;  y1 = 3;  z1 = 0;
x2 = 5;  y2 = 3;  z2 = 0;
x3 = 5;  y3 = 5;  z3 = 0;
xS = 4;  yS = 4;  zS = H;

k = 1;

F = [x1 x2 x3 k*xS
     y1 y2 y3 k*yS
     z1 z2 z3 k*H
      1  1  1 k];

P = C * F^(-1);
X = @(x,y,z)(P(1,:)*[x,y,z,1]');
Y = @(x,y,z)(P(2,:)*[x,y,z,1]');
Z = @(x,y,z)(P(3,:)*[x,y,z,1]');
T = @(x,y,z)(P(4,:)*[x,y,z,1]');

% k = Z(xS,yS,h) ./ T(xS,yS,h);
% x=4;y=4;z=0; %is send onto point 0.5 0.5 0
% x=5;y=5;z=0; %is send onto point 1 1 0
% x=4;y=4;z=3; %is send onto point 0.5 0.5 1

Xp = X(x,y,z) ./ T(x,y,z);
Yp = Y(x,y,z) ./ T(x,y,z);
Zp = Z(x,y,z) ./ T(x,y,z);

fprintf('The image of %s %s %s is %s %s %s\n',...
num2str(x),num2str(y),num2str(z),num2str(Xp),num2str(Yp),num2str(Zp));