Approssimare funzioni seno e coseno

In un precedente articolo (si veda [7]) si è cercato di fare luce su una delle funzioni più famose e discusse nell'ambito della programmazione di videogiochi. Sarebbe interessante vedere se esistono implementazioni, altrettanto veloci, per approssimare le funzioni seno e coseno dato che queste sono le funzioni che in assoluto, insieme alla radice quadrata, vengono usate più spesso. Un esempio di implementazione efficiente è quello della libreria DirectXMath di Microsoft che fornisce supporto ad applicazioni DirectX:

inline void XMScalarSinCos(float* pSin, float* pCos, float  Value)
{
 assert(pSin);
 assert(pCos);
 
 // Map Value to y in [-pi,pi], x = 2*pi*quotient + remainder.
 float quotient = XM_1DIV2PI * Value;
 
 if (Value >= 0.0f)
 {
  quotient = static_cast<float>(static_cast<int>(quotient + 0.5f));
 }
 else
 {
  quotient = static_cast<float>(static_cast<int>(quotient - 0.5f));
 }
 
 float y = Value - XM_2PI * quotient;
 
 
 // Map y to [-pi/2,pi/2] with sin(y) = sin(Value).
 float sign;
 
 if (y > XM_PIDIV2)
 {
  y = XM_PI - y;
  sign = -1.0f;
 }
 else if (y < -XM_PIDIV2)
 {
  y = -XM_PI - y;
  sign = -1.0f;
 }
 else
 {
  sign = +1.0f;
 }
 
 float y2 = y * y;
 
 // 11-degree minimax approximation
 *pSin = (((((-2.3889859e-08f * y2 + 2.7525562e-06f) * 
    y2 - 0.00019840874f) * y2 + 0.0083333310f) * 
    y2 - 0.16666667f) * y2 + 1.0f) * y;
 
 // 10-degree minimax approximation 
 float p = ((((-2.6051615e-07f * y2 + 2.4760495e-05f) * 
    y2 - 0.0013888378f) * y2 + 0.041666638f) * y2 - 0.5f) * y2 + 1.0f;

 *pCos = sign * p;
}

Rappresentazione Floating Point

Continuo vs Discreto

Su un computer, l'insieme continuo ed infinito dei numeri reali può essere, per ovvi motivi, solo approssimato con un insieme finito e discreto. Questo vuol dire che solo un sottoinsieme molto ristretto di numeri reali è esattamente rappresentabile e memorizzabile in memoria. La stragrande maggioranza dei numeri reali, invece, deve accontentarsi di una buona approssimazione. Quindi, se ad esempio definiamo una variabile e le assegniamo il valore $12.53$ non è detto che questo coincida proprio con uno dei valori esattamente rappresentabili nel sistema in uso. Nella maggior parte dei casi sarà approssimato al valore più vicino esattamente rappresentabile.

Circocentro in 2D e 3D

Il circocentro è il centro del circumcerchio (la circonferenza circoscritta ad un poligono). Se il poligono è un triangolo, questo è il punto di intersezione o degli assi dei suoi lati (figura sotto a sinistra) oppure delle bisettrici dei suoi angoli (figura sotto a destra).

Riferendosi all'immagine a destra e al caso 2D, si definiscano

$A = (x_a, y_a)$,  $B = (x_b, y_b)$,  $C = (x_c, y_c)$, $P = (x_p, y_p)$
$\overrightarrow{AB} = (B - A) = ((x_b - x_a), (y_b - y_a)) = (X_{AB}, Y_{AB})$
$\overrightarrow{AC} = (C - A) = ((x_c - x_a), (y_c - y_a)) = (X_{AC}, Y_{AC})$
$\overrightarrow{AP} = (P - A) = ((x_p - x_a), (y_p - y_a)) = (X_{AP}, Y_{AP})$
$\overrightarrow{BP} = (P - B) = ((x_p - x_b), (y_p - y_b)) = (X_{BP}, Y_{BP})$
$\overrightarrow{CP} = (P - C) = ((x_p - x_c), (y_p - y_c)) = (X_{CP}, Y_{CP})$

Radice quadrata inversa veloce

Oltre che per la sua indiscussa utilità, il calcolo del reciproco della radice quadrata è diventata famosa anche per aver innescato per anni accesi dibattiti in rete sulla paternità, sulle origini e sul funzionamento di una sua particolare implementazione (si veda [1] per una breve ricostruzione storica). L'implementazione del metodo in questione è utilizzata prevalentemente in ambiti (come quello grafico) dove non è richiesta una accuratezza elevata ma solo una buona approssimazione. La particolarità che l'ha resa nota, però, è stata quella di essere tanto efficiente quanto criptica.

1  float FastInvSqrt(float x)
2  {
3     float xhalf = 0.5f*x;          // xhalf = x/2
4     int i = *(int*)&x;             // converte x nella sua rappresentazione intera
5     i = 0x5f3759df - (i >> 1);     // prima approssimazione di 1/sqrt(x) (nella sua rappresentazione intera)
6     x = *(float*)&i;               // ritorna alla rappresentazione floating-point
7     x = x * (1.5f - xhalf * x*x);  // singolo passo di Newton-Raphson per ottenere una approssim. di 1/sqrt(x) migliore
8     return x;                      // ritorna tale valore
9  }

