Se l'area è nulla, infatti, 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$. Per questo $A_p = |\mathbf{v} \times \mathbf{w}|$
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.
\[\mathbf{n_1} \times \mathbf{n_2} =\begin{vmatrix} i & j & k\\ x_b - x_a & y_b - y_a & 0\\ x_c - x_a & y_c - y_a & 0\\ \end{vmatrix} = [(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a)] \mathbf{k}\]
A questo punto si potrebbe obiettare che il risultato del prodotto vettoriale è un vettore mentre a noi interessa il suo modulo. In effetti il risultato del determinante visto sopra è un vettore dipendente esclusivamente dalla sua componente $z$ (poiché ortogonale ai due vettori usati come operandi del prodotto vettoriale e che giacciono sul piano $xy$) ed il suo modulo, come sappiamo, indica l'area di un parallelogramma. Ma tale modulo altro non è che il modulo della sua componente $z$, e cioè $|(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a)|$.
Infatti, come è noto, per ogni vettore $\mathbf{v} = a\mathbf{i} + b\mathbf{j} + c\mathbf{k}$, $a, b$ e $c$ sono grandezze scalari che indicando la lunghezza, con segno, delle proiezioni di $\mathbf{v}$ lungo i versori $\mathbf{i}, \mathbf{j}$ e $\mathbf{k}$. Se $a$ e $b$ sono entrambi nulli, il modulo di $\mathbf{v}$ è semplicemente la lunghezza, senza segno, della proiezione di $\mathbf{v}$ su $\mathbf{k}$, e cioè $|\mathbf{v}| = |c|$. Quindi, quando questo vettore $\mathbf{v}$ è il risultato di un prodotto vettoriale i cui vettori operandi giacciono sul piano $xy$, $|c|$ (il modulo di $\mathbf{v}$) rappresenta l'area del parallelogramma e, di conseguenza, $c$ la sua area con segno (se i punti che descrivono i due vettori operandi del prodotto vettoriale sono in ordine antiorario l'area avrà segno positivo, altrimenti sarà negativo). 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 quindi se tre punti giacciono sulla stessa retta equivale a verificare che la componente $z$ del prodotto vettoriale dei due vettori individuati dai tre punti sia nulla.
\begin{equation}\label{eq:1}(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a) = 0\end{equation}
Lo stesso risultato può essere ottenuto calcolando il determinante i cui elementi sono le coordinate omogenee dei tre punti. Da notare che l'espressione ottenuta in questo modo è, come in precedenza, la lunghezza con segno della proiezione del vettore risultante dal prodotto vettoriale $(B - A) \times (C - A)$ su $\mathbf{k}$ e, quindi, l'area con segno del parallelogramma individuato dai due vettori operandi del prodotto vettoriale. Basta poco per verificarlo:
\[\begin{vmatrix} x_a & y_a & 1\\ x_b & y_b & 1\\ x_c & y_c & 1\\ \end{vmatrix} =\\ 1(x_{b}y_c - x_{c}y_b) - 1(x_{a}y_c - x_{c}y_a) +1(x_{a}y_b - x_{b}y_a) =\\ x_{b}y_c - x_{c}y_b - x_{a}y_c + x_{c}y_a + x_{a}y_b - x_{b}y_a \color{red}{+ x_{a}y_a - x_{a}y_a} =\\ y_{c}(x_b - y_a) - y_{a}(x_b - x_a) + y_{b}(x_a - x_c) - y_{a}(x_a - x_c) =\\ (x_b - x_a)(y_c - y_a) + (x_a - x_c)(y_b - y_a) =\\ (x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a)\]
In ultimo, allo stesso risultato si perviene anche partendo dalla classica equazione della retta per un punto $A = (x_a, y_a)$
\[y - y_a = m(x - x_a)\]
Tra le infinite rette passanti per $A$, si impone il passaggio per il punto $B = (x_b, y_b)$
\[y_b - y_a = m(x_b - x_a)\]
\[m = \frac{y_b - y_a}{x_b - x_a}\]
Sostituendo $m$ (coefficiente angolare della retta passante per due punti, $A$ e $B$) nell'equazione della retta per $A$
\[y - y_a = \frac{y_b - y_a}{x_b - x_a}(x - x_a)\]
\[\frac{y - y_a}{y_b - y_a} = \frac{x - x_a}{x_b - x_a}\]
Un altro punto $C = (x_c, y_c)$, che passa per la retta passante per $A$ e $B$, deve soddisfare l'equazione sopra e quindi
\[\frac{y_c - y_a}{y_b - y_a} = \frac{x_c - x_a}{x_b - x_a}\]
\[(y_c - y_a)(x_b - x_a) = (x_c - x_a)(y_b - y_a)\]
\[(x_b - x_a)(y_c - y_a) - (x_c - x_a)(y_b - y_a) = 0\]
L'implementazione in codice C/C++ è quindi abbastanza semplice. L'equazione d'interesse è \eqref{eq:1}. 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 [1] e [2].
// se si lavora con interi bool collinear_i(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c) { int t = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); if (t == 0) return true; else return false; } // se si lavora con numeri reali bool collinear_f(const Vertex_Coord a, const Vertex_Coord b, const Vertex_Coord c) { float t = (b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1]); if (AlmostEqual(t, 0.0f)) return true; else return false; } 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; }
Riferimenti:
[1] Real-time Collision Detection - Ericson
[2] Comparing Floating Point Numbers, 2012 Edition
Nessun commento:
Posta un commento