Collinearità e Complanarità in 3D

Collinearità di tre punti in 3D

Verificare se tre punti $A, B, C$ in 3D giacciono sulla stessa retta a livello teorico è del tutto equivalente a quanto accade in 2D (si veda [1]). Si tratta cioè di calcolare l'area del parallelogramma individuato dai due vettori $\,\, \mathbf{n_1} = B - A \,\,$ e $\,\, \mathbf{n_2} = C - A \,$. Se l'area è nulla i tre punti sono collineari.

Com'è noto dall'algebra lineare, l'area di questo parallelogramma può essere calcolata tramite il modulo del prodotto vettoriale $|\mathbf{n_1} \times \mathbf{n_2}|$. Infatti l'area $A_p$ di un parallelogramma è uguale a $A_p = b * h$ e il prodotto vettoriale è definito come $\mathbf{v} \times \mathbf{w} = \hat{\mathbf{n}} |\mathbf{v}||\mathbf{w}|sin\theta$ (con $\hat{\mathbf{n}}$ versore normale al piano in cui giacciono $\mathbf{v}$ e $\mathbf{w}$). Quindi $|\mathbf{v} \times \mathbf{w}| = |\mathbf{v}||\mathbf{w}|sin\theta$ (si ricordi che $|\hat{\mathbf{n}}|=1$). Ma come si può vedere dalla figura sotto $|\mathbf{w}| = b$ e $|\mathbf{v}| sin\theta = h$.

Sempre dall'algebra lineare si sa che il prodotto vettoriale, espresso componente per componente, può essere ottenuto calcolando il determinante della matrice quadrata che ha come elementi sulla prima riga (o colonna) i versori $\mathbf{i}, \mathbf{j}$, $\mathbf{k}$ e sulle ultime due righe (o colonne) le componenti dei due vettori.

\begin{equation}\label{eq:1}\mathbf{n_1} \times \mathbf{n_2} = \begin{vmatrix} i & j & k\\ x_b - x_a & y_b - y_a & z_b - z_a\\ x_c - x_a & y_c - y_a & z_c - z_a\\ \end{vmatrix} = \\ \\ [(y_b - y_a)(z_c - z_a) - (y_c - y_a)(z_b - z_a)] \mathbf{i} - [(x_b - x_a)(z_c - z_a) - (x_c - x_a)(z_b - z_a)] \mathbf{j} + [(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a)] \mathbf{k}\end{equation}

Per ogni vettore $\mathbf{v} = a\mathbf{i} + b\mathbf{j} + c\mathbf{k}$ si ha che il modulo di $\mathbf{v}$ (indicato con $\mathbf{|v|}$) è uguale a $\mathbf{|v|} = \sqrt{a^2 + b^2 + c^2}$, con $a, b$ e $c$ grandezze scalari che indicando la lunghezza, con segno, delle proiezioni di $\mathbf{v}$ lungo i versori $\mathbf{i}, \mathbf{j}$ e $\mathbf{k}$.
Quando $\mathbf{v}$ è il vettore risultante di un prodotto vettoriale il suo modulo (o lunghezza senza segno) è uguale all'area del parallelogramma individuato dai vettori operandi del prodotto vettoriale. Perché l'area sia nulla, quindi, la sua lunghezza deve essere nulla. Per questo, se $\mathbf{|v|} = 0$, $\mathbf{|v|} = \sqrt{a^2 + b^2 + c^2} = 0$ e, di conseguenza, $a = b = c = 0$. Se non fosse così, infatti, la proiezione di $\mathbf{|v|}$ su almeno uno dei versori sarebbe non nulla portando così ad una contraddizione. Infine, quando ci si ritrova in un sistema cartesiano ogni vettore $\mathbf{v} = a\mathbf{i} + b\mathbf{j} + c\mathbf{k}$ può essere identificato univocamente con $\mathbf{v} = (a, b, c)$ con $a$ componente $x$, $b$ componente $y$ e $c$ componente $z$ delle coordinate di $\mathbf{v}$.
Verificare se tre punti giacciono sulla stessa retta equivale quindi a verificare che le componenti $x, y$ e $z$ del prodotto vettoriale dei due vettori individuati dai tre punti siano nulle. Per comprendere l'ordine dei fattori e dei termini della equazione \eqref{eq:2} sotto si noti il segno meno prima del coefficiente di $\mathbf{j}$ sul lato destro in \eqref{eq:1})

