3.4 L'algoritmo di Powell-Sabin

L'algoritmo che presentiamo in questo paragrafo, dati $(\mathbf{x}_i, f_i, \nabla \mathbf{f}_i) i=1,\ldots,N$ restituisce una superficie che è globalmente $C^1$.

Figura 3.2: Dati dell'algoritmo di Powell-Sabin

Sia assegnata una triangolazione $T$ dei dati $\mathbf{x}_i$ e consideriamo un generico triangolo $T_i$ (rappresentanto nella figura 3.2): per tale triangolo si calcola l'incentro, e poi si scelgono i punti $M_1, M_2, M_3$ che sono l'intersezione fra un lato del triangolo e il segmento che unisce gli incentri dei triangoli che condividono quel lato. In formule abbiamo:

$$\mathbf{o}=\left(\frac{l_1}{l_1+l_2+l_3}, \frac{l_2}{l_1+l_2+l_3}, \frac{l_3}{l_1+l_2+l_3},\right)$$

$$M_1 = (1 - \alpha_1) \mathbf{A_2} + \alpha_1 \mathbf{A_3}$$

$$M_2 = (1 - \alpha_2) \mathbf{A_3} + \alpha_2 \mathbf{A_1}$$

$$M_3 = (1 - \alpha_3) \mathbf{A_1} + \alpha_3 \mathbf{A_2}$$

con $0 \le \alpha_i \le 1, i=1,2,3$. Indichiamo poi le direzioni dei lati con

$$\mathbf{d}_1 = \mathbf{A}_2 - \mathbf{A}_1$$

$$\mathbf{d}_2 = \mathbf{A}_3 - \mathbf{A}_2$$

$$\mathbf{d}_3 = \mathbf{A}_1 - \mathbf{A}_3$$

Quindi possiamo esprimere le derivate direzionali in ciascun vertice con

$$\left\{ \begin{array}{ccc} p_1 & = & D_{\mathbf{d}_1} f(\mathbf{A... ...A}_1) = - \nabla \mathbf{f} (\mathbf{A}_1) \mathbf{d}_3 \\ \end{array} \right.$$

$$\left\{ \begin{array}{ccc} p_2 & = & D_{\mathbf{d}_2} f(\mathbf{A... ...{A}_2) = - \nabla \mathbf{f} (\mathbf{A}_2) \mathbf{d}_2 \\ \end{array}\right.$$

$$\left\{ \begin{array}{ccc} p_3 & = & D_{\mathbf{d}_3} f(\mathbf{A... ...{A}_3) = - \nabla \mathbf{f} (\mathbf{A}_3) \mathbf{d}_2 \\ \end{array}\right.$$

Figura 3.3: Macroelemento di Powell-Sabin

La tecnica che utilizza l'algoritmo di Powell-Sabin è quella di fare una sottotriangolazione di ogni triangolo e su ciascuno di essi definire un patch quadratico. In riferimento alla figura (3.3) è possibile individuare $6$ sottotriangoli per ciascun triangolo, congiungento l'incetro con gli $M_i$ ed i vertici. Per definire lo schema interpolatorio allora si fissano le condizioni su ciascun sottotriangolo. Se indichiamo con $b_i^{(j)}$ l'$i$-esimo punto di controllo del triangolo $j$-esimo, allora le condizioni per l'interpolazione dei vertici del triangolo sono

$$\mathbf{b}_1^{(1)} = \mathbf{b}_3^{(6)} = f(\mathbf{A}_1)$$

$$\mathbf{b}_1^{(3)} = \mathbf{b}_3^{(2)} = f(\mathbf{A}_2)$$

$$\mathbf{b}_1^{(5)} = \mathbf{b}_3^{(4)} = f(\mathbf{A}_3)$$

Analizziamo come fissare i punti per raccordare le derivate sui vertici del triangolo, deve risultare:

$$\mathbf{b}_2^{(1)} = \mathbf{b}_1^{(1)} - \frac{\alpha_3}{2} \mathbf{p}_1$$

$$\mathbf{b}_2^{(3)} = \mathbf{b}_1^{(3)} - \frac{\alpha_1}{2} \mathbf{p}_2$$

$$\mathbf{b}_2^{(5)} = \mathbf{b}_1^{(5)} - \frac{\alpha_2}{2} \mathbf{p}_3$$

ed

$$\mathbf{b}_2^{(6)} = \mathbf{b}_3^{(6)} - \frac{1-\alpha_2}{2} \mathbf{q}_1$$

$$\mathbf{b}_2^{(2)} = \mathbf{b}_3^{(2)} - \frac{1-\alpha_3}{2} \mathbf{q}_2$$

$$\mathbf{b}_2^{(4)} = \mathbf{b}_3^{(4)} - \frac{1-\alpha_1}{2} \mathbf{q}_3$$

Vanno poi fissate le condizioni sulle derivate direzionali dell'incentro ( $\mathbf{o}=(w_1, w_2, w_3)$):

$$\mathbf{b}_4^{(1)} = \mathbf{b}_5^{(6)} = \mathbf{b}_1^{(1)} - \frac{1}{2} \left( w_1 \mathbf{p}_1 + w_3 \mathbf{q}_1 \right)$$

$$\mathbf{b}_4^{(3)} = \mathbf{b}_5^{(2)} = \mathbf{b}_1^{(3)} - \frac{1}{2} \left( w_3 \mathbf{p}_2 + w_1 \mathbf{q}_2 \right)$$

$$\mathbf{b}_4^{(5)} = \mathbf{b}_5^{(4)} = \mathbf{b}_1^{(5)} - \frac{1}{2} \left( w_1 \mathbf{p}_1 + w_2 \mathbf{q}_3 \right)$$

Rimangono da fissare le condizioni per il raccordo dei sottotriangoli:

$$\mathbf{b}_1^{(2)} = \mathbf{b}_3^{(1)} = (1-\alpha_3)\mathbf{b}_2^{(1)} - \alpha_3 \mathbf{b}_2^{(2)}$$

$$\mathbf{b}_1^{(4)} = \mathbf{b}_3^{(3)} = (1-\alpha_1)\mathbf{b}_2^{(3)} - \alpha_1 \mathbf{b}_2^{(4)}$$

$$\mathbf{b}_1^{(6)} = \mathbf{b}_3^{(5)} = (1-\alpha_2)\mathbf{b}_2^{(5)} - \alpha_2 \mathbf{b}_2^{(6)}$$

e

$$\mathbf{b}_5^{(1)} = \mathbf{b}_4^{(2)} = (1-\alpha_3)\mathbf{b}_4^{(1)} - \alpha_3 \mathbf{b}_5^{(2)}$$

$$\mathbf{b}_5^{(3)} = \mathbf{b}_4^{(4)} = (1-\alpha_1)\mathbf{b}_4^{(3)} - \alpha_1 \mathbf{b}_5^{(4)}$$

$$\mathbf{b}_5^{(5)} = \mathbf{b}_4^{(6)} = (1-\alpha_2)\mathbf{b}_4^{(5)} - \alpha_2 \mathbf{b}_5^{(6)}$$

Infine la condizione per il raccordo sull'incentro è:

$$\mathbf{b}_6^{(j)} = w_1 \mathbf{b}_4^{(1)} + w_2 \mathbf{b}_4^{(3)} + w_3 \mathbf{b}_4^{(5)} \qquad\qquad j=1,\ldots,6$$