Appunti dal corso di Metodi Numerici per la Grafica

Carli Samuele
E-mail: carlisamuele@csspace.net
www.csspace.net

9 ottobre 2010

Indice

1 Nozioni di base
 Definizioni e richiami
  Notazione
  Rappresentazione in Scilab
   Bulk operations
 Operazioni tra vettori
  Prodotto scalare
  Prodotto vettoriale
 Combinazioni di punti
  Spazi di lavoro
 Trasformazioni affini
  Movimenti rigidi
   ProprietÓ
 Rappresentazioni parametriche
  Riparametrizzazioni
   Lunghezza dell’arco
 Frenet Frame
 Curvatura e Torsione
2 Curve parametriche
 Polinomi di Bernstein
  ProprietÓ
   Non negativitÓ
   Partizionamento dell’unitÓ
   Simmetria
   Indipendenza lineare
   Precisione lineare
   Variation diminishing
   Valori agli estremi
   Derivata
 Curve di BÚzier
  ProprietÓ
   Convex Hull
   Interpolazione agli estremi
   Invarianza alle trasformazioni affini
   Simmetria
   Precisione lineare
   Controllo pseudo-locale
   Variation diminishing
   Derivata
   Tangenza agli estremi
 Algoritmi per curve di BÚzier
  Disegno della curva: de Casteljau
  Degree elevation
  Suddivisione
 Curve di BÚzier razionali
  Coordinate omogenee
   Interpretazione geometrica in d = 3
   Curve di BÚzier razionali
  ProprietÓ
   Simmetria
   Convex Hull
   Invarianza per trasformazioni affini
   Linear precision
   Valori agli estremi
   Variation diminishing
   Controllo pseudolocale
   Derivata
   Tangenza agli estremi
  Curve di BÚzier razionali in forma standard
 Sezioni coniche
  Coordinate baricentriche bidimensionali
  Rappresentazione di curve coniche
  Archi di coniche con curve di BÚzier
3 Curve SPLINE
 BÚzier-spline
  ContinuitÓ dei raccordi
  ContinuitÓ geometrica
   Implementazioni
 B-SPLINE
  Curve SPLINE monovariate
  Base delle potenze troncate
  Base delle B-SPLINE
  ProprietÓ
   Supporto compatto
   Non negativitÓ
   Partizione dell’unitÓ
   Derivata
   Indipendenza lineare
  Generalizzazione della base delle B-SPLINE
   Derivata
  Curve B-SPLINE
   ProprietÓ
   Invarianza alle trasformazioni affini
   Controllo locale
   Strong ConvexHull
   Linear Precision
   Derivata
   End point interpolation
   Derivata agli estrimi: raccordo
  Algoritmi per curve BSPLINE
   Disegno della curva: de Boor
   Knot Insertion
 Curve NURBS
  ProprietÓ
4 Interpolazione
 Parametrizzazioni
 Interpolazione mediante SPLINE cubica C1
   Schemi di approssimazione delle derivate
   Schema FMILL
   Tangenti di Bessel
 Interpolazione mediante SPLINE cubica C2
  Scelta delle tangenti wi
 Spline cubiche C1 G2
5 Superfici parametriche
 Rappresentazione parametrica
  Esempi di superfici parametriche
   Sfera
   Cono
 Riparametrizzazione ammissibile
 Tipologie di superfici
  Superfici di rivoluzione
  Superfici Rigate
  Superficie bilineare
 Superfici tensor-product
  Superfici tensor-product come patch BÚzier
  ProprietÓ
   Curve di contorno
   Posizione negli angoli
   Altre proprietÓ
   Forma compatta
  Algoritmi
   de Casteljau
   Degree-elevation
   Subdivision
  Interpolazione di superfici con patch di BÚzier
   Interpolazione su punti sparsi
   Interpolazione su griglia rettangolare
 Patch triangolari di BÚzeir
   Coordinate baricentriche
  Polinomi di Bernstein generalizzati
   ProprietÓ
   Non negativitÓ
   Partizione dell’unitÓ
   Linear Precision
  Patch triangolari
   ProprietÓ
   Comportamento ai bordi
   Controllo pseudolocale
  End-point interpolation
   Forma ricorsiva
  Algoritmi
   de Casteljau
   Degree-elevation
A Codici
 A.1 BÚzier SPLINE
 A.2 Curve di BÚzier
 A.3 Curve BSPLINE
 A.4 Interpolazione
 A.5 Superfici parametriche
 A.6 Mouse input
 A.7 Lista di import per tutti gli script

Capitolo 1
Nozioni di base

Definizioni e richiami

Notazione

dove non diversamente specificato.

Rappresentazione in Scilab

I vettori sono convenzionalmente intesi come vettori colonna .

Bulk operations

Gli oggetti grafici sono usualmente rappresentati come liste di punti; in generale una curva sarÓ rappresentata con una matrice d×n In questo contesto d=2 o 3 .
Per migliorare l’efficienza le operazioni tra vettori sono implementate in modo da eseguire la stessa operazione in blocco su liste di vettori di lunghezza arbitraria (ovvero punto a punto lungo una curva).

Operazioni tra vettori

Prodotto scalare

Definizione 1.1 (Prodotto scalare).

