\begin{equation}\label{eq:1}P - A = v(B - A) + w(C - A)\end{equation}\begin{equation}\label{eq:2}P = A + v(B - A) + w(C - A) = (1 - v - w)A + vB + wC = uA + vB + wC\end{equation}
Con $u = 1 - v - w$ (e quindi $u + v + w = 1$). Se, oltre a ciò, si aggiunge anche il vincolo $0\le v, w \le 1$, allora $P$ risiede all'interno di $\triangle ABC$.
L'ultimo membro a destra di \eqref{eq:2}, combinazione lineare dei tre punti $A, B$ e $C$, definisce le coordinate baricentriche ($u, v$ e $w$) di $P$ rispetto a $\triangle ABC$. Indicando le coordinate di $A, B$ e $C$ con $A = (a_x, a_y), B = (b_x, b_y), C = (c_x, c_y)$, la \eqref{eq:2} in formato matriciale diventa
\begin{align}
a_xu + b_xv + c_xw &= p_x \notag \\
a_yu + b_yv + c_yw &= p_y \notag \\
u + v + w &= 1 \notag
\end{align}
Quindi per trovare le coordinate baricentriche di $P$ basta risolvere tale sistema. Si potrebbe, ad esempio, usare la regola di Cramer che afferma che, per trovare l'$i$-esima incognita, basta calcolare il rapporto tra il determinante della matrice $\mathbf{A_i}$ (formata dai coefficienti del sistema tranne che per l'$i$-esima colonna, i cui valori sono presi dal vettore dei termini noti) e il determinante della matrice dei coefficienti $\mathbf{A}$.
\[x_i = \frac{det(\mathbf{A_i})}{det(\mathbf{A})}\]
Quindi, per trovare $u, v$ e $w$
\[u = \frac{\begin{vmatrix} p_x & b_x & c_x\\ p_y & b_y & c_y\\ 1 & 1 & 1\\ \end{vmatrix}}{\begin{vmatrix} a_x & b_x & c_x\\ a_y & b_y & c_y\\ 1 & 1 & 1\\ \end{vmatrix}} \quad v = \frac{\begin{vmatrix} a_x & p_x & c_x\\ a_y & p_y & c_y\\ 1 & 1 & 1\\ \end{vmatrix}}{\begin{vmatrix} a_x & b_x & c_x\\ a_y & b_y & c_y\\ 1 & 1 & 1\\ \end{vmatrix}} \quad w = \frac{\begin{vmatrix} a_x & b_x & p_x\\ a_y & b_y & p_y\\ 1 & 1 & 1\\ \end{vmatrix}}{\begin{vmatrix} a_x & b_x & c_x\\ a_y & b_y & c_y\\ 1 & 1 & 1\\ \end{vmatrix}} = 1 - u - v\]
Considerando quanto visto in [2], il determinante di matrici che hanno per elementi le coordinate omogenee di tre punti indica l'area con segno del parallelogramma individuato dai vettori individuati, a loro volta, dai tre punti. Dall'aritmetica si sa che moltiplicando numeratore e denominatore per una stessa costante il rapporto non cambia. Essendo l'area del parallelogramma il doppio di quella del triangolo, si può affermare che, in generale, le coordinate baricentriche in 2D sono date dal rapporto tra funzioni proporzionali alle aree di triangoli. Nello specifico si ha che
\begin{equation}\label{eq:3}
u = \frac{x*\!Area(\triangle PBC)}{x*\!Area(\triangle ABC)} \quad
v = \frac{x*\!Area(\triangle PCA)}{x*\!Area(\triangle ABC)} \quad
w = \frac{x*\!Area(\triangle PAB)}{x*\!Area(\triangle ABC)} = 1 - u - v
\end{equation}
Dunque, è relativamente semplice trovare le coordinate di un punto $P$ in 2D (si veda [5] per formula e codice).
Spostandosi nello spazio 3D le cose non cambiano di molto ma, per giustificare questa affermazione, è oppurtuna una breve digressione sul rapporto tra area di un parallelogramma individuata da due vettori, prodotto vettoriale e prodotto scalare ortogonale (perp dot product), quest'ultimo definito come $\mathbf{u^{\perp}}\cdot \mathbf{v} = (-u_y, u_x) \cdot (v_x, v_y) = -u_yv_x + u_xv_y = u_xv_y - u_yv_x$ (con $\mathbf{u^{\perp}}$ vettore ortogonale a $\mathbf{u}$). Come è noto (si veda ancora [2] per maggiori dettagli), il modulo del prodotto vettoriale tra due vettori $\mathbf{u}$ e $\mathbf{v}$ ($|\mathbf{u}\times \mathbf{v}|$) indica l'area del parallelogramma individuato dai due vettori. La forma componente per componente del prodotto vettoriale si può ottenere tramite determinante della matrice delle componenti di $\mathbf{u}$ e $\mathbf{v}$ (si veda \eqref{eq:5}). Proiettando ortogonalmente questi vettori sul piano $xy$ le rispettive componenti z saranno nulle e il calcolo del determinante, di fatto, si rivela essere proprio il prodotto scalare ortogonale tra $\mathbf{u}$ e $\mathbf{v}$ (si vedano \eqref{eq:4} e \eqref{eq:6}). Quindi anche il prodotto scalare ortogonale, come il prodotto vettoriale, indica l'area di un parallelogramma. Guardano più attentamente \eqref{eq:5} infatti, si può notare come i coefficienti di $\mathbf{i}, \mathbf{j}$ e $\mathbf{k}$ siano proprio i prodotti scalari ortogonali delle proiezioni ortogonali dei vettori 3D $\mathbf{u}$ e $\mathbf{v}$ sui piani $zy$, $zx$ e $xy$, rispettivamente, ed indicano quindi le aree dei parallelogrammi individuati dalle proiezioni ortogonali di $\mathbf{u}$ e $\mathbf{v}$ su tali piani.
\begin{equation}\label{eq:4}\mathbf{u^{\perp}}\cdot \mathbf{v} = (-u_y, u_x) \cdot (v_x, v_y) = -u_yv_x + u_xv_y = u_xv_y - u_yv_x\end{equation}
\begin{equation}\label{eq:5}\mathbf{u} \times \mathbf{v} =\begin{vmatrix} i & j & k\\ u_x & u_y & u_z\\ v_x & v_y & v_z\\ \end{vmatrix} = (u_yv_z - u_zv_y)\mathbf{i} - (u_xv_z - u_zv_x)\mathbf{j} + (u_xv_y - u_yv_x)\mathbf{k}\end{equation}\begin{equation}\label{eq:6}\mathbf{u} \times \mathbf{v} =\begin{vmatrix} i & j & k\\ u_x & u_y & 0\\ v_x & v_y & 0\\ \end{vmatrix} = (u_xv_y - u_yv_x)\mathbf{k}\end{equation}
Capire quanto tutto ciò possa tornare utile ai fini del calcolo delle coordinate baricentriche in 3D risulta evidente quando si considera il fatto che le coordinate baricentriche sono invarianti rispetto alla proiezione ortogonale. Se fosse così si potrebbe ridurre il caso 3D a quello 2D effettuando la proiezione di $\triangle ABC$ su uno dei piani tra $zy$, $zx$ e $xy$ e calcolando le coordinate baricentriche con \eqref{eq:3}. In effetti, come è noto (si veda [4]), la proiezione ortogonale preserva non solo il parallelismo ma anche il rapporto tra le lunghezze. Guardando la figura sotto, le basi e le altezze proiettate di $\triangle ABC$ e degli altri triangoli con vertice in $P$ sul piano $zy$ possono cambiare ma il loro rapporto non cambia. Per questo, anche le aree proiettate potrebbero essere diverse ma, essendo legate alle misure di basi e altezze, i loro rapporti non cambiano. Dunque, usare la \eqref{eq:3} anche per il caso 3D è possibile semplicemente prendendo in considerazione la proiezione ortogonale del triangolo su un piano.
Per quanto riguarda l'implementazione in codice C/C++, per evitare situazioni degeneri è consigliabile proiettare il triangolo sul piano che offre aree più grandi. Per questo basta trovare la componente più grande (in valore assoluto) della normale al triangolo. Se, ad esempio, la componente $x$ è quella maggiore in valore assoluto, allora si proietta sul piano $zy$. Se è maggiore la componente $y$ allora si proietta sul piano $zx$. Altrimenti sul piano $xy$. La normale al piano non è normalizzata dato che, per l'utilizzo che se ne fa, non è necessario. Per il suo calcolo si usa quanto visto in [3] (si vedano (1) e (2)). Per il calcolo delle (doppie) aree dei triangoli proiettati che hanno vertice in $P$ si usa quanto visto in [2] e [5] mentre per l'area (doppia) della proiezione di $\triangle ABC$ si sfrutta quanto visto in questo articolo riguardo la relazione tra componenti del prodotto vettoriale ed area del parallelogramma proiettato sui piani $zy, zx$ e $xy$. A questo proposito si noti come l'area (doppia) della proiezione sul piano $zx$ venga invertita di segno a causa del segno meno prima del coefficiente di $\mathbf{j}$ nel prodotto vettoriale (vedi (1) in [3]). Infine, si fornisce una funzione (di facile implementazione per quanto detto all'inizio dell'articolo) per verificare se un punto è all'interno o meno di un triangolo.
#include <iostream> #include <cmath> typedef float Coord2D[2]; typedef float Coord3D[3]; struct Vertex3D { Coord3D pos; }; struct Vector3D { Coord3D pos; }; void barycentric(const Vertex3D&, const Vertex3D&, const Vertex3D&, const Vertex3D&, float&, float&, float&); float area2(const Coord2D, const Coord2D, const Coord2D); bool in_triangle(const Vertex3D&, const Vertex3D&, const Vertex3D&, const Vertex3D&); int main(int argc, char* argv[], char* env[]) { // A, B, C Vertex3D coords[] = { {-1.0, -1.0, 0.0}, {1.0, -1.0, 0.0}, {0.0, 1.0, 0.0} }; std::cout << (in_triangle(coords[0], coords[1], coords[2], { 0.0, 0.0, 0.0 }) ? "punto si trova in triangolo" : "punto non si trova in triangolo") << std::endl; std::cout << (in_triangle(coords[0], coords[1], coords[2], { 0.0, 2.0, 0.0 }) ? "punto si trova in triangolo" : "punto non si trova in triangolo") << std::endl; return 0; } void barycentric( const Vertex3D &a, const Vertex3D &b, const Vertex3D &c, const Vertex3D &p, float &u, float &v, float &w) { // normale al piano (vedi (1) e (2) in [3]) Vector3D norm = { (b.pos[1] - a.pos[1]) * (c.pos[2] - a.pos[2]) - (c.pos[1] - a.pos[1]) * (b.pos[2] - a.pos[2]), (c.pos[0] - a.pos[0]) * (b.pos[2] - a.pos[2]) - (b.pos[0] - a.pos[0]) * (c.pos[2] - a.pos[2]), (b.pos[0] - a.pos[0]) * (c.pos[1] - a.pos[1]) - (c.pos[0] - a.pos[0]) * (b.pos[1] - a.pos[1]) }; float nu, nv, duv; float x = fabs(norm.pos[0]); float y = fabs(norm.pos[1]); float z = fabs(norm.pos[2]); if (x >= y && x >= z) { Coord2D p2 = { p.pos[1], p.pos[2] }; Coord2D a2 = { a.pos[1], a.pos[2] }; Coord2D b2 = { b.pos[1], b.pos[2] }; Coord2D c2 = { c.pos[1], c.pos[2] }; nu = area2(p2, b2, c2); nv = area2(p2, c2, a2); duv = 1.0f / norm.pos[0]; } else if (y >= x && y >= z) { Coord2D p2 = { p.pos[0], p.pos[2] }; Coord2D a2 = { a.pos[0], a.pos[2] }; Coord2D b2 = { b.pos[0], b.pos[2] }; Coord2D c2 = { c.pos[0], c.pos[2] }; nu = area2(p2, b2, c2); nv = area2(p2, c2, a2); // inverte segno meno prima di coefficiente versore j (vedi (1) in [3]) duv = 1.0f / -norm.pos[1]; } else { Coord2D p2 = { p.pos[0], p.pos[1] }; Coord2D a2 = { a.pos[0], a.pos[1] }; Coord2D b2 = { b.pos[0], b.pos[1] }; Coord2D c2 = { c.pos[0], c.pos[1] }; nu = area2(p2, b2, c2); nv = area2(p2, c2, a2); duv = 1.0f / norm.pos[2]; } u = nu * duv; v = nv * duv; w = 1.0f - u - v; } float area2(const Coord2D a, const Coord2D b, const Coord2D c) { return (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); } bool in_triangle(const Vertex3D &a, const Vertex3D &b, const Vertex3D &c, const Vertex3D &p) { float u, v, w; barycentric(a, b, c, p, u, v, w); // non è necessario verificare anche u dato che è calcolata come 1-v-w return v >= 0.0f && w >= 0.0f && (v + w) <= 1.0f; }
Riferimenti
[1] Real-time Collision Detection - Ericson
[2] Collinearità di tre punti in 2D
[3] Collinearità e Complanarità in 3D
[4] http://jwilson.coe.uga.edu/MATH7200/Sect4.1.html
[5] Area di un poligono
Nessun commento:
Posta un commento