Coordinate baricentriche definite da un triangolo in 2D e 3D

Dati tre punti non collineari $A, B$ e $C$ nel piano, le coordinate di un punto $P$ qualsiasi nello stesso piano possono essere calcolate attraverso la combinazione lineare dei due vettori ($B - A$) e ($C - A$)

\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$.

Verificare posizione di un punto rispetto a circumcerchio e circumsfera

L'idea alla base dell'algoritmo per il calcolo della triangolazione di Delaunay (vista in [1]) può essere sfruttata per verificare (in 2D) se un punto $D$ è interno, esterno o appartiene al cerchio circoscritto al triangolo individuato da tre punti $A$, $B$ e $C$.

Triangolazione di Delaunay in 2D

Problema: Dato un insieme $\mathbf{P}$ di $n$ punti nel piano, trovare la triangolazione di Delaunay $\mathbf{D}(\mathbf{P})$ di $\mathbf{P}$. Trovare cioè quella triangolazione tale che nessun punto di $\mathbf{P}$ sia interno al cerchio circoscritto di ogni triangolo di $\mathbf{D}(\mathbf{P})$.

Soluzione: Convex Hull in 3D
Complessità: $O(n^2)$

Convex Hull in 3D

Problema: Dato un insieme di $n$ punti nello spazio, trovare il più piccolo involucro convesso che contiene gli $n$ punti. Trovare cioè i punti, tra gli $n$ dell'insieme dato, che si trovano sui bordi del più piccolo poliedro convesso che contiene tutti gli altri.

Soluzione: Incrementale
Complessità: $O(n^2)$

Convex Hull in 2D

Problema: Dato un insieme di $n$ punti nel piano, trovare il più piccolo involucro (o inviluppo) convesso che contiene gli $n$ punti. Trovare cioè i punti, tra gli $n$ dell'insieme dato, che si trovano sui lati del più piccolo poligono convesso che contiene tutti gli altri.

Soluzione: Graham Scan
Complessità: $O(n\log n)$

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.

Collinearità di tre punti in 2D

Verificare se tre punti $A, B, C$ nel piano giacciono sulla stessa retta risulta abbastanza semplice se si riduce il problema al calcolo dell'area del parallelogramma individuato dai due vettori $\quad \mathbf{n_1} = B - A \,\,$ e $\,\, \mathbf{n_2} = C - A$.
Se l'area è nulla, infatti, i tre punti sono collineari.

Triangolare un poligono con ear clipping

Problema: Dato un poligono $P$ descritto da $n$ punti nel piano trovare una triangolazione di $P$ in modo che le diagonali del poligono siano segmenti non intersecanti ed interni al poligono e che ogni regione interna sia un triangolo.

Soluzione: Ear clipping
Complessità: $O(n^2)$

Per triangolare un poligono la soluzione più semplice, ma anche più costosa, è quella di prendere tutti i segmenti che uniscono due vertici del poligono (ci sono $O(n^2)$ possibili candidati) e, per ogni segmento, controllare che sia interno al poligono e che non ci siano intersezioni con i lati del poligono (costo $O(n)$). Queste operazioni hanno quindi costo $O(n^3)$ e devono essere ripetute sino a quando non si trovano tutte le diagonali (che sono $n-3$). Il costo totale dunque arriva a $O(n^4)$. L'algoritmo è molto semplice ma con un piccolo sforzo in più si può ridurre il costo computazionale a $O(n^2)$ (in realtà si potrebbe arrivare fino a $O(n)$ ma il livello di complessità inizierebbe a farsi piuttosto elevato; per metodi più "semplici" con costo $O(n\log n)$ si vedano invece [2] e [3]).

Verificare l'intersezione di due segmenti in 2D

Verificare se due segmenti nel piano si intersecano risulta relativamente semplice se si considera che, preso uno qualsiasi dei due segmenti (ad esempio il segmento $AB$ della figura sotto), un vertice dell'altro segmento ($CD$) si troverà a sinistra del primo segmento ($D$ in questo caso) e l'altro a destra ($C$).

Area di un poligono

Per calcolare l'area di un poligono la soluzione più logica e semplice sarebbe quella di triangolare il poligono e sommare le aree dei triangoli che compongono il poligono.

Questa soluzione risulta abbastanza intuitiva per poligoni convessi, meno per quelli concavi. Utilizzando però il giusto metodo si può sfruttare questa stessa soluzione a prescindere dal tipo di poligono.
Il classico calcolo dell'area di un triangolo $A=(b \times h)/2$, però, si rivela insufficiente in questo caso. Invece, dall'algebra lineare sappiamo che dati tre vertici $a, b, c$, sia il modulo del prodotto vettoriale ($(b - a)\times (c - a)$) che il determinante della matrice composta dalle coordinate (omogenee) dei tre vertici

Matrice inversa in Computer Graphics

La decomposizione $LU$ è, almeno in teoria, il metodo più efficiente attraverso il quale invertire una matrice. In pratica, però, quando si ha a che fare con matrici dalle dimensioni contenute (massimo $4\times4$) ci si può affidare alla classica formula analitica

\begin{equation}\label{eq:1}A^{-1} = \frac {1}{det(A)}\cdot C_{A}^T\end{equation}