\begin{equation}\label{eq:2}(y_b - y_a)(z_c - z_a) - (y_c - y_a)(z_b - z_a) = 0 \\ (x_c - x_a)(z_b - z_a) - (x_b - x_a)(z_c - z_a) = 0 \\ (x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a) = 0\end{equation}



Complanarità di quattro punti in 3D

Verificare se quattro punti $A, B, C, D$ in 3D giacciono sullo stesso piano equivalente a verificare se è nullo il volume del parallelepipedo individuato dai tre vettori

$\mathbf{n_1} = B - A$
$\mathbf{n_2} = C - A$
$\mathbf{n_3} = D - A$

Com'è noto dall'algebra lineare il volume di questo parallelepipedo può essere calcolato tramite il modulo del prodotto misto $|(\mathbf{n_1} \times \mathbf{n_2}) \cdot \mathbf{n_3}|$.
Sempre dall'algebra lineare si sa che il prodotto misto può essere ottenuto calcolando il determinante della matrice quadrata che ha come elementi le componenti dei vettori operandi del prodotto misto.

\[(\mathbf{n_1} \times \mathbf{n_2})\cdot \mathbf{n_3} =\begin{vmatrix} x_b - x_a & y_b - y_a & z_b - z_a\\ x_c - x_a & y_c - y_a & z_c - z_a\\ x_d - x_a & y_d - y_a & z_d - z_a\\ \end{vmatrix} = \\ (x_d - x_a) [(y_b - y_a)(z_c - z_a) - (y_c - y_a)(z_b - z_a)] - \\ (y_d - y_a) [(x_b - x_a)(z_c - z_a) - (x_c - x_a)(z_b - z_a)] + \\ (z_d - z_a) [(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a)]\]

Il risultato del prodotto misto è uno scalare ed il suo modulo rappresenta il volume di un parallelepipedo. Di conseguenza il risultato del prodotto misto rappresenta lo stesso volume ma con segno. Questo sarà positivo se $D$ vede $A, B$ e $C$ in senso orario. Altrimenti sarà negativo. Quindi per verificare se i quattro punti sono complanari basta verificare che il risultato del calcolo del determinante della matrice definita in precedenza è nullo.