∙,⋅†   †Per brevit`a talvolta si omette completamente l’operatore: vw                               : ℝd × ℝd → ℝ

ProprietÓ

Su una base canonica: vw = viwi = vT w
Dati v(t),w(t):

-d            dv(t)       dw-(t)
dt v(t)w (t) =  dt  w (t) +   dt  v(t)

Prodotto vettoriale

Definizione 1.2 (Prodotto vettoriale).

        d     d     d
×, ∧ : ℝ  × ℝ  →  ℝ

con d = 2, 3

ProprietÓ

Dati v(t),w(t):

d-               dv-(t)          dw-(t)
dt(v(t) ∧ w (t)) =  dt  ∧ w (t) +   dt  ∧ v(t)

Casi particolari:

(    )   (    )
   a        c
   b   ∧    d   = (ad − cb)

(    )    (     )   (               )
| v1 |    | w1  |   |  v2w3 − v3w2  |
( v2 )  ∧ ( w2  ) = (  v3w1 − v1w3  )
  v3        w3         v1w2 − v2w1


: Cross product 2D/3D su array di vettori
 
1function cp = crossprod(x,y) 
2  [rx,cx] = size(x); 
3  [ry,cy] = size(y); 
4  if or([rx˜=ry,cx˜=cy]) then printf(Cannot croosprod different vectors!); cp=0; end; 
5  if rx==2 then 
6    cp = diag( x * [y(2,:);-y(1,:)] ) 
7  end; 
8  if rx==3 then 
9    cp = [] 
10    for n=1:cx 
11      cp = [cp,[x(2,n)*y(3,n)-x(3,n)*y(2,n); x(3,n)*y(1,n)-x(1,n)*y(3,n); x(1,n)*y(2,n)-x(2,n)*y(1,n)]] 
12    end; 
13   end; 
14endfunction

Combinazioni di punti

Una combinazione di punti:

P = ∑  αiPi,   αi ∈ ℝ,Pi ∈ Ed
     i

si dice baricentrica se iαi = 1, convessa se anche αi > 0i
Se la combinazione Ŕ baricentrica, pu˛ essere scritta come la somma di uno qualsiasi degli addendi e un vettore:

∑              ∑          ∑              ∑
   αiPi = Pk +     αiPi −    αiPk = Pk +    αi (Pi − Pk )
 i               i         i              i

Una combinazione baricentrica convessa Ŕ utile per rappresentare parametrizzazioni di oggetti:

Esempio 1.1 (Punto medio di un segmento).  
s = b0b1; b0,b1 2;
M = i=1,2αibi = b0 + α0(b0 b0) + α1(b1 b0) = b0 + α1(b1 b0); α1 = α2 = 12

Esempio 1.2 (Punti di un segmento). Generalizzando l’esempio precedente:
s = b0b1; b0,b1 2; t [0, 1];
P(t) = tb0 + (1 t)b1 = b0 + (1 t)(b1 b0)

Spazi di lavoro

Definizione 1.3 (Spazio di lavoro). Si definisce spazio affine di ordine d, a partire da un punto O d, l’insieme:

 d   {                 d}
E  =  P  = O  + v|v ∈ ℝ

Definizione 1.4 (Convex Hull).

           {             ∑      }
CHb0...bn =   P ∈ Ed |P  =    αibi
                          i

in cui gli αi formano una combinazione baricentrica convessa.
CH Ŕ il pi¨ piccolo insieme convesso che contiene i punti b0...bn

Trasformazioni affini

Definizione 1.5 (Trasformazione affine).

        d     d
Φ (x) : E →  E

lineare rispetto ad x.

Generalmente si usano funzioni del tipo

Φ (x) = Ax  + v
(1.1)

in cui A Md×d e v d che sono lineari e permettono di definire facilmente le trasformazioni pi¨ comuni in grafica.

Movimenti rigidi

Traslazione La traslazione si ottiene con una funzione del tipo (1.1) in cui A I, per cui Φ(x) = x + v. Alternativamente Ŕ possibile modificare la matrice A:

(        ) (  )    (      )
|1   0  a| | x|    | x + a|
(0   1  b) ( y) =  ( y + b)
 0   0  1    1         1

Rotazione La rotazione intorno all’origine (del sistema di riferimento ortonormale corrente) si ottiene ancora con una funzione del tipo (1.1) in cui v 0 (vettore nullo, non c’Ŕ traslazione) ed in cui la matrice, nel caso tridimensionale, assume la forma :

(                  )
   cosθ  − sinθ  0
|(  sin θ   cosθ   0 |)
    0       0    1

Scalatura La scalatura si ottiene imponendo un fattore di scala ad ogni dimensione della matrice identitÓ:

(          )
   α  0  0
|(  0  β  0 |)
   0  0  γ

ProprietÓ delle trasformazioni affini della famiglia (1.1)

Le trasformazioni affini sono lineari rispetto alle combinazioni baricentriche. Infatti:

           (       )      (        )     ◜= ◞1◟-◝
Φ (P ) = Φ  ∑  α P    = A  ∑   α P   + v ∑  α =  ∑  α (AP  + v ) = ∑  α Φ(P )
             i   i i         i  i i       i  i    i  i    i         i  i    i

Inoltre, se A Ŕ ortogonale ATA = AAT = I , una trasformazione della famiglia (1.1) lascia invariati angoli e lunghezze.
Infatti, dati x,y vettori, per le lunghezze:

∥Φ (x)− Φ (y )∥2 =  ∥Ax + v− Ay − v ∥2 = ∥A(x − y)∥2 =    (A(x − y))T(A(x − y)) = ∥x− y ∥2
                                                    def

Per gli angoli, se x,y sono tali che xT y = |x||y| cos θ ed, in quanto vettori, possono essere letti come x = a b; y = c d (con a,b,c,d d):

Φ(x )T Φ(y) = (Φ (a) − Φ (b))T(Φ(c) − Φ(d)) = (Ax )TAy = xT y

Rappresentazioni parametriche

Definizione 1.6 (Curva). Si definisce curva l’immagine di un’applicazione X : I Ed continua e localmente iniettiva.

Definizione 1.7 (Rappresentazione Parametrica). L’applicazione X : I Ed continua e localmente iniettiva si dice rappresentazione a parametri in I della curva; in generale una stessa curva pu˛ ammettere pi¨ di una rappresentazione parametrica. Ad esempio, un segmento di retta (curva degenere) tra l’origine e il punto (1,1) pu˛ essere rappresentata in forma implicita (ax + by + c = 0), esplicita (y = mx + q) o parametrica su un’intervallo I(= [0,1] per esempio):
X(t) = {
   x = t
   y = t o equivalentemente X(t) = {         k
  x = 1− tk
  y = 1− t

Definizione 1.8 (Rappresentazione parametrica differenziabile). La rappresentazione X(t) di dice differenziabile se esiste la funzione:

        dX (t)
X˙(t) = ------
          dt

Se 0 su tutto l’intervallo di definizione I, X si dice differenziabile regolare

Riparametrizzazioni

La non unicitÓ della rappresetazione parametrica suggerisce che si possa passare da una rappresentazione all’altra di una stessa curva, utilizzando di volta in volta la parametrizzazione con le caratteristiche pi¨ adatte al contesto.

Definizione 1.9 (Funzioni di riparametrizzazione). Sia τ : I I2. Si definisce riparametrizzazione di X rispetto a τ una funzione  ˜
X tale che  ˜
X(τ(t)) = X(t).

Si ha che:

dX˜-  dX---1-   dX--dt
dτ  =  dt dτ =  dt d τ
          dt

in cui dX∕dt Ŕ la tangente alla curva e dt∕dτ Ŕ uno scalare > 0.
Se τ C1 e ˙τ > 0 su tutto I la riparametrizzazione si dice ammissible; in questo caso la curva Ŕ percorsa tutta nello stesso verso e la sua tangente non cambia verso nÚ direzione, ma solo il modulo.

Osservazione 1.1 (Tangente alla curva). La tangente alla curva non pu˛ essere usata per caratterizzarla, in quanto potenzialmente cambia modulo al cambiare della rappresentazione.
Per caratterizzare una curva occorre definire il versore tangente:

        X˙(t)-    ˙
e1(t) =  |X ˙(t)|, X (t) ⁄= 0

Definizione 1.10 (Rappresentazione parametrica regolare). Se (t)0 t I la rappresentazione parametrica si dice regolare.

Lunghezza dell’arco

Definizione 1.11 (Lunghezza dell’arco). Si definisce lunghezza dell’arco la funzione:

        ∫  |      |
         t ||dX--(z-)||
S (t) =  a || dz   ||dz

La lunghezza dell’arco pu˛ essere usata come funzione di riparametrizzazione ammissibile, che ha anche la caratteristica di essere lineare rispetto alle lunghezze sulla curva. In particolare S : [a,b] [0,l], e dS(dtt) = ||    ||
|ddX((t)t)| implica che se X(t) Ŕ regolare anche S(t) lo Ŕ, e vale anche:

dX--  dX--dt-   dX---1---
dS  =  dt dS =   dt ||dX-|| = 1
                    |dt|

Osservazione 1.2 (Invarianza rispetto alla rappresentazione). La lunghezza dell’arco Ŕ invariante rispetto alla rappresentazione parametrica della curva, infatti, cambiando parametrizzazione:

        ∫  |   |      ∫  |      |
          z||dX-||       k ||dX--dt||dy-
S (z) =  a ||dy ||dy =   a ||dt dy ||dt dt = S(k)

Frenet Frame

Definizione 1.12 (Fernet Frame). Il Frenet Frame Ŕ un sistema di coordinate ortonormali definito localmente in un punto di una curva X(t).

╚ possibile definire un Frenet frame L’insieme di tutti i Frenet frame Ŕ rappresentato tramite tre vettori ortonormali e1,e2,e3; in realtÓ ogniuno di essi Ŕ una funzione del parametro della rappresentazione della curva (C’Ŕ un Frenet frame diverso per ogni punto della curva). a partire dalle derivate di una qualsiasi rappresentazione parametrica X(t) di una curva:

(
||||  e (t) = -X˙(t)                                                  versore tangente
||||   1     |X˙(t)|
||||
{          de1(t)     ˙       Ę      ˙         ˙    ˙     Ę
||  e2(t) = ---dt-- =  X-(t) ×-(X-(t) ×-X-(t))-=-X-(t)(X-(t)-⋅X-(t))    versore normale
||||         |de1(t)|    |X ˙(t)||X Ę(t) × X˙(t)|  |X˙(t)(X ˙(t) ⋅ ĘX (t))|
||||            dt
||(                                                                                ‡
   e3(t) = e1(t) × e2(t)                                          versore binormale

IdentitÓ “BAC-CAB”: A × (B × C) = B(A C) C(A B).

Curvatura e Torsione

Definizione 1.13 (Curvatura). Si definisce curvatura la quantitÓ:

        | ˙X (t) × XĘ(t)|
K (t) = --------------
           |X˙(t)|3

La curvatura Ŕ una grandezza caratteristica della curva indipendente dalla sua rappresentazione parametrica.
L’espressione della curvatura si semplifica notevolmente quando la rappresentazione Ŕ lunghezza dell’arco:

K (s) = |X Ę(s)|

in quanto (s) 1.


: Calcolo della curvatura
 
1function K = curvature(dx,ddx) 
2  K=(sum(crossprod(dx,ddx).^2,r)^0.5)./(sum(dx.^2,r)^0.5)^3 
3endfunction

Definizione 1.14 (Torsione). Si definisce torsione la quantitÓ:

           [           ...   ]
       det( X˙(t);XĘ(t);X (t))
τ(t) = ----||------------||2-----
           |X˙(t) × XĘ(t)|

La torsione si annulla solo quando il prodotto misto Ŕ nullo:

                ...         [           ...   ]
(X˙(t) × XĘ(t)) ⋅X (t) = det( X˙(t);XĘ(t);X (t)) = 0

La torsione Ŕ non nulla solo nello spazio e non su un piano.


: Calcolo della torsione
 
1function T = torsione(dx,ddx,dddx) 
2  [d,n] = size(dx) 
3  den = (sum(crossprod(dx,ddx).^2,r)); 
4  T = []; 
5  for i=1:n 
6    T = [T,det([dx(:,i),ddx(:,i),dddx(:,i)])/den(i)]; 
7    end; 
8endfunction

Capitolo 2
Curve parametriche

Polinomi di Bernstein

Definizione 2.1 (Base di Bernstein). I polinomi di Bersntein sono una famiglia di polinomi definiti, per ogni n +, come:

( B0 (t)† = 1
||||   0
||{          ( )
  Bni (t) = ni ti(1 − t)n−i ∈ Πn, i ≤ n
||||
||(   n                    ‡
  B i (t) = 0,i > n,i < 0


PIC (a) Terzo grado
PIC (b) Quarto grado

Figura 2.1: Basi di Bernstein Bin(t)


Definizione 2.2 (Polinomi di Bernstein in forma ricorsiva).

B00 = 1
Bjn = 0,  j0, 1,...,n
Bin = (1 t)B in1 + tB i1n1
Le due definizioni sono equivalenti, infatti:

Caso (i = 0 e i = n).

  n            n     n    n
B 0(t) = (1 − t) , B n = t

Caso (0 < i < n). in quanto (n)
 i = (n−1)
  i + (n−1)
 i−1:

Bin(t) = (  )
  n
  iti(1 t)ni = ((      )   (      ) )
   n − 1  +   n − 1
     i        i − 1ti(1 t)ni
= (1 t)(      )
 n − 1
    iti(1 t)ni1 + t(     )
 n − 1
 i − 1ti1(1 t)ni
= (1 t)Bin1 + tB i1n1

Definizione 2.3 (Polinomi di Berstein di grado n). Si definisce Polinomio di Bernstein di grado n la serie:

       ∑n     n
p(t) =    ckB k(t)
       k=0


PIC

Figura 2.2: Basi di Bernstein di ventesimo grado


ProprietÓ

I polinomi di Bernstein sono definiti su tutto , ma ci si limita a studiarne le proprietÓ nell’intervallo [0, 1] (sufficiente in questo contesto).

Osservazione 2.1 (MappabilitÓ della base). I polinomi di Bernstein possono essere definiti, in maniera equivalente, su qualsiasi intervallo [a,b] mappabile a [0, 1]:

         (  )
  n       n  (t-−-1)i(b-−-t)n−i
B i (t) =  i      (b − a)n    ,  t ∈ [a, b]

Non negativitÓ

Ogni polinomio della base Ŕ non negativo in quanto prodotto di componenti non negative ((n)
 i,ti, (1 t)ni).

Partizionamento dell’unitÓ
∑n   n
    Bi (t) = 1
i=0

infatti, t :

                       ∑n (  )              ∑n
1 =  (t + (1 − t))n =def      n  ti(1 − t)n− i ≡    Bn(t)
                       i=0  i                i=0  i


PIC

Figura 2.3: Partizionamento dell’unitÓ, la somma dei polinomi Ŕ evidenziata in rosso


Simmetria
Bn (t) = Bn  (1 − t)
  i       n−i

infatti:

        (  )              (      )
 n        n  i      n− i      n          n−ii      n
Bi (t) =     t (1 − t)    =         (1 − t)   t = †B n− i(1 − t)
          i                 n − i

Per il binomio di Newton vale la proprietÓ ( )
 ni = (   )
 nn−i

Indipendenza lineare

I polinomi di Bernstein sono una base per Πn; infatti, Bin pu˛ essere scritto come:

        (  )                (  )    (     )
  n       n   i      n−i   i  n  n∑−i n −  i   j    n−i−j n−i− j
B i (t) =    t(1 − t)   =  t                [1 ] − 1     t
          i                   i  j=0    j

da cui:

(   )    ( 1  ...  ...  ...  1)  (  )
 Bn0     |                   .|   t0
||   ||    || 0  ...            ..||  ||  ||
||  ..||    || ..       (n)       ..||  || ..||
||  .||  = || .        i        .||  || .||
|(   |)    |( ...            ...  ...|)  |(  |)
 Bnn                              tn
           0                 1

Se la matrice del cambio di base Ŕ non singolare i polinomi di Bernstein sono linearmente indipendenti; la matrice Ŕ non singolare in quanto triangolare superiore con tutti elementi non nulli sulla diagonale (determinante non nullo).

Precisione lineare
∑n i- n
   n Bi (t) = t
i=0

Infatti, per induzione su n:

  1
B 1(t) = t

e supponendo che valga per n 1:

i=0n i
nBin(t) = i=0ni
--
n[        n−1     n− 1]
 (1 − t)B i  +  tBi−1
= n-−-1-
  n(1 t) i=0n1---i--
n − 1Bin1(t) + tn-−-1-
  n i=1n--i---
n − 1Bi1n1(t)
= n-−-1-
  n(1 t)t + n-−--1
  nt j=0n1j +-1-
n − 1Bjn1(t) =
= n-−-1-
  n(1 t)t + n-−--1
  nt⌊n− 1                      n−1              ⌋
⌈ ∑  --j---Bn− 1(t) + --1---∑   j +-1Bn −1(t)⌉
  j=0 n − 1  j        n − 1 j=0 n − 1  j
= n-−-1-
  n(1 t)t + n-−--1
  nt[         ]
 t + --1---
     n − 1 =
= n − 1
--n---(1 t)t + n −  1
--n---t2 + t
n- = t

PIC

Figura 2.4: Precisione lineare: in verde i=0ni
--
nBin(t).


Variation diminishing

Un polinomio di Bernstein p(t) = k=0nc kBkn(t) a coefficienti c 0,..,cn, dato V il numero di variazioni di segno della sequenza di costanti e N numero di radici reali di p(t) contate con la loro molteplicitÓ, gode della proprietÓ:

N ≤  V,   e   N  = V −  2k,  t ∈ [0, 1].

Infatti,

        n   (  )
p(t) = ∑  c  n  tk(1 − t)n− k
       k=0  k k

parametrizzato su u [0,], ovvero considerando t = -u--
1+u e (1 t) = -1--
1+u, diventa:

           n   (  )                   n    ( )     k                n   (  )
          ∑      n   -u--k -1--n−k   ∑      n  ---u----    ---1----∑      n   k
p (t(u)) =    ck  k  (1+u) (1+u)    =     ck k  (1 + u)n =  (1 + u )n   ck k  u
          k=0                        k=0                           k=0

che Ŕ un polinomio sulla base canonica per il quale vale la regola di Cartesio Il numero di radici positive (contato con molteplicitÓ) Ŕ dato dal numero di cambi di segno (variazioni) fra due coefficienti consecutivi. .

Valori agli estremi

Bin(0) = ({
  1,  i = 0
( 0,  i ⁄= 0
Bin(1) = (
{ 1,  i = n
(
  0,  i ⁄= n

Derivata

La derivata dei polinomi di Bernstein pu˛ essere espressa in termini della base di Bernstein stessa:

dBni-
 dt = -d
dt(  )
  n
  iti(1 t)ni =
= ---in!---
i!(n − i)!ti1(1 t)n1 (n-−-i)n!
i!(n − i)!ti(1 t)n1i
= n   (n − 1)!
---------------
(i − 1)!(n − i)!ti1(1 t)ni n   (n − i)!
-------------
i!(n − i − 1)!ti(1 t)ni1
= n((      )                (      )              )
   n − 1  i−1       n− 1    n − 1   i      n−i− 1
          t  (1 − t)   −          t(1 − t)
   i − 1                     i
= n(                  )
 Bni−−11(t) − Bni−1(t)
La stessa derivata per˛ pu˛ anche essere rappresentata in modo da evidenziare il comportamento della funzione:
dBn
--i--
 dt =  d
--
dt(n )

  iti(1 t)ni
= (  )
  n

  iti1(1 t)ni1[i(1 t) t(n i)]
= (  )
  n
  iti1(1 t)ni1[i nt]
da cui si nota che i punti non banali in cui la derivata si annulla sono facilmente determinabili:  ^
ti = i-
n, con i = 0,,n

Curve di BÚzier


PIC

Figura 2.5: Esempio di curve di BÚzier; in rosso il poligono di controllo


Definizione 2.4 (Curva di BÚzier). Dati b0,...,bn Ed insieme di punti di controllo, si definisce curva di BÚzier il polinomio:

        ∑n    n
X (t) =    biBi (t), t ∈ [0, 1]
        i=0

Definizione 2.5 (Poligono di controllo). Si definisce poligono di controllo il poligono che si ottiene unendo ordinatamente i punti di controllo.

ProprietÓ

Convex Hull

La curva giace completamente nella convex hull dei punti di controllo, in quanto ogni punto della curva Ŕ una combinazione baricentrica di questi I polinomi di Bernstein partizionano l’unitÓ .

Interpolazione agli estremi

Segue direttamente dai valori agli estremi dei polinomi di Bernstein che X(0) b0 e X(1) bn.

Invarianza alle trasformazioni affini

Data Φ : Ed Ed di tipo Φ(p) = Ap + v, si ha:

                        ∑n             ∑n           ∑n                  ∑n
Φ (X (t)) = AX  (t) + v =   AbiBn  (t) +    vBn (t) =   (Abi + v)Bn (t) =    Φ (bi)Bn (t)
                        i=0     i      i=0   i      i=0           i      i=0        i

Simmetria

Le curve di BÚzier sono simmetriche rispetto all’ordine dei punti di controllo, infatti, dati ci = bni:

        ∑n           ∑n                 ∑n
X˜(t) =    ciBni (t) =   ciBnn− i(1 − t) =    biBni (1 − t) = X (1 − t)
        i=0           i=0                i=0

Precisione lineare

Se i punti di controllo sono allineati ed equidistanti la convex hull Ŕ il segmento b0bn e la curva giace su una retta:

        ∑n                           ∑n                   n∑
X (t) =    (b0 + i(bn − b0))Bni (t) = b0    Bni (t) + (bn − b0)  iBni (t) = †b0 + t(bn − b0)
        i=0      n                     i=0                  i=0 n

Precisione lineare dei polinomi di Bernstein


PIC
PIC

Figura 2.6: Precisione lineare: poligono di controllo e curva che giace su di esso; si nota che i punti della curva sono distribuiti uniformemente rispetto al parametro t.


Controllo pseudo-locale

Modificando uno dei punti di controllo b0,...,bn:

     (
     {bi,  i ⁄= j
ci = (
      bi + v  i = j

si ottiene:

˜      ∑n     n     ∑n     n        n
X (t) =    ciB i (t) =    biB i (t) + vB j (t)
       i=0          i=0

ovvero tutti i punti della curva si spostano nella direzione di v con un’incidenza che dipende da t e dalla forma del jesimo polinomio di Bernstein. In generale, in quanto i polinomi hanno una forma a campana, i punti della curva sono interessati in maniera proporzionale alla distanza dal punto di controllo modificato.


PIC
PIC

Figura 2.7: Controllo pseudolocale: la curva (verde) pu˛ essere fatta coincidere con uno qualsiasi dei polinomi della base


Variation diminishing

La proprietÓ di variation diminishing aiuta a capire l’andamento della curva; si pu˛ infatti sapere il numero massimo di intersezione che la curva pu˛ avere con una retta (N) a partire da quello del poligono di controllo (Np). In particolare, N Np e N = Np 2k.
Sia r = αx + βy + γ, le intersezioni si possono ricavare, separando le dimensioni, con l’equazione:

   n              n               n
  ∑      n       ∑       n       ∑    n
α    bixB i (t) + β   biyB i (t) + γ  B i (t) = 0
  i=0             i=0             i=0

da cui:

∑n                    n
   (◟αbix-+◝β◜biy +-γ-)◞B i (t) = 0
i=0        ci

e la proprietÓ segue direttamente da quella di variation diminishing dei polinomi di Bernstein.

Derivata (o curva Hodograph)

(t) = i=0nb iBin(t) = n i=0nb i(  n−1       n− 1  )
 B i−1 (t) − Bi  (t)
= n[          n−1                      n−1]
 (b1 − b0)B 0  (t) + ...+  (bn − bn− 1)B n−1
= n i=0n1(b i+1 bi)Bin1(t) = n i=0n1Δb iBin1(t)

dove Δbi = def(bi+1 bi).
La derivata di una curva di BÚzier Ŕ quindi ancora una curva dello stesso tipo ma di grado minore e con punti di controllo differenti.


: Calcolo della hodograph
 
1function bh=Bezier˙hodograph(b) 
2[d,n] = size(b) 
3for i=1:d 
4       bh(i,:) = diff(b(i,:)); 
5end 
6bh=bh*n; 
7endfunction

Derivate successive

Dalla forma della derivata seconda:

                n∑−2
XĘ(t) = n(n − 1)    Δ2biBni−2(t)
                i=0

si pu˛ generalizzare Falling factorial: xk = x(x 1)...(x k + 1), ovvero i primi k termini del fattoriale. Nota: xk = --x!-
(x−k)! :

  (r)       rn∑− r  r   n−r
X   (t) = n-    Δ biB i  (t)
             i=0

Tangenza agli estremi

In quanto:

X˙(0) = nΔb0,   X˙(1) = nΔbn − 1

la curva di BÚzier Ŕ sempre tangente ai segmenti che uniscono i primi e gli ultimi due punti del poligono di controllo.


PIC

Figura 2.8: Tangenza agli estremi



PIC

Figura 2.9: Esempio di hodograph prima e seconda e curvatura



PIC

Figura 2.10: Altro esempio di hodograph prima e seconda e curvatura


Algoritmi per curve di BÚzier

Disegno della curva: de Casteljau

I punti della curva possono essere calcolati ricorsivamente definendo:

(
{ b0i = bi,                     i = 0, ...,n
( br=  (1 − t)br−1 + tbr− 1,    i = 0, ...,n − r;r = 1,...,n
   i          i       i+1

Si dimostra per induzione che questa formula ricorsiva, in r passi, calcola i punti della curva. Per r = 1:

                        1
b1 = (1 − t)b  + tb   = ∑  b   B1 (t)
 i          i     i+1   j=0 i+j  j

Supponendo che valga per r 1, si ha:

bir(t) = (1 t)b ir1(t) + tb i+1r1(t)
= (1 t) j=0r1b i+jBjr1(t) + t j=0r1b i+j+1Bjr1(t)
= (1 t) k=ii+r1b kBkir1(t) + t k=i+1i+rb kBki1r1(t)
= k=ii+r[                               ]
 (1 − t)bkBrk−− 1i(t) + tbkBrk−−1i−1(t)
= k=ii+rb kBkir(t) = j=0rb j+iBjr(t)
ovvero variando t si possono ottenere tutti i punti della curva.

: Algoritmo di de Casteljau
 
1function Curve = Bezier˙casteljau(b,t) //b poligono, t punti di valutazione 
2lt = length(t) 
3[d,n] = size(b) 
4f = t 
5e = 1-t 
6b = b 
7 
8for i=1:d 
9     C(:,:,i) = b(1:n-1,i)*e+b(2:n,i)*f 
10end 
11for r=2:n-1 
12     for row=1:n-r 
13          for i=1:d 
14               C(row,:,i) = C(row,:,i).*e+C(row+1,:,i).*f 
15          end 
16     end 
17end 
18for i=1:d 
19     Curve(i,:)=C(1,:,i) 
20end 
21endfunction

Degree elevation

Una curva di BÚzeier Ŕ rappresentabile con un’altra di grado pi¨ elevato, ovvero esiste soluzione all’equazione:

        ∑n    n      n∑+1    n+1
X (t) =    biB i (t) =   ciB i  (t)
        i=0           i=0

Per la proprietÓ di interpolazione agli estremi si deve avere c0 = b0 e cn+1 = bn.
Moltiplicando a sinistra per (t + (1 t)) e usando la proprietÓ di partizionamento dell’unitÓ dei polinomi di Bernstein, si ha:

             n∑            ∑n  (  ) [                             ]   n∑+1
(t + (1 − t))    biBn (t) =    bi n   ti(1 − t)n−i+1 + ti+1(1 − t)n−i =     ciBn+1 (t)
             i=0    i      i=0    i                                    i=0    i

da cui, considerando b1 = bn+1 = 0:

i=0nb i(n )

  i[                            ]
 ti+1(1 − t)n−i + ti(1 − t)n− i+1
= i=0nb i(  )
 n

  iti+1(1 t)ni + i=0nb i(  )
  n

  iti(1 t)ni+1
= i=0n+1b i1(     )
   n
 i − 1ti(1 t)ni+1 + i=0n+1b i(  )
 n
  iti(1 t)ni+1
= i=0n+1[(   )       ( )  ]
  ni−1 bi−1 +  ni bi
------(n+1)-------
        iBin+1(t)
Si sono trovati i punti di controllo ci:
     [( n )       (n) ]
       i−1 bi−1 +  i bi      i         (      i   )
ci = -------(n+1)------ =  n-+-1bi−1 +  1 − n-+-1- bi
              i

Osservazione 2.2. I punti di controllo ci sono ottenuti come combinazione lineare dei due vertici adiacenti del poligono di controllo di partenza, quindi il nuovo poligono Ŕ inscritto nel vecchio.


: Degree elevation
 
1function B = Bezier˙degelev(b) 
2[d,n] = size(b) 
3B(:,1) = b(:,1) 
4B(:,n+1) = b(:,n) 
5for i=2:n 
6     k = (i-1)/(n); 
7     B(:,i) = k*b(:,i-1) + (1-k)*b(:,i); 
8end 
9endfunction


PIC
PIC

Figura 2.11: Esempi di elevazione del grado


Osservazione 2.3 (Repeated degree elevation). Per disegnare una curva si pu˛ pensare di innalzare il grado del poligono di controllo finnchŔ questo non sia sufficientemente vicino alla curva; i coefficienti dell’r-esima iterazione sono:

                (   )
      ∑n   (n )   r
c(ir)=     bi    (i−j)-
      j=0    j   n+ir

Fissato t, se r →∞ il poligono tende effettivamente alla curva, ma il costo computazionale Ŕ molto pi¨ elevato dell’algoritmo di de Casteljau in quanto servono molte elevazioni per raggiungere una precisione paragonabile.

Suddivisione

Una curva pu˛ essere suddivisa in due curve dello stesso grado giunte in un punto.
Dati ^t]0, 1[ e:

        n
       ∑      n
X (t) =    biB i (t)
        i=0

si cercano:

(         ∑n                                 t
||||CL (τ) =     ckBnk(τ )  = X (t),t ∈ [0,^t],τ = - ∈ [0,1]
{         k=0                                ^t
||          n∑      n                          t − ^t
||(CR (u ) =    dkB k(τ)  = X (t),t ∈ [^t,1],u = 1-−-t ∈ [0,1]
          k=0

Definizione 2.6 (Operatori identitÓ e shift).  
Si definiscono gli operatori I,E tali che Ibi = bi,Ebi = bi+1.

Osservazione 2.4 (Rappresentazione alternativa delle curve di BÚzier). I punti della curva di BÚzier (algoritmo di de Casteljau) si possono rappresentare a partire dagli operatori shift e identitÓ: bir = [(1 − t)I + tE ] r1b i.

I coefficienti ck e dk si ottengono dall’algoritmo di de Casteljau come:

      k ^          n−k ^
ck = b0(t),  dk = bk  (t)

infatti Il caso CR si dimostra in maniera analoga :

CL(τ) = k=0nb 0k(^t )(n )

  kτk(1 τ)nk
= k=0n[             ]
 (1 − ^t)I + ^tEkb 0(  )
  n
  kτk(1 τ)nk
= k=0n(  )
 n
 kτk[    ^     ^ ]
 (1 − t)I + tEk(1 τ)nkInkb 0
= k=0n(  )
 n
 k[τ 1 − τ ^tI + τ ^tE ]k[(1 − τ )I ] nkb 0
= [                        ]
 τI − τ^tI + τ^tE + I − τ Inb 0
= [     ^    ^  ]
 (1 − t)I + tEnb 0

: Subdivision
 
1function [b1,b2] = Bezier˙subdiv(b,t) 
2[d,n] = size(b) 
3f = t 
4e = 1-t 
5b1(:,1) = b(:,1) 
6b2(:,n) = b(:,n) 
7b = b 
8for i=1:d 
9     C(:,:,i) = b(1:n-1,i)*e+b(2:n,i)*f 
10end 
11for i=1:d 
12     b1(i,2) = C(1,1,i) 
13     b2(i,n-1) = C(n-1,1,i) 
14end 
15for r=2:n-1 
16     for row=1:n-r 
17          for i=1:d 
18               C(row,:,i) = C(row,:,i).*e+C(row+1,:,i).*f 
19          end 
20     end 
21for i=1:d 
22     b1(i,r+1) = C(1,1,i) 
23     b2(i,n-r) = C(n-r,1,i) 
24end 
25end 
26endfunction


PIC

Figura 2.12: Esempio di subdivision



PIC

Figura 2.13: Altro esempio di subdivision



PIC

Figura 2.14: Esempio di curva tridimensionale con curvatura e torsione



PIC

Figura 2.15: Altro sempio di curva tridimensionale con curvatura e torsione


Curve di BÚzier razionali

Le curve di BÚzieri razionali sono un’estensione alle curve di BÚzier che mira a migliorarne la controllabilitÓ.

Coordinate omogenee

Un insieme Ŕ rappresentabile mediante coordinate omogenee se Ŕ possibile associare un vettore numerico (x1,...,x2) ad ogni elemento dell’insieme, in modo tale che due vettori distinti siano associati a due oggetti diversi solo se linermente indipendenti.

Definizione 2.7 (Spazio proiettivo). Lo spazio proiettivo P(K) associato ad uno spazio vettoriale K Ŕ definito come l’insieme dei suoi sottospazi vettoriali di dimensione uno (rette).

Se ogni vettore di uno spazio vettoriale di dimensione finita Ŕ rappresentabile tramite le sue coordinate (x1,...,xn), ogni retta Ŕ descrivile come lo span lineare di un vettore non nullo: λ(x1,...,xn).

Interpretazione geometrica in d = 3 Si immagini un piano E2 immerso in uno spazio E3 ed allineato col sistema di riferimento, ad esempio perpendicolare all’asse z. I punti di E3 sono rappresentabili tramite l’insieme delle triple (x,y,z) oppure, a partire dai punti di E2, come unione dei punti al finito ed i punti all’infinito:

{                 2}   {                    2 }
 [x,y,1]|(x,y) ∈ E   ∪  [x,y,0 ]|(x,y) ∈ P (E  )

Analogamente Ŕ possible rappresentare i punti del piano a partire dai punti dello spazio attravero un’operazione di proiezione: (x,y,z) (x
z,y
z, 1).

Curve di BÚzier razionali Una curva di BÚzier razionale in uno spazio Ed Ŕ la proiezione di una curva di BÚzier classica definita in uno spazio Ed+1. In particolare si definisce Per comoditÓ di notazione si intende lo spazio proiettivo perpendicolare alla prima componente del vettore, ovvero (x,y,z,...,cn) (1,cxn,ycn,...) :

        (         )
        |  y0(t)  |    n∑                    (  1 )
X (t) = |(    ...    |) =     BiBni (t),  Bi = βi       ,  bi ∈ Ed
          y   (t)      i=0                      bi
           d+1

in cui i Bi sono vettori in Ed+1. La prima componente della curva (d + 1)-dimensionale Ŕ pertanto:

       ∑n
y0(t) =    βiBn (t)
       i=0    i

ed applicando direttamente il concetto di proiezione si ottiene la curva proiettiva in Ed:

       (       )    ∑n       n
          y0(t)        βibiB i (t)
X (t) = ||    ..  ||  = i=0----------,  βi > 0
       (    .  )     ∑n     n
          yd(t)         βiB i (t)
                     i=0

Osservazione 2.5. Se β0 = β1 = ... = βn, X(t) i=0nb iBin(t) ovvero rimane una BÚzier definita dai punti di controllo bi.


: Algoritmo di Casteljau per curve di BÚzeir razionali
 
1function Curve = RATBezier˙casteljau(b,beta,t) 
2[d,n] = size(b) 
3n2 = length(beta) 
4B=[beta] 
5for i=1:d //Aumenta di una dimensione 
6     B = [B;beta.*b(i,:)] 
7end 
8C = Bezier˙casteljau(B,t) 
9for i=1:d 
10Curve(i,:) = C(i+1,:)./C(1,:) //Proietta 
11end 
12endfunction

ProprietÓ

Simmetria

Sia s = 1 t. La proprietÓ segue direttamente da quella dei polinomi di Bernstein:

        ∑n                   ∑n                   ∑n
            βn−ibn−iBni (s)      βn−ibn−iBnn− i(t)       βjbjBnj (t)
        i=0--------------    i=0-----------------  j=0----------
X (s) =   ∑n              =    n∑                =  ∑n
             βn −iBni (s)           βn−iBnn−i(t)          βjBnj (t)
          i=0                  i=0                  j=0

Convex Hull

Definendo la curva come:

       ∑n                  βiBni (t)         ∑n     n
X (t) =    bivi(t),   vi(t) =   w(t)  ,  w (t) =    βiB i (t)
       i=0                                   i=0

si ha che:

 n              n
∑  v (t) = --1--∑  β Bn (t)=  1
i=0  i     w (t)i=0  i i
               ◟----◝◜----◞
                   w(t)

ovvero la curva Ŕ una combinazione convessa dei punti di controllo ed appartiene alla convex hull di questi.

Invarianza per trasformazioni affini

Deriva direttamente dal fatto che la curva Ŕ una combinazione baricentrica di punti; la dimostrazione Ŕ analoga a quella per le curve di bÚzier classiche.

Linear precision

In generale non vale, infatti:

       ∑n                      ∑n                   1  ∑n
X (t) =    (b0 + inv)vi(t) = b0 +    (invi(t))v = b0 + -----   inβiBni (t)v
       i=0                     i=0                 w (t) i=0

da cui:

X (t) = b0 + --t--v
            w (t)

Se i βi sono omogenei per˛ la proprietÓ si mandiene valida.


PIC

Figura 2.16: Linear precision: la distribuzione dei pesi (ciano) varia la densitÓ dei punti della curva (blu) rispetto al parametro t.


Valori agli estremi

La curva interpola i punti di controllo agli estremi:

X(0) = β0b0-
 β0 = b0
X(1) = βnbn-
 βn = bn

Variation diminishing

I punti della curva X(t) = i=0nb ivi(t), con vi(t) = βiBni (t)
--------
  w (t),w(t) = i=0nβ iBin(t), che giacciono sulla retta ax + by + c = 0 sono:

a i=0nb ixvi(t) + b i=0nb iyvi(t) + c i=0nv i(t) = 0
i=0n(ab ix + bbiy + c)vi(t) = 0
i=0nβi(abix +-bbiy +-c)
       w(t)Bin(t) = 0
quindi, ponendo di = βi(abix-+-bbiy +-c)
      w (t) si rientra nel caso della curve di BÚzier classica, per cui la proprietÓ Ŕ dimostrata.

Controllo pseudolocale

Modificando uno dei pesi:

^
βj = βj + Δ βj

si ottiene la curva:

        ∑n        n             n
^X (t) =  --i=0-βibiB-i (t) +-Δ-βjbjB-j (t)
          ∑n     n           n
             βiB i (t) + Δ βjB j (t)
          i◟=0--------◝◜---------◞
                    l(t)

Confrontando la curva modificata con quella originale si ha:

(X^(t) X(t))l(t) = i=0nβ ibiBin(t) + Δβ jbjBjn(t) i=0nβ ibiBin(t) X(tβ jBjn(t)
= ΔβiBin(t)[b j X(t)]
La differenza tra le curve vale quindi:
                        Δ βiBni (t)
( ^X (t) − X (t)) = ∑n------n------------n--[bj − X (t)]
                   i=0 βiBi (t) + Δ βjBj (t)

Modificando un peso si ha quindi un effetto su ogni punto della curva, direttamente proporzionale alla distanza bj X(t): la curva tende ad avvicinarsi al punto bj se Δβj Ŕ positivo, allontanarsi altrimenti.


PIC

Figura 2.17: Esempio di differenti distribuzioni dei pesi


Derivata

Sia:

               ∑n         n
X (t) = c(t)-= -∑i=0βibiB-i (t)
        w(t)      ni=0βiBni (t)

Si pu˛ calolare la derivata di X(t):

 ′        ′             ′        c′ − Xw ′
c = (Xw  ) = X˙w   + Xw   ⇒  ˙X =  ---------
                                    w

Le derivate di c(t) e w(t) sono note:

c(r) = nr i=0nrΔrβ ibiBinr
w(r) = nr i=0nrΔrβ iBinr
quindi, sostituendo in (t) si ottiene:
          [∑                    (∑                )     ]
        n    ni=−0rΔ βibiBni −r(t) −    ni−=r0 ΔβiBni−r(t) X (t)
X˙(t) = ------------------∑n-------n---------------------
                            i=0βiB i (t)


PIC

Figura 2.18: Un altro esempio di differenti distribuzioni dei pesi


Tangenza agli estremi

La tangenza agli estremi Ŕ mantenuta, infatti:

(0) = nβ1
β--
 0Δb0
(1) = nβn−1-
 βnΔbn1
Se si volesse ottenere una curva chiusa si deve quindi imporre:
(βn− 1 + β1)b0 = βn− 1bn− 1 + β1b1

Curve di BÚzier razionali in forma standard

Definizione 2.8 (Curva di BÚzier razionale in forma standard). Una curva di BÚzier razionale si dice in forma standard se i pesi associati agli estremi del poligono di controllo sono unitari.

Per trasformare una curva razionale in una in forma standard, si utilizza la riparametrizzazione:

t = (1-+-α-)u- ≡ ϕ(u),   α ∈ ℝ, u ∈ [0,1 ]
     1 + αu

La derivata della funziona di riparametrizzazione ϕ(u) vale:

             1                                      1 + α
ϕ ′(u) =  --------2 [(1 + α)(1 + α) − α(1 + α )u ] =--------2-
         (1 + αu )                                 (1 + αu )

quindi la funzione Ŕ monotona per α > 1, e la riparametrizzazione Ŕ ammissibile in quel range.
Ha senso quindi applicarla:

                                                                 Bn
                  (  )                                ◜(--)-------◞i◟----------◝
           ∑n      n  (1-+-α)iui-(1-−-u)n−i-   ∑n       n   i      n− i      i
              βibi i   (1 + αu )i(1 + αu )n−i      βibi  i  u (1 − u )  (1 + α )
X (t(u)) = i=0n---(--)--------i-i--------n−i- = i=0n---(--)--------------------- =
            ∑      n  (1 +-α-)u-(1-−-αu-)--     ∑      n   i      n− i      i
               βi  i       (1 + αu )n              βi  i  u (1 − u )  (1 + α )
            i=0                                 i=0

   ∑n                   n∑
      βibi(1 + α )iBni       ˜βibiBni
=  i=0----------------=  i=0--------
   ∑n           i n     ∑n  ˜  n
       β◟i(1-+◝◜-α-)◞ Bi        βiB i
   i=0     ˜βi           i=0

e richiedere che siano verificate:

(
|||  ˜βn ≡ βn (1 + α )n = 1,  α  = √-1--−  1
{                               nβn
||
|(  ˜β0 ≡ β0(1 + α)0 = 1,    se β0 = 1
(2.1)

La forma standard si pu˛ ottenere solo nel caso in cui la condizione β0 = 1 sia giÓ verificata inizialmente.


: Algorito per il calcolo dei beta normalizzati
 
1function beta=RATBezier˙scalebeta(beta) 
2n = length(beta); 
3beta=beta/beta(1) 
4alpha=(1/beta(n))^(1/n-1); 
5for i=1:n 
6     beta(i) = (alpha^(i-1))*beta(i); 
7end 
8endfunction

Sezioni coniche

Le sezioni coniche sono rappresentabili utilizzando curve di BÚzier di secondo grado.

Coordinate baricentriche bidimensionali

Un punto p che giace nella convex hull di altri tre (b0,b1,b2) ogniuno ha componenti xi,yi che delimitano un triangolo Ŕ identificabile univocamente con un sistema di coordinate baricentriche (u0,u1,u2) a riferimenti b0,b1,b2:

(
{ p  =  u0b0 + u1b1 + u2b2
(
  1  =  u0 + u1 + u2

Il sistema di coordinate sopra descritto Ŕ univocamente determinato, si dimostra che in:

 (           ) (   )    (  )
  x0   x1  x2    u0       x
 |( y0  y1  y2|) |( u1|) =  |( y|)
   1   1   1     u        1
◟------◝◜------◞   2     ◟ ◝◜-◞
       A                  p

la matrice A Ŕ nonsingolare e quindi il sistema ammette soluzione unica. Il determinante di A Ŕ legato all’area del triangolo b0b1b2; fissando una base per esso si ha:

           1      1
area(A ) = 2|Δ | = 2|(b1 − b0) × (b2 − b0)|

Applicando il metodo ci Cramer:

                  (            )
                    x  x1  x2
detA0 (x,y) = det |( y  y1   y2 |) ;
                    1   1   1

                  (            )
                  | x0  x  x2  |
detA1 (x,y) = det ( y0  y   y2 ) ;
                     1  1   1

                  ( x   x   x  )
                  |  0    1    |
detA2 (x,y) = det ( y0  y1  y  ) ;
                     1   1  1

e quindi:

     detAi    Δi (x, y)
ui = det-A-≡  ---Δ----

dove Δi = det Ai(x,y) e Δ = det A. Il determinante di A Ŕ non nullo quindi le ui sono tutte definite. Il sistema perde per˛ l’unicitÓ nel caso in cui tutti e tre i punti siano allineati.

Rappresentazione di curve coniche

In generale una curva conica Ŕ rappresentabile come:

ax2 + 2bxy +  cy2 + dx + ey + f = 0

o in forma matriciale, dato il punto p = (  )
  x

  y :

  (a   b)
pT        p + (d,e)p + f = 0
    b  c

Il determinante della matrice determina il carattere della conica:

           (
    (    ) || > 0   ellisse
     a  b  {
det  b  c  || = 0   parabola
           ( < 0   iperbole

Si ricordano alcune definizioni:

Archi di coniche con curve di BÚzier

Siviluppando la formula di una curva di BÚzier razionale in forma standard di secondo grado e etichettandone in modo opportuno i termini:

          u0(t)      -u1(t)--     u2(t)
          ◜◞1◟ ◝    ◜  ◞◟2  ◝    ◜ ◞2◟◝
X (t) = b0-B0(t)+b1-β1B-1(t)+b2-B-2(t) ,  v0(t) =  u0(t),v1(t) = u1(t),v2(t) = u2(t)
          B20(t) + β1B21(t) + B22(t)             w(t)          w(t)         w (t)
          ◟----------◝◜----------◞
                    w (t)

i punti della curva possono essere espressi in maniera semplice come:

X (t) = b0v0(t) + b1v1(t) + b2v2(t)

Dalle proprietÓ dei polinomi di Bernstein segue direttamente che vi(t) = 1 ovvero ogni punto della curva Ŕ una combinazione baricentrica dei tre punti di controllo, leggibile anche come coordinate baricentriche rispetto agli stessi.
Calcolando il quadrato della componente v1 si ottiene:

           t2(1 − t)2
v21(t) = 4β1----------=  4β1v0(t)v2(t) ⇒ v21(t) − 4β1v0(t)v2(t) = 0
             w2 (t)

sostituendo ai vi la relativa notazione risultante dal teorema di Cramer vista nella sezione precedente (vi = Δi(x,y)Δ), si ha che punti della curva rispettano:

  2
Δ 1(x, y) − 4 β1Δ0 (x, y)Δ2(x, y) = 0

che un’equazione di secondo grado, ovvero una curva conica che cambia comportamento in base al valore di β1:

   (|
   |{ < 1  ellisse
β1 | = 1  parabola
   |( > 1  iperbole


: Algoritmo per il calcolo di curve coniche
 
1function C = RATBezier˙coniche(b,s,ttab) 
2// beta = s/(s-1) 
3C = RATBezier˙casteljau(b,[1,s,1],ttab) 
4endfunction

: Calcolo delle coniche mostrate in fig (2.19)
1pc = [1,0,1;1,1,1]; 
2ellipse1 = RATBezier_coniche(pc,0.5,xcol); 
3ellipse2 = RATBezier_coniche(pc,0.5,xcol); 
4parabola = RATBezier_coniche(pc,1,xcol); 
5iperbole = RATBezier_coniche(pc,3,xcol); 
6pc2 = [.5,0,1;1,0,.5]; 
7circ1 = RATBezier_coniche(pc2,0.5,xcol); 
8circ2 = RATBezier_coniche(pc2,0.5,xcol); 
9}


PIC

Figura 2.19: Rappresentazione di curve coniche


Capitolo 3
Curve SPLINE

BÚzier-spline

Uno dei possibili modi per migliorare la controllabilitÓ e localitÓ di una curva di BÚzier Ŕ quello di considerare sottointervalli della curva come curve separate, preoccupandosi di gestire i raccordi. Dato un intervallo di definizione della curva [a,b] ed una sua partizione in m sottointervalli τ : a u1 < u2 < ... < um < um+1 b disgiunti ed adiacenti, si definisce la curva polinomiale a tratti:

                d            ⋃
X (u) : [a, b] → E ,  X (u ) =   Xi (t),t ∈ [ui,uu+1 )

in cui:

         ∑n      n (u − ui)
Xi (u) =    bk,iB k  ------  ,  hi = ui+1 − ui
         k=0          hi

e bk,i Ŕ il k-esimo punto di controllo della curva i-esima. I punti della curva X(ui) sono detti nodi (knot) o punti di giunzione. La condizione di continuitÓ si ottiene imponendo:

Xi(ui+1) = Xi+1 (ui+1) ⇒  bn,i−1 = b0,i

grazie alla proprietÓ di end-point interpolation; ad esempio date due curve disconnesse si pu˛ aggiungere un punto a tutte e due che stia spazialmente tra l’ultimo della prima e il primo della seconda.


PIC

Figura 3.1: BeziÚr-spline con continuitÓ C0: anche visivamente la curva Ŕ continua ma ha dei punti in cui non Ŕ derivabile. A destra la hodograph, visibilmente disconnessa.



PIC

Figura 3.2: BeziÚr-spline con continuitÓ C1: la curva appare continua. A destra la hodograph, di continuitÓ C0.


ContinuitÓ dei raccordi

Per aumentare la regolaritÓ della curva si pu˛ voler imporre la continuitÓ delle derivate; ad esempio per ottenere una curva in C1 si impone:

(
||| Xi(ui+1) = Xi+1 (ui+1)
|{
|
|||( -d-X (u   ) = -d-X   (u   )
  du  i  i+1    du   i+1  i+1

Si ha che:

d--
duXi(ui+1) = -d-
duXi+1(ui+1)
d
---
duXi(1) d
---
duu − u
-----i
  hi = d
---
duXi+1(0)u − u
-----i+1-
  hi+1
Δbn −1,i
--h----
    i = Δb0,i+1
--h----
   i+1
ovvero:
                hi+1-
b1,i+1 = b0,i+1 +  h  (bn,i − bn−1,i).
                  i

Analogamente per avere continuitÓ C2 si impone anche:

  2             2
d--Xi(ui+1) = -d-Xi+1 (ui+1 )
du            du

da cui:

Δ2bn-−2,i    Δ2b0,i+1-
   h2    =   h2
    i          i+1

ovvero:

Δb      = Δb     +  hi+1(Δb   −  Δb     ).
   1,i+1      0,i+1    hi     n,i      n−1,i

Quanto appena descritto Ŕ un metodo per ottenere una curva di BÚzier definita a tratti parametrizzata in modo tale che rispetto al parametro u la velocitÓ di percorrenza della curva sia proporzionale alla dimensione dei diversi intervalli assegnati ad ogni curva. In pratica Ŕ comodo considerare l’intervallo [a,b] suddiviso equamente tra i vari tratti di curva (curva polinomiale a tratti uniformi); in questo caso la funzione di riparametrizzazione Ŕ la stessa per tutte (a meno dell’intervallo) e non c’Ŕ bisogno di riparametrizzarle, perchŔ si pu˛ dare per scontato che ogni curva sarÓ valutata in [0, 1] quando u Ŕ nell’intervallo corrispondente. Si ha che:
Le derivate delle due curve sono facili da derivare:

d-
dtXi(t) = -d-
duXi+1(u)
n k=0n1Δb k,iBkn(t) = m k=0m1Δb k,i+1Bkn(u)
Sappiamo di voler allineare l’ultimo punto della prima curva con il primo della seconda, quindi:
n k=0n1Δb k,iBkn1(1) = m k=0m1Δb k,i+1Bkn1(0)
nΔbn1,i = mΔb0,i+1
in quanto in 0 e 1 non si annullano soltanto rispettivamente il primo e l’ultimo polinomio di Bersntein. Da qui, tenendo conto che per avere C0 si impone b n.i = b0,i+1 segue direttamente la relazione tra i quattro vertici:
b1,i+1 = (m  + n)bn,i − mbn −1,i

Analogamente per avere C2

n(n 1) k=0n2Δ2b k,iBkn2(1) = m(m 1) k=0m2Δ2b k,i+1Bkn2(0)
n(n 1)Δ2b n2,i = m(m 1)Δ2b 0,i+1
da cui, imponendo le condizioni per C0 e C1:
b2,i+1 = 2b1,i+1 − bn,i + m-(m-−--1)(bn,i − 2bn−1,i + bn−2,i)
                        n(n − 1)


PIC

Figura 3.3: BeziÚr-spline con continuitÓ C2: visivamente la curva Ŕ continua e regolare. A destra la prima e la seconda hodograph: la prima Ŕ C1 e la seconda C0.


ContinuitÓ geometrica

Due curve si raccordano con continuitÓ geometrica Gs nel punto p se esiste una riparametrizzazione ammissibile per cui le curve hanno continuitÓ Cs in p.
Per avere continuitÓ G1, date due curve X(u),Y (t) e la funzione u(t) : [t 0,t1] [u0,u1]:

d
--
dtX(t1) =  d
--
dtY (t1)
ω1-d-
duX(u1) = d-
dtY (t1), ω1 = -d
dtu(t1) > 0
ovvero, formalmente due curve si raccordano G1 se:
                  d            d
∃ ω1 >  0 tale che--Y (t1) = ω1 --X (u1)
                  dt           du

Il caso G2 Ŕ analogo:

 2
d--
dtX(t1) =  2
d--
d2tY (t1)
d--
du(                )
  d--     -d
  duX (u1)dt u(t1) = d2-
d2tY (t1)
d
---
duX(u1)d2
-2-
d tu(t1) +  d2
--2-
d uX(u1)( d      )
  --u (t1)
  dt2 = d2
-2-
d tY (t1)
e come prima le due curve hanno raccordo G2 se:
                          d2            d2             d
∃ ω1 >  0,ω2 > 0, tali che -2-Y (t1) = ω2 -2-X (u1) + ω1---X (u1)
                          d t           d u           du


PIC

Figura 3.4: BeziÚr-spline con continuitÓ G2: visivamente la curva Ŕ continua e regolare. A destra la prima e la seconda hodograph: la prima Ŕ C1 e la seconda C0 dopo che sia stata applicata la scalatura relativa al segmento di spline.


Come prima, nel caso pratico considerando equiparametrizzate le curve si ottiene, per le continuitÓ G1 e G2:

nΔbn1,i = ω1mΔb0,i+1 b1,i+1 = (n +-m-ω1)bn,i −-nbn−1,i-
         m ω1
n(n 1)Δ2b n2,i = m(m 1)Δ2b 0,i+1 b2,i+1 =
(n(n-−-1)-−-m-(m-−-1-)ω2-)bn,i −-2n-(n-−-1)bn−1,i +-n(n-−-1)bn−2,i-+-2ω2m-(m-−-1)b1,i+1
                                   m (m − 1 )ω2


PIC

Figura 3.5: BeziÚr-spline con continuitÓ G2: visivamente la curva Ŕ continua e regolare. A destra la prima e la seconda hodograph: la prima Ŕ C1 e la seconda C0 dopo che sia stata applicata la scalatura relativa al segmento di spline.


Implementazioni

Vedi appendice A.1.

B-SPLINE

Curve SPLINE monovariate

In generale, una curva SPLINE monovariata Ŕ una funzione S : [a,b] d polinomiale a tratti. In particolare, definita una partizione τ dell’intervallo [a,b]:

       ⋃
[a,b] =   Ii,   Ii = [τi,τi+1), τ : a ≡ τ1 < τ2 < ...< τn < τn+1 ≡  b

Ii Ij = con ij, Ii Ii+1 [ti,ti+2) la funzione S Ŕ rappresentabile come:

     ⋃                       d
S =    Pi,  Pi : [τi,τi+1] → ℝ ,  Pi ∈ Πm

Si indica con Sm,τ una generica funzione SPLINE sulla partizione τ e composta di polinomi tutti di grado m. Considerando l’intervallo [a,b] suddiviso in n sottointervalli, si Ŕ definito uno spazio di dimensione:

dim (Sm,τ) = n × (m  + 1)

in quanto ogni polinomio di grado m Ŕ identificato da m+1 coefficienti. In pratica, volendo rappresentare una curva, si vuole che S sia continua (ovvero lim tti+1Pi(t) = Pi+1(ti+1)) e che anche le derivate si raccordino (S Cm1) ; questo lascia uno spazio di dimensione:

dim (Sm, τ) = n × (m + 1) − (n − 1) × m =  n + m

infatti ci sono n 1 punti di raccordo per ogniuno dei quali si impogono m condizioni di di raccordo (una di continuitÓ e le derivate).

Base delle potenze troncate

Si definiscono:

             (
        m    {(x − τi)m+,     x < τi
(x − τi)+ =  (
              0,             x ≤ τi

Definizione 3.1 (Base delle potenze troncate). Sulla partizione τ : a τ1 < ... < τl b si definisce base delle potenze troncate l’insieme di polinomi:

     0  1     m         m              m
{1,x ,x  ,...,x  ,(x − τ1)+ ,...,(x − τl−1)+}

Teorema 3.1. L’insieme delle potenze troncate Ŕ una baste per lo spazio delle spline Sm,τ.

Dimostrazione. Si deve dimostrare che i polinomi della base sono l + m linearmente indipendenti, ovvero che:

 m        l−1
∑  a xi + ∑  b(x − τ )m = 0 ⇔  a =  b = 0
i=0 i     i=1  i     i +         i    i

Se x I1 [τ12] si ha:

m∑     i  ∑l−1          m
   aix +     bi(x◟-−◝◜τi)+◞=  0
i=0        i=1      =0

ovvero la tesi Ŕ dimostrata in quanto la base Ŕ quella canonica. Per x I2 si ha:

∑m    i            m   l∑−1          m
   aix +  b1(x − τ1)+ +     bi(◟x-−◝τ◜i)+◞= 0
i=0                    i=2      =0

ma la propritÓ vale per gli ai, il che implica che debba valere anche b1. Procedendo allo stesso modo si dimostra la tesi per tutti i bi __

Osservazione 3.1 (Derivata). (x τ)+m Cm e dm-
dx(x τ)+m = m! se x > τ.

Base delle B-SPLINE

La base delle potenze troncate presenta problemi di malcondizionamento; questi possono essere risolti definendo, sullo stesso concetto, una base a supporto compatto.
Si definisce una partizione di [a,b] vettore dei nodi pi¨ rilassata della precedente:

ν : a ≡ ν0 ≤ ν1 ≤ ...≤  νl ≡ b,   l intervalli

e, data la base:

          (
          { 1,  νi ≤ x ≤ νi+1
Ni,1(x) = (
            0,   altrimenti

si definisce ricorsivamente:

          --x-−-νi---           -νi+k-−-x--
Ni,k(x ) = νi+k− 1 − νiNi,k−1(x) + νi+k − νi+1Ni+1,k−1(x)

Nota (Abuso di notazione). Per comoditÓ, se uno dei denominatori di un termine in Ni,k(x) si annulla, questo viene considerato nullo.

Esempio 3.1 (Caso k = 2).

          (
          || -x-−-νi--
          |||| νi+1 − νi,    x ∈ [νi,νi+1)
          ||{
Ni,2(x) =
          |||| -νi+2 −-x-,  x ∈ [νi+1, νi+2)
          |||| νi+2 − νi+1
          (0              altrimenti

Esempio 3.2 (Caso k = 3).

          (  x − νi
          |||| ---------Ni,2(x),                         x ∈ [νi,νi+1)
          |||| νi+2 − νi
          ||||
          |||| -x-−-νi--         --νi+3-−-x--
          ||{ νi+2 − νiNi,2(x) + νi+3 − νi+2Ni+1,2(x),  x ∈ [νi+1,νi+2)
Ni,3(x) =
          ||||
          |||| -νi+3-−-x--Ni+1,2,                        x ∈ [ti+2,ti+3)
          |||| νi+3 − νi+1
          ||||
          ||(
            0                                        altrimenti


PIC (a) k=2 PIC (b) k=3
PIC (c) k=5 PIC (d) k=10

Figura 3.6: Base delle BSPLINE, in rosso Ŕ evidenziata la somma delle componenti.


ProprietÓ

Supporto compatto

Ni,k(x) Ŕ a supporto compatto, infatti Ni,k = 0 se t∕∈[ti,ti+k].

Non negativitÓ

Ni,k(x) 0t [ti,ti+k], in quanto somma di soli termini positivi in quell’intervallo.

Partizione dell’unitÓ
∑n
   Ni,k(t) = 1∀t ∈ [νk− 1,νn+1 ]
i=0

Per induzione, la tesi Ŕ vera banalmente per Ni,1(t) e t [t0,tn+1]; supponendola vera per Ni,k1(t) , t [νk2n+1] per Ni,k(t) e t [νk1n+1] vale:

i=0nN i,k(t) = i=0n[   t − νi               νi+k − t            ]
 -----------Ni,k−1(t) + -----------Ni+1,k−1(t)
 νi+k−1 − νi           νi+k − νi+1
= t − ν0
--------
νk−1−ν0 N◟0,k◝−◜1(t)◞0t[ν0k1] + i=1n[                        ]
   t − νi      νi+k−1 − t
 -----------+  -----------
 νi+k−1 − νi   νi+k− 1 − νiNi,k1(t)
+ -νn+k-−--t--
νn+k − νn+1 N       (t)
◟-n+1,k◝−◜-1-◞0t[νn+1n+k1]
= 0 + i=1nN i,k1(t) + 0 + N◟0,k◝−◜1(t)◞=0,t[νk1n+1]
= i=0nN i,k1(t)
e in quanto [νk1n+1] [νk2n+1], la proprietÓ Ŕ dimostrata.


PIC

Figura 3.7: Partizione dell’unitÓ: la somma delle basi Ŕ unitaria per t [νk1n+1]. (in questo caso k = 3,T = {0, 1, 2, 3, 4, 5, 6},n = 3)


Derivata

La derivata della base delle B-SPLINE Ŕ definita come segue:

                  [                         ]
dNi,k(x) = (k − 1)  Ni,k−1(x)--− Ni+1,k−1(x)
   dx               νi+k− 1 − νi   νi+k − νi+1

Indipendenza lineare
∑n
   aiNi,1(x) = 0 ⇔  ai = 0, ∀i
i=0

si verifica banalmente, t [a,b]. Analogamente, se k = 2:

∑n
   aiNi,2(x) = 0 ⇔  ai = 0, ∀i
i=0

infatti nei nodi (x = νj) solo Nk1,2(νj)0 e aj deve essere diverso da zero, per ogni j. Per dimostrare la proprietÓ per k > 2 occorre procedere per induzione, utilizzando i precedenti come casi base. Tenendo conto che vale:

∑n                         ∑n   [ N     (x)    N       (x)]
   aiNi,k(x) = 0 ⇒ (k − 1)    ai ---i,k−1---- − --i+1,k−1---  = 0
i=0                        i=0    νi+k−1 − νi   νi+k − νi+1

espandendo e risommando i termini della sommatoria nella derivata si ottiene:

        n∑  [  a −  a   ]
(k − 1)     ---i----0-- Ni,k−1(x) = 0
        i=1  νi+k−1 − νi

da cui si deve avere ai = ai+1 ovvero tutti i coefficienti coincidenti. Riportando questo risultato alla sommatoria di partenza:

 ∑n
a   Ni,k(t) = 0 ⇔ a =  0
 i=0

in quanto gli Ni,k(x) sono una combinazione baricentrica.

Generalizzazione della base delle B-SPLINE

La definizione di SPLINE data nella sezione precedente richiede che la curva abbia lo stesso tipo di continuitÓ su tutti i nodi. Questo vincolo nella pratica pu˛ essere troppo stringente, e talvolta conviele usare uno spazio delle spline allargato.
Data la partizione:

τ : a ≡ τ0 < τ1 < ...< τn ≡ b,  Ii = [τi,τi+1)

ed un vettore:

M  = m1, ...,mn −1,  0 ≤  mi ≤ k

Questo vettore sta ad indicare il tipo di raccordo che si vuole nei punti τ2,..,τn, se mi = 0 non si chiede neanche la continuitÓ, se mi = s si chiede la continuitÓ delle derivate dalla zeresima alla s. Si definisce:

Sm,τ,M ={f : [a,b] d per cui  p 1(t),...,pn(t) Πm,f(t)|Ii = pi(t), e
pi1(s)(τ i) = pi(s)(τ i),s = {0,...,mi 1} se mi > 0,
nessuna condizione imposta in τi se mi = 0}
per cui se mi = k Sm,τ,M = Sm,τ mentre se mi = 0, Sm,τ,M = Pm,τ ovvero non sussiste neanche la condizione di raccordo semplice. Si ha quindi che Πm Sm,τ m,τ,M Pm,τ, e lo spazio appena definito ha dimensione Come in precedenza, n intervalli, m grado :
                             n∑−1          n∑−1
dim(Sm, τ,M ) = (m  + 1) × n −    mi  = n +     nm--−-nmi-+--mi-= n +  μ
                             i=1          i=1      n − 1
                                          ◟--------◝μ◜--------◞

Per comoditÓ, usualmente il vettore M non si riferisce direttamente al grado di derivabilitÓ nei nodi, ma indica direttamente la molteplicitÓ del nodo i-esimo nel vettore dei nodi, e prende il nome di vettore delle molteplicitÓ. In una curva di grado k la molteplicitÓ del nodo (m) e il grado di derivabilitÓ della curva (d) sono legate dalla relazione:

m =  k − d + 1
(3.1)

ovvero per ottenere una curva a derivabilitÓ massima (che coincide con le B-SPLINE classiche) la molteplicitÓ dei nodi deve essere minima (1) e viceversa, per ottenere la derivabilitÓ minima (nel caso pratico si richiede sempre la continuitÓ C0) la molteplicitÓ Ŕ massima (k).

Definizione 3.2 (Vettore dei nodi esteso). A partire dall’insieme M e dal vettore dei nodi ν (che ha n intervalli) si definisce il vettore dei nodi esteso :

T : t0 ≤ ... ≤ tk−2≤ tk−1 ≤ ...≤  tn+ μ+1 ≤ tn+μ+2 ≤ ... ≤ tn+μ+k
   ◟-----◝◜-----◞   ◟-------◝◜-------◞   ◟--------◝◜--------◞
         (a)                 (b)                    (c)

in cui:

Come prima, si definisce ricorsivamente la base per le spline:

Ni,1(t) = (|
|{1   t ∈ [ti,ti+1)
|
|(0   altrimenti o se t = t
                    i    i+1
Ni,k(t) = (|
||||0                                            ti = ...= ti+k
||||
||||    --ti+k-−-t-
||||0 + t   −  t  Ni+1,k−1(t)                    ti = ...= ti+k −1 < ti+k
|{     i+k   i+1
|
||||---t −-ti- N    (t) + 0                      t < ... < t     =  t
||||ti+k−1 − ti i,k− 1                             i        i+k −1    i+k
||||
||||   t − ti              ti+k − t
|(---------- Ni,k− 1(t) + ----------Ni+1,k−1(t)  ti < ... < ti+k
 ti+k−1 − ti           ti+k − ti+1

Osservazione 3.2 (Caso particolare: polinomi di Bernstein).  
Ponendo [τ01] = [0, 1], il vettore dei nodi estesi per una spline di grado k diventa:
t : τ0 = t0 = ... = tk1 < tk = ... = t2k1 = τ1.
Dalla definizione di Ni,k, per i < k segue:

           t − t            t    − t
N1,k(t) = -----1-N1,k−1(t)+ -1+k-----N2,k−1(t)
          tk − t1 ◟--◝◜---◞  t1+k − t2
                    =0

quindi, ricorsivamente:

              t − t              t   − t
N2,k−1(t) = ------2--N2,k−2(t)+ -2+k-----N3,k−2(t)
            tk+1 − t2◟---◝◜---◞  t2+k − t3
                        =0

Si pu˛ scrivere la formula generica:

Ni,k(t) = t1+k −-t-t2+k −-t...Ni+k− 1,1(t)
          t1+k − t2t2+k − t3

ma t0 = ... = tk1 e tk = ... = t2k1 da cui:

          (       )k− 1
Ni,k(t) =  -τ1 −-t     Ni+k− 1,1(t)
           τ1 − τ0

analogamente, per i k:

                                                            (        )k−1
N   (t) = --t-−-ti--N     (t) + -ti+k-−-t--N       (t)=  ...=   -t −-τ0     N      (t)
  i,k      ti+k− 1 − ti i,k−1      ti+k − ti+1 ◟-i+1,◝k◜−1--◞         τ1 − τ0      2k−1,1
                                              =0

Da qui, sostituendo τ0 = 0 e τ1 = 1 si ottiene:

N   (t) = (1 − t)k−1 ≡ Bk −1(t),  N   (t) = tk− 1 ≡ Bk− 1,   t ∈ [τ ,τ ]
  1,k                    1          k,k              k−1         0  1

In generale si pu˛ dimostrare che:

             k−1
Ni+j,k(t) = B k  (s),  t ∈ [ti+j,ti+j+1],s ∈ [τ0,τ1]


: Algoritmo per il calcolo della base delle BSPLINE
 
1function base = Bspline˙base(i,k,T,t) 
2 
3 
4if k==1 then 
5   if tíT(i+1) & T(i)í=t then base=1; else base=0; end; 
6else 
7      if  T(i)==T(i+k) then 
8         base=0; 
9      else 
10         if (T(i+k-1)˜=T(i)) then 
11            lp = ((t-T(i))/(T(i+k-1)-T(i)))* Bspline˙base(i,k-1,T,t) 
12         else 
13            lp = 0 
14         end 
15 
16         if (T(i+k)˜=T(i+1)) then 
17            rp = ((T(i+k)-t)/(T(i+k)-T(i+1)))* Bspline˙base(i+1,k-1,T,t) 
18         else 
19            rp = 0 
20         end 
21 
22         base=lp+rp; 
23      end 
24end 
25endfunction

Derivata

La definizione di derivata Ŕ l’estensione diretta di quella della base B data in precedenza:

                  [                        ]
dNi,k(t)             Ni,k− 1(t)    Ni+1,k−1(t)
--------= (k − 1)  ---------- − -----------
  dt               ti+k− 1 − ti   ti+k−ti+1


PIC (a) k=3
PIC (b) k=8

Figura 3.8: Due esempi di curve BSPLINE, una di 3 e una di 8 grado. Si nota chiaramente la mancanza di interpolazione agli estremi.


Curve B-SPLINE

La proprietÓ di somma unitaria della B-base estesa in [tk1,tn+1] permette di utilizzarla in modo analogo ai polinomi di Bernstein per definire una curva:

       ∑n
X (t) =    diNi,k(t),  t ∈ [tk−1,tn+1]
        i=0

in cui di sono punti in Ed.


: Algoritmo per il disegno delle BSPLINE: ricorsivo
 
1function Curve = Bspline˙curve(k,T,b,t) 
2lt = length(t) 
3[d,n] = size(b) 
4lT = length(T) 
5Curve=[] 
6for p=1:lt-1 
7      c=[] 
8      count = 1 
9      for i=1:lT-k 
10      c=c+b(:,count).*Bspline˙base(i,k,T,t(p)) 
11      count = count+1 
12      end 
13      Curve=[Curve,c] 
14end 
15endfunction


PIC (a) k=3
PIC (b) k=8

Figura 3.9: MolteplicitÓ dei nodi: i due esempi precedenti di curve, con un nodo a molteplicitÓ massima (derivabilitÓ minima).


ProprietÓ

Invarianza alle trasformazioni affini

                         n∑
Φ (X (t)) = AX  (t) + v =    Φ(di)Ni,k(t)
                         i=0

Controllo locale Se t [tr,tr+1]:

          ∑r
X (t) =        diNi,k(t)
        i=r− k+1

in quanto tutte le altre Ni,k(t) sono nulle, ne consegue che la modifica di di ha effetti sulla curva solo nel corrispondente intervallo ti,ti+k.

Strong ConvexHull Come nel caso dei polinomi di Bernstein, la somma i=0nd iNi,k(t), t [tk1,tn+1] Ŕ una combinazione baricentrica dei di, quindi la curva giace nella convex hull di questi; in pi¨ la proprietÓ di localitÓ garantisce che la curva, per t [tr,tr+1] sta anche nella convex hull dei punti [ti,ti+k], ovvero si ha la convex hull forte.

Linear precision (come conseguenza della strong convex hull) Se si hanno (k 1) punti di controllo allineati dj+i = dj + αiv,r = 0,..,k 1 l’equazione della curva diventa:

                                         (                     )
             j+∑k−1                        j+k∑− 1
X (tj+k−1) =      diNi,k(tj+k− 1) = dj + v (     αi−jNi,k(tj+k−1))
              i=j                           i=j

ovvero se almeno k 1 punti sono allineati la curva tocca il segmento su cui giacciono i punti di controllo, e se almeno k lo sono la curva si schiaccia sul segmento stesso.


PIC (a) k=4, 3 nodi allineati
PIC (b) k=4, 4 nodi allineati

Figura 3.10: Linear precision: esempi con k 1(= 3) e k(= 4) nodi allineati.


Derivata

dX-(t)
  dt = (k 1) i=0nd i[                         ]
  -Ni,k−-1(t)-−  Ni+1,k−1(t)
  ti+k− 1 − ti  ti+k − ti+1
= (k 1) i=1nd  − d
-i----i−1-
ti+k−1 − tiNi,k1(t)
scrivendo il coefficiente -di −-di−1
ti+k− 1 − ti come di[1] (e ricorsivamente d i[r] =  [r−1]   [r− 1]
di----−-di−1-
  ti+k−1 − ti) si pu˛ generalizzare:
 r                        ∑n
d-X-(t) = (k − 1)...(k − r)   d[r]Ni,k−1(t),  t ∈ [tk−1,tn+1]
  drt                     i=1 i

End point interpolation In generale la proprietÓ di end point interpolation non vale per le B-SPLINE, ma la si pu˛ forzare.
Dato il vettore ausiliario t : t0,...,ti,...,ti+k2,...,tk+μ+k, se i k-1 nodi ti,...,ti+k2 coincidono si ha che X(ti) = di1. Quindi, forzando t0 = ... = tk1 = τ0 e tk+μ+2 = ... = tk+μ+k = τn, N0,k(t) = B0k1(        )
  t-−-τ0-
  τ1 − τ0 con t [τ01] e Nn,k(t) = Bk1k1(          )
 -t −-τn−1
 τn − τn−1 con t [τn1n], ovvero la proprietÓ Ŕ verificata.


PIC (a) k=3
PIC (b) k=8

Figura 3.11: End point interpolation: i due esempi precedenti di curve a cui Ŕ stata forzata l’end point interpolation.


Derivata agli estrimi: raccordo Per disegnare una curva chiusa a derivate continue si deve imporre:

 r          r
d-X-(τ0)-=  d-X-(τn),  r = 0,...,k − 2
  drt        drt

ovvero, per ogni r:

                ∑n  [r]                            ∑n  [r]
(k − 1)...(k − r)    di Ni,k− 1(τ0) = (k − 1)...(k − r)   di Ni,k−1(τn)
                i=1                                i=1

Per le proprietÓ degli Ni,k:

                k− 2                                  n
(k − 1)...(k − r) ∑  d[r]N     (τ ) = (k − 1)...(k − r)  ∑     d[r]N     (τ )
                 i=0  i  i,k− 1 0                    i=n−k+2  i  i,k−1  n

da cui derivano direttamente le condizioni da imporre, ovvero:


PIC (a) k=3
PIC (b) k=8

Figura 3.12: Derivata agli estremi: i due esempi precedenti di curve, con T modificato per avere una curva chiusa.


Algoritmi per curve BSPLINE

Disegno della curva: de Boor

L’algoritmo di de Boor calcola i punti della curva in maniera analoga a quello di de Casteljau per le curve di BÚzier.
Considerando t [tk1,tn+1]

X(t) = i=0nd i(0)N i,k(t)
= i=0nd i(0)[                                          ]
 ---t −-ti-            -ti+k-−-t--
 ti+k−1 − tiNi,k−1(t) + ti+k − ti+1 Ni+1,k−1(t)
= i=1n ⌊                                ⌋
|                                |
||   t − ti   (0)   ti+k− 1 − t (0)||
|| ----------di  +  ----------di−1||
⌈ t◟i+k−◝1◜ −-ti◞     ◟ti+k−1◝◜−-ti◞    ⌉
    α(1)(t)           1−α(1)(t)
◟----i----------◝◜-----i---------◞di(1)(t)Ni,k1(t)
a questo punto, definendo αi(j)(t) =   t − t
-------i--
ti+k−j − ti e di conseguenza di(j)(t) = α i(j)(t)d i(j1)(t)+(1α i(j)(t))d i1(j1)(t):
X(t) = i=jnd i(j)(t)N i,kj(t)
= i=k1nd i(k1)(t)N i,1(t)
ovvero, se t [tr,tr+1], X(t) = dr(k1)(t).

Per t [tr,tr+1] si possono escludere dalla sommatoria iniziale i termini nulli:

          ∑r    (k−1)
X (t) =        di   (t)Ni,1(t)
        i=r− k+1

ed Ŕ possibile ottimizzare l’algoritmo.


: Algoritmo per il disegno delle BSPLINE: deBoor
 
1function Curve = Bspline˙deboor(k,T,b,t) 
2lt = length(t) 
3[d,n] = size(b) 
4lT = length(T) 
5 
6 
7intervals=[] 
8for i=1:(lT-1) 
9      if T(i)íT(i+1) 
10         intervals=[intervals,[T(i);i]] 
11      end 
12end 
13intervals=[intervals,[T(lT);min(find(T==T(lT)))]] 
14 
15Curve = [] 
16for p=1:lt-1 
17      r = intervals(2,max(find(intervals(1,:)í=t(p)))) 
18      d=b 
19      for j=1:k-1 
20         for i=r:-1:r-k+j+1 
21             a = (t(p)-T(i))/(T(i+k-j)-T(i)) 
22          d(:,i)=a*d(:,i)+(1-a)* d(:,i-1) 
23         end 
24      end 
25      Curve=[Curve,d(:,r)] 
26end 
27endfunction

Knot Insertion

Si vuole inserire un nodo t∕∈T tale che t [t r,tr+1] per creare il vettore T = T ∪{t}, e definire una curva su T che sia identica a quella precedentemente definita su T. In particolare si definisce la curva:

  ∗     n∑+1  ∗  ∗      ∑n
X  (t) =     diNi,k(t) ≡    diNi,k(t) = X (t)
        i=0            i=0

in cui:

T : t i = (
||{ti,    i = 0,...,r
 t∗,    i = r + 1
||(
 ti−1,  i = r + 2, ...,n + μ + k
Definendo:
αL,i = (|   t∗ − ti
{ ----------,  i = r − k + 2,...,r
|( ti+k−1 − ti
  1,           i = r − k + 1
αR,i = (| ti+k−1 − t∗
{ -----------,  i = r − k + 1,...,r − 1
|( ti+k−1 − ti
  1,            i = r
si posso scrivere i nuovi Ni,k a partire da quelli vecchi:
Ni,k(t) = (
||{ Ni,k(t)                       i = 0, ...,r − k
  αL,iNi,k(t) + αR,i+1Ni+1,k(t) i = r − k + 1,...,r
||(
  Ni+1,k(t)                     i = r + 1,...,n + 1
Espletando parte dell’equazione della nuova curva si ottiene:
X(t) = i=0nd iNi,k(t) = i=0n+1d iN i,k(t)
= i=0rkd iN i,k(t) + i=rk+1rd i(α L,iNi,k(t) + αR,i+1Ni+1,k(t)) + i=r+1n+1d iN i+1,k(t)
= i=0rk+1d iNi,k(t) + i=rk+2r(d iαL,i + di1αR,i)Ni,k(t) + i=r+1n+1d i1Ni,k(t)
da cui si possono ricavare direttamente i di:
     (
     ||{ di,                 i = 0,...,r − k + 1
d∗i =   αL,idi + αR,i+1di−1 i = r − k + 2,...,r
     ||(
       di−1                i = r + 1,...,n + 1


: Algoritmo di knot insertion
 
1function [nb,nT] = Bspline˙knotins(k,T,b,p) 
2[d,n] = size(b); 
3lT = length(T); 
4 
5 
6intervals=[] 
7for i=1:(lT-1) 
8      if T(i)íT(i+1) 
9         intervals=[intervals,[T(i);i]] 
10      end 
11end 
12intervals=[intervals,[T(lT);min(find(T==T(lT)))]]; 
13 
14r = intervals(2,max(find(intervals(1,:)í=p))); 
15 
16nT = [T(1:r),p,T(r+1:lT)]; 
17 
18nb = b(:,1:r-k+1); 
19 
20for i=r-k+2:r 
21      a = (p-T(i))/(T(i+k-1)-T(i)); 
22      nb = [nb,a*b(:,i)+(1-a)*b(:,i-1)]; 
23end 
24 
25nb = [nb,b(:,r:n)]; 
26 
27 
28endfunction


PIC (a) k=3
PIC (b) k=8

Figura 3.13: Knot insertion:i due esempi precedenti di curve, con una operazione di knot insertion eseguita per ogni intervallo.


Curve NURBS

Le curve NURBS (Non Uniform Rational B-SPLINEs) sono un’estensione delle B-SPLINEs analoga alla razionalizzazione delle curve di BÚzier.
La curva NURBS Ŕ definita come:

        ∑n
           βidiNi,k(t)
X (t) = i=0n----------,  t ∈ [tk− 1,tn+1 ],  T  = {t0,...,tk− 1,...,tn+1,...,tn+k}
         ∑  β N   (t)
         i=0  i  i,k

Per comoditÓ si indica con w(t) = i=0nβ iNi,k(t) e di conseguenza con Ri,k(t) = βiNi,k(t)
  w (t); cosý l’espressione della curva diventa:

        ∑n
X (t) =    diRi,k(t)
        i=0


: Algoritmo per il calcolo delle curve NURBS
 
1function Curve = Bspline˙rational(k,T,b,beta,t) 
2[d,n] = size(b) 
3n2 = length(beta) 
4B=[beta] 
5for i=1:d //Aumenta di una dimensione 
6     B = [B;beta.*b(i,:)] 
7end 
8C = Bspline˙deboor(k,T,B,t) 
9for i=1:d 
10Curve(i,:) = C(i+1,:)./C(1,:) //Proietta 
11end 
12endfunction

ProprietÓ

Le proprietÓ delle curve NURBS derivano direttamente da quelle delle B-SPLINEs, se βi > 0  i.

Influenza dei coefficenti

Per studiare l’influenza dei coefficienti si pu˛ calcolare cosa succede quando uno di questi cambia. Si definisce:

      (
      {βj + Δ βj,  i = j
β»i =  (
       βi,          altrimenti

e di conseguenza:

        n
w»(t) = ∑  »β N  (t) = w (t) + Δ β N   (t)
        i=0  i i,k                j  j,k

La nuova curva differisce dalla nuova di:

(            )    Δ β N   (t) − X (t)Δβ  N  (t)   Δ β N   (t)(d − X  (t))
 X» (t) − X (t) =  ---j--j,k-------------j-j,k----= ----j-j,k-----j--------
                              »w(t)                        w»(t)

ovvero l’effetto Ŕ limitato dalla proprietÓ di localitÓ di Nj,k(t), ha verso (dj X(t)) (cioŔ Ŕ in direzione del nodo interessato dal cambio di peso) e ha ampiezza Δβj.


PIC (a) k=3
PIC (b) k=8

Figura 3.14: NURBS: i due esempi precedenti di curve, con il peso relativo ad un punto modificato in direzioni opposte.



PIC

Figura 3.15: Disegno di una circonferenza con una curva NURBS


Capitolo 4
Interpolazione
usando curve parametriche

Dato un’insieme di coppie di valori del tipo:
(xi,fi),   x0 < ...<  xn,  i = 0,...,n

tali che xi k e f i Ed, il problema dell’interpolazione consiste nel trovare una curva f(t) (tipicamente una funzione polinomiale: f(t) Πn(t)) ed una parametrizzazione t = g(xi) per cui la funzione interpoli i dati forniti, ovvero:

f(g(xi)) = fi,   i = 0, ...,n

e scelte in modo che siano compatibili:

(
{f : [a, b] →  Ed
(
 g : [x0,xn ] → [a, b]

Parametrizzazioni

Supponendo di vole usare la base di Lagrange Li,n(x) = j⁄=i
j=1n x− xj
x-−-x-
 i   j per definire un polinomio interpolante ai dati, si avrebbe un’espressione del tipo:

       ∑n             ∑n   ∏n
p(x) =    fiLi,n(x) =    fi   -x −-xj
       i=0            i=0   j⁄=ixi − xj
                           j=1
(4.1)

Da qui, per passare dal polinomio (4.1), strettamente legato agli xi, a una curva parametrica in t su un intervallo [a,b] arbitrario che sia interpolante gli stessi punti, si definisce:

         n
X (t) = ∑  PiLi,n(t)
        i=0

per la quale i punti di controllo Pi coincidono con i punti da interpolare fi, e i parametri ti corrispondenti dipendono dalla parametrizzazione utilizzata.
Alcuni schemi di parametrizzazione sono:

In generale questi stessi schemi di parametrizzazione si possono applicare nel dominio delle curve SPLINEs per segmentare il dominio della curva nei suoi sottodomini.


PIC

Figura 4.1: Parametrizzazioni a confronto


Interpolazione mediante SPLINE cubica C1 alla Hermite

Si supponga di avere a disposizione i dati da interpolare in forma di triple:

       ′
(ti,fi,fi),   i = 0,..., n

in cui fi rappresenta il valora da interpolare e fila tangente alla curva nel punto di interpolazione; si cerca una funzione spline f(t) Πn che interpoli i dati in questo modo:

(| f(t)   = f
|{    i      i
| f′(ti)  = fi′
|(

Ricapitolando, dati:

(ui,Pi,wi),  a = u0 < ...< un =  b

si pu˛ costruire una curva spline di terzo grado con continuitÓ C1 interpolante a partire dalla base di Bernstein. In particolare, si cerca una curva:

X(u) = Xi(u), u [ui,ui+1]
Xi(u) = k=03b k,iBk3(  u − ui )
 ---------
 ui+1 − ui
per la quale valgano:
(
||| Xi (ui)      =  Pi
||||
||||
|||| dXi-(ui)    =  wi
||{    du

|||| X  (u   )    =  P
||||   i  i+1         i+1
||||
|||| dXi (ui+1)
|( ----du----  =  wi+1

Dalle proprietÓ dei polinomi di Bernstein segue direttamente che per le condizioni agli estremi si devono imporre:

(
|| Xi(ui)    = Xi (0) = b0,i = Pi
{
||
( Xi(ui+1)  = Xi (1) = b3,i = Pi+1
mentre per le derivate:
(| dXi (ui)           1    dXi (0 )       3
|||| --------    =  ----------------=  ---------(b1,i − b0,i) = wi
|{    du          ui+1 − ui  du      ui+1 − ui
|
|||| dXi-(ui+1)     ----1----dXi-(1-)   ----3----
|(     du      =  u   −  u   du   =  u   −  u (b3,i − b2,i) = wi+1
                  i+1    i           i+1    i
da cui si ricavano direttamente:
(
||| b1,i = ui+1-−-uiwi + b0,i
|{           3
|
|||(         ui+1 −-ui
  b2,i = −     3    wi+1 + b3,i


: Spline cubiche
 
1function [Curve,hodo,hodo2]=Interp˙c1(p,w,tparam,npoints) 
2      [d,n] = size(p) 
3      Curve=[] 
4      hodo=[] 
5      hodo2=[] 
6      dp=diff(tparam) 
7      for i=1:n-1 
8         b=[] 
9         b(:,1)=p(:,i) 
10         b(:,2)=dp(i)*w(:,i)/3+b(:,1) 
11         b(:,4)=p(:,i+1) 
12         b(:,3)=-dp(i)*w(:,i+1)/3+b(:,4) 
13         tab = [0:dp(i)/npoints:1] 
14         Curve=[Curve,Bezier˙casteljau(b,tab)] 
15         hodo=[hodo,1/dp(i)*Bezier˙casteljau(Bezier˙hodograph(b),tab)] 
16         hodo2=[hodo2,1/dp(i)^2*Bezier˙casteljau(Bezier˙hodograph(Bezier˙hodograph(b)),tab)] 
17      end 
18endfunction

Schemi di approssimazione delle derivate

Se non si hanno a disposizione triplette di dati ma soltanto coppie (ui,fi), si possono definire degli schemi per approssimare i valori fied utilizzare comunque un approccio alla Hermite per l’interpolazione.

Schema FMILL Scegliendo di approssimare le derivate in questo modo:

wi = Pi+1-−-Pi−-1
  hi− 1 + hi, i = 1,...,n
w0 = P  − P
--1----0
   h0,wn = P  − P
-n-----n−1
   hn−1
in cui hi = ui+1 ui, si ha un’approssimazione al prim’ordine della curva. Infatti, se F(u) fosse la funzione che ha originato le coppie di interpolazione (ui,fi), si avrebbe:
Pi+1 = F(ui+1) = F(ui) + hiF(ui) + 1
2-hi2F′′(ξ)
Pi1 = F(ui1) = F(ui) + hi1F(ui) + 1
--
2hi12F′′(ζ)
Pi+1 Pi1 = (hi1 + hi)F(ui) + O(h2)
Pi+1 − Pi−1
------------
 hi−1 + hi = F(ui) + O(h2)

: Schema FMILL
 
1function w=I˙tang˙fmill(p,u) 
2      [d,n] = size(p) 
3      h=diff(u) 
4      w = [(p(:,2)-p(:,1))/h(1)] 
5      for i=2:n-1 
6         w = [w,(p(:,i+1)-p(:,i-1))/(h(i-1)+h(i))] 
7      end 
8      w = [w,(p(:,n)-p(:,n-1))/h(n-1)] 
9endfunction


PIC

Figura 4.2: Schema FMILL


Tangenti di Bessel L’idea alla base delle tangenti di Bessel Ŕ quella di considerare la derivata di un polinomio di secondo grado interpolante tre punti consecutivi come terzo valore per il polinomio interpolante di terzo grado.
Siano qi(u) curve polinomiali a tratti in Π2 tali che:

q (u ) = P ,  j = i − 1,i,i + 1  u ∈ [u  ,u   ]
 i  j     j                            i− 1  i+1

definite in questo modo:

                 (            )
        ∑2      2  -u-−-ui−1--
qi(u ) =    bk,iB k  u   − u
        k=0          i+1    i− 1

definendo t(u) = (            )
  -u-−-ui−1--
  ui+1 − ui− 1 la precente pu˛ essere scritta:

qi(ui) = k=03b k,iBk2(t ui) = Pi1B02(t ui) + b1,1B12(t ui) + Pi+1B22(t ui)
= Pi1(1 tui)2 + 2b 1,itui(1 tui) + Pi+1tui2 = defPi
da cui:
      Pi − Pi−1(1 − tu)2 − Pi+1t2
b1,i = ----------------i---------ui
              2tui(1 − tui)

Mantenendo la notazione hi = ui+1 ui, si ha che:

t  = ---hi−-1--,  1 − t  =  ---hi----
ui   hi−1 + hi        ui   hi−1 + hi

e si pu˛ riscrivere l’equazione precedente in questa forma:

        (                                   )
      1- Pi(hi−1-+-hi)2   Pi−1hi-  Pi+1hi−-1
b1,i = 2      h   h      −  h     −     h
              i−1 i         i−1         i

Ora che i polinomi qi sono definiti Ŕ possibile ricavare i valori wi per i punti di interpolazione in questo modo:

wi = dqi(ui)-
  du = ---1-----
h   + h
 i− 1    idqi(ti)
  dt = h Pi-−-Pi−-1+ h    Pi+1 −-Pi-
--i--hi−1-------i−-1---hi----
          h   +  h
           i−1    i
Agli estremi si fissano:
      2(P1-−-P0)-       2-(Pn--−-Pn−1)
w0 =   h  + h      wn =   h  + h
         0    1            n     n− 1

Questo schema fornisce un’approssimazione della derivata della curva, infatti:

Pi = F(ui)
Pi1 = Pi F(ui)hi1 + 1
--
2hi12F′′(u i) + O(hi13)
Pi+1 = Pi + F(ui)hi + 1
--
2hi2F′′(u i) + O(hi3)
da cui:
Pi+1 −-Pi-
   h
    i = F(ui) + 1-
2hiF′′(ui) + O(hi2)
Pi −-Pi−-1
  hi−1 = F(ui) 1-
2hi1F′′(ui) + O(hi12)
ed infine:
w  = F ′(u ) + O (h2)
  i       i


: Schema di Bessel
 
1function w=I˙tang˙bessel(p,u) 
2      [d,n] = size(p) 
3      h=diff(u) 
4      w = [2*(p(:,2)-p(:,1))/h(1)+h(2)] 
5      for i=2:n-1 
6         w = [w,((p(:,i)-p(:,i-1))*h(i)/h(i-1)+(p(:,i+1)-p(:,i))*h(i-1)/h(i))/(h(i-1)+h(i))] 
7      end 
8      w = [w,2*(p(:,n)-p(:,n-1))/h(n-1)+h(n-2)] 
9endfunction


PIC

Figura 4.3: Schema di Bessel


Interpolazione mediante SPLINE cubica C2

Dati:

a =  u <  ...<  u  = b,  P  ∈ Ed,   i = 0,...,n
      0         n         i

si vuole definire una spline C2 composta dalle curve:

                                  (          )
          3∑                          u − ui
Xi (u )‡ =    bk,iB3k(t(u)),  t(u) =   ---------  †
         k=0                        ui+1 − ui

per brevitÓ si pone ui+1 ui = hi per le quali valgono le condizioni:

Xi(ui) = Xi1(ui) = Pi
dXi-(ui)-
  du = dXi−-1(ui)-
   du
d2Xi(ui)-
  d2u = d2Xi−-1(ui)-
   d2u
Dalle condizioni per le continuitÓ C0 e C1 deriva direttamente che i cofficienti devono valere, come ricavato in precedenza:
b0,i = Pi
b1,i = Pi + hi
--
3wi
b2,i = Pi+1 hi
3wi+1
b3,i = Pi+1

Scelta delle tangenti wi

Dalle condizioni di continuitÓ della derivata seconda d2Xi(ui)-
 d2u = d2Xi−1(ui)
   d2u, segue che:

    2              2
-1-d-Xi(0)-=  -1--d-Xi−-1(1)
h2i   d2t      h2i− 1   d2t

ovvero:

1--
h2i(b  − 2b   + b  )
 2,i    1,i    0,i = -1---
h2i−1(b     − 2b     + b    )
  3,i−1     2,i−1    1,i−1
da cui:
 1
-2-
hi[(        hi    )     (     hi   )       ]
  Pi+1 −  --wi+1  − 2  Pi + -- wi  + (Pi)
          3                  3
=  1
-2---
hi−1[        (            )   (                ) ]
               hi− 1               hi−1
 (Pi) − 2 Pi − -----wi  +   Pi−1 + ----wi−1
                 3                  3
la quale pu˛ essere manipolata, moltiplicando per 3hihi1, fino a questa forma:
                                       [    P    − P      P  − P    ]
hiwi− 1 + 2 (hi− 1 + hi)wi + hi−1wi+1 = 3 hi−1-i+1-----i+ hi--i----i−1-
                                     ◟----------hi--◝◜-------hi−1---◞
                                                    Ti

I Tn sono tutti termini noti, e si ha un sistema di n1 considerando n intervalli equazioni in n + 1 tutti i wi incognite, che pu˛ essere scritto nella forma:

Aw  = T
(4.2)

in cui:

w  = (w  ,...,w )T,   T = (T  ,...,T    )T
        0     n             1     n−1

ed A Ŕ una matrice in M(n1)×(n1) ad elementi:

     (                                                        )
     | h1  2(h0 + h1)     h0        0          ⋅⋅⋅         0  |
     || 0       h2      2(h1 + h2)   h1         ⋅⋅⋅         0  ||
A =  || ...      ...         ...      ...         ...         ...  ||
     ||                                                        ||
     ( 0      ⋅⋅⋅          0      hn −1  2(hn−2 + hn− 1)  hn− 2)

Rimangono da fissare i valori delle derivate agli estrimi w0 e wn. Alcune tra le scelte pi¨ comuni sono le seguenti:


PIC

Figura 4.10: Spline C2 naturale 3d



PIC

Figura 4.11: Spline C2 not a knot 3d



PIC

Figura 4.12: Spline C2 periodica 3d


Considerazioni

Sia s(u) una spline cubica C2 tale che:

s(ui) = fi,  i = 0,...,n

e y(u) una funzione qualsiasi C2 a sua volta tale che:

y(ui) = fi,  i = 0,...,n

La quantitÓ:

∫ un
    Ęs2(u)du
 u0

Ŕ una possibile misura del grado di oscillazione della curva.
Si dimostra che vale la relazione:

∫ un           ∫ un
     Ęs2(u)du ≤      Ęy2(u )du
  u0             u0

Infatti:

ř(u) = (ř(u) Ęs (u)) + Ęs (u)
ř2(u) = Ęs 2(u) + (ř(u) Ęs (u))2 + 2Ęs (u)(ř(u) Ęs (u))
u0un ř2(u)du = u0un Ęs 2(u)du + u0un (ř(u) Ęs (u))2du + 2 u0un Ęs (u)(ř(u) Ęs (u))du
integrando per parti l’ultimo integrale della sommatoria si ha:
∫                                                ∫
  unĘs(u )(Ęy(u ) − Ęs(u))du = Ęs(u)(˙y(u) − ˙s(u ))|un −   un.s..(u)(˙y(u) − ˙s(u))du
 u0                        ◟--------◝◜-------u0◞    u0
                                    α            ◟----------◝β◜----------◞

in cui α = 0 ed in quanto ... s(u) Ŕ costante a tratti, β pu˛ essere trasformato nella somma:

n−1
∑  si[y(u) − s(u)]ui+1
i=0               ui

che a sua volta Ŕ nulla perchŔ valgono le condizioni y(ui) = s(ui).
Un’altra misura della stessa quantitÓ pu˛ essere ottenuta dalla curvatura, considerando come parametro la lunghezza dell’arco:

∫ L          ∫ un 2                    ∫ un
   k(s)ds =      k ∥X˙(u)∥du,   s(t) =     ∥X˙(t)∥dt
 0            u0                        u0

Rappresentando la funzione vettorialmente:

         (     )            (     )            (     )
            u                  1                  0
X  (u ) =  y (u ) ,   X (u) =   ˙y(u)  ,  X (u) =   Ęy(u)

la curvatura assume espressione:

       X˙(u)-×-XĘ(u)    ---yĘ(u)----
k(u) =    ∥X˙(u)∥3   =  (1 + y˙2(u )) 32

da cui:

∫                            ∫
  un k2(u)∥X˙(u)∥du,   s(t) =   un ---Ęy(u)----du
 u0                           u0  (1 + ˙y2(u))

Spline cubiche C1 G2

Data una spline in forma di BÚzier:

X(u) = Xi(u), u [ui,ui+1]
Xi(t) = i=03b k,iBk3(t), t = ( u − ui)
  ------
    hi
che ha coefficienti:
(
|||| b0,i  = Pi
||||             hi
|||{ b1,i  = Pi +  3 wi
                 hi
||| b2,i  = Pi+1 −  3-wi+1
|||| b    = P
||||  3,i      i+1
(

Imponendo le condizioni di raccordo, considerando ω1,i = 1 si ottiene:

 2          2
d-Xi(u ) − d-Xi-−1(u ) = νw  ,  w  = dXi-(u ),  ν =  ω
 du    i     du     i     i i    i    du   i     i    2,i

Le condizioni di raccordo permettono di scrivere infine le equazioni:

6--
h2i[                                   ]
 (P   −  hiw   ) − 2(P +  hiw ) + P
   i+1   3  i+1       i   3   i    i
-6---
h2i+1[                                        ]
 Pi − 2(Pi − hi−-1wi) + (Pi− 1 + hi−-1wi−1)
               3                 3νiwi = 0
manipolando e moltiplicando per hihi−-1
   2 si ottiene infine:
       [                      ]                  [                            ]
                      νi                               Pi+1 −-Pi-    Pi-−-Pi−1-
hiwi +  2(hi−1 + hi) + 2 hihi−1 wi + hi− 1wi+1 = 3  hi− 1   hi     + hi  hi− 1

ovvero si ottiene un sistema analogo a quello delle spline C2 con un termine in pi¨ sulla diagonale:

                       ( d         )
                       |  0  .     |
(A + D )w =  T,   D =  |(      ..   |)
                                 dn

in cui d0,dn sono assegnati e di = νi
2hi1hi.
Per ottenre una spline chiusa, supponendo di avere il primo punto uguale al primo (P0 = Pn) si ha invece un sistema in n incognite in cui le prime e le ultime equazioni di A sono analoghe a quelle per la spline periodica precedentemente ricavate, ovvero del tipo:

                                       (                            )
h0wn −1 + 2(hn−1 + h0)w0 + hn− 1w1 = 3  hn− 1P1-−-P0-+ h0 P0-−-Pn−-1
                                               h0           hn −1

e inoltre

d0 = dn =  ν0hn−1h0
           2

╚ interessante osservare cosa succede per νmin →∞. Innanzitutto serve calcolare:

(A + D)w = T
(AD1 + I)(Dw) = T
(I + E)w˜ = T, E = AD, ˜w = Dw
Se νmin →∞ anche D →∞; AD1 Ŕ non singolare e per D sufficientemente grande E< 1. In quanto:
         −1    ----1---
∥(I + E )  ∥ ≤ 1 − ∥E ∥

si ha che:

                       1
(I + E )˜w ⇒ ∥ ˜w∥ ≤  1 −-∥E∥-∥T ∥

ovvero se D →∞, w 0.


: ν-spline
 
1function [w,A,T]=I˙tang˙nu˙nat(p,u,nu) 
2      [d,n] = size(p) 
3      h=diff(u) 
4      A = [] 
5      T = [] 
6      A(1,:) = [2*h(1),h(1),zeros(1,n-2)] 
7      T(:,1) = 3*(p(:,2)-p(:,1)) 
8      for i=2:n-1 
9      A(i,:)=[zeros(1,i-2),h(i),2*(h(i)+h(i-1)),h(i-1),zeros(1,n-i-1)] 
10      T(:,i) = 3*((h(i-1)/h(i))*(p(:,i+1)-p(:,i))+(h(i)/h(i-1))*(p(:,i)-p(:,i-1))) 
11      end 
12      A(n,:) = [zeros(1,n-2),h(n-1),2*h(n-1)] 
13      T(:,n) = 3*(p(:,n)-p(:,n-1)) 
14      A(1,1)=A(1,1)+nu(1) 
15      for i=2:n-1 
16      A(i,i)=A(i,i)+(nu(i)*h(i-1)*h(i)/2) 
17      end 
18      A(n,n)=A(1,1)+nu(n) 
19      w=(AT’) 
20endfunction


PIC
Figura 4.13: ν-spline naturale



: ν-spline chiusa
 
1function [w,A,T]=I˙tang˙nu˙closed(p,u,nu) 
2      [d,n] = size(p) 
3      h=diff(u) 
4      A = [] 
5      T = [] 
6      A(1,:) = [2*(h(1)+h(n-1)),h(n-1),zeros(1,n-4),h(1),0] 
7      T(:,1) = 3*(h(n-1)*((p(:,2)-p(:,1))/h(1))+h(1)*((p(:,n)-p(:,n-1))/h(n-1))) 
8      for i=2:n-1 
9      A(i,:)=[zeros(1,i-2),h(i),2*(h(i)+h(i-1)),h(i-1),zeros(1,n-i-1)] 
10      T(:,i) = 3*((h(i-1)/h(i))*(p(:,i+1)-p(:,i))+(h(i)/h(i-1))*(p(:,i)-p(:,i-1))) 
11      end 
12      A(n,:) = [0,h(n-1),zeros(1,n-4),h(1),2*(h(n-1)+h(1))] 
13      T(:,n)=T(:,1) 
14      A(1,1)=A(1,1)+(nu(1)/2)*h(1)*h(n-1) 
15      for i=2:n-1 
16      A(i,i)=A(i,i)+(nu(i)*h(i-1)*h(i)/2) 
17      end 
18      A(n,n)=A(1,1)+(nu(1)/2)*h(1)*h(n-1) 
19      w=(AT’) 
20endfunction


PIC
Figura 4.14: ν-spline chiusa


Capitolo 5
Superfici parametriche

Rappresentazione parametrica

Una superficie, in maniera analoga ad una curva, pu˛ essere rappresentata in modo parametrico come l’immagine di una funzione continua e iniettiva di due variabili reali X : (A 2) E3 tale che:

           (       )
             x(u,v)
X (u, v) = |( y(u,v)|) ,  (u, v) ∈ ℝ2

             z(u,v)

Definizione 5.1 (Superficie parametrica regolare).

La superficie parametrica si dice regolare se valgono le condizioni:

Definizione 5.2 (Vettori associati alla superficie: tangenti e normale). Per comoditÓ nello studio delle proprietÓ della superficie si definiscono i vettori:

Xu = ( ∂x (u,v))
| --------|
||    ∂u   ||
||         ||
|| ∂y-(u,v)||
||    ∂u   ||
||         ||
( ∂z-(u,v))
     ∂u, Xv = ( ∂x(u, v))
| --------|
||   ∂v    ||
||         ||
|| ∂y(u,v-)||
||   ∂v    ||
||         ||
( ∂z(u,v-))
    ∂v
delle derivate parziali rispetto ai parametri, ed il vettore normale alla superficie di conseguenza:
          -Xu-∧-Xv--
N (u,v) = |Xu ∧ Xv |

Osservazione 5.1 (Caratteristiche dei vettori Xu,Xv,N(u,v)). Se la superficie Ŕ regolare si ha che i vettori Xu e Xv sono linearmente indipendenti e non paralleli, ovvero Xu Xv0; di conseguenza anche N(u,v)0 (u,v) A.

Esempi di superfici parametriche

Sfera La sfera centrata nell’origine degli assi Ŕ definibile come il luogo geometrico dei punti equidistanti dall’origine:

 2    2    2
x +  y +  z =  r

Una possibile parametrizzazione per i punti della superficie sferica si pu˛ dare a partire dai parametri u [π
2,π
2] e v [0, 2π]:

          (          )
            cosu cosv
S (u, v) = |( cosu sin v|)

              sin u

I vettori tangenti alla superficie sono quindi:

      (            )          (            )
        − sin u cosv             − cosu sin v
Su =  |( − sin u sin v|) ,  Sv =  |(  cosu cosv |)
           cosu                      0


PIC (a) PIC (b)
PIC (c) PIC (d)

Figura 5.1: Esempio di sfera parametrizzata come descritto nel testo. Nella prima e seconda figura si evidenziano in risoluzione le coordinate verticali od orizzontali, in basso l’unione delle due figure con e senza linee ad unire i punti nelle due direzioni.


Cono Un cono Ŕ lo spazio geometrico dei punti che soddisfano l’equazione:

 2    2    2
z =  x +  y

Una parametrizzazione in due parametri u [L,L], v [0, 2π] Ŕ:

          ( cu cosv)
          |        |
C (u,v) = ( cu sin v)
               cu

I vettori tangenti sono, di conseguenza:

      (       )          (         )
        ccos v             − cu sinv
Cu =  |( csinv |) ,  Cv =  |( cu cosv |)
          c                   0


PIC

Figura 5.2: Esempio di cono parametrizzato come descritto nel testo.


Riparametrizzazione ammissibile

Si definiscono due nuovi parametri in funzione di quelli vecchi:

({
  ˜u = ˜u (u, v)
( ˜v = ˜v(u,v )

Data la matrice:

          ( ∂˜u   ∂˜u )
          | ---  ---|
J(u,v ) = || ∂u   ∂v ||
          |( ∂˜v   ∂˜v |)
            ---  ---
            ∂u   ∂v

la riparametrizzazione si definisce ammissibile se det J(u,v)0. In particolare si ha che:

Xu Xv = (    ∂˜u      ∂ ˜v)
  X ˜u---+ X ˜v---
     ∂u      ∂u(   ∂ ˜u      ∂˜v )
 X ˜u--- + X ˜v---
    ∂v       ∂v = (X ˜u ∧ X ˜v) ( ∂˜u ∂˜u    ∂˜v ∂˜v)
  ------−  ------
◟-∂u-∂v-◝◜-∂u-∂v-◞det J(u,v)

Tipologie di superfici

Superfici di rivoluzione

Una superficie di rivoluzione Ŕ una superficie ottenuta ruotando una curva (chiamata curva generatrice o profilo) intorno ad un asse (detto asse di rotazione). La circonferenza ottenuta interesecando un piano perpendicolare all’asse di rotazione con la superficie Ŕ detta parallelo, la curva ottenuta intersecando la superficie con un piano passate per l’asse di rotazione Ŕ detta meridiano.
In generale una superficie di rotazione Ŕ rappresentabile in modo parametrico. Sia una curva C definita parametricamente come:

        (
        { C (u)
C (u ) =    x      ,  u ∈ [a,b]
        ( Cy(u)

La curva C pu˛ essere fatta ruotare intorno all’asse z definendo la superficie:

          (
          ||{ Cx (u )cos(v)
X (u,v) =   Cy (u )sin (v)   ,  v ∈ [0,2π]
          ||(
            Cx (u )


PIC

Figura 5.3: Esempio di curva di BÚzier utilizzata come curva parametrica di rotazione.


Superfici Rigate

Una superficie rigata Ŕ una superficie ottenuta da un’unione di rette.

Definizione 5.3 (Superficie rigata). Una superficie S si dice rigata se esiste una famiglia di rette rα tale che S Ŕ l’unione delle rette di tale famiglia: S = αrα. Una definizione equivalente Ŕ che S Ŕ rigata se per ogni punto di S passa una retta rs che sia tutta contenuta in essa.

Date due curve C0(u),C1(u) e due parametri u [a,b],v [0, 1] una superficie rigata pu˛ essere difinita parametricamente come:

S (u, v) = (1 − v )C0(u) + vC1(u )

ovvero si sfrutta l’esistenza di una corrispondenza univoca tra i punti di C0 e quelli di C1.

Superficie bilineare

Se C0 e C1 sono segmenti, ovvero sono parametrizzabili come:

(
{ C0(u) = (1 − u)p00 + up01
(
  C1(u) = (1 − u)p10 + up11

la superficie rigata S pu˛ essere scritta come:

                                            (        ) (      )
                                             p00  p01    1 − v
S(u,v ) = (1 − v)C0 (u) + vC1(u) = (1 − u,u ) p   p        v
                                               10   11


PIC

Figura 5.4: Esempio di curva di BÚzier utilizzata come bordi per una superficie rigata.


Superfici tensor-product

Definizione 5.4 (Tensor product). Date due funzioni monovariate S1,S2 tali che:

S1 =< f0,...,fn >, f0,...,fn funzioni monovariate in u [a,b]
S2 =< g0,...,gm >, g0,...,gm funzioni monovariate in v [c,d]
si definisce tensor product S1 S2 la funzione:
S1 ⊗ S2 =<  f0g0,...,f0gm,...,fng0,...,fngm  >

S1 S2 Ŕ pertanto una funzione bivariata, e vale:

                        ∑n  m∑
S (u, v) ∈ S ⇔  S(u,v) =        bijfi(u)gj(v),  (u,v) ∈ [a, b] × [c,d]
                        i=0j=0

Il tensor product pu˛ quindi essere utilizzato per definire superfici, considerando funzioni del tipo:

           n  m
          ∑   ∑
X (u,v) =        bijfi(u)gj(v),  (u,v) ∈ [a, b] × [c,d]
          i=0j=0

Superfici tensor-product come patch BÚzier

La patch di BÚzier Ŕ definita a partire dalla base omonima nel modo seguente:

           n  m
X (u,v ) = ∑  ∑  b Bn (u)Bm (v)
           i=0 j=0  ij  i     j

ProprietÓ

Curve di contorno La superficie Ŕ delimitata da quattro curve di contorno boundary :

(                            ∑m        m
|||| u = 0,v ∈ [0,1],  X (0,v) = ∑ j=0 b0jB j (v )
|||| u = 1,v ∈ [0,1],  X (1,v) =   mj=0 b1jBmj (v )
{                            ∑n        n
|| v = 0,u ∈ [0,1],  X (u,0 ) = ∑ i=0 bi0B i (u)
|||| v = 1,u ∈ [0,1],  X (u,1 ) =  ni=0 bi1Bni (u)
||(

Posizione negli angoli Negli angoli la supericie Ŕ parallela al piano individuato dai punti di controllo. Infatti, considerando le derivate prime parziali della superficie lungo le due direzioni negli angoli, si ha:

u(0, 0) = nΔ(1,0)b 0,0 = n(b1,0 b0,0)
u(1, 0) = n i=0n1Δ(1,0)b i,0Bin1(1) = nΔ(1,0)b n1,0 = n(bn,0 bn1,0)
v(0, 0) = nΔ(0,1)b 0,0 = m(b0,1 b0,0)
v(0, 1) = m j=0m1Δ(0,1)b 0,jBjm1(1) = mΔ(0,1)b 0,m1 = m(b0,m b0,m1)

Altre proprietÓ In quanto somme di polinomi di Bernstein valgono le seguenti proprietÓ, dimostrabili in modo analogo a quanto giÓ mostrato in modo esplicito per le curve di BÚzier:

Forma compatta La superficie X(u,v) pu˛ essere espressa in forma compatta come:

                             (               ) (       )
                               b0,0  ⋅⋅⋅  b0,m     Bm0 (v)
X (u, v) = (Bn (u),...,Bn(u )) ||  ..         ..  || ||    ..  ||
             0         n     (  .         .  ) (    .  )
                               bn,0  ⋅⋅⋅  bn,m     Bmm (v)

Algoritmi

de Casteljau

L’algoritmo di de Casteljau per le superfici funziona in maniera del tutto analoga a quello descritto precedentemente per le curve.

Dati bij00 = b ij come punti iniziali per l’algoritmo (i = 0,..,n; j = 0,...,m), si definisce il passo dell’algoritmo unidirezionale (nelle due direzioni, rispettivamente) in questo modo:

bijk0(u) = (1 u)b ijk1,0 + ub i+1,jk1,0, i = 0,...,n k
bij0k(v) = (1 v)b ij0,k1 + vb i,j+10,k1, j = 0,...,m k

Svolgendo contemporaneamente nelle due direzioni si ha:

bijkl(u,v) = (1 v)b ijk,l1 + vb i,j+1k,l1
= (1 v)[       k−1,l−1     k−1,l−1]
(1 − u)bij     + ubi+1,j + v[       k− 1,l− 1    k− 1,l− 1]
 (1 − u )bi,j+1   + ubi+1,j+1
e se m = n,
bijkl(u,v) = ((1 u),u)( bk−1,l−1  bk−1,l−1)
|  ij        i+1,j  |
(  k−1,l−1   k−1,l−1)
  bi,j+1    bi+1,j+1( 1 − v)
|      |
(      )
    v, i = 0,...,n k; j = 0,...,m l

: Algoritmo per il calcolo delle patch di BÚzeir
 
1function Curve = Bezier˙patch˙casteljau(b,u,v) //b poligono, u,v punti di valutazione 
2Lu = length(u) 
3Lv = length(v) 
4[d,n,m] = size(b) 
5Curve =[] 
6tcurve=[] 
7for im=1:m 
8      tcurve(:,:,im) = Bezier˙casteljau(b(:,:,im),v); 
9end 
10for iv=1:Lv 
11      Curve(:,:,iv)=Bezier˙casteljau(matrix(tcurve(:,iv,:),3,-1),u) 
12end 
13endfunction


PIC
PIC

Figura 5.5: Esempio di Patch di BÚzier da due angolazioni di verse.



PIC
PIC

Figura 5.6: Esempio di Patch di BÚzier con soli quattro punti di controllo: risalta la proprietÓ di tangenza agli angoli.


Degree-elevation

L’algoritmo di degree elevation consente di elevare il grado della patch di BÚzier lungo una delle due direzioni aggiungendo una riga di punti di controllo e modificando quelli esistenti in modo opportuno. In particolare, data la superficie:

             ∑n ∑m     n     m
Xn,m (u,v) =       bijBi (u )B j (v)
             i=0 j=0

si vuole aumentarne il grado lungo una direzione, ovvero ottenere:

Xn+1,m(u,v) = i=0n+1 j=0mc ijBin+1(u)B jm(v)
oppure
Xn,m+1(u,v) = i=0n j=0m+1d ijBin(u)B jm+1(v)
Per le proprietÓ dei polinomi di Bernstein segue direttamente che:
Xn,m(u,v) = i=0n j=0mb ijBin(u)B jm(v)
= Xn+1,m(u,v) = i=0n+1 j=0mc ijBin+1(u)B jm(v) = j=0m(n+1           )
 ∑  c  Bn+1 (u)
 i=0 ij  iBjm(v)
da cui:
 m ( n          )           m (n+1            )
∑    ∑      n       m      ∑    ∑      n+1       m
        bijB i (u) B j (v) =        cijBi   (u ) B j (v)
j=0  i=0                    j=0  i=0

che implica:

 n            n+1
∑      n      ∑       n+1
   bijB i (u) =    cijB i  (u),  ∀j
i=0           i=0

I nuovi punti di controllo sono pertanto:

cij = n +-1 −-ibij + --i--bi−1,j,   ∀j, i = 0,...,n + 1
       n +  1       n + 1

Gli stessi ragionamenti si applicano in modo perfettamente analogo anche all’altro caso.

In generale si preferisce aumentare il grado contemporaneamente in entrambe le direzioni. Applicando lo stesso ragionamento anche nell’altra direzione si ottiene:

   (             )               (                )
n∑+1  ∑m                       n∑+1  m∑+1
   (    cijBmj (v)) Bni+1 (u) =    (     dijBm+j1  (v)) Bni+1 (u)
i=0  j=0                      i=0  j=0

da cui, come prima:

∑m     m       m∑+1     m+1
   cijBj (v) =     dijB j   (v )
j=0            j=0

che implica:

      m-+-1-−-j-     --j---
dij =   m + 1   cij + m +  1ci,j−1,  ∀i, j = 0,...,m +  1

Esplicitando i cij si ottiene:

dij = m  + 1 − j
----------
  m  + 1(n + 1 − i        i       )
 ---------bij + -----bi−1,j
   n + 1        n + 1 +    j
------
m  + 1(n + 1 − i          i          )
 ---------bi,j−1 + ------bi−1,j− 1
   n + 1          n + 1
= (  i    n + 1 − i)
 n-+-1-,--n +-1--(              )
| bi−1,j− 1  bi−1,j|
(              )
   bi,j− 1    bij(     j     )
    ------
||   m + 1   ||
||           ||
( m-+-1-−-j-)
    m + 1

: Degree elevation nelle due direzioni
 
1function B = Bezier˙patch˙degelev(b) 
2[d,n,m] = size(b) 
3B=b 
4 
5for i=1:n 
6 B(:,i,1:m+1)=matrix(Bezier˙degelev(matrix(b(:,i,:),3,-1)),3,1,m+1) 
7end 
8b=B 
9for i=1:m+1 
10 B(:,1:n+1,i)=Bezier˙degelev(b(:,:,i)) 
11end 
12 
13endfunction


PIC
PIC

Figura 5.7: Esempio di applicazione della degree elevation, da due angolazioni di verse. In rosso il poligono di controllo originale, in verde quello di grado elevato, in grigio la sovrapposizione esatta della vecchia superficie (in blu) e della nuova (grigia).



PIC
PIC

Figura 5.8: Esempio di applicazione della degree elevation, da due angolazioni di verse. In rosso il poligono di controllo originale, in verde quello di grado elevato, in ciano una seconda elevazione. In grigio la sovrapposizione esatta della vecchia superficie (in blu) e delle nuova (grigia).


Subdivision

L’algoritmo di subdivision lungo una delle due direzioni Ŕ analogo a quello per le curve. Data la superficie:

          ∑n  m∑
X (u,v ) =       bijBni (u)Bmj (v)
           i=0 j=0

la si vuole suddividere in due superfici lungo una delle due direzioni:

(
|||           ∑n    ∑m       n     m                     ^       t
{XL (τ,v ) =  i=0   j=0 cijB i (τ)B j (v) = X (t,v), t ∈ [0,t], τ = ^t ∈ [0,1 ]
||            ∑n   ∑m                                           t − ^t
|(XR  (τ, v) =   i=0   j=0 dijBni (τ )Bmj (v) = X (t,v ), t ∈ [^t,1], τ =----^ ∈ [0,1]
                                                               1 − t

oppure lungo l’altra:

(
||            ∑n    ∑m       n     m                             t
|{ XL (u,τ ) =  i=0   j=0 cijB i (u)B j (τ ) = X (u, t), t ∈ [0,^t], τ = ^ ∈ [0, 1]
|             ∑    ∑                                            tt − ^t
||( XR (u,τ ) =   ni=0   mj=0 dijBni (u )Bmj (τ) = X (u,t), t ∈ [^t,1 ], τ =---- ∈ [0,1]
                                                                1 − ^t

In entrambi i casi, in modo analogo al caso bidimensionale i coefficienti si ottengono a partire da quelli dell’algoritmo di de Casteljau:

       k,0             n−k,0
ckj = b0j (τ), dkj = bk,j  (τ )

per la direzione u e:

cik = b0i,0k(τ),  dik = b0,i,mk− k(τ )

lungo v.


: Subdivision nelle due direzioni
 
1function [B1,B2,B3,B4] = Bezier˙patch˙subdiv(B,t) 
2[d,n,m] = size(B) 
3 
4for i=1:n 
5      [a,b] = Bezier˙subdiv(matrix(B(:,i,:),3,-1),t) 
6      B5(1:3,i,1:m)=matrix(a,3,1,m) 
7      B6(1:3,i,1:m)=matrix(b,3,1,m) 
8end 
9 
10for i=1:m 
11      [a,b] = Bezier˙subdiv(B5(:,:,i),t) 
12      B1(1:3,1:n,i)=a 
13      B2(1:3,1:n,i)=b 
14      [a,b] = Bezier˙subdiv(B6(:,:,i),t) 
15      B3(1:3,1:n,i)=a 
16      B4(1:3,1:n,i)=b 
17end 
18endfunction


PIC
PIC

Figura 5.9: Esempio di subdivision applicata a t=0.5, vista da due angolazioni diverse. In magenta i punti di controllo originali; i nuovi punti di controllo sono colorati e le relative superfici hanno colori concordi.



PIC

Figura 5.10: Esempio di subdivision applicata a t=0.5, vista dall’alto si nota bene come ogni quadrante sia diviso in quattro quadranti uguali.



PIC
PIC

Figura 5.11: Esempio di subdivision applicata a t=0.2, vista da due angolazioni diverse. In magenta i punti di controllo originali, in grigio la superficie originale; i nuovi punti di controllo sono colorati e le relative superfici hanno colori concordi.



PIC

Figura 5.12: Esempio di subdivision applicata a t=0.2, vista dall’alto si nota la distribuzione della suddivisione; le griglie delle nuove superfici, campionate uniformemente come in tutte gli esempi precedenti, enfatizzano la diversa scalatura delle quattro suddivisioni.


Interpolazione di superfici con patch di BÚzier

Una superficie tensor-product di BÚzier, come definita in precedenza, Ŕ una funzione del tipo:

          ∑n ∑m     n     m
X (u,v) =       bijBi (u )Bj (v),  (u,v) ∈ [0,1] × [0,1]
          i=0j=0

In questo caso, definiti:

S1 = < Bn0,...,Bnn >,  S2 = < Bm0 ,...,Bmm  >,   S = S1 ⊗ S2

si ha che si sta lavorando su uno spazio di dimensione:

N =  dim S = (n + 1)(m  + 1)

Il problema dell’interpolazione di superfici quindi Ŕ completamente determinanto se si hanno a disposizione N campionamenti della superficie da interpolare; i dati del problema sono rappresentabili come:

Xin (ui,vj) = Pk,   k = 0,...,N  − 1,  Pk  ∈ E3

Interpolazione su punti sparsi

Se si hanno a disposizione campionamenti sparsi della superficie che si vuole interpolare:

Xin(ui,vi) = Pk,  i = 0,...,N ;

si va a risolvere il sistema:

           ∑n ∑m      n     m
X (ui,vi) =        bijB i (ui)Bj (vi) = Xin (ui,vi) = Pk
           i=0j=0

il quale pu˛ essere scritto nella forma:

                                       (     )
                                          b00
(                                    ) ||   .. ||    (   )
  Bn0(u1)Bm0 (v1)  ⋅⋅⋅  Bnn(u1)Bmm(v1)  ||   . ||     P1
||        ..                    ..      || || b0m || =  || .. ||
(        .                    .      ) ||  b10 ||    ( . )
  Bn0 (un)Bm0 (vn)  ⋅⋅⋅  Bnn(un)Bmm(vn)  ||   .. ||     Pk
                                       (   . )
                                         bnm

Interpolazione su griglia rettangolare

Se i campionamenti sono distribuiti su una griglia rettangolare:

X (ui,vi) = Pij,  i = 0,...,n; j = 0,...,m

si vuole risolvere, come prima, il sistema:

            ∑n ∑m
X (ui,vj) =       bijBni (ui)Bmj (vj) = Xin (ui,vj) = Pij
            i=0j=0

Per ogni coppia di punti si ha:

                                 ( b00  ⋅⋅⋅  b0m ) ( Bm (vr))
            (  n          n    ) ||  .         . || ||   0.   ||
X  (uk, vr) =  B0 (uk ),...,B n(uk)  (  ..         .. ) (    ..   ) =  Pkr
                                   bn0  ⋅⋅⋅  bnm     Bmm(vr)
                                 ◟------◝◜-------◞
                                         B

Definita la matrice dei punti noti:

     (               )
       P00  ⋅⋅⋅  P0m
P  = ||  ..          .. ||
     (  .          . )
       Pn0  ⋅⋅⋅  Pnm

e le matrici:

     (                      )         (                      )
       Bn0(u0)  ⋅⋅⋅  Bn0 (un )            Bm0 (v0)  ⋅⋅⋅  Bm0 (vm)
F  = ||    ..             ..   || ,  G =  ||    ..             ..   ||
     (   n.           n .   )         (   m.           m .   )
       B n(u0)  ⋅⋅⋅  Bn (un )            Bm (v0)  ⋅⋅⋅  Bm (vm)

il problema pu˛ essere rappresentato in forma matriciale come

  T
F  BG  =  P

Patch triangolari di BÚzeir

L’idea alla base delle patch triangolari Ŕ quella di usare una combinazioni baricentrica di coordinate come parametri per muoversi su una superficie. In particolare in questo caso si considera come dominio un triangolo non degenere T,di vertici ABC; i punti dell’area sono rappresentabili come combinazione baricentrica dei vertici:

P ∈ T,   P  = uA +  vB + wC,    u + v + w =  1

e la superficie come una funzione:

X  : T → E3

Coordinate baricentriche

Le coordinate baricentriche descritte sopra sono in relazione biunivoca con quelle cartesiane. Si ha:

(  )    (             ) (  )
  x       xA   xB  xC     u
|( y|)  = |( yA   yB  yC |) |( v|)
  1        1   1    1     w
        ◟------◝◜------◞
               T

AffinchŔ la relazione sia biunivoca (ovvero il sistema ha soluzione unica) se e solo se la matrice T Ŕ non singolare, ovvero il triangolo ABC non Ŕ degenere.

Polinomi di Bernstein generalizzati

I polinomi di Bernstein possono essere generalizzati alle coordinate baricentriche. Per comoditÓ si utilizza la seguente notazione:

Definizione 5.5 (polinomio di Bernstein generalizzato). Si definisce polinomio di Bernstein generalizzato il polinomio:

  n      n!  I
B I(μ) = I!μ

Tale polinomio Ŕ bivariato e di grado n, infatti μI = uivjwk e i + j + k = n.
Fissare |I| = n riduce il numero di indici distinti a (n + 1)(n + 2)2.

ProprietÓ

Non negativitÓ Tutti i fattori sono positivi o nulli, e la propritÓ Ŕ banalmente soddisfatta.

Partizione dell’unitÓ Vale che:

 ∑
     Bn (μ ) ≡ 1,  ∀μ ∈  T
|I|=n  I

Dimostrazione. si ha infatti che:

1 = 1n = (u + v + w)n = (u + (v + w))n
= i=0n(  )
  n
  iui(v + w)ni
= i=0n(  )
  n
  iui j=0ni(     )
 n −  i
   jvjw◜--◞◟---◝
n − i − k  k
= i=0n j=0ni(  )
 n
 i(     )
  n − i
   juivjwk
= i=0n j=0ni n!
-----
i!j!k!uivjwk
= |I|=nn!
I!μI
= BIn(μ)
__

Linear Precision Si pu˛ dimostrare che vale:

μ = ∑   -IBn (μ)
    |I|=nn   I

la dimostrazione Ŕ componente per componente e si riconduce direttamente alla proprietÓ analoga giÓ dimostrata nel caso monovariato.

Patch triangolari

Ora Ŕ possibile dare una definizione pi¨ rigorosa delle patch triangolari di BÚzier.

Definizione 5.6 (Patch triangolare di BÚzier). Dato un dominio triangolare T, una patch triangolare di BÚzeire di grado n Ŕ una funzione X : T E3 del tipo:

        ∑
X (μ ) =     bIBn (μ)
        |I|=n     I

in cui bI E3 sono chiamati punti di controllo.

ProprietÓ

Comportamento ai bordi Lungo uno dei bordi (ovvero quando uno dei tre parametri in μ Ŕ nullo si ha:

X(  )
 u
|(v |)
 0 = |I|=n
k=0bIBIn(  )
 u
|(v |)
  0
= i=0nb (i,ni,0)----n!-----
(n − i)!i!0!uivniw0
= i=0nb (i,ni,0) (n)
    uivn−iw0
◟i---◝◜-----◞Bin(n), v=1u
ovvero lungo i bordi si ha una curva di BÚzier. Gli altri casi sono analoghi.

Controllo pseudolocale Spostando un punto di controllo:

˜b =  b + δb
 I    I

si ha:

˜X (μ) = X (μ) + δbBnI (μ)

ovvero come nel caso delle curve di BÚzier la modifica di un punto di controllo ha effetti su tutta la curva, in maniera proporzionale al valore del polinomio di Bernstein corrispondente.

End-point interpolation

Ai vertici del triangolo si ha:

X(1 )
|  |
(0 )
 0 = |Ik|==n0bIBIn i=0nb (i,0,0)(1)
| |
(0)
 0 = |Ij|==n0
k=0bIBIn(1)
| |
(0)
 0
= i=nnb (i,0,0)B(i,0,0)n(1 )
|  |
(0 )
 0
= b(n,0,0) n!
n!0!0!-1nv0w0 = b (n,0,0)
La dimostrazione per gli altri due casi Ŕ analoga.

Forma ricorsiva

I polinomi di BÚrnstein generalizzati posso essere espressi in forma ricorsiva:

  n          n−1          n−1          n−1
B I(μ)† = uB I−e1(μ) + vBI− e2(μ) + wB I−e3(μ)
(5.1)

in cui con e1,e2,e3 si indicano i versori unitari lungo le direzioni u,v,w Infatti:

BIe1n1(μ) =  (n − 1)!
-----------
(i − 1)!j!k!ui1vjwk
BIe2n1(μ) =  (n − 1)!
-----------
(j − 1)!i!k!uivj1wk
BIe3n1(μ) = -(n-−-1)!--
(k − 1)!i!j!uivjwk1
e svolgendo la somma (5.1) si ottiene:
[ 1     1    1]        (n − 1)!                  n!
  ---+ ---+ --  ----------------------uivjwk =  -----uivjwk =def BnI(μ)
  jk   ik   ij  (i − 1 )!(j − 1)!(k − 1)!         i!j!k!

Algoritmi

de Casteljau

Come nel caso delle patch quadrate di BÚzeir, l’algoritmo di de Casteljau si definisce a partire dalla definizione ricorsiva del polinomi di Bernstein.
Si definiscono i punti di controllo al passo zeresimo:

b0=  bI  ∀ |I | = n
 I

La superficie Ŕ definita dai punti:

X(μ) = |I|=nbI0B In(μ)
= |I|=nbI0[uBn −1(μ ) + vBn −1 (μ) + wBn −1 (μ)]
    I−e1        I−e2         I−e3
= |J|=n1(ubJ+e10 + vb J+e20 + wb J+e30)B Jn1(μ)

e definendo:

bn = (ubn− 1 + vbn−1  + wbn− 1)
 J      J+e1     J+e2     J+e3

in generale la superficie pu˛ essere rappresentata come:

          ∑
X (μ ) =       brJBnJ−r(μ)
        |J|=n −r

i cui singoli punti sono quindi in:

         n
X (μ) = b(0,0,0)


: de Casteljau per patch triangolari
 
1function Curve = Bezier˙tripatch˙casteljau(b,u,v) 
2Lu = length(u) 
3Lv = length(v) 
4 
5for cu=1:Lu 
6      tcurve = [] 
7      for cv=1:Lv 
8         tcurve = [tcurve, tripatch˙point(b,u(cu),v(cv))] 
9      end 
10      Curve(:,:,cu) = tcurve 
11end 
12endfunction 
13 
14function P = tripatch˙point(b,cu,cv) 
15[d,n,m] = size(b) //n,m must be the same! 
16tb1=b 
17tb2=b 
18cw = 1-cu-cv 
19for r=1:n-1 
20      for cn=1:n-1-r+1 
21         for cm=1:m-cn 
22            tb2(:,cn,cm) = cw*tb1(:,cn,cm)+cu*tb1(:,cn+1,cm)+cv*tb1(:,cn,cm+1) 
23         end 
24      tb1=tb2 
25      end 
26end 
27P = tb1(:,1,1) 
28endfunction


PIC
PIC

Figura 5.13: Esempio di patch di bezier triangolare, da due angolazioni diverse.



PIC
PIC

Figura 5.14: Esempio di patch di bezier triangolare, da due angolazioni diverse.


Degree-elevation

Data la superficie:

         ∑      n
X (μ) =      bIB I(μ ),   μ ∈ T
        |I|=n

la si vuole rappresentare con una patch di bezier di grado pi¨ elevato:

          ∑        n+1
Xe(μ ) =       cIBI   (μ) = X (μ),  μ ∈ T
        |I|=n+1

In quanto u + v + w = 1 si pu˛ scrivere:

X(μ) = |I|=nbIBIn(μ) = (u + v + w) |I|=nbIBIn(μ)
= |I|=nbI-n!--
i!j!k!ui+1vjwk + |I|=nbI--n!-
i!j!k!uivj+1wk + |I|=nbI--n!-
i!j!k!uivjwk+1
= |I|=nbIi-+-1-
n + 1(n-+-1)!-
(I + e1)!UI+e1 + |I|=nbIi +-1-
n + 1-(n-+--1)!
(I + e2)!UI+e2 + |I|=nbI-i +-1
n + 1(n-+-1)!-
(I + e3)!UI+e3
= |J|=n+1(  i             j            k       )
 ------bJ−e1 + -----bJ−e2 + -----bJ− e3
 n + 1         n + 1        n + 1BJn+1(μ)
= |J|=n+1cIBJn+1(μ)asd

Chiamando uivjwk = UI Considerato che (n+-1)!-
(I + e1)!UI+e1 = BI+e 1n+1 e cosý via

ovvero i coefficienti cI possono essere ricavati direttamente:

     (                                      )
         i            j            k
cI =   -----bJ−e1 + -----bJ− e2 + ------bJ−e3
       n + 1        n + 1        n + 1


: Degree elevation per patch triangolari
 
1function d = Bezier˙tripatch˙degelev(b) 
2[dim,n,m] = size(b) 
3d = hypermat([dim,n+1,m+1]) 
4 
5 
6 
7for i=2:n 
8      for j=2:n+2-i 
9         k = (n-j+1-i+1) 
10         d(1:3,i,j)=(1/(n))*((i-1)*b(1:3,i-1,j)+(j-1)*b(1:3,i,j-1)+(k)*b(1:3,i,j)) 
11      end 
12end 
13 
14for j=2:n 
15      k=n+1-1-j+1 
16      d(1:3,1,j) = (1/(n)) * ((j-1)*b(1:3,1,j-1)+(k)*b(1:3,1,j)) 
17end 
18 
19for i=2:n 
20      k=n+1-1-i+1 
21      d(1:3,i,1) = (1/(n)) * ((i-1)*b(1:3,i-1,1)+(k)*b(1:3,i,1)) 
22end 
23 
24d(1:3,1,1) = b(1:3,1,1) 
25d(1:3,n+1,1) = b(1:3,n,1) 
26d(1:3,1,m+1) = b(1:3,1,m) 
27 
28endfunction


PIC
PIC

Figura 5.15: Esempio di degree elevation applicata ad una patch di bezier triangolare, da due angolazioni diverse.



PIC

Figura 5.16: Esempio di degree elevation applicata ad una patch di bezier triangolare, su un griglia piana.


Appendice A
Codici

Qui sono riportati i codici troppo lunghi o non abbastanza importanti per essere interlacciati al testo.
Per provare i vari codici e sufficiente avviare scilab nella directory che contiene tutti gli script e lanciare uno qualsiasi degli script di esempio con il comando
  exec script.sci

All’inizio di ogni script di esempio, tutti denominati esempi˙*.sci, c’Ŕ una lista di selettori per selezionare quali figure si vogliono prodotte.

A.1 BÚzier SPLINE: acquisizione punti condizionale e calcolo


: Acquisizione e calcolo BÚzier spline C0
 
1exec glib.sci; 
2xcol = [0:0.001:1]; 
3gcount=1; 
4g = input(numero e grado curve? ) 
5 
6po(:,:,1) = mouse˙input˙2d(g(2)); 
7plot(po(1,:,1),po(2,:,1),Color, [1/g(1),0.5,1-1/g(1)],LineStyle,--,Marker,o); 
8handler=get(current˙axes); 
9handler.isoview=on;handler.grid=[12,12]; 
10for i=1:g(1)-1 
11po(:,:,i+1) = [po(:,size(po,2),i),mouse˙input˙2d(g(2)-1)]; 
12plot(po(1,:,i+1),po(2,:,i+1),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
13end 
14 
15scf(gcount);gcount=gcount+1; 
16subplot(1,2,1) 
17for i=1:g(1) 
18 c(:,:,i) = Bezier˙casteljau(po(:,:,i),xcol); 
19 plot(po(1,:,i),po(2,:,i),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
20plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
21end 
22handler=get(current˙axes); 
23handler.isoview=on;handler.grid=[12,12]; 
24subplot(1,2,2) 
25for i=1:g(1) 
26rpo = Bezier˙hodograph(po(:,:,i)); 
27c(:,:,i) = Bezier˙casteljau(rpo,xcol); 
28plot(rpo(1,:),rpo(2,:),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
29plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
30end 
31handler=get(current˙axes); 
32handler.isoview=on;handler.grid=[12,12];

: Acquisizione e calcolo BÚzier spline C1
 
1exec glib.sci; 
2xcol = [0:0.001:1]; 
3gcount=1; 
4g = input(numero e grado curve? ) 
5 
6po(:,:,1) = mouse˙input˙2d(g(2)); 
7plot(po(1,:,1),po(2,:,1),Color,[1/g(1),0.5,1-1/g(1)],LineStyle,--,Marker,o); 
8handler=get(current˙axes); 
9handler.isoview=on;handler.grid=[12,12]; 
10for i=1:g(1)-1 
11//tutte le curve hanno lo stesso grado, i coefficienti si semplificano 
12tpo = [po(:,size(po,2),i),(2*po(:,size(po,2),i)-po(:,size(po,2)-1,i))]; 
13plot(tpo(1,:),tpo(2,:),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
14po(:,:,i+1)=[tpo,mouse˙input˙2d(g(2)-2)] 
15plot(po(1,:,i+1),po(2,:,i+1),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
16end 
17scf(gcount);gcount=gcount+1; 
18subplot(1,2,1) 
19for i=1:g(1) 
20 c(:,:,i) = Bezier˙casteljau(po(:,:,i),xcol); 
21 plot(po(1,:,i),po(2,:,i),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
22plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
23end 
24handler=get(current˙axes); 
25handler.isoview=on;handler.grid=[12,12]; 
26subplot(1,2,2) 
27for i=1:g(1) 
28rpo = Bezier˙hodograph(po(:,:,i)); 
29c(:,:,i) = Bezier˙casteljau(rpo,xcol); 
30plot(rpo(1,:),rpo(2,:),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
31plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
32end 
33handler=get(current˙axes); 
34handler.isoview=on;handler.grid=[12,12];

: Acquisizione e calcolo BÚzier spline C2
 
1exec glib.sci; 
2xcol = [0:0.001:1]; 
3gcount=1; 
4g = input(numero e grado curve? ) 
5 
6po(:,:,1) = mouse˙input˙2d(g(2)); 
7plot(po(1,:,1),po(2,:,1),Color,[1/g(1),0.5,1-1/g(1)],LineStyle,--,Marker,o); 
8handler=get(current˙axes); 
9handler.isoview=on;handler.grid=[12,12]; 
10for i=1:g(1)-1 
11//i coefficienti si semplificano se le curve hanno tutte lo stesso grado 
12tpo = [po(:,size(po,2),i)] 
13tpo = [tpo(:,1),(2*tpo(:,1)-po(:,size(po,2)-1,i))] 
14tpo1 = -2*po(:,size(po,2)-1,i) + po(:,size(po,2)-2,i) +2*tpo(:,2) 
15tpo = [tpo,tpo1] 
16 
17 
18plot(tpo(1,:),tpo(2,:),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
19po(:,:,i+1)=[tpo,mouse˙input˙2d(g(2)-3)] 
20plot(po(1,:,i+1),po(2,:,i+1),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
21end 
22 
23scf(gcount);gcount=gcount+1; 
24subplot(1,3,1) 
25for i=1:g(1) 
26 c(:,:,i) = Bezier˙casteljau(po(:,:,i),xcol); 
27 plot(po(1,:,i),po(2,:,i),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
28plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
29end 
30handler=get(current˙axes); 
31handler.isoview=on;handler.grid=[12,12]; 
32subplot(1,3,2) 
33for i=1:g(1) 
34rpo = Bezier˙hodograph(po(:,:,i)); 
35c(:,:,i) = Bezier˙casteljau(rpo,xcol); 
36plot(rpo(1,:),rpo(2,:),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
37plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
38end 
39handler=get(current˙axes); 
40handler.isoview=on;handler.grid=[12,12]; 
41subplot(1,3,3) 
42for i=1:g(1) 
43rrpo = Bezier˙hodograph(Bezier˙hodograph(po(:,:,i))); 
44c(:,:,i) = Bezier˙casteljau(rrpo,xcol); 
45plot(rrpo(1,:),rrpo(2,:),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
46plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
47end 
48handler=get(current˙axes); 
49handler.isoview=on;handler.grid=[12,12];

: Acquisizione e calcolo BÚzier spline G1
 
1exec glib.sci; 
2xcol = [0:0.001:1]; 
3gcount=1; 
4g = input(numero, grado, scala curve? ) 
5k = g(3) 
6 
7po(:,:,1) = mouse˙input˙2d(g(2)); 
8plot(po(1,:,1),po(2,:,1),Color,[1/g(1),0.5,1-1/g(1)],LineStyle,--,Marker,o); 
9handler=get(current˙axes); 
10handler.isoview=on;handler.grid=[12,12]; 
11for i=1:g(1)-1 
12//tutte le curve hanno lo stesso grado, i coefficienti si semplificano 
13tpo = [po(:,size(po,2),i),((k+1)*po(:,size(po,2),i)-po(:,size(po,2)-1,i))/k]; 
14 
15plot(tpo(1,:),tpo(2,:),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
16po(:,:,i+1)=[tpo,mouse˙input˙2d(g(2)-2)] 
17plot(po(1,:,i+1),po(2,:,i+1),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
18end 
19 
20scf(gcount);gcount=gcount+1; 
21subplot(1,2,1) 
22for i=1:g(1) 
23 c(:,:,i) = Bezier˙casteljau(po(:,:,i),xcol); 
24 plot(po(1,:,i),po(2,:,i),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
25plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
26end 
27handler=get(current˙axes); 
28handler.isoview=on;handler.grid=[12,12]; 
29subplot(1,2,2) 
30for i=1:g(1) 
31rpo = Bezier˙hodograph(po(:,:,i)); 
32c(:,:,i) = Bezier˙casteljau(rpo,xcol); 
33plot(rpo(1,:),rpo(2,:),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
34plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
35if i┐1 then 
36plot(k*c(1,:,i),k*c(2,:,i),Color,[i/g(1),0.7,(1-i/g(1))]); 
37end 
38end 
39handler=get(current˙axes); 
40handler.isoview=on;handler.grid=[12,12];

: Acquisizione e calcolo BÚzier spline G2
 
1exec glib.sci; 
2xcol = [0:0.001:1]; 
3gcount=1; 
4g = input(numero, grado, scala1, scala2 curve? ) 
5k = g(3) 
6s = g(4) 
7 
8po(:,:,1) = mouse˙input˙2d(g(2)); 
9plot(po(1,:,1),po(2,:,1),Color,[1/g(1),0.5,1-1/g(1)],LineStyle,--,Marker,o); 
10handler=get(current˙axes); 
11handler.isoview=on;handler.grid=[12,12]; 
12for i=1:g(1)-1 
13//tutte le curve hanno lo stesso grado, i coefficienti si semplificano 
14b1 = ((k+1)*po(:,size(po,2),i)-po(:,size(po,2)-1,i))/k 
15b2=((-s+1)*po(:,size(po,2),i)-2*po(:,size(po,2)-1,i)+po(:,size(po,2)-2,i)+2*s*b1)/s 
16 
17tpo = [po(:,size(po,2),i),b1,b2]; 
18 
19plot(tpo(1,:),tpo(2,:),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
20po(:,:,i+1)=[tpo,mouse˙input˙2d(g(2)-3)] 
21plot(po(1,:,i+1),po(2,:,i+1),Color,[(i+1)/g(1),0.5,1-(i+1)/g(1)],LineStyle,--,Marker,o); 
22end 
23 
24scf(gcount);gcount=gcount+1; 
25subplot(1,3,1) 
26for i=1:g(1) 
27 c(:,:,i) = Bezier˙casteljau(po(:,:,i),xcol); 
28 plot(po(1,:,i),po(2,:,i),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
29plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
30end 
31handler=get(current˙axes); 
32handler.isoview=on;handler.grid=[12,12]; 
33subplot(1,3,2) 
34for i=1:g(1) 
35rpo = Bezier˙hodograph(po(:,:,i)); 
36c(:,:,i) = Bezier˙casteljau(rpo,xcol); 
37plot(rpo(1,:),rpo(2,:),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
38plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
39if i┐1 then 
40plot(k*c(1,:,i),k*c(2,:,i),Color,[i/g(1),0.7,i/g(1)]); 
41end 
42end 
43handler=get(current˙axes); 
44handler.isoview=on;handler.grid=[12,12]; 
45subplot(1,3,3) 
46for i=1:g(1) 
47rrpo = Bezier˙hodograph(Bezier˙hodograph(po(:,:,i))); 
48c(:,:,i) = Bezier˙casteljau(rrpo,xcol); 
49plot(rrpo(1,:),rrpo(2,:),Color,[(i)/g(1),0.5,1-(i)/g(1)],LineStyle,--,Marker,o); 
50plot(c(1,:,i),c(2,:,i),Color,[i/g(1),0.5,1-i/g(1)]); 
51if i┐1 then 
52plot(s*c(1,:,i),s*c(2,:,i),Color,[i/g(1),0.7,i/g(1)]); 
53end 
54end 
55handler=get(current˙axes); 
56handler.isoview=on;handler.grid=[12,12];

A.2 Curve di BÚzeir: esempi del testo


Listing: