1

Il calcolo di radiosity

Daniele Marini

2

Radiosity di superfici finite

BAiEA i

A iBA j

j FA i Aj

Sistema lineare di N equazioni in N incognite:

1-1F11 -1F12 -1F1n

-2F21 1-2F22 -2F2n

-nFn1 -nFn2 1-nFnn

*

B1

B2

Bn

=

E1

E2

En

3

Il sistema lineare

• riscriviamo il sistema lineare da risolvere:Mb = e

• dove: M è la matrice dei fattori di forma, b è il vettore incognito delle radiosity di ogni superficie, e è il vettore dei termini noti, emissività delle superfici

Mi,j = i,j - i Fi,j

• i,j è il delta di Kronecker, =1 per i=j; =0 altrimenti

4

• in forma esplicita il sistema si scrive come:

M i, jB j E i con 1i Nj1

N

M i,i 1 iFi,i

• 0 ≤ i < 1 per ragioni fisiche• 0 ≤ Fi,i ≤ 1 ; in superfici piane Fi,i = 0 (la faccia non

può vedere sè stessa)• quindi Mi,i > 0 e il sistema si può riscrivere:

Bi M i, j

M i,i

B j E i

M i,ij1; ji

N

con 1i N

5

Metodo di Jacobi

• metodo iterativo; si sceglie una stima iniziale del vettore b, a ogni passo k successivo si valuta l’equazione risolvente, usando il vettore già calcolato al passo k-1:

Bi(k )

M i, j

M i,i

B j(k 1)

E i

M i,ij1; ji

N

con 1i N

6

Metodo di Gauss-Seidel

• Modifica di quello di Jacobi, ad ogni passo k si usa sia la valutazione del passo k-1 sia quella del passo k stesso:

Bi(k )

M i, j

M i,i

B j(k )

j1; ji

N

M i, j

M i,i

B j(k 1)

E i

M i,ij1; ji

N

con 1i N

7

Gauss-Seidel

• Vantaggi: più veloce, converge sempre se la matrice è a diagonale dominante (e questo è proprio il caso!):

M i,i 1 iFi,i ma 1 Fi, j

j1

N

quindi :

M i,i Fi, j

j1

N

iFi,i inoltre : i i 1 quindi Fi, j iFi, j

quindi : M i,i > iFi, j

j1

N

iFi,i = iFi, j

j1; ji

N

M i, j

j1, ji

N

8

“guess” iniziale

• I metodi iterativi si basano su una ipotesi iniziale del valore della soluzione, chiamata guess

• in generale si può usare come valore di b quello del vettore di emissività e

9

Gathering vs. Shooting

• L’iterazione con G-S o Jacobi produce un risultato solo alla fine

• il metodo Shooting premette di vedere il risultato in volo

• la complessità è però analoga: si può vedere il risultato parziale in anticipo, poichè tutte le patch hanno fin dall’inizio una radiosità associata

10

Shooting: metodo di Southwell

• Riscriviamo l’equazione come: K = dove:

i AiBi

i AiE i

K i, j Ai

A j

M i, j i, j iF j,i

• K dipende a Fj,i invece che Fi,j come in M (reciprocità)

• a ogni passo k di iterazione k) approssima la radiosità con un errore residuo

11

Shooting (cont)

• il residuo al passo k è:

r(k ) K (k )

• al crescere delle iterazioni il residuo tende a 0

• il metodo procede generando una approssimazione di al passo k+1 modificando una sola componente del vettore

12

Shooting (cont)

• a ogni passo si sceglie la componente di radiosità cui corrisponde il massimo residuo, altri residui possono crescere ma alla fine ciascuno tende a zero.

• l’ordine di aggiornamento non è noto a priori

13

Shooting (cont)

• al passo k si supponga che la componente i dia luogo al residuo massimo ri ; imponendo che esso tenda a zero richiede di risolvere la:

0 r(k1) i K i, j

j1

N

j(k )

• solo la i viene modificata, quindi

j i, j(k1) j

(k )

14

Shooting (cont)• dalla:

• si ricava:

• oppure:

0 r(k1) i K i, j

j1

N

j(k )

i(k1)

1

K i,i

i K i, j

j1

N

j(k )

i(k1) i

(k ) ri

(k )

K i,i

15

Shooting (cont)

• la forma del residuo diventa:

r(k1) K (k1) r(k ) K( (k1) (k ))

• poiché una sola componente di k vien modificata, allora il termine

( (k1) (k ))ha solo un componente diverso da zero (quello di indice i) e basta

una colonna di K per aggiornare il residuo

16

Shooting (cont)

• in conclusione il residuo si calcola:

rj(k1) rj

(k ) K j,i

K i,i

ri(k )

17

Shooting (cont)

• Date le due equazioni:

i(k1) i

(k ) ri

(k )

K i,i

rj(k1) rj

(k ) K j,i

K i,i

ri(k )

• la prima permette il calcolo di radiosità per la patch i, e i residui sono aggiornati con la seconda fino a una soglia minima. Ad ogni passo si usa una sola colonna della matrice K

18

Shooting (cont)

• il guess per può essere 0; il guess per il residuo può essere il vettore di emissività

• da qui il significato fisico: si minimizza il residuo massimo, cioè la maggiore sorgente di luce, e la si “spara” su tutte le patch; quando il residuo è inferiore alla soglia si passa alla sorgente successiva:

19

Shooting (cont)• infatti:

• ma:

• e assumendo Fi,i=0:

• quindi l’energia emessa dalla sorgente si ripartisce su tutte le altre patch

rj(k1) rj

(k ) K j,i

K i,i

ri(k )

K i, j Ai

A j

M i, j i, j iF j,i

rj(k1) rj

(k ) jFi, jri(k )

20

Calcolo del form factor

• trascurando occlusioni si può riformulare:

FAi A j

1

Ai

vij

cos i cos j

ri, j2

Aj

Ai

dAidAj

FA i A j

1

Ai

cosi cos j

ri, j2

yA j

xA i

dxdy

21

Analogia di Nusselt

22

Soluzione analitica

• possibile solo per forme semplici, esistono formule per rettangoli, dischi, cilindri• l’integrale d’area può essere sostituito da un integrale di contorno (teorema di Stokes):

• dove Ci, Cj sono i contorni delle patch Ai, Aj e i due differenziali sono i vettorilungo i contorni– attenzione alle singolarità per r->0

FA i A j

1

2Ai

ln ri, jdxdyC j

C i

23

Integrazione per contorno

• Nel caso di patch poligonali l’integrale si può calcolare esplicitamente, stimando il fattore di forma tra un punto e la patch:

Fx,P 1

2n

gG

.G