\[(x_d - x_a) [(y_b - y_a)(z_c - z_a) - (y_c - y_a)(z_b - z_a)] - \\ (y_d - y_a) [(x_b - x_a)(z_c - z_a) - (x_c - x_a)(z_b - z_a)] + \\ (z_d - z_a) [(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a) =\]
\begin{equation}\label{eq:3}(x_d - x_a) [(y_b - y_a)(z_c - z_a) - (y_c - y_a)(z_b - z_a)] + \\ (y_d - y_a) [(x_c - x_a)(z_b - z_a) - (x_b - x_a)(z_c - z_a)] + \\ (z_d - z_a) [(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a) = 0\end{equation}

Allo stesso risultato si perviene (dopo tediosi calcoli, sviluppando e riarrangiando i termini) calcolando il determinante della matrice che ha come elementi le coordinate omogenee dei quattro punti.

\[\begin{vmatrix} x_a & y_a & z_a & 1\\ x_b & y_b & z_b & 1\\ x_c & y_c & z_c & 1\\ x_d & y_d & z_d & 1 \\ \end{vmatrix} = \\ (x_d - x_a) [(y_b - y_a)(z_c - z_a) - (y_c - y_a)(z_b - z_a)] - \\ (y_d - y_a) [(x_b - x_a)(z_c - z_a) - (x_c - x_a)(z_b - z_a)] + \\ (z_d - z_a) [(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a)]\]

Se invece si fosse interessati al volume del tetraedro individuato dai quattro punti $A, B, C, D$ (piuttosto che al volume del parallelepipedo individuato dai 3 vettori $\mathbf{n_1}, \mathbf{n_2}$ e $\mathbf{n_3}$) si può sfruttare il fatto che il volume del parallelepipedo è 6 volte quello del tetraedro. Infatti, indicando con $A_{\diamond}$ l'area della base del parallelepipedo e con $h$ la sua altezza (definita come proiezione sul versore normale alla base) si ha

$A_{\diamond} = |\mathbf{n_1} \times \mathbf{n_2}|$

$h = \mathbf{n_3}\cdot \displaystyle\frac{\mathbf{n_1} \times \mathbf{n_2}}{|\mathbf{n_1} \times \mathbf{n_2}|}\quad$ (si ricordi che $\mathbf{v} \cdot \mathbf{w} = |\mathbf{v}||\mathbf{w}|cos\theta$ e se $\mathbf{w}$ è normalizzato si ha $\mathbf{|proj_w v|} = \mathbf{|v|}cos\theta = \mathbf{v} \cdot \mathbf{w}$)

e quindi, indicando con $V_p$ e $V_t$ i volumi rispettivamente del parallelepipedo e del tetraedro e ricordando che la piramide a base triangolare è un tetraedro il cui volume è $V = (S_b * h)/3$ (con $S_b$ superficie di base della piramide), si ha

\[V_p = A_{\diamond} * h = |\mathbf{n_1} \times \mathbf{n_2}| * \mathbf{n_3}\cdot \frac{\mathbf{n_1} \times \mathbf{n_2}}{|\mathbf{n_1} \times \mathbf{n_2}|} = \mathbf{n_3}\cdot (\mathbf{n_1} \times \mathbf{n_2})\]
\[V_t = \frac{1}{3}h * \frac{1}{2}A_{\diamond} = \frac{1}{6} (\mathbf{n_1} \times \mathbf{n_2}) \cdot \mathbf{n_3}\]

Quindi basta moltiplicare il prodotto misto per $\frac{1}{6}$ per ottenere il volume del tetraedro.

Per quanto riguarda l'implementazione in codice C/C++ le equazioni di interesse sono \eqref{eq:2} e \eqref{eq:3}. L'unica accortezza è quella richiesta quando si opera con numeri floating point. A causa dei possibili errori di approssimazione nella rappresentazione dei numeri è consigliabile adottare una strategia difensiva e operare su intervalli di tolleranza quando si usano operatori di comparazione. Per maggiori informazioni si vedano [2] e [3].

// se si lavora con interi
bool collinear3D_i(const Vertex_Coord *a, const Vertex_Coord *b, const Vertex_Coord *c)
{
 return (b[1] - a[1]) * (c[2] - a[2]) - (c[1] - a[1]) * (b[2] - a[2]) == 0 &&
  (c[0] - a[0]) * (b[2] - a[2]) - (b[0] - a[0]) * (c[2] - a[2]) == 0 &&
  (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]) == 0;
}
 
bool coplanar(const Vertex_Coord *a, const Vertex_Coord *b, const Vertex_Coord *c, const Vertex_Coord *d)
{
 int vol;
 int dx, dy, dz, bx, by, bz, cx, cy, cz;
 
 dx = d[0] - a[0];
 dy = d[1] - a[1];
 dz = d[2] - a[2];
 bx = b[0] - a[0];
 by = b[1] - a[1];
 bz = b[2] - a[2];
 cx = c[0] - a[0];
 cy = c[1] - a[1];
 cz = c[2] - a[2];
 
 vol = dx * (by*cz - bz * cy) + dy * (bz*cx - bx * cz) + dz * (bx*cy - by * cx);
 
 return vol == 0;
}
 
// se si lavora con numeri reali
bool collinear3D_f(const Vertex_Coord *a, const Vertex_Coord *b, const Vertex_Coord *c)
{
 return AlmostEqual((b[1] - a[1]) * (c[2] - a[2]) - (c[1] - a[1]) * (b[2] - a[2]), 0.0f) &&
  AlmostEqual((c[0] - a[0]) * (b[2] - a[2]) - (b[0] - a[0]) * (c[2] - a[2]), 0.0f) &&
  AlmostEqual((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]), 0.0f);
}
 
inline float max_3(float a, float b, float c)
{
 return ((a > b) ? ((a > c) ? a : c) : (b > c) ? b : c);
}
 
bool AlmostEqual(float a, float b, float epsilon = FLT_EPSILON)
{
 if (fabs(a - b) <= epsilon * max_3(1.0f, fabs(a), fabs(b)))
  return true;
 
 return false;
}
 
bool coplanar(const Vertex_Coord *a, const Vertex_Coord *b, const Vertex_Coord *c, const Vertex_Coord *d)
{
 float vol;
 float dx, dy, dz, bx, by, bz, cx, cy, cz;
 
 dx = d[0] - a[0];
 dy = d[1] - a[1];
 dz = d[2] - a[2];
 bx = b[0] - a[0];
 by = b[1] - a[1];
 bz = b[2] - a[2];
 cx = c[0] - a[0];
 cy = c[1] - a[1];
 cz = c[2] - a[2];
 
 vol = dx * (by*cz - bz * cy) + dy * (bz*cx - bx * cz) + dz * (bx*cy - by * cx);
 
 return AlmostEqual(vol, 0.0f);
}

Riferimenti:
[1] Collinearità di tre punti in 2D
[2] Real-time Collision Detection - Ericson
[3] Comparing Floating Point Numbers, 2012 Edition