Tesi - Simulazione e Rendering di materiali deformabili in ... · Questa caratteristica rende...

POLITECNICO DI BARI

DIPARTIMENTO DI INGEGNERIA MECCANICA E GESTIONALE

Tesi per il conseguimento della

Laurea di Primo Livello in Ingegneria Meccanica

Simulazione e Rendering di materiali

deformabili in Realtà Virtuale

Laureando: Relatore:

Stefano Cristiano Chiar.mo Prof. Antonio Uva

CoRelatore:

Chiar.mo Prof. Giuseppe Monno

A Babbo e Mamma

Sommario

Il Virtual Sculpturing è ancora un campo di ricerca aperto, che si compone simulazioni

basate su modelli fisici ad approcci semplificati. Le tecniche di modellazione

volumetrica sono più intuitive di quelle basate sulle superfici, ma sono limitate dal

consumo di memoria e dalla complessità del calcolo, che risulta generalmente molto

elevata. In questa tesi è stato sviluppato un approccio innovativo di simulazione basato

su particelle che supporta la manipolazione interattiva di materiale deformabile

internamente rappresentato mediante sfere rigide. Inoltre è supportata l’addizione o la

rimozione di materiale in tempo reale, poiché essa avviene senza alcun processo di

ricalcolo della topologia. Questa caratteristica rende possibili manipolazioni arbitrarie

dei volumi di materiali simulati con fori, rotture e cambiamenti di topologia.

Il rendering è fatto visualizzando la isosuperficie generata dalle particelle. L’estrazione

di questa isosuperficie è stata fatta in un primo momento con una versione ottimizzata

del classico algoritmo dei Marching Cubes. Sono in corso esperimenti di rendering

sfruttando le ultimissime potenzialità degli acceleratori grafici presenti sul mercato

(GPU) con risultati molto incoraggianti, tanto da dedicare un ampio capitolo per

descriverle.

Il modello ipotizzato in questa tesi è stato realizzato praticamente in linguaggio C++

sotto piattaforma Windows, utilizzando le librerie grafiche DirectX. Il programma è in

fase di continua evoluzione e può essere utilizzato anche da utenti inesperti per

deformare interattivamente materiale virtuale su un normale PC.

I risultati sia in termini di performance sia in termini di feedback dell’utente sono stati

molto interessanti.

Ringraziamenti

Ringrazio prima di tutti la mia famiglia, mio Padre Antonio, mia Madre Adriana e mio

Fratello Domenico perché è grazie ai loro sacrifici e alla loro unità che sono potuto

giungere a questo importante traguardo.

Un caloroso abbraccio va a Federica, per il suo inestimabile supporto nei momenti più

difficili e per la sua cucina deliziosa, senza la quale sarebbe stato molto più pesante

scrivere questa tesi.

Un ringraziamento speciale va a Dario “Pellicu$” Pelella, per la sincera amicizia e

straordinaria intelligenza che lo contraddistingue.

Ringrazio anche tutti i miei amici e le mie amiche: Erica, Tommaso, Faebio, Sauè,

Giuseppe,Marco,Andrea,Alessandra,Luigi,Danao,Tano,Mimmo,Ruggiero,Daniela,Auro

ra, i componenti del gruppo, Alessandro Special, Tonino Special, Walter Special e

Marco Special.

Un ringraziamento d’oltre Manica a Maurizio “Cthulhu” Sciglio per la sua amicizia e

per le continue e stimolanti discussioni tecniche che facciamo in continazione.

La mia speciale gratitudine va al Prof. Antonio Uva per l’idea di base di questo lavoro

e per le continue discussioni teoriche e tecniche che hanno portato alla sua

realizzazione. Un rigraziamento è dovuto al Prof. Giuseppe Monno per la sua attenta

supervisione ed il suo incoraggiamento. Ringrazio infine il dott. Michele Fiorentino,

che ha visionato questa tesi.

Indice Generale

Capitolo 1: Introduzione ......................................................................................................................1

1.1 REAL-TIME VIRTUAL SCULPTURING ..............................................................................................1 1.1.1 Classificazione .....................................................................................................................2

Modelli superficiali ..................................................................................................................................2 Modelli volumetrici..................................................................................................................................2

Modelli Basati sulla fisica ...................................................................................................................3 Modelli Semplificati............................................................................................................................3 Sistemi Particellari ..............................................................................................................................3

1.2 SISTEMI PARTICELLARI..................................................................................................................4 1.2.1 Sistemi indipendenti e sistemi collegati .................................................................................5 1.2.2 Sistemi particellari nel campo della fisica .............................................................................5 1.2.3 Sistemi particellari in Computer Graphics ............................................................................6

1.3 CONSIDERAZIONI ..........................................................................................................................7

Capitolo 2: Simulazione .......................................................................................................................8

2.1 PRIMITIVE ....................................................................................................................................8 2.1.1 Ellissoidi e sfere ................................................................................................................ 10 2.1.2 Ellissoidi ............................................................................................................................ 10 2.1.3 Sfere................................................................................................................................... 11

2.2 REGOLE DEL SISTEMA ................................................................................................................. 11 2.2.1 Dinamica Molecolare ......................................................................................................... 11 2.2.2 Potenziale di Lennard-Jones (LJ)........................................................................................ 12 2.2.3 Tensione superficiale .......................................................................................................... 13 2.2.4 Dalle forze agli spostamenti................................................................................................ 14 2.2.5 Funzione di risposta ........................................................................................................... 16

2.3 RISOLUZIONE ITERATIVA DEI VINCOLI. RELAXATION .................................................................... 19 2.3.1 Campi di utilizzo del “relaxation” ...................................................................................... 19 2.3.2 Metodi verlet ...................................................................................................................... 20

2.4 IMPACCHETTAMENTO.................................................................................................................. 21 2.4.1 Densità............................................................................................................................... 21 2.4.2 Impacchettamenti CFC e Esagonale.................................................................................... 22

2.5 STRUTTURE DI PARTIZIONAMENTO SPAZIALE ............................................................................... 23 2.5.1 Introduzione alle strutture di partizionamento..................................................................... 24 2.5.2 Bounding Volume Hierarchies ............................................................................................ 25 2.5.3 Sphere Tree ........................................................................................................................ 27 2.5.4 Octree ................................................................................................................................ 29 2.5.5 Loose Octree ...................................................................................................................... 31

2.6 POLITICA DEGLI AGGIORNAMENTI................................................................................................ 32 2.7 STRUMENTI DI MODELLAZIONE .................................................................................................... 34

2.7.1 Strumenti composti di sfere ................................................................................................. 34 2.7.2 Strumenti poligonali ........................................................................................................... 36 2.7.3 Performance....................................................................................................................... 36

Capitolo 3: Rendering in CPU ........................................................................................................... 38

3.1 INTRODUZIONE ........................................................................................................................... 39 3.1.1 Superfici Implicite .............................................................................................................. 40 3.1.2 Funzione D1....................................................................................................................... 40 3.1.3 Funzione D2....................................................................................................................... 43 3.1.4 Funzione D3....................................................................................................................... 43 3.1.5 Funzione D4....................................................................................................................... 43 3.1.6 Funzioni a paragone........................................................................................................... 44

3.2 MARCHING CUBES ...................................................................................................................... 45 3.2.1 Marching Squares .............................................................................................................. 45 3.2.2 L’algoritmo ........................................................................................................................ 47 3.2.3 Ottimizzazioni..................................................................................................................... 50

Ricalcolo zone modificate ......................................................................................................................50 Marcia Guidata ......................................................................................................................................50 Utilizzo della struttura di partizionamento...............................................................................................51

3.2.4 Modalità di calcolo della normale....................................................................................... 51 3.2.5 Rendering adattativo .......................................................................................................... 53 3.2.6 Requisiti di memoria........................................................................................................... 54 3.2.7 Considerazioni Finali ......................................................................................................... 54

Capitolo 4: Rendering in GPU ........................................................................................................... 55

4.1 CLASSIFICAZIONE DELLE MODERNE TECNICHE DI VISUALIZZAZIONE.............................................. 57 4.1.1 Campi vettoriali.................................................................................................................. 58 4.1.2 Simboli matematici ............................................................................................................. 59 4.1.3 Campi scalari 2D ............................................................................................................... 59 4.1.4 Campi scalari 3D ............................................................................................................... 60

4.2 TEORIA DEL RENDERING VOLUMETRICO ...................................................................................... 61 4.3 TECNICHE DI RENDERING VOLUMETRICO..................................................................................... 63

4.3.1 Texture 3D ......................................................................................................................... 63 Definizione............................................................................................................................................63 Filtering.................................................................................................................................................64 Memoria................................................................................................................................................65 Applicazioni ..........................................................................................................................................65

4.3.2 Texture slices...................................................................................................................... 66 Texture Slices 2D...................................................................................................................................66 Texture Slices 3D...................................................................................................................................67

4.3.3 Pre integrated volume rendering......................................................................................... 69 4.3.4 GPU Raytraced volumetric Rendering ................................................................................ 71

Introduzione al raytracing.......................................................................................................................71 L’algoritmo............................................................................................................................................72 Problemi dello Step Vector.....................................................................................................................73 Raytracing con Transfer Function...........................................................................................................73 Raytracing ‘First Hit’ .............................................................................................................................74 Calcolo della normale.............................................................................................................................75 Migliorie ...............................................................................................................................................76 Z-Buffering............................................................................................................................................78

Capitolo 5: L’applicativo.................................................................. Errore. Il segnalibro non è definito.

Capitolo 6: Appendice A .................................................................................................................... 91

6.1 TEST DI INTERSEZIONE SFERA-CUBO ........................................................................................... 91 6.2 FUNZIONE DI RISPOSTA................................................................................................................ 93 6.3 CALCOLO RAGGIO FUNZIONE DI INFLUENZA DEL MC.................................................................... 93 6.4 RICERCA DELL’OTTANTE GIUSTO DEL LOOSEOCTREE .................................................................... 94 6.5 CODICE DI CRAZIONE TEXTURE PRE-INTEGRATA ........................................................................... 95 6.6 DETTAGLI TECNICI SUL GPU RENDERING ATTUALE ...................................................................... 96

6.6.1 La pipeline ......................................................................................................................... 96 6.6.2 Vertex Program.................................................................................................................. 97 6.6.3 Pixel (o fragment) program................................................................................................. 99 6.6.4 Novità dei modelli VS e PS 3.0.......................................................................................... 101

Dichiarazioni Input/Output...................................................................................................................101 Predicazione ........................................................................................................................................102 Controllo di flusso statico e dinamico ...................................................................................................103 Altre novità..........................................................................................................................................105

Conclusioni ....................................................................................................................................... 106

Indice Figure

Figura 1. Sistemi di riferimento ....................................................................................................4

Figura 2. Il programma Paricle Juice ...........................................................................................7

Figura 3. Un modello superficiale a particelle orientate [18]........................................................9

Figura 4. Rottura a trazione e a taglio simulati in MD [21]........................................................ 12

Figura 5. Potenziale LJ (rosso) and forza (blu) funzione della distanza (r). .............................. 13

Figura 6. Forze repulsive tra particelle....................................................................................... 15

Figura 7. Forze repulsive tra particelle....................................................................................... 18

Figura 8. Effetti del potenziale LJ sulle particelle ...................................................................... 18

Figura 9. Varie tipologie di impacchettamento [27] ................................................................... 22

Figura 10. Costruzione geometrica dell’impacchettamento ......................................................... 22

Figura 11. Impacchettamento esagonale centrato ...................................................................... 262

Figura 12. Una Bounding Volume Hierarchy e la sua rappresentazione ad albero..................... 26

Figura 13. Una Bounding Sphere ‘obesa’ su un oggetto lungo e stretto....................................... 27

Figura 14. La costruzione di uno Sphere Tree e ........................................................................... 30

Figura 15. La costruzione di un Octree ........................................................................................ 30

Figura 16. Comparazione tra un Octree standard ed un Loose Octree ....................................... 31

Figura 17. Fasi successive dello schema di propagazione della deformazione ............................. 32

Figura 18. Generazione di sample points (Monte Carlo) [31] ...................................................... 33

Figura 19. Un tool bidimensionale composto da sfere .................................................................. 34

Figura 20. Flood-Fill per riempire di particelle una forma chiusa............................................... 35

Figura 21. Grafico delle perfomance del sistema.......................................................................... 37

Figura 22. Flood-Fill per riempire di particelle una forma chiusa............................................... 39

Figura 23. Una metaball................................................................................................................ 41

Figura 24. Due Metaball ............................................................................................................... 41

Figura 25. Allontanamento delle metaball.................................................................................... 42

Figura 26. Variazione del fattore di iso-livello.............................................................................. 42

Figura 27. Metaball usando linee come primitive......................................................................... 42

Figura 28. Grafici sovrapposti delle funzioni D1,D2,D3,D4 ......................................................... 44

Figura 29. Ingrandimento di una cella di confine con Marching Square..................................... 45

Figura 30. L’approssimazione fatta dal Marching Square........................................................... 46

Figura 31. Tutti i possibili casi del marching square.................................................................... 47

Figura 32. Caso ambiguo di Marching Square............................................................................. 47

Figura 33. Casi ‘base’ dei Marching Cubes.................................................................................. 48

Figura 34. Metaball usando linee come primitive......................................................................... 49

Figura 35. Risoluzione casi ambigui Marching Cubes ................................................................. 49

Figura 36. Si cammina esaminando le celle di superficie contigue. .............................................. 50

Figura 37. Normale calcolata come risultante normali ai poligoni vicini..................................... 52

Figura 38. Rendering adattativo a diversi livelli di LOD ............................................................. 53

Figura 39. Rappresentazione di pezzi di una forma ai diversi livelli LOD .................................. 53

Figura 40. Curva del numero di transistor installati sulle GPU recenti....................................... 55

Figura 41. Confronto qualitativo tra performance GPU/CPU..................................................... 56

Figura 42. Il tipo di dati implica un determinato algoritmo [34].................................................. 57

Figura 43. Classificazione delle varie tipologie di rendering [34]................................................. 57

Figura 44. Visualizzazione di Campi vettoriali fluidodinamici e streaklines [34]........................ 58

Figura 45. Visualizzazione di un simbolo matematico [34] .......................................................... 59

Figura 46. Visualizzazione di un HeightField [34]........................................................................ 59

Figura 47. Campo scalare 3D ‘Voxel’ [38].................................................................................... 60

Figura 48. Filtri comunemente usati nella ricostruzione [38]...................................................... 61

Figura 49. Applicazione di una TF per evidenziare l’interno di una scansione[37]..................... 61

Figura 50. Integrale di Volume Rendering (parzialmente preso da [38]) .................................... 62

Figura 51. Comparazione tra una Texture 2D e una Texture 3D [42] ......................................... 63

Figura 52. Comparazione Filtering Anisotropico su Texture 2D e Texture 3D [42].................... 64

Figura 53. Composizione di ‘fette’ 2D (slice) per il rendering 3D [38]......................................... 66

Figura 54. Composizione di slice allineate all’oggetto [38]........................................................... 66

Figura 55. Problema del ‘cambio’ del set di slices [38]................................................................. 67

Figura 56. Blending di una slice intermedia interpolata tra due diverse slices [38]..................... 67

Figura 57. Slices Orientate alla camera [39]................................................................................. 68

Figura 58. A destra una singola slice. A sinistra una composizione di slice ................................. 68

Figura 59. Due Transfer Function diverse sulo stesso oggetto ..................................................... 69

Figura 60. Un raggio che passa attraverso due slice..................................................................... 70

Figura 61. Funzionamento basilare di un raytracer classico........................................................ 71

Figura 62. Cubo con vertici mappati sulla texture 3D.................................................................. 72

Figura 63. Lo step vector interpolato dal PS si ‘accorcia’ nei pixel centrali ................................ 73

Figura 64. Step successivi nel pixel shader per trovare la superficie ........................................... 74

Figura 65. Algoritmo approssimato per il calcolo della normale nel PS ...................................... 76

Figura 66. Caso ‘difficile’ da trattare nel calcolo della normale .................................................. 77

Figura 67. Un oggetto esterno dovrebbe ‘bloccare’ i raggi........................................................... 78

Figura 68. L’applicazione realizzata per questa tesi .................................................................... 79

Figura 69. La rappresentazione interna del blocco di materiale.................................................. 80

Figura 70. Un teschio proveniente da una scansione MRI ........................................................... 80

Figura 71. Modalità di deformazione 2D...................................................................................... 81

Figura 72. Stacchiamo un pezzo di materiale per usarlo come tool ......... Errore. Il segnalibro non è

definito.

Figura 73. Risultati della deformazione........................................................................................ 82

Figura 74. Seleziono particelle del tool ......................................................................................... 82

Figura 75. Le particelle vicine al tool rimangono ‘attaccate ........................................................ 83

Figura 76. Il blocco di materiale e la punta pronta per ‘incidere’ qualcosa ................................ 83

Figura 77. Trasciniamo la punta e incidiamo il materiale............................................................ 84

Figura 78. Continuiamo il percorso .............................................................................................. 84

Figura 79. Chiudiamo il percorso ................................................................................................. 85

Figura 80. Un blocco composto da poche particelle...................................................................... 85

Figura 81. Iniziamo a premere sul materiale. La superficie non si rompe................................... 86

Figura 82. La superficie si è rotta. ................................................................................................ 86

Figura 83. Il buco .......................................................................................................................... 87

Figura 84. Il buco visto da lontano................................................................................................ 87

Figura 85. Ecco la ‘vittima’ da rompere....................................................................................... 88

Figura 86. La pallina tocca la superficie ....................................................................................... 88

Figura 87. La pallina entra nella superficie.................................................................................. 89

Figura 88. Oramai il danno è fatto................................................................................................ 89

Figura 89. Alcune immagini random............................................................................................ 90

Figura 90. Trasformazioni geometriche nel Test AABB-Sfera .................................................... 92

Figura 91. Numerazione dei nodi figli dell’octree......................................................................... 94

Figura 92. Diagramma schematico di una Pipeline di GPU moderna.......................................... 96

Figura 93. Il Vertex Program all’internod della pipeline............................................................. 97

Figura 94. Comunicazione unidirezionale da CPU a GPU........................................................... 99

Figura 95. Posizionamento del Pixel Program nella Pipeline ....................................................... 99

Figura 96. Rasterizing di un triangolo........................................................................................ 100

Figura 97. Il Pixel Program viene applicato ad ogni pixel disegnato ......................................... 100

Tesi per il conseguimento della Laurea in Ingegneria Meccanica L3

______________________________________________________________________________________________

______________________________________________________________________________________________

Stefano Cristiano; Simulazione e Rendering di materiali deformabili in Realtà Virtuale Pag. 1

Capitolo 1: Introduzione Nel campo della Computer Graphics (CG), la modellazione della forma di oggetti con

proprietà tipiche dei solidi deformabili è un campo di ricerca molto attivo ed aperto. Lo

scopo finale di tale ricerca è quello di creare delle tecniche di modellazione di forme

che potrebbero essere utilizzate nei campi più disparati. Le applicazioni nei campi

dell’istruzione, arte, entertainment, realtà virtuale sono immediate, ma anche il campo

del disegno industriale potrebbe beneficiare di tali strumenti. Molti designer e artisti

preferiscono provare le loro idee in forma “reale”, con carta e penna oppure utilizzando

materiale deformabile tipo argilla, invece di farlo direttamente su un computer. Allo

stato attuale non esistono tecniche di modellazione volumetrica real-time del tutto

generali e convincenti.

Questo lavoro di tesi si propone di descrivere un nuovo modello basato sulla

modellazione di particelle, che permetta di deformare interattivamente dei materiali,

usando come primitive delle sferette di raggio unitario. Tale sistema è stato

implementato praticamente e può essere facilmente utilizzato da utenti non esperti per

deformare blocchi di materiale virtuali con un certo grado di realismo, su un normale

PC. Il rendering viene effettuato al momento utilizzando una versione ottimizzata del

famoso algoritmo di Marching Cubes (MC), che sfrutta le informazioni di

partizionamento spaziale del simulatore. I risultati sono molto incoraggianti: le

performance sono ampiamente sopra la soglia minima del real-time e la sensazione di

manipolare materiale deformabile è fortemente presente.

1.1 Real-Time Virtual Sculpturing Il termine “virtual sculpturing” (VS) indica quel settore della modellazione volumetrica

CAD che simula il processo di interazione tra un utente ed un blocco di materiale

“virtuale”.

Il termine “real-time” (in tempo reale) invece vuole sottolineare la scelta di base che è

stata fatta. La scelta è stata quella di ipotizzare un sistema che fosse realizzabile

praticamente sfruttando le tecnologie odierne, sfruttando l’odierna potenzialità di

computazione dei normali PC.


______________________________________________________________________________________________

______________________________________________________________________________________________


1.1.1 Classificazione Il campo della modellazione volumetrica in tempo reale è in una fase di costante

crescita, in svariate direzioni. I diversi “filoni” di ricerca che sono nati vengono portati

avanti parallelamente, e molto spesso è possibile assistere a degli scambi proficui tra di

essi.

Le varie tecniche di VS possono essere classificate in base a quattro attributi principali:

modelli basati sulla fisica, modelli semplificati, volumetrici e di superfici. Ognuno di

questi può essere simulato in modi diversi, ma la gran parte di essi utilizza sistemi

particellari (approccio Lagrangiano) o griglie di densità di una qualche sostanza

(approccio Euleriano).

Esistono modelli fisici superficiali, modelli fisici volumetrici, ma anche modelli

semplificati volumetrici o superficiali. Esistono addirittura modelli che sono

contemporaneamente volumetrici e superficiali ma basati sulla fisica reale, come ad

esempio [18].

Se si concentra l’attenzione sul real-time, le tecniche di modellazione volumetrica sono

più intuitive di quelle superficiali, ma fortemente limitate dalle potenzialità dei moderni

calcolatori, sia in termini di memoria che di tempo computazionale. La complessità

degli algoritmi associati a tali tecniche spesso esplode in modo incontrollabile. Di

seguito verranno brevemente esposti i principali ‘filoni’ di ricerca che sono presenti in

letteratura, indicando in quale di essi intende collocarsi questo lavoro.

Modelli superficiali I modelli basati su superfici poligonali ([10][11][12][13]) sono facili da implementare e

molto veloci, ma richiedono trucchi complessi per assicurare la possibilità di eseguire

grandi deformazioni. La gran parte di questi modelli non assicura la conservazione del

materiale., anche se esistono alcuni modelli [19] che sfruttano leggi fisiche, ad esempio

l’equazione dei gas perfetti per ovviare a questo inconveniente.

I modelli superficiali vengono spesso accoppiati anche con lo splatting [20] quando si

utilizzano particelle orientate [18]

Modelli volumetrici Le tecniche di modellazione volumetrica sono sicuramente le più intuitive, giacché

l’utente può attualmente deformare il materiale come farebbe nel mondo reale, con vera


______________________________________________________________________________________________

______________________________________________________________________________________________


‘creta’ o argilla e veri strumenti, per effettuare operazioni fisiche come premere, curvare

e stirare il materiale.Esistono soluzioni che modellano materiali deformabili, materiali

elastici ed anche modelli ibridi,attraverso dei ‘trucchi’ opinabili.

Modelli Basati sulla fisica La ricerca corrente sui modelli phisically-based [7] si è evoluta nella creazione di grandi

sistemi del tipo molla-smorzatore, da risolvere ad ogni frame. Questi sistemi sono molto

costosi a livello computazionale, e spesso non è realistico poterli applicare in tempo

reale. Questo problema aumenta ancor di più quando bisogna effettuare deformazioni

che richiedono grandi cambiamenti topologici dell’oggetto in tempo reale. Tali

deformazioni infatti richiederebbero il ricalcolo completo dei ‘collegamenti’ tra le

‘particelle’ che compongono il sistema.

Modelli Semplificati Alcuni ricercatori [1][2][6] hanno deciso di abbandonare la correttezza fisica del

modello, concentrando i loro sforzi nello sviluppo di approcci semplificati che possono

effettivamente essere simulati in tempo reale. La maggior parte dei modelli volumetrici

semplificati vengono definiti utilizzando campi scalari, che vengono salvati in una

griglia o controllati attraverso primitive. La superficie di un oggetto è nella maggior

parte dei casi una isosuperficie di questo campo, che gioca un ruolo simile alla densità

di materiale. Questa rappresentazione è molto utile per operazioni di deformazione

locale, come quelle descritte precedentemente,ma anche l’addizione di materiale risulta

semplificata. Quest’ultima infatti viene effettuata semplicemente attraverso una

modifica locale dei valori del campo. Queste rappresentazioni volumetriche quindi

possono gestire qualunque tipo di deformazione che implichi cambiamenti di topologia

senza il bisogno di codificare meccanismi specifici.

Una soluzione utilizzando i cellular automata è stata proposta per calcolare le

deformazioni locali su materiale di argilla virtuale. [3][4]

Sistemi Particellari Gran parte dei modelli qui esposti utilizzano sistemi particellari. Essi sono trattati in

maggior dettaglio nel prossimo paragrafo poiché sono la base del sistema simulativo

realizzato.


______________________________________________________________________________________________

______________________________________________________________________________________________


1.2 Sistemi particellari I sistemi particellari [18] sono formati da un insieme di masse puntiformi con delle

forze ad essere associate che si muovono seguendo determinate leggi della fisica. La

letteratura sull’argomento è vastissima e facilmente reperibile (Tonnesen 1991; Szeliski

and Tonnesen 1992; Szeliski, Tonnesen and Terzopoulos, 1993; SZeliski, Tonnesen and

Terzopoulos,1993)

Figura 1. Sistemi di riferimento

Ad ogni particella vengono assegnati un insieme di attributi, come

massa,posizione,velocità e accelerazione, che ne descrivono lo stato.

I potenziali sono generalmente utilizzati per generare forze che agiscano sulle particelle

stesse, e successivamente il movimento viene imposto dalle leggi della fisica

Newtoniana e cioè:

)(1 ;)()(

dttd

mt i

i

i vf= )(2 ;

)()(

dttd

t ii

xv =

Dove le fi sono le forze, le vi sono le velocità, le xi sono le posizioni e le mi sono le

masse dell’i-esima particella. Partendo da determinate condizioni iniziali, questi sistemi

possono essere simulati nel tempo integrando le equazioni del moto. A seconda delle

forze applicate, questi sistemi possono modellare una verità di sistemi complessi, che si

evolvono nel tempo.

Questo è un approccio che potremmo definire “Lagrangiano” a differenza di quello

euleriano dove lo stato del sistema è definito in punti fissi di una griglia.


______________________________________________________________________________________________

______________________________________________________________________________________________


Il sistema presentato in questa tesi è fondamentalmente Lagrangiano, ma la politica

degli aggiornamenti considera, come vedremo più avanti, una serie di “celle” nello

spazio in cui sono contenute delle particelle da aggiornare tutte contemporaneamente,

nel quale si può leggere, forse con un po’ di forzatura, uno schema Euleriano.

1.2.1 Sistemi indipendenti e sistemi collegati Un’altra distinzione che si fa all’interno dei sistemi particellari è quella che li

categorizza in base all’interazione tra le particelle.

Un primo sistema è quello chiamato delle particelle indipendenti, dove le forze su ogni

particella sono indipendenti dalle restanti particelle del sistema. Le particelle

indipendenti hanno il vantaggio di poter essere trattate separatamente l’una dall’altra ma

allo stesso tempo, mancando di informazioni topologiche si necessita di grandi numeri

di essere per “riempire” tutto lo spazio rinchiuso da una determinata forma.

Il secondo sistema è detto delle particelle collegate, dove ogni particella è collegata ad

altre particelle attraverso dei connettori. Un sistema rigido di punti materiali cade ad

esempio in questa categoria; i vincoli di rigidità assolvono alla funzione di “connettori”

delle particelle.

Ovviamente cambiando il tipo di connettore si possono simulare sistemi anche molto

diversi dai corpi rigidi. Inserendo delle molle ad esempio come connettori, si possono

simulare tessuti oppure solidi con proprietà elastiche. Questo sistema ha il vantaggio di

avere una rappresentazione più compatta nello spazio, poiché per rinchiudere un volume

è necessario definire e collegare solo le particelle “di superficie” o, più

matematicamente, definite sulla frontiera del nostro oggetto stesso.

1.2.2 Sistemi particellari nel campo della fisica Nel campo della fisica i sistemi particellari sono stati usati per modellare una varietà di

fenomeni, tra cui l’evoluzione delle galassie, il plasma, le proprietà dei campi

magnetici, dei gas comprimibili, dei cambiamenti di fase della materia etc.

In tutti questi sistemi ogni particella rappresenta un elemento primitivo del fenomeno

che si sta studiando, ad esempio una stella o una molecola di gas o un elettrone. Il

problema fondamentale di questi approcci è che le simulazioni richiedono delle

interazioni complesse tra le particelle con alto grado di precisione numerica. Se poi si

considera il fatto che per ottenere risultati degni di una certa validità scientifica, bisogna

utilizzare numeri di particelle dell’ordine dei miliardi, si capisce come allo stesso


______________________________________________________________________________________________

______________________________________________________________________________________________


tempo questo potente approccio sia limitato dalle possibilità di calcolo e dalla memoria

dei moderni elaboratori.

E’ per questo motivo che in tutti i campi la ricerca si spinge nel cercare delle opportune

ipotesi semplificative per evitare di risolvere tutte le possibili interazioni tra particelle e

resto del sistema e tra le particelle stesse. Ad esempio nel campo dell’astrofisica

(Hockey e Eastwoord, 1988, Hegie, 1987) per modellare l’evoluzione delle stelle ci si

basò sulla densità dei corpi e sui campi gravitazionali. Ebbene, per modellare piccoli

numeri di stelle, nell’ordine delle migliaia, il sistema poteva semplicemente considerare

ogni particella una stella e calcolare le forze direttamente.

Tale sistema è ovviamente di complessità )( 2NO per N particelle. Per sistemi più

grandi si tende ad organizzare tutto in “clusters”, in regioni di influenza, considerando

ad esempio che le forze di attrazione tra corpi vanno quasi a zero dopo un determinato

range, rendendo inutile l’interazione tra corpi lontani.

La dinamica molecolare allo stesso modo distingue delle interazioni a piccola distanza,

proporzionali ad nr − , modellando invece le forze a larga distanza come mr − . Allo

stesso modo si simulano fenomeni di temperatura.

Il concetto fondamentale alla base di tutti questi esempi è quello di eliminare le

interazioni non necessarie oltre una certa distanza, ed è questo il ragionamento alla base

del nostro schema di risoluzione delle equazioni, a livello della singola particella, con

l’ambiente esterno.

1.2.3 Sistemi particellari in Computer Graphics Nella grafica al calcolatore, i sistemi particellari sono la norma per quanto riguarda la

simulazione di fenomeni naturali visivamente complessi, quali ad esempio fumi, nuvole,

fluidi ma anche nell’emulazione di solidi deformabili,elastici o rigidi. E’ da sottolineare

che in questo campo ciò che conta è l’effetto visivo finale, e non l’accuratezza della

simulazione come in campo scientifico.

Questi sistemi sono in genere sistemi di particelle indipendenti, così come sono state

definite precedentemente. Per creare effetti complessi, queste tecniche usano un gran

numero di particelle che sono soggette a forze, come ad esempio la gravità, a collisioni

col mondo esterno e a campi di forze esterne, quali ad esempio un vento e la turbolenza.

Le particelle vengono create e distrutte nel tempo in base alle regole stabilite dal

particolare fenomeno che si sta modellando.


______________________________________________________________________________________________

______________________________________________________________________________________________


Il Laureando ha avuto diverse esperienze riguardo i sistemi particellari in computer

grafica, tanto da mettere a punto un programma che permette di controllare e cambiare

in tempo reale una moltitudine di parametri che governano il sistema. La potenza del

sistema deriva anche dal fatto che tali parametri possono anche essere funzioni del

tempo, cioè è possibile “disegnarne” l’andamento con il mouse. Questo programma è

stato già utilizzato per la realizzazione di effetti grafici in un prodotto multimediale

commerciale attualmente in fase di sviluppo da parte di un azienda di Bari.

Figura 2. Il programma Paricle Juice

1.3 Considerazioni E’ stato ideato in questa tesi un modello volumetrico particellare, non basato sulla

fisica. Sono queste le condizioni necessarie che ne permettono una implementazione in

real-time che sia rispondente alle finalità del lavoro stesso. L’utilizzo di un modello

volumetrico, come sarà spiegato nel capitolo della simulazione, non impedirà di

simulare effetti come quelli di tensione superficiali che in genere si ritiene simulabili

solamente con modelli prettamente superficiali.


______________________________________________________________________________________________

______________________________________________________________________________________________


Capitolo 2: Simulazione La simulazione del sistema è effettuata attraverso un sistema particellare volumetrico.

Di conseguenza è necessario definire quali primitive si intende utilizzare e quali siano le

regole che governano il sistema stesso. Successivamente si vedrà come possano essere

soddisfatte queste regole e quali sono gli schemi algoritmici che si possono

implementare per accelerare i calcoli. Infine si descriverà il modello ideato per la

gestione dei tool.

2.1 Primitive La scelta fondamentale che bisogna fare nel progettare un sistema particellare è quello

di definire la o le primitive associate alle nostre particelle. Questa scelta è fondamentale

perché molte delle proprietà del sistema,nonché la sua complessità, dipendono dalle

proprietà delle singole particelle.

Una prima distinzione può essere quella di scegliere particelle orientate oppure

particelle normali.

La configurazione geometrica di una particella normale nello spazio è definita

semplicemente da una terna di valori che ne individua la posizione rispetto ad un fissato

sistema di assi cartesiano. Questo genere di particelle può anche avere un volume, ma si

ignora l’orientamento della superficie che lo racchiude. Tale superficie si muoverà nel

tempo rimanendo sempre parallela a se stessa, ed è per questo che in genere le primitive

utilizzate in questo caso hanno delle forti simmetrie si più assi.

Le particelle orientate invece hanno bisogno anche di una terna di coseni direttori

rispetto ad un sistema di riferimento, tale da fissarne l’orientamento nello spazio.

Questo genere di particelle è in generale più costoso da utilizzare, in quanto bisogna

tirare in gioco anche le equazioni di equilibrio a momento, considerando spesso le

singole particelle come veri e propri corpi rigidi. Questo genere di particelle si può

trovare in [18].


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 3. Un modello superficiale a particelle orientate [18]

Come si può vedere nella figura, nel paper citato si usano, nell’ambito del free-form

modelling, dei dischi circolari orientati come primitive. Queste primitive devono essere

tangenti ad altri dischi orientati (con un vincolo di attrazione che funge da connettore) e

insieme quindi definiscono un volume racchiuso che rappresenta l’oggetto in questione.

Nell’esempio possiamo vedere una sfera che viene trasformata in un toroide da due

sfere rigide che fungono da “scalpello”.


______________________________________________________________________________________________

______________________________________________________________________________________________


Questo lavoro si propone di creare un sistema che funzioni in modo analogo a quello

presentato, ma che sia molto più performante. Consultando la letteratura

sull’argomento, in molti casi si può vedere come vengano usate le particelle orientate

per definire delle superfici, ma in tutti questi lavori è presente il problema irrisolto della

performance. Nessuno di essi riesce a gestire un numero elevato di particelle in tempo

reale sfruttando la potenza attuale dei calcolatori.

2.1.1 Ellissoidi e sfere Tra le varie possibilità che sono state esplorate c’è stata quella di utilizzare delle

primitive orientate, primi tra tutti gli ellissoidi. L’idea alla base era quella di poter

orientare questi ellissoidi in modo da farli allineare lungo piani preferenziali che

potessero definire una superficie, anziché riempire un volume. In questo frangente si è

valutato cosa avrebbe comportato il fatto di utilizzare ellissoidi di raggi variabili al

posto di utilizzare un raggio fisso.

2.1.2 Ellissoidi Il primo problema che si è posto, e che è risultato insormontabile visti i presupposti di

funzionamento in real-time di questo lavoro, è stato quello dei test di intersezione

ellissoide/ellissoide. Se entrambi gli ellissoidi fossero stati di raggi identici, sarebbe

basta una semplice trasformazione di scala per portarli in uno spazio nel quale essi

assumano la forma di circonferenze, delle quali è banale trovare l’intersezione.

Il secondo problema, ancora più insormontabile del primo, consiste nel fatto che per lo

schema del sistema particellare che si vuole progettare c’è bisogno di una informazione

aggiuntiva nei test di intersezione, la distanza di penetrazione.

Come in molti metodi Impulse-Based, nell’ambito della fisica dei rigidi al calcolatore,

la conoscenza della distanza di penetrazione tra un timestep ed un altro permette una

semplice risoluzione della funzione di risposta all’urto, se così lo vogliamo chiamare.

Un test di distanza di penetrazione tra ellissoidi orientati a raggio variabile è troppo

costoso per poter essere implementato su un gran numero di particelle in tempo reale,

per tale motivo questa opzione è stata scartata.

Una volta scartata l’idea di poter avere raggio variabile e di poter orientare gli ellissoidi,

si è deciso di abbandonare questa strada perché non portava nessun vantaggio

sostanziale in nessun ambito.


______________________________________________________________________________________________

______________________________________________________________________________________________


2.1.3 Sfere La sfera è stata la seconda forma geometrica della quale sono state studiate le

potenzialità.

Il problema ancora una volta però, è quello di capire se ci siano dei reali vantaggi

nell’utilizzo di un raggio variabile. Ancora una volta però si è dovuta constatare

l’inutilità di tale caratteristica per il sistema ideato, che preferisce approssimare le forme

con un grande numero di piccole particelle, nelle quali non sarebbe possibile notare le

differenze di raggio tra una particella ed un’altra. Il test di distanza di penetrazione tra

sfere di raggi unitari è oltremodo leggero sotto il profilo computazionale.

2.2 Regole del sistema Come è stato spiegato nella breve introduzione ai sistemi particellari, c’è bisogno di

stabilire delle regole che governino il sistema ideato. Le regole consistono

sostanzialmente in un insieme di equazioni da soddisfare. Queste equazioni, assieme ad

i metodi di integrazione, devono però essere ben studiate, poiché è molto semplice

portare il sistema in condizioni di “esplosione”, dove in sostanza il sistema perde di

significato fisico a causa delle pesanti discrepanze della precisione limitata dei

calcolatori. Analizziamo in dettaglio questi argomenti, partendo dalla dinamica

molecolare, che ha studiato in dettaglio questi argomenti e le problematiche connesse.

2.2.1 Dinamica Molecolare La Dinamica Molecolare (MD) è la scienza che usa i computer per calcolare il moto

delle singole molecole di un solido, liquido o gas nel tempo.

Il moto dei singoli atomi dipende da certe forze che ogni atomo esercita sugli altri. Tali

forze dipendono dall’energia potenziale, la quale a sua volta è funzione della distanza

tra due atomi.

Esistono diversi modelli che definiscono queste funzioni potenziali, ognuna con i suoi

limiti,ma probabilmente la più diffusa è quella del potenziale di Lennard Jones LJ, detto

anche potenziale “6-12”. Ad esempio in [21] si utilizza la formula di LJ per simulare le

proprietà dei materiali elastici. In questa simulazione, che si avvicina concettualmente

molto a quella proposta in questa tesi, è possibile vedere come questo genere di

potenziale, nonostante le sue ipotesi semplificative, conserva a livello macroscopico

delle proprietà tipiche dei materiali reali. Negli esperimenti effettuati infatti è possibile


______________________________________________________________________________________________

______________________________________________________________________________________________


applicare degli sforzi di trazione o taglio su una “provetta” fatta di queste molecole

virtuali, e riscontrare successivamente il fenomeno della strizione con rottura per carico

eccessivo e anche la rottura per sforzo di taglio.

Figura 4. Rottura a trazione e a taglio simulati in MD [21]

2.2.2 Potenziale di Lennard-Jones (LJ) Il potenziale LJ [16] descrive l’interazione tra due molecole prive di carica o tra atomi.

Questo potenziale è poco attrattivo quando le due molecole prive di carica o gli atomi si

avvicinano da distanza, ma è fortemente repulsivo quando la loro distanza diventa

minima. Il potenziale risultante è mostrato in figura. All’equilibrio il paio di atomi o

molecole tende ad andare in separazione lungo il minimo del potenziale di LJ [17].

Il potenziale risultante da queste interazioni di attrazione e repulsione è descritto dalla

seguente equazione:

)(3

−

=

6124)(

rrrVij

σσε

Dove σ e ε sono I parametric specifici di LJ, caratterisciti del materiale. Il termine

121~ r è dominante alle brevi distanze, e modella la repulsione tra atomi quando essi


______________________________________________________________________________________________

______________________________________________________________________________________________


vengono a trovarsi molto vicini. Il termine 61~ r invece è il più forte alle grandi

distanze, e costituisce la parte attrattiva dell’espressione.

La forza di LJ tra le due molecole è quindi data dall’equazione:

)(4

−

−= 7

6

13

12224)(

rrrFij

σσε

Figura 5. Potenziale LJ (rosso) and forza (blu) funzione della distanza (r).

2.2.3 Tensione superficiale Molti dei sistemi analizzati dalla letteratura fanno una distinzione tra forze di volume e

forze di superficie. Questo viene fatto per definire matematicamente il concetto fisico di

tensione superficiale. In letteratura sono presenti diversi modi per capire quali siano le

particelle di confine, ma buona parte dei lavori converge quando si arriva al punto di

come esse debbano essere trattate. Una volta individuate tali particelle infatti nella

maggior parte dei casi si impongono dei connettori governati da equazioni di vario tipo,

per lo più di tipo molla-smorzatore. Nel processo di risoluzione delle intersezioni dei

tool con il materiale, si controlla che le deformazioni prodotte dagli spostamenti delle

particelle non siano superiori ad una certa soglia, oltre la quale il connettore viene

“rotto” e la particella inserita in una lista che raccoglie tutte quelle “da aggiornare”.

Lo smorzatore assolve alla funzione di stabilizzare il sistema, soprattutto perché gli

errori numerici del calcolatore nella risoluzione dell’equazione del connettore diventano

molto grandi all’aumentare del numero di particelle se si vuole mantenere costante il


______________________________________________________________________________________________

______________________________________________________________________________________________


frame rate. Questo modo di considerare le particelle superficiali quindi presenta non

pochi problemi per quanto riguarda le performance in un ambito real-time.

In una prima formulazione di questo lavoro si è pensato di sfruttare le informazioni

derivanti da una struttura dati particolari per capire velocemente quali fossero le

particelle “superficiali”. Se si immagina di inserire tutto il sistema in una struttura dati,

dove ogni particella è dotata di una “valenza” che rappresenta nel nostro caso il numero

massimo di particelle tangenti a quella data (in analogia con la chimica), si può dire che

una particella è di confine quando la sua valenza non è piena. Una volta individuate le

particelle superficiali si sarebbero potute formulare delle equazioni da applicare solo in

quel caso, per rendere più difficile l’operazione di “rottura” della superficie, simulando

di fatto una tensione superficiale. L’analisi dei costi computazionale di questa

operazione, associata al gravoso vincolo del real-time è risultata negativa, e quindi si è

pensato che il potenziale LJ potesse assurgere anche alla funzione di creatore delle forze

superficiali con opportune considerazioni.

2.2.4 Dalle forze agli spostamenti Gran parte della letteratura sul Virtual Sculpturing fa largo uso di equazioni differenziali

del moto per la risoluzione del sistema ad ogni frame, aggiungendo problematiche

enormi di complessità e stabilità numerica.

Se si riflette un poco sulle modalità di modellazione plastica di blocchi di materiale

virtuale argilloso, ci si rende conto di come questa sia in buona parte un insieme di

azioni di tipo statico. In altre parole le azioni dinamiche che si instaurano tra gli

elementi del sistema particellare hanno ben poco di dinamico, ma vengono tutte

modellate per il raggiungimento di un quanto più celere possibile, stato di equilibrio. Ci

si è chiesto a questo punto se fosse possibile tagliare in qualche modo parte delle

equazioni del moto delle particelle, cercando esclusivamente di arrivare direttamente

allo stato di equilibrio che soddisfi dei vincoli fondamentali di un sistema di

modellazione di Virtual Sculpturing.

I vincoli considerati sono stati: la conservazione del volume di materiale su cui si sta

operando e la non compenetrazione del materiale con l’ambiente esterno in generale, e

con i tool di modellazione in particolare. Entrambi questi vincoli non necessitano di

nessuna informazione o risoluzione “dinamica” del sistema.


______________________________________________________________________________________________

______________________________________________________________________________________________


Si è pensato dunque di abbandonare la dinamica e provare a calcolare solo degli

spostamenti funzione della distanza tra due particelle. Questi spostamenti però non

vengono calcolati a caso, ma con una funzione di repulsione e attrazione modellata

esattamente come quella di LJ.

Figura 6. Forze repulsive tra particelle

Nella figura possiamo vedere le particelle 1 e 2 che si intersecano,quindi bisogna

imporre uno spostamento che le faccia allontanare, contemporaneamente però le

particelle 1,2,3 e 4 devono attrarsi tra di loro vicine in un punto che è sostanzialmente il

baricentro del sistema da esse formato. Un esempio così semplice probabilmente non

rende evidente la complessità del problema.

Sono molte le domande alle quali bisogna trovare risposta per poter praticamente

implementare un sistema del genere.

Ad esempio:

§ Come possiamo riuscire a soddisfare tutti i vincoli di non compenetrazione tra

particelle e contemporaneamente “compattare” il blocco di materiale?

§ Come possiamo garantire la costanza del volume di materiale simulato?

§ Come possiamo soprattutto garantire questi vincoli con prestazioni real-time

processando solo 2 particelle per volta?

La risposta a questi quesiti è l’oggetto delle prossime due sezioni, che descrivono la

‘funzione di risposta’ ed il fondamentale approccio iterativo, con relativa metodologia

di aggiornamento.

.


______________________________________________________________________________________________

______________________________________________________________________________________________


2.2.5 Funzione di risposta La funzione di risposta è definita come una funzione da R3 ad R3.

)(5 33:)( ℜ→ℜxf

)(6 x = C1-C2

La variabile in utilizzata è il vettore differenza tra i centri delle due particelle in

considerazione. La funzione associa a questo vettore libero, un vettore applicato in C2.

Questa funzione deve anche assicurarsi che la sua applicazione trasformi uno stato

invalido di configurazione tra due particelle in uno stato valido, cioè uno stato di

equilibrio. Questo stato di equilibrio però deve anche avere la proprietà di essere stabile,

inteso non tanto nell’accezione classica di ritornare alla sua configurazione iniziale

dopo un’azione di disturbo quanto nel senso di non innescare autonomamente fenomeni

vibratori incontrollati che portino all’esplosione del sistema.

L’equazione deve essere sostanzialmente stabilizzare il sistema, indipendentemente

dalla configurazione alla quale viene applicata.

Come vedremo avanti è possibile predisporre le particelle in uno stato di equilibrio

stabile prima di iniziare la simulazione, stando attenti al cosiddetto “impacchettamento”

iniziale, che sia costruito in modo da costituire stato di energia globale già minimo.

La funzione di risposta è definita formalmente come:

)(7x

xx

x

x)(

)(

:)( 33

kf

f

−=

ℜ→ℜ

La funzione k(y) è la caratteristica del materiale. Come è ovvio dalla notazione, essa è

una funzione scalare, che associa ad x un numero reale. In poche parole questa funzione

indica quanto si debbano attrarre o respingere le due particelle. Precedentemente

abbiamo analizzato la funzione di Lennard-Jones come una buona funzione con delle

proprietà interessanti quindi ricordando che:

)(8 ( )

−

−= 7

6

13

12224

rrrFLJ

σσε

Si definisce la k(y) come:


______________________________________________________________________________________________

______________________________________________________________________________________________


)(97137

6

13

12 21224)(

2)()(

ya

ya

yyyLJ

yLJyk

−=

−

=

=

σσε

I primi test eseguiti in due dimensioni ci hanno subito mostrato come questa funzione

faccia raggiungere rapidamente l’equilibrio ad un sistema di particelle sparse. La

variabilità della risposta che si riesce ad avere variando a1 ed a2 è troppo poco

controllabile, minime variazioni di uno dei due parametri causano esplosioni o

implosioni anomale del materiale. Si è così deciso di discretizzare questa funzione e di

approssimarla collegando i punti con delle rette. Successivamente è stato reso possibile

modificare questi punti graficamente, alterando di fatto la funzione di risposta. Nei test

fatti si è visto che variando di poco opportunamente questi punti si riusciva a simulare

materiali con proprietà leggermente differenti tra loro. Ad esempio permettendo una

piccola compenetrazione tra particelle, spostando il lato sinistro della funzione di LJ

poco più a sinistra si notava che il materiale assorbiva le deformazioni imposte

“comprimendosi” seppur di poco.

Bisogna ora chiarire perché k(x) è stata definita come LJ(x)/2.

Si è visto che per portare il sistema in equilibrio stabile entrambe le particelle che

definiscono il vettore spostamento devono essere allontanate o avvicinate di una

quantità uguale alla metà della funzione LJ per fare in modo che la distanza alla fine

dell’applicazione della LJ sia esattamente quella prevista. Si sarebbe potuto anche

suddividere lo spostamento in percentuali diverse tra le due particelle, per esempio

spostando una del 20% della LJ e l’altra dell’80% della LJ, ma si è notato che questo è

motivo di forte instabilità negli step successivi, poiché non garantisce l’equilibrio del

tipo di “relaxation” spiegato nel prossimo paragrafo.


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 7. Forze repulsive tra particelle

Da notare che la funzione di risposta ha il suo zero in 2r, che ovviamente è il caso nel

quale le due palline sono tangenti. Come vedremo nei dettagli implementativi 6.2 per

evitare che errori di precisione numerica portino due particelle a vibrare continuamente

al di sopra e al disotto di questo zero, saranno imposte delle tolleranze al di sotto delle

quali la funzione non verrà applicata.

Praticamente l’effetto di questa funzione è quello di portare ad equilibrio dopo un certo

numero di iterazioni, il sistema formato da due particelle. La figura di seguito

esemplifica tale comportamento.

Figura 8. Effetti del potenziale LJ sulle particelle

La stabilità del sistema non dipende esclusivamente dalla funzione di risposta. Vedremo

più avanti che un ruolo molto importante è tenuto dagli impacchettamenti iniziali, che

minimizzano l’energia globale del materiale.


______________________________________________________________________________________________

______________________________________________________________________________________________


2.3 Risoluzione iterativa dei vincoli. Relaxation Data la natura della simulazione, nella quale due particelle vengono prese e la loro

interazione risolta per poi passare ad altre particelle, si capisce come non sia sufficiente

all’equilibrio eseguire tale procedura per tutte le possibili coppie di particelle.

Quando si risolve un vincolo tra due sferette, molto probabilmente le loro nuove

posizioni avranno invalidato il rapporto con altre particelle, con le quali magari sono

state già risolte. L’idea è quindi quella di ripetere iterativamente questa procedura

“sperando” di convergere ad una soluzione di equilibrio globale generato a partire

dall’equilibrio locale di ogni coppia di due particelle. Sebbene questa procedura possa

sembrare alquanto empirica, la letteratura ci mostra che essa ha basi matematiche molto

forti e giustificate.

In [23] è presentato un metodo chiamato di “relaxation” (o di Jacobi o delle iterazioni di

Gauss-Seidele in base a come lo si fa esattamente). In questo metodo si prova a

risolvere un sistema con molte incognite, cercando di soddisfare consecutivamente i vari

vincoli locali e iterando. Se le condizioni sono giuste, questa procedura convergerà ad

una configurazione globale che soddisfa tutti i vincoli contemporaneamente.

2.3.1 Campi di utilizzo del “relaxation” L’approccio relaxation è utile in tutte quelle situazioni nelle quali diversi vincoli tra di

loro inter-dipendenti devono essere soddisfati simultaneamente. Il numero di iterazioni

necessarie dipende dalla natura del sistema fisico simulato. Questo sistema può anche

essere reso adattivo, misurando lo scarto dall’ultima iterazione. Se si ferma l’iterazione

prematuramente, il risultato potrebbe non essere propriamente valido, ma poiché la

nostra risoluzione di vincoli partirà da uno stato che è “maggiormente equilibrato”

rispetto al precedente, nelle iterazioni dei frame successivi della simulazione, continuerà

a convergere partendo da tale stato.

I metodi di relaxation sono in genere la norma per buona parte delle simulazioni

nell’ambito della dinamica molecolare, in quanto non è possibile pretendere di risolvere

simultaneamente con esattezza tutte le equazioni che legano migliaia o milioni di

molecole nella simulazione di un gas, un liquido o un solido. Questi metodi sono inoltre

utilizzati in computer grafica per la simulazione di tessuti, a volte in accoppiata con i

metodi di integrazione Verlet che si sposano perfettamente con i sistemi particellari a

vincoli fissi tra particelle.


______________________________________________________________________________________________

______________________________________________________________________________________________


2.3.2 Metodi verlet I metodi Verlet danno una formulazione delle equazioni del moto dove la velocità è

espressa implicitamente dallo scarto di posizione tra una iterazione e la precedente.

In sostanza per capire l’idea alla base del metodo Verlet, si può considerare il semplice

schema di integrazione di Eulero:

)(10 ,''

tt

∆⋅+=∆⋅+=

avvvxx

dove ∆t è il time step, e a è l’accelerazione calcolata usando la legge di Newton f=ma

(dove f è la risultante di tutte le forze esterne e interne al sistema applicate alla

particella). Invece di “ricordare” la velocità e la posizione di ogni particella all’istante

considerato, si può utilizzare la posizione corrente x e la precedente posizione x*. Se si

tiene fisso il passo di integrazione, si ha allora:

)(11 .2'

*

2*

xxaxxx

=

∆⋅+−= t

Questo è quello che viene chiamato “Metodo di integrazione Verlet”, presente in molti

lavori sperimentali sulle proprietà molecolari dei fluidi usando i potenziali di Lennard-

Jones (vedere ad esempio il documento originale di Verlet [25] per ulteriori

informazioni). Questo metodo è decisamente stabile, poiché la velocità viene data

implicitamente, e di conseguenza è più difficile per la posizione e la velocità perdere di

sincronizzazione. Questo sistema si basa sul fatto che 2x-x*=x+(x-x*) e x-x* è

un’approssimazione della velocità corrente (in realtà sarebbe la distanza percorsa

dall’ultimo passo di integrazione). I risultati non sono sempre molto accurati, poiché è

presente una intrinseca dissipazione di energia, ma lo schema di integrazione è veloce e

stabile. Se si diminuisce il coefficiente 2, ad esempio ponendolo ad esempio a 2-ε con

ε>0 e ε<0,1 si può introdurre un leggero smorzamento sul sistema.

Questo metodo è stato anche usato nel campo dell’entertainment per simulare corpi

rigidi vincolati utilizzati come bambole di pezza (rag-doll) per i personaggi di un noto

videogioco [24].

Questo lavoro non ha utilizzato i metodi Verlet per la scelta iniziale di trascurare

l’aspetto dinamico della situazione, ma sono in corso esperimenti per valutare se sia

possibile abilitarlo o disabilitarlo a quando se ne ha necessità.


______________________________________________________________________________________________

______________________________________________________________________________________________


2.4 Impacchettamento Se si considera il potenziale di Lennard-Jones, definito in termini di distanza euclidea da

una partice, è una funzione di energia potenziale simmetrica nello spazio. Questo vuol

dire che prese due particelle i e j, l’insieme degli stati a minima energia (che nel nostro

sistema sono posizioni) per una particella j relativa a i p il luogo dei punto che forma

una sfera centrata sul centro della particella i. Quando le forze esterne sono trascurabili,

le particelle si organizzano in strutture fortemente impacchettate, per minimizzare la

loro energia totale. E’ stata proprio questa l’intuizione che ha permesso al simulatore di

particelle realizzato di partire già da una condizione di equilibrio stabile, senza aspettare

lunghi periodi prima di potervi pervenire.

Per funzioni energetiche potenziali che presentano simmetria circolare in 2D, come il

potenziale LJ, le particelle si dispongono in gruppi esagonali. In 3D invece si

dispongono lungo livelli successivi di piani 2D, ognuno dei quali presenta un

ordinamento esagonale, rendendo di fatto il potenziale LJ ottima per la modellazione di

volumi di materiale.

2.4.1 Densità Il denso impacchettamento delle particelle pone immediatamente una questione, e cioè

quante particelle saranno necessarie per riempire un dato volume La risposta a questo

quesito è di grande aiuto per capire l’efficienza di una tecnica di ricerca delle particelle

“vicine” ad una data. Possiamo riformulare la domanda chiedendoci quante sfere di

eguale dimensione saranno necessarie per riempire un dato volume.

Il fattore di impacchettamento volumetrico è definito come il rapporto del numero di

atomi presenti in un volume unitario. Per un dato impacchettamento questo è

equivalente al rapporto del volume di una sfera alla regione associata del diagramma di

Voronoi tridimensionale, dove però il diagramma di Voronoi è calcolato sui centri delle

sfere.Consultando [27] e alcuni libri di tecnologia dei materiali si è potuto calcolare la

densità degli impacchettamento scelti


______________________________________________________________________________________________

______________________________________________________________________________________________


2.4.2 Impacchettamenti CFC e Esagonale

Figura 9. Varie tipologie di impacchettamento [27]

Per capire come funziona questo particolare tipo di impacchettamento, si immagini di

disporre 6 sfere come in figura per formare un triangolo equilatero, e successivamente si

piazzi un’altra sfera sopra questo livello per creare una primaide triangolare. Ora se

immaginiamo di creare un altro gruppo di sette sfere, e piazziamo le due piramidi in

modo tale da farle affacciare in direzioni opposte, ci rendiamo conto di aver costruito un

cubo. Collegando i centri di queste 14 sfere si ottiene quella che viene chiamata stella

octangula.

Figura 10. Costruzione geometrica dell’impacchettamento

Considerando il cubo definito dai centri delle 14 sfere, questa cella “unitaria” contiene

otto ottavi di sfere (uno per ogni verice del cubo) e sei emisferi. Il volume totale di sfere

nella cella unitaria è quindi:

)(12 333/ 3

163

443

4216

818 rrrV cellasfere π

ππ=⋅=

⋅+⋅=

La diagonale della faccia è 4r cosi ogni lato è r22 . Il volume di tale cella è quindi:


______________________________________________________________________________________________

______________________________________________________________________________________________


)(13 33 216)22( rrVcella ==

E quindi la densità di impacchettamento è:

)(1423216

316

3

2ππ

η ==r

rcfc =0,74048

La densità del CFC è uguale a quella dell’impacchettamento esagonale.

Figura 11. Impacchettamento esagonale centrato

Si dimostra che questi impacchettamenti sono quelli con la massima densità

teoricamente possibile.

2.5 Strutture di partizionamento spaziale Come abbiamo precedentemente visto, la funzione di LJ che guida l’interazione tra le

particelle va quasi a zero oltre un certo range. Per questo motivo si è pensato fosse

completamente inutile formulare un approccio brute force, che eseguisse ad ogni

iterazione l’intero step di relaxation su tutte le possibili coppie di particelle. A questo

proposito si è pensato di inserire il sistema in una struttura di partizionamento spaziale

che fosse la più efficiente possibile, considerando il tipo di primitive che stiamo

utilizzando, e cioè delle sfere di raggio unitario. L’aspetto qui presentato è quello che ha

subito le maggiori evoluzioni nella fase di progettazione del sistema stesso. Sono state

provate diverse soluzioni che sembravano avere vantaggi teorici di performance ma che

una volta implementate hanno mostrato subito grandi problemi. In generale però tutte le

strutture esaminate appartengono alla categoria di strutture di partizionamento

gerarchiche dove lo spazio è ricorsivamente partizionato in zone più piccole che sono

contenute all’interno delle zone “generatrici”, le quali a loro volta sono contenute in

altre zone ancora. Tali “zone” vengono chiamate con terminologia più tecnica Bounding

Volumes (BV) [28] e sono i nodi della struttura ad albero. Questo genere di strutture


______________________________________________________________________________________________

______________________________________________________________________________________________


permette di scartare con pochi test geometrici, grosse porzioni di spazio dove non

interessa eseguire l’aggiornamento.

Tutte queste strutture sono inoltre conservative. Questo termine significa che alla

domanda “quali sono le particelle che si trovano in questa determinato volume di test?”

verranno sicuramente indicate tutte quelle che lo sono realmente, ma potranno essere

indicate alcune particelle che sono esterne al volume di test.

I parametri di performance fondamentali che sono stati presi in considerazione sono

quindi la complessità computazionale, il numero di sfere riportate rapportato al numero

di sfere che dovrebbero essere riportate e la memoria occupata dalla struttura, molto

importante nel nostro caso di simulazione di milioni di particelle. Tutte questi vincoli

devono essere rispettati tenendo bene a mente che non si sta trattando di oggetti statici,

ma di migliaia di elementi dinamici che si muovono ad ogni iterazione.

2.5.1 Introduzione alle strutture di partizionamento In generale, una struttura dati di partizionamento organizza della geometria in uno

spazio n-dimesionale. In generale ci si limita a casi bi- e tri-dimensionali, ma i concetti

possono facilmente essere generalizzati per dimensioni arbitrarie. Queste strutture dati

possono essere usate per accelerare dei test che vogliono sapere se determinate

geometrie si sovrappongono. Questi test, chiamati in genere “query”, possono essere

usati in una gran varietà di operazioni, come ad esempio algoritmi di culling oppure nei

test di intersezione e raytracing o per il rilevamento delle collisioni.

L’organizzazioni di strutture dati di partizionamento spaziali è genericamente

gerarchica. Questo significa, in parole povere, che il livello superiore racchiude il livello

sotto di esso, che racchiude a sua volta un ‘altro livello e così via. Quindi questa

struttura è di tipo ricorsivo, come dicono gli informatici. La ragione principale per usare

una gerarchia è che vari tipi di query vengono fortemente accelerate, da complessità di

tipo O(n) a O(log n). E’ da sottolineare tra l’altro che la costruzione della maggior parte

delle strutture dati di partizionamento è abbastanza costosa . equindi viene fatta prima di

cominciare la simulazione vera e propria.

Alcuni tipi di strutture di partizionamento sono le Bounding Volume Hierarchies (BVH

s) , gli alberi Binary Space Partition (BSP) e gli octrees.

Gli alberi BSP e gli octree sono strutture dati basate sulla suddivisione spaziale. Questo

significa che l’intero spazio della scena è suddiviso e catturato nella struttura dati. Per

esempio, l’unione dello spazio di tutti i nodi “finali” della struttura (in inglese leaf


______________________________________________________________________________________________

______________________________________________________________________________________________


nodes, cioè i nodi privi di figli) equivale all’intero spazio della scena. Questi nodi finali

non si sovrappongono mai, ad eccezione dei loose octree che hanno delle proprietà

molto interessanti però come analizzeremo in seguito. Le varianti dei BSP partizionano

lo spazio in modo irregolare, a differenza degli octree e derivati, in cui lo spazio è

suddiviso uniformemente. Anche se a volte può sembrare restrittiva, questa uniformità

può essere fonte di efficienza computazionale.

Una Bounding Volume Hierarchy invece,non è propriamente classificabile come una

struttura di partizionamento spaziale. Piuttosto racchiude regioni dello spazio che

circondano oggetti geometrici e quindi le BVHs non necessariamente racchiudono tutto

lo spazio.

2.5.2 Bounding Volume Hierarchies Un Bounding volume (BV) è un volume che racchiude un insieme di oggetti. Il punto di

forza di un BV consiste nel fatto che esso sia una forma geometricamente più semplice

dell’oggetto contenuto, in modo tale che fare test su un BV può essere fatto più

velocemente che usando gli ogetti da soli. Esempi di BV sono sfere, parallelepipedi

allineati agli assi (Axis-Aligned Bounding Boxes AABBs), parallelepipedi orientati

(Oriented Bounding Boxes OBBs) e k-DOP (k Discretely Oriented Polytope). Una BV è

comunemente usata per velocizzare le query di intersezione.

Quando si usa una BVH, la scena viene organizzata in un albero gerarchico, che

consiste in un nodo radice, nodi interni e nodi finali (leaves). Questi nodi, come ogni

albero, sono legati da relazioni di “parentela” del tipo padre/figlio. Il nodo radice è

quello che non ha padre, ed è anche un nodo interno, tranne nel caso sia l’unico nodo

nell’albero. Ogni nodo, e anche i nodi finali, nell’albero hanno un BV associato che

racchiude la geometria dell’intero sotto-albero che diparte da esso. Questo significa che

il nodo radice ha un VB che racchiude l’intera scena. Un esempio è mostrato in figura.


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 12. Una Bounding Volume Hierarchy e la sua rappresentazione ad albero

La struttura sottostante della BVH è un albero, e nell’ambito dell’informatica la

letteratura sulle strutture dati ad albero è vastissima. Qui verranno presentati solo alcuni

risultati notevoli. Per maggiori informazioni si consulti ad esempio il libro Introduction

to Algorithms [29]

Si consideri un albero k-nario, cioè un albero dove ogni nodo interno ha k figli. In

albero con solo un nodo (la radice) è detto a livello 0. Un nodo finale del nodo radice è

al livello 1 e così via. Un albero bilanciato è un albero dove tutti i nodi sono al livello h

o h-1. In generare il livell, h, di un albero bilanciato è logk n, dove n è il numero totale

di nodi (interni e finali) dell’albero.

Un albero si dice completo se tutti i nodi finali sono allo stesso livello h. Alcune

proprietà degli alberi completi sono esposte di seguito. Il numero totale di nodi può

essere calcolato come somma geometrica di:

)(151

1...1

110−

−=+++=

+−

kkkkkkn

hhh

Di conseguenza il numero di nodi finali, l è l=kh, e il numero di nodi interni i è i = n-l =

11

−−

kk h

. Assumendo che sia noto saolametne il numero dei nodi finali, allora il numero

totale dei nodi è n= i+l = 11

−−

kkl . Per un albero binario, dove k=2, questo ci da n = 2l1.

Da notare che un k più alto genera un albero con un minor numero di livelli, che

significa si impiegherà meno tempo per percorrere l’albero,ma richiederà un maggior

carico di lavoro ad ogni nodo. L’albero binario è in genere la scelta più semplice, e da


______________________________________________________________________________________________

______________________________________________________________________________________________


buone performance. Ad ogni modo per applicazioni specifiche si usano altri valori di k,

in genere 2,4 e 8.

Il BVH sono eccellenti per l’accelerazione di vari tipi di query e sono potenti anche in

scene fortemente dinamiche. Se per esempio un oggetto contenuto in un VB si è mosso,

basta controllare se esso è ancora contenuto nel suo BV padre. Se è così allora il VB è

ancora valido, altrimenti il nodo che contiene quell’oggetto è rimosso, ed il BV del

padre viene ricalcolato. Il nodo rimosso verrà successivamente reinserito ricorsivamente

nell’albero partendo dal nodo radice. Un altro metodo è di far crescere il BV del nodo

padre per mantenere i nodi figli al suo interno fino ad una certa soglia. Ad ogni modo,

con entrambi questi metodi l’albero può diventare non bilanciato e inefficiente man

mano che avvengono queste rimozioni ed inserzioni.

2.5.3 Sphere Tree Gli sphere tree, o più semplicemente alberi di sfere, sono una BVH il cui Bounding

Volume è una sfera, come il nome stesso suggerisce.

In letteratura sono disponibili i dettagli implementativi [30] e le possibilità ma anche le

limitazioni di una struttura di questo tipo. Inizialmente si è provato ad implementare

questa struttura poiché sembrava avere tutti i requisiti posti in fase di progetto.

In primo luogo, considerando una sfera come Bounding Volume significa in generale

sovrastimare l’effettivo volume occupato da oggetti lunghi e stretti (vedere la Figura

11).

Figura 13. Una Bounding Sphere ‘obesa’ su un oggetto lungo e stretto


______________________________________________________________________________________________

______________________________________________________________________________________________


Questo non era un problema per il sistema ideato, dato che i nostri “oggetti” sono sfere

di raggio unitario, cioè quanto di meglio si possa offrire in pasto ad un albero di sfere.

Innanzitutto cerchiamo di capire come funziona lo sphere tree di fronte ad una

moltitudine di oggetti in movimento.

Quando una particella cambia posizione, si controlla se essa sia ancora contenuta nella

sfera del nodo padre. In questo caso l’algoritmo si ferma subito, poiché la posizione del

nodo associato al BV di quell’oggetto nell’albero è ancora valida. Per gli oggetti statici

o che non si sono mossi (e nel nostro sistema è una condizione frequente), questa è

addirittura l’unica operazione da eseguire per mantenere aggiornato l’albero, al di là

dell’inserimento iniziale.

Quando uno spostamento invece,fa in modo che il BV dell’oggetto in questione

oltrepassi quello del suo BV padre,allora il nodo finale viene rimosso e reinserito

nell’albero partendo dal nodo radice. In realtà per efficienze di implementazione,

quando i nodi figli vengono staccati dai padri, vanno a finire in una coda di

reintegrazione di tipo FIFO (First In First Out). Il padre è invece inserito in una coda di

ricalcolo FIFO per mantenere bilanciato l’albero. Ad ogni step di simulazione, prima di

eseguire le query, vengono svuotate queste code, reintegrando e ricalcolando il

necessario in un unico loop.

Un problema delle strutture di tipo quadtree od octree, come vedremo, consiste nel fatto

che è possibile per un singolo nodo finale avere molteplici nodi padre. Se un oggetto

attraversa un confine di nodo, necessiterà di essere rappresentato in tutti i nodi che ha

“toccato”. Ebbene gli alberi di sfere non hanno questo problema, poiché nessun nodo

finale può essere contenuto in più di una sfera padre. Qualunque proprietà di un nodo

padre, ricade esattamente su tutti i nodi figli. Se ad esempio scopriamo che un

determinato nodo padre non interseca un nostro volume di test, possiamo stare sicuri

che anche i nodi figli del nodo padre non intersecheranno il volume.

In figura possiamo vedere un esempio di questa tecnica in azione in un caso

bidimensionale. Le circonferenze grandi e medie dal bordo sottile sono i nodi interni, il

cui unico scopo è di contenere altri nodi. Le circonferenze piccole sono i nodi finali e

rappresentano fisicamente gli oggetti di cui si sta tenendo traccia in questa struttura di

partizionamento. La circonferenza di bordo spesso è quella di test. Come si può vedere

si è prima cercato di capire quali fossero le circonferenze “grossolane”, di livello più

basso,toccate da quella di test. Successivamente si è provveduto a vedere quali nodi figli


______________________________________________________________________________________________

______________________________________________________________________________________________


di questi grandi nodi intersecassero effettivamente la circonferenza di test, andando così

via via fino ai nodi finali.

Figura 14. La costruzione di uno Sphere Tree

Una volta implementata questa struttura dati, si è avuta conferma delle buone

prestazioni della struttura. L’unica pecca è forse l’elevato numero di test di intersezione,

dovuta alla conservativa in senso stretto. Questa struttura infatti riporta solo i nodi che

sono esattamente contenuti nel volume di test, piuttosto che riportare blocchi di nodi

che “potrebbero” essere contenuti nel volume di test. Questo effetto potrebbe

ovviamente ottenersi fermandosi a livelli fissati dell’albero invece di andare fino in

profondità, ma in tal caso la forma “sovradimensionata” della sfera, rispetto al

contenuto ritorna un numero di possibili oggetti intersecanti troppo alto.

Si è provveduto quindi ad esplorare altre strutture di partizionamento.

2.5.4 Octree Un octree può essere definito come un albero BSP allineato agi assi.. Un cubo viene

diviso simultaneamente lungo tutti e tre gli assi, con dei piani passanti per il centro del

cubo. Questa procedura crea otto nuovi cubi, da cui il nome octree. Questa costruzione

rendere la struttura regolare, ed alcune query possono diventare più efficienti grazie a

questo fatto.

Un octree è costruito racchiudendo l’intera scena nel più piccolo cubo allineato agli assi.

Il resto della procedura è naturalmente ricorsiva, e si ferma quando si raggiunge un

determinato criterio di stop. Questi criteri possono essere determinati da un numero

massimo di livelli di ricorsione, oppure dal fatto che si vada sotto una determinata

soglia di numero di primitive in un cubo. Se si ferma la ricorsione,l’algoritmo associa le

primitive al cubo e termina la ricorsione, altrimenti suddivide il cubo lungo i suoi assi


______________________________________________________________________________________________

______________________________________________________________________________________________


coordinati, usando tre piani,e quindi formando otto cubi di eguale dimensione. Ogni

cubo è testato e possibilmente suddiviso ancora in 2x2x2 cubi più piccoli.

Figura 15. La costruzione di un Octree

In figura si può vedere come funziona questa struttura in due dimensioni, con il nome di

quadtree in questo caso caso.

L’octree presenta notevoli vantaggi di velocità in tutte le operazioni, sia di

aggiornamento sia di query di intersezioni e sembra il perfetto candidato per

partizionare lo spazio della simulazione sin qui ideata. Purtroppo però, come si può

vedere dalla figura, esistono oggetti completamente contenuti all’interno di un nodo

finale ed oggetti, come la mela o la stella ad esempio, che appartengono a più di un

nodo.

Una soluzione a questo problema potrebbe esser quella di suddividere gli oggetti in

questione, introducendo quindi un numero maggiore di primitive. Questa soluzione però

non è praticabile per il sistema particellare ideato, in quanto la primitiva base, la sfera

unitaria, è per noi unità indivisibile.

Un’altra soluzione potrebbe invece essere quella di inserire praticamente questi oggetti

in più di un nodo, ma ciò renderebbe difficoltose le operazioni di aggiornamento, in

quanto si dovrebbe tener traccia di 8 ‘possibili’ nodi padre per ogni particella. Questo

causerebbe problematiche di memoria e di cache miss in CPU enormi, poiché si

dovrebbero inserire svariati test condizionali da eseguire su ogni particella.

Pur sembrando questo un problema insormontabile, nel paragrafo successivo

scopriremo che una semplice modifica di questa struttura riuscirà a soddisfare tutti i

vincoli e le necessità del sistema particellare considerato.


______________________________________________________________________________________________

______________________________________________________________________________________________


2.5.5 Loose Octree L’idea base dei loose octree [14] è la stessa degli octree ordinari, ma la scelta della

dimensione del singolo cubo viene trattata differentemente. Se la dimensione del lato di

un normale cubo è l, allora si usa kl con k > 1.

Figura 16. Comparazione tra un Octree standard ed un Loose Octree

Come si può vedere dalla figura, l’octree standard avrebbe dovuto inserire l’oggetto in

nei due nodi “toccati”. Se invece allarghiamo ogni nodo di un fattore k opportuno,

vediamo che l’oggetto è completamente contenuto in almeno un nodo, anche se

ovviamente potrebbe anche essere contenuto in più di uno. In questo caso vediamo che

l’oggetto è completamente contenuto nel nodo in alto a sinistra.

E’ importante notare che i centri dei cubi sono gli stessi. Essendo ogni oggetto inserito

in un solo nodo, l’eliminazione dall’octree è semplicissima. Alcuni vantaggi teorici ci

sono usando k = 2. In primo luogo l’inserzione e l’eliminazione delgli oggetti è O(1),

cioè costante. Questo perché la conoscenza della dimensione dell’oggetto significa

conoscere immediatamente il livello dell’octree in cui può essere inserito, essendo

completamente contenuto in un cubo loos. Il centroide dell’oggetto determina in quale

cubo loose dell’octree esso debba essere messo. Vedremo nella Appendice A al

paragrafo 6.4 dettagli di codice di queste operazioni. In poche parole si può sapere in

quale “ottante” è andata a finire la nostra particella,esaminando i segni del vettore

differenza tra il centro del nodo padre e la particella stessa.


______________________________________________________________________________________________

______________________________________________________________________________________________


Per via di queste proprietà dunque, questa struttura si presta molto bene a gestire oggetti

dinamici, a costo ovviamente di una minima efficienza e la perdita di un rigoroso ordine

quando si cammina l’albero. Inoltre se un oggetto si muove di poco da una iterazione

alla successiva, il precedente cubo è in generale valido. Di conseguenza solo una

frazione di oggetti nell’octree ha bisogno di essere aggiornata nel corso di un frame.

2.6 Politica degli aggiornamenti Quando il numero di particelle simulate cresce molto, diventa impossibile aggiornarle

tutte ad ogni frame in tempo reale. Di conseguenza bisogna istituire una politica degli

aggiornamenti, per aggiornare solo una porzione dell’intero volume. Gli step di

simulazione non vengono quindi eseguiti ad ogni frame su tutte le particelle ma solo su

alcune di esse. Queste considerazioni scaturiscono dall’osservazione del comportamento

dei materiali clay-like: le maggiori deformazioni infatti sono localizzate “vicino” il tool.

Inoltre nell’arco di pochi frames, gli strumenti spesso si spostano molto poco, se non

rimangono addirittura fermi. La politica degli aggiornamenti è molto importante quindi

per questa simulazione, ed è stata studiata per sfruttare le seguenti due considerazioni.

Immaginiamo una situazione nella quale l’utente sta avvicinando alla superficie del

materiale con uno strumento a forma di cubo.

Figura 17. Fasi successive dello schema di propagazione della deformazione

Lo schema di aggiornamento dovrebbe prima cercare di aggiornare le particelle “vicine”

al tool, e poi “propagare” questi spostamenti alle altre particelle nella simulazione. Tra

le varie strategie ideate quella che è risultata più efficace è stata quella alla “Monte-

Carlo”.


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 18. Generazione di sample points (Monte Carlo) [31]

Il nome di questa tecnica prende spunto dalla teoria del calcolo delle probabilità, dove si

usa il metodo di Monte-Carlo sfruttando il risultato notevole:

)(16 ∫ ∫ ∑=

≈=N

i xpxf

Ndxxp

xpxfdxxf

1 )()(1)(

)()()(

Con f(x) funzione generica e p(x) funzione densità di probabilità della variabile in

considerazione. Questa tecnica viene utilizzate per eseguire una integrazione di una

funzione definita su una sfera, prendendo dei punti non proprio completamente “a

caso”, ma scalandoli di un certo fattore che da un peso minore ai punti che “più

probabilmente” saranno scelti. Per scegliere questi punti casuali di integrazione, si

genera una moltitudine di coppie di numeri casuali indipendenti, denotati come numeri

casuali canonici xξ ed yξ entrambi appartenenti all’intervallo [0,1[ , che possiamo

vedere nella figura [ref di figura] distribuiti in un quadrato. Successivamente questi

numeri vengono “mappati” in coordinate sferiche usando la trasformazione:

)(17 ),()2),1arccos(2( ϕϑπξξ →− yx

Ebbene in questo lavoro viene utilizzata questa distribuzione per generare dei punti

casuali nei quali eseguire la nostra operazione di “integrazione” delle equazioni del

moto delle singole particelle. Ovviamente non è una integrazione in senso stretto,

matematico in quanto, come già detto precedentemente, sono stati abbandonati gli

aspetti dinamici del problema. Utilizzando questo convergenza metodo però si è visto

che il sistema tende rapidamente all’equilibrio, uniformemente in tutte le direzioni.

Il sistema dunque funziona in questo modo: gli aggiornamenti seguono la superficie di

una sfera di base centrata sul tool. Questa sfera ha un raggio chiamato ‘raggio di


______________________________________________________________________________________________

______________________________________________________________________________________________


espansione’ che aumenta fino ad arrivare ad un massimo per tornare poi a zero. Tale

parametro è stato collegato al numero di iterazioni che viene eseguita ad un determinato

livello. Questo vuol dire che le zone che sono più vicine al tool vengono aggiornate più

frequentemente delle zone lontane. Il ragionamento ha perfettamente senso se si

considera quanto detto precedentemente, e cioè che le deformazioni sono inizialmente

‘localizzate’ vicino al tool, e successivamente si propagano al resto del materiale.

La propagazione ha sostanzialmente il compito di conservare del volume totale del

materiale simulato. Il risultato finale è abbastanza piacevole, l’esperienza interattiva

dell’utente risulta notevolmente migliorata.

2.7 Strumenti di modellazione Gli strumenti utilizzati per modellare il materiale deformabile che è stato messo a punto

in questo lavoro possono essere di due tipi: a sfere o poligonali.

2.7.1 Strumenti composti di sfere Gli strumenti a sfere vengono creati approssimando la forma del ‘tool’ in

considerazione con delle particelle. Tali particelle sono in tutto e per tutto uguali a

quelle della simulazione, a parte un flag che le marca. Il risolutore iterativo provvede a

considerare le particelle del tool come costitute da una massa infinita durante la

risoluzione del vincolo. Questa fa sì che esse siano capaci solo di “spingere” altre

particelle,ma non possano mai essere “spinte”.

Figura 19. Un tool bidimensionale composto da sfere

Questa formulazione dei tool è molto vantaggiosa in quanto evita di eseguire pesanti

test di intersezione sfera/poligono, sostituendoli con semplici test sfera unitaria/sfera

unitaria, che è il più semplice che si possa immaginare. Tra l’altro, poiché questo

sistema riesce a simulare un discreto numero di particelle, l’approssimazione del tool


______________________________________________________________________________________________

______________________________________________________________________________________________


poligonale è “ben fatta”, nel senso che le impronte lasciate da questo tipo di tool sono

quasi identiche a quelle dei tool poligonali. Rimane il problema di capire come si possa

approssimare una forma arbitraria con un insieme di sfere.

Innanzitutto se la forma che bisogna approssimare è una forma convessa chiusa, nulla è

più facile poiché basta prendere il bounding volume che racchiude la forma, inserirlo in

pezzo di materiale già impacchettato e scartare tutte le particelle che non sono contenute

all’interno. Il test di contenimento all’interno per una forma convessa chiusa il test è

semplicissimo,ma si utilizzerà quello descritto di seguito,valido per le forme concave

chiuse e quindi anche per quelle convesse chiuse.

Per le forme concave chiuse si potrebbero chiamare in causa diversi teoremi della

geometria computazionale, ma una soluzione semplice potrebbe essere quella di

utilizzare un algoritmo di flood-fill. In poche parole si parte da punto che si è certi sia

all’interno della forma in questione e si creano altre sfere nelle tre direzioni cardinali,

nonché nelle loro possibili combinazioni. Prima di creare fisicamente una nuova sfera

però, si valuta se il segmento che va dal centro della sfera corrente, al centro della nuova

sfera da creare oltrepassi un qualunque poligono dell’oggetto da approssimare.

Figura 20. Flood-Fill per riempire di particelle una forma chiusa

I test di intersezione raggio-poligono sono anche trattati nei paragrafi riguardante i

dettagli implementativi.


______________________________________________________________________________________________

______________________________________________________________________________________________


Per le forme aperte purtroppo, data l’impossibilità di stabilire quale sia il “dentro” e

quale sia il “fuori” non è possibile usare questi algoritmi. Si pensi al nodo di Moebius

per esempio per il quale non è nemmeno possibile definire il concetto di normale

uscente. Ebbene in questi casi l’unica possibilità potrebbe essere quella di fare un

sampling dei triangoli del tool, ottenendo così dei punti sui quali poter piazzare le sfere.

Con questo schema però risulta difficile garantire il perfetto impacchettamento delle

particelle stesse. In realtà questo non dovrebbe essere un grosso problema per le

particelle del tool in quanto esse non possono essere modificate dalla funzione di

risposta, quindi non si hanno problemi di “vibrazione” attorno allo zero dell’energia

poiché non si parte dallo stato di energia minimo. Il vero problema però è il fatto che

questi tool siano “vuoti” al loro interno, e questo è un problema perché per la natura

della simulazione, alcune particelle potrebbero “oltrepassare” la sottile superficie del

tool e venir inglobate dal tool stesso, con un effetto totalmente irrealistico.

Come ultimo aspetto si vuole sottolineare che è possibile inserire in questo contesto

strumenti di tipo force feedback. Basterebbe accumulare la quantità di energia

“assorbita” dal tool (semplicemente calcolando il lavoro virtuale dello spostamento

impedito) e passarlo al device di Force Feedback. Allo stesso tempo è possibile

ovviamente fare in modo che il tool sia comandato direttamente con le mani da un

sistema di VR, e non attraverso l’innaturale mouse e tastiera.

2.7.2 Strumenti poligonali Gli strumenti poligonali, come indica il nome, sono delle mesh arbitrarie che possono

essere modellate in un qualunque programma CAD o di grafica 3D in generale.

Questi strumenti hanno il vantaggio di rappresentare con grande esattezza la forma del

tool che si vuole utilizzare nella simulazione, ma allo stesso tempo lo svantaggio di

dover eseguire dei pesanti test di intersezione con distanza di penetrazione di sfera con

poligono. Da alcuni test svolti si è avuta la certezza che questa rappresentazione dei tool

non da apprezzabili vantaggi nella modellazione di Virtual Sculpturing e allo stesso

tempo causa una notevole perdita di performance. Per tali motivi questa soluzione è

stata scartata.

2.7.3 Performance La dimensione del tool e la forma è molto importante per quanto riguarda le

performance. Il raggio di espansione, descritto nel paragrafo sulla politica degli


______________________________________________________________________________________________

______________________________________________________________________________________________


aggiornamenti è proporzionale alla dimensione più grande dello strumento. Gli oggetti

stretti ed allungati hanno prestazioni peggiori di quelli più ‘cubici’, a parità di nmero di

particelle. Ovviamente seguendo questo ragionamento, gli strumenti piccoli sono più

leggeri computazionalmente di quelli grossi. La misura del tempo di simulazione è

molto difficoltosa, poiché dipende dalla forma del tool, dalla dimensione ma anche dalla

velocità con la quale lo strumento stesso viene mosso dall’utente. Una maggior velocità

di movimento dello strumento causa un numero maggiore di particelle toccate, e quindi

un maggior numero di cicli della CPU che devono essere impiegati per aggiornarle. I

tempi di simulazione per diversi numeri di particelle sono stati rilevati su un Athlon XP

2000+, utilizzando un tool a forma di cubo composto da 10*10*10 particelle.

I risultati sono mostrati in tabella:

N° of Particles Simulation time [ms]

60*60*60=216000 16

80*80*80=512000 34

100*100*100=1000000 54

I risultati sono molto interessanti poiché i tempi di simulazione sembrano essere quasi

lineari col numero di particelle, e assolutamente accettabile anche per un milione di

essere

Figura 21. Grafico delle perfomance del sistema


______________________________________________________________________________________________

______________________________________________________________________________________________


Capitolo 3: Rendering in CPU Il rendering [28] è il più grande collo di bottiglia di tutto il sistema fin qui ideato. Il

problema fondamentale è che la simulazione viene eseguita in CPU e

indipendentemente dalla tecnica di rendering che si intende utilizzare ci si trova di

fronte ad una quantità enorme di dati da passare dalla CPU all’acceleratore Video.

Ovviamente questo problema deriva dall’architettura stessa dei PC che non è progettata

per questo genere di applicazioni che sfruttano pesantemente la banda passante.

Esistono altri tipi di architetture, come ad esempio quelle delle moderne console o dei

sistemi SGI che utilizzano architetture memoria di tipo UMA (Unified Memory

Architecture). In poche parole in queste architetture CPU e chip grafici possono leggere

e scrivere nella stessa memoria, riducendo enormemente i tempi di latenza nel

‘passaggio’. Prima dell’avvento delle porte AGP sulle architetture PC, la situazione era

drammatica, in quanto il bus PCI era insufficiente a far passare le decine di giga di dati

per secondo che sono la norma per i rendering odierni. Giusto per fare un esempio di

come sia possibile raggiungere queste moli di dati, facciamo un paio di conti

approssimati. Questi conti ovviamente non tengono in considerazione fattori tipo il

caching ed il DMA (direct memory access) che sono presenti nelle architetture

moderne, ma ad ogni modo l’ordine di grandezza dei numeri che verranno presentati è

corretta.

Immaginiamo di scrivere un singolo pixel. Prima di scrivere il pxiel, b isogna leggere

dallo Z-buffer il pixel stesso (ZR), e se il test viene passato, allora questa z deve essere

riscritta nel buffer (ZW) e ovviamente anche il color buffer dovrebbe essere aggiornato

(CW), ma per ora lo trascuriamo. Se si utilizzasse qualche forma di blending (come

vedremo nella parte sui pixel shader) sarà necessario anche leggere dal color buffer

(CR). Se immaginiamo poi che le nostre primitive abbiano una texture applicata per

aumentarne il realismo, sono necessari uno o più letture da texture (TR). Immaginiamo

il caso comune di due TR con trilinear mipmapping per pixel, con costo totale di 8 x 3 =

24 bytes. Immaginiamo anche che ZR,ZW,CW e CR usino 32 bits (cioè 4 bttres) di

accesso in memoria. Un pixel quindi ‘costa’:

ZR+ZW+CW+2TR= 60 bytes/pixel

Sembrano pochi questi 60 bytes per un pixel, ma se si immagina un contesto reale, con

60 frame al secondo alla risoluzione di 1280 x 1024 abbiamo


______________________________________________________________________________________________

______________________________________________________________________________________________


60x1280x1024x60 bytes/sec à4,5 Gbytes/s

Questo conto ovviamente ha tenuto conto di scrivere ogni singolo pixel una sola volta.

Questo ovviamente non è un caso realistico, in media si ha un fattore di depth

complexity che è circa 4, quadruplicando il numero di tests.

In questo contesto è ben poco quello che si può fare poiché l’ostacolo non può essere

aggirato. Le due tipologie di rendering che è possibile utilizzare in questo campo sono

due: il rendering poligonale e il rendering per pixel.

3.1 Introduzione La prima cosa da dire riguardo al rendering in CPU è che esso è sostanzialmente di tipo

poligonale, in quanto genera delle primitive poligonali (in genere triangoli).

Successivamente questi triangoli vengono mandati alla GPU per essere esclusivamente

disegnati, magari con qualche effetto di illuminazione particolare come vedremo.

Questa è stata la prima soluzione adottata nel rendering del sistema particellare descritto

nel precedente capitolo poiché è più semplice da realizzare e ha permesso di vedere

praticamente se i risultati della simulazione fossero coerenti con quanto previsto. Per

capire come si sia potuto passare dalle singole particelle ai poligoni si osservi la

seguente immagine.

Figura 22. Flood-Fill per riempire di particelle una forma chiusa

A destra vediamo le particelle disegnate nella loro consistenza di sfere a raggio unitario.

Al centro e a sinistra vediamo invece la iso-superficie che esse generano. Questa iso-

superficie è anche detta superficie implicita generata dalla funzione f, e sarà analizzata

dal seguente paragrafo.


______________________________________________________________________________________________

______________________________________________________________________________________________


3.1.1 Superfici Implicite La modellazione geometrica 3D classica si basa sulla combinazione di oggetti base,

detti anche ‘primitive’, che possono essere ad esempio linee, piani, cubi etc. Esiste però

una serie di oggetti curvi, magari deformabili che sono difficili o poco efficienti da

rappresentare con queste primitive, anche utilizzando superfici parametriche come sfere

o superfici tipo bezier/spline. E’ per questo motivo che sono stati inventati metodi per

poter disegnare delle iso-superfici di funzioni scalari in 3D.

Non è stato usato a caso il termine iso-superficie poiché si sta parlando di una funzione

del tipo:

ℜ→ℜ3:f

)(18 kf =)(x

dove f è un campo vettoriale, funzione di tutte le particelle. In generale la nostra

funzione f sarà di questo tipo:

)(19 ( )

ℜ→ℜ

−= ∑−

=:

)(1

0D

Df in

icxx

Ovviamente è la funzione D che fa da padrone in questa formula. Consideriamo una

serie di possibili formulazione per la D.

3.1.2 Funzione D1 Si consideri la funzione:

)(202

1)(x

xD =

E una griglia di punti di controllo nello spazio 3D. x è la distanza di un punto nello

spazio da un particolare punto di controllo. Questi punti di controllo poi sono nel nostro

caso specifico i centri delle sfere,ma per ora continuiamo a chiamarli così, in modo che

il discorso sia più generale. Il campo totale generato dalla funzione f definita

precedentemente, è la somma delle singole D calcolate su tutti i punti di controllo

rispetto ad un fissato punto dello spazio. Ripetendo tale procedura per tutti i punti dello

spazio riusciamo a definire il nostro campo vettoriale.


______________________________________________________________________________________________

______________________________________________________________________________________________


Se proviamo a disegnare f(x)=k con k un valore arbitrario otteniamo la isosuperficie.

Proviamo a vedere cosa succede nel caso il numero di punti di controllo sia uguale ad

uno:

)(21 ( ) ( ) ( )111

1)(1cxcx

cxx−⋅−

=−= Df

Quindi f1(x)=k è il luogo dei punti il cui inverso della distanza al quadrato dal centro

del punto di controllo è uguale a k, che è ovviamente una sfera. Immaginiamo per ora di

“sapere” come si possa fare per disegnare i poligoni ‘giusti’ per approssimare questa

superficie (lo vedremo nel paragrafo sull’algoritmo di Marching Cubes).

Figura 23. Una metaball

Proviamo ora a vedere cosa succede se abbiamo due punti di controllo distinti:

)(22 ( ) ( ) ( ) ( ) ( ) ( )221121

11)(2cxcxcxcx

cxcxx−⋅−

+−⋅−

=−+−= DDf

La nostra iso-superficie diventa di questo tipo:

Figura 24. Due Metaball

Questa superficie assomiglia abbastanza una ellissoide. Proviamo a variare la distanza

tra i due punti di controllo:


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 25. Allontanamento delle metaball

La forma che assume la iso-superficie proiettata su un piano ricorda molto le curve ad

iso-potenziale in un campo elettrico generato da due cariche puntiformi. Proviamo a

variare k, e quindi a cambiare il ‘livello’ di iso-superficie scelto, a parità di distanza:

Figura 26. Variazione del fattore di iso-livello

Come vediamo dalla figura, abbassando il fattore k diminuisce l’influenza ‘attrattiva’ tra

le due particelle. E’ lecito immaginare quindi che non sia necessario calcolare

l’interazione in ogni punto come somma dei contributi di tutte le particelle, ma basterà

limitarsi solamente a quelle che sono più ‘vicine’ al punto considerato.

In realtà noi abbiamo considerato solo il caso di oggetti puntiformi, ma è possibile

generalizzare il concetto di distanza |x-c| con una distanza di x da un ente geometrico

qualsiasi, ad esempio la distanza di un punto da un segmento, utilizzando note formule

della geometria. In tal caso avremmo isosuperfici del tipo (per vari valori di k):

Figura 27. Metaball usando linee come primitive


______________________________________________________________________________________________

______________________________________________________________________________________________


3.1.3 Funzione D2 La funzione D2 è probabilmente una delle prime utilizzate, modellata dal famoso Jim

Blinn (e da lui chiamata Blobby Function) sulle formule della densità di carica in

elettromagnetismo:

)(232

)( bxaexD −=

Il fattore ‘b’ è collegato alla deviazione standard della curva, e il fattore ‘a’ all’altezza.

Questa funzione è riconoscibile come una curva gaussiana centrata all’origine (in realtà

mezza curva in quanto per noi r>0).

3.1.4 Funzione D3 Uno dei problemi della funzione D2 delle consiste nel fatto che la funzione si estende

all’infinito. Questo significa che per calcolare il campo in un punto dello spazio,

bisogna eseguire la sommatoria di tutti i contributi delle primitive di controllo. Questo

ovviamente diventa molto pesante a livello computazionale, al crescere del numero di

primitive di controllo. Questa funzione non presenta questo problema in quanto è

definita come:

)(24

−

−

=

0

12

3

31

)(2

2

2

bxa

bxa

xD

xb

bxb

bx

≤

≤≤

≤≤

3

30

Con ‘a’ fattore di scala, e ‘b’ la massima distanza oltre la quale una primitiva di

controllo contribuisce al calcolo del campo

3.1.5 Funzione D4 Questa funzione è attribuita al lavoro dei fratelli Wyvill, ed è in pratica composta dai

primi termini di una espansione di Taylor della funzione esponenziale, troncata per

restringerla all’intervallo di influenza:


______________________________________________________________________________________________

______________________________________________________________________________________________


)(25

−+−

=

09

229

17941

)( 2

2

4

4

6

6

bx

bx

bxa

xDbx

bx

>

≤≤0

Il fattore ‘a’ scala la funzione, e il fattore ‘b’ come negli altri casi taglia l’influenza oltre

una certa soglia. Questa funzione è più leggera delle metaballs perché usa solo potenze

di 2 della distanza, e quindi non è necessario calcolare radici quadrate per ottenere la

norma prima di elevarla ad una potenza arbitraria.

3.1.6 Funzioni a paragone

Figura 28. Grafici sovrapposti delle funzioni D1,D2,D3,D4

Questo grafico ci mostra le diverse funzioni a paragone. Come vediamo a parte la 1/r^2

vanno tutte a zero oltre una certa soglia e quindi il loro calcolo può venir enormemente

semplificato.

La funzione utilizzata in questo lavoro è la D2. Essa è stata scelta poiché appartiene al

gruppo di quello che è possibile “tagliare” oltre un certo range e poiché tende a

‘gonfiare’ molto di meno la zona compresa tra due control points.


______________________________________________________________________________________________

______________________________________________________________________________________________


3.2 Marching Cubes L’algoritmo dei Marching Cubes (MC) è stato ideato con lo scopo approssimare la

isosuperficie generata da una funzione implicita con dei poligoni. Il concetto di base

dietro ai MC è di suddividere lo spazio in una griglia di piccoli cubi. L’algoritmo poi

provvederà a ‘marciare’ attraverso questi cubi, controllando i vertici di tali cubi e

sostituendo ai cubi stessi dei poligoni. L’insieme di tutti i poligoni generati sarà una

superficie che approssima quella descritta dalla funzione stessa, calcolata però solo nei

vertici di questa griglia regolare. I MC sono anche molto usati per la ricostruzione di

superfici provenienti da scansioni volumetriche in ambito medico. Per esempio le

scansioni MRI descrivono un volume 3D di campioni disposti sui vertici di una mesh

regolare.

Per capire meglio di cosa si tratta è forse meglio guardare al corrispondente algoritmo

bidimensionale (che alcuni chiamano ‘marching squares’).

3.2.1 Marching Squares

Figura 29. Ingrandimento di una cella di confine con Marching Square

L’algoritmo dei marching-squares parte dal presupposto che sia possibile definire un

“dentro” ed un “fuori” per la forma considerata da approssimare. Se la funzione è

analitica, questo si fa semplicemente controllando che sia rispettivamente f(x)<k o

f(x)>k. Successivamente si crea una griglia attorno alla forma stessa, dividendone i

punti che cadono dentro alla forma dagli altri. Preso un singolo quadrato, sapendo quali

vertici sono dentro e quali sono fuori è facile definire per quali lati passerà la linea ‘di


______________________________________________________________________________________________

______________________________________________________________________________________________


confine’ della forma. Se si prende in esame il quadrato in questione infatti:

Figura 30. L’approssimazione fatta dal Marching Square

Ci si rende facilmente conto che collegando tutti i lati del quadrato nei punti in cui la

funzione va a zero, si ottiene una “approssimazione” del contorno ‘vero’ della funzione.

Trovare i punti e ed f dove va a zero la funzione, può venir fatto calcolando

analiticamente la derivata della funzione e ponendola a zero se e possibile farlo, oppure

per interpolazione lineare (caso comune quando si deve approssimare un dataset che

non ha rappresentazione analitica). In tal caso useremo una formula del tipo:

12

1211

))((VV

PPVkPP−

−−+=

Dove i vari Pi sono i punti nella griglia, i Vi sono i valori della funzione nella griglia e

K è la costante della isosuperficie, la stessa usata per f(x)=k.

Per far funzionare correttamente tutti i casi del Marching Cubes, è ora solo questione di

enumerare correttamente tutti i possibili casi. In 2D abbiamo 24 possibili combinazioni

poiché ogni vertice può avere due stati, dentro o fuori. Ecco le 16 combinazioni:


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 31. Tutti i possibili casi del marching square

I casi 5 e 10 sono tratteggiati poiché appartengono alla categoria dei casi ambigui, basti

guardare il seguente esempio:

Figura 32. Caso ambiguo di Marching Square

Ovviamente non si è in grado di dire quale delle due soluzioni sia ‘corretta’.

Ovviamente in una situazione come questa la cosa migliore da fare è diminuire il

‘passo’ della maglia del marching squares.

3.2.2 L’algoritmo L’algoritmo vero e proprio dei marching cubes è la semplice trasposizione in 3D dei

marching squares. In questo caso abbiamo però 28 possibili combinazioni, per un totale

di 256. Considerando però simmetrie e rotazioni, ci si può ricondurre a 15 casi ‘base’.


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 33. Casi ‘base’ dei Marching Cubes

Ovviamente anche in 3D esiste il problema dei casi ambigui, ed è possibile vederlo

nella Figura 32.


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 34. Metaball usando linee come primitive

Essi vengono però risolti aggiungendo famiglie addizionali di casi combinati:

Figura 35. Risoluzione casi ambigui Marching Cubes


______________________________________________________________________________________________

______________________________________________________________________________________________


3.2.3 Ottimizzazioni L’algoritmo dei marching cubes ‘standard’ non è molto veloce e vanno quindi adottati

alcuni espedienti per poterlo utilizzare in real-time su un gran numero di sfere da

tassellare. Si elencano di seguito le ottimizzazioni che sono state realmente

implementate nell’applicativo a corredo di questa tesi.

Ricalcolo zone modificate La parte simulativa del sistema particellare aggiorna ad ogni frame solo un piccolo

subset dell’intero sistema. Di conseguenza non è necessario ricalcolare la funzione che

definisce la iso-superficie su tutti i punti della griglia ogni volta, ma solo su quelli

‘vicini’ alle particelle spostate. Il valore effettivo da attribuire alla distanza limite oltre

la quale si evita di aggiornare, dipende ovviamente dalla funzione utilizzata. Abbiamo

visto per esempio che per la funzione D1 questa distanza è infinita, in quanto la

funzione stessa tende asintoticamente a zero. Per altre funzioni fortunatamente non è

così, ed è quindi stata impostata una procedura generica per il calcolo della zona di

influenza.

Marcia Guidata Jonsson [32] suggerisce un ottimo metodo per ottimizzare la ‘marcia’ dell’algoritmo.

Figura 36. Si cammina esaminando le celle di superficie contigue.

In sostanza suggerisce di partire da un punto di controllo, e sparare un raggio verso

l’esterno fino ad incontrare la superficie. Una volta trovato un cubo di MC che contiene


______________________________________________________________________________________________

______________________________________________________________________________________________


un pezzo di superficie, si ‘cammina’ nelle direzioni che contengono altri pezzi di

superficie. Ottenere questa informazione è banale perché si tratta semplicemente di una

tabella che associa ad un particolare caso del MC le possibili direzioni sulle quali

potrebbe continuare la superficie, che poi sono quelle su cui giacciono i lati dei triangoli

del caso corrente.

Utilizzo della struttura di partizionamento Il loose octree risulta ancora utile anche in fase di rendering. Infatti per calcolare la

funzione di isosuperficie in una determinata ‘zona’, viene prima eseguita una query sul

loose octree per ottenere solo le sfere ‘vicine’, che potrebbero dare una qualche

influenza sul risultato finale. Questa è in assoluto l’ottimizzazione più pesante, dato che

la gran parte del tempo CPU richiesto dall’algoritmo viene speso per calcolare la iso-

funzione nei punti di griglia.

3.2.4 Modalità di calcolo della normale Come si può vedere dal programma in azione, la normale viene calcolata in tre diversi

modalità, ognuna con dei vantaggi e svantaggi propri.

La modalità mesh-normals provvede a calcolare la normale in ogni vertice, prendendo

la risultante delle normali di tutti i triangoli che lo condividono e normalizzandola.

Ovviamente questa procedura non è matematicamente corretta, ma il risultato è

indistinguibile dalle normali calcolate analiticamente per superfici con forme dolci ed

arrotondate. In presenza di cuspidi o curvature molto strette questo metodo fallirebbe

miseramente. Inoltre il vero problema di questa modalità esce fuori quando essa si

accoppia con il ricalcolo della superficie soltanto nelle zone modificate consiste nel

fatto che bisogna pre-allocare un vettore di possibili normali su tutti i punti della griglia.

Nel paragrafo riguardante i requisiti di memoria ci si renderà conto che questo è un

problema non da poco.


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 37. Normale calcolata come risultante normali ai poligoni vicini.

In figura si vede come la normale al vertice evidenziato dal pallino è la risultante

normalizzata (resa versore) delle normali di tutti i poligoni che in esso convergono.

La modalità analitic normals invece calcola la normale differenziando la funzione f:

ℜ→ℜ3:f kf =)(x

)(26 ( ) ℜ→ℜ−= ∑−

=:)(

1

0DDf i

n

icxx

)(27

∂∂

∂∂

∂∂

=∇zf

yf

xff ,,)(x

( )( ) ( ) ( )

( )∑∑∑∑∑

∑∑∑∑

−

=

−

=

−

=

−

=

−

=

−

=

−

=

−

=−

=

∇=

∂

∂∂

∂∂

∂=

∂

∂

∂

∂

∂

∂

=

∂

∂

∂

∂

∂

∂

=∇=∇

1

0

1

0

1

0

1

0

1

0

1

0

1

0

1

01

0

,)(

,)(

,)(

)(,

)(,

)(

,,)(

n

ii

n

i

iii

n

ii

n

ii

n

ii

in

ii

n

ii

n

ii

n

i

Dz

Dy

Dx

Dz

D

y

D

x

D

z

D

y

D

x

D

Df

yyyyyyy

yyy

yx

L’ultima modalità è quella delle flat normals che calcola la normale per ogni triangolo e

quindi permette solo un rendering in flat shading, inaccettabile per un sistema di

modellazione.


______________________________________________________________________________________________

______________________________________________________________________________________________


3.2.5 Rendering adattativo Variando la densità della griglia si possono ottenere diverse approssimazioni della

funzione di isospuerficie, ad esempio si consideri la seguente figura:

Figura 38. Rendering adattativo a diversi livelli di LOD

Il rendering adattativo consiste nel variare la risoluzione della griglia nelle zone a

maggior curvatura, in modo da dare più dettaglio dove effettivamente ce n’è bisogno. In

pratica questa approssimazione viene eseguita creando diverse griglie sovrapposte l’una

all’altra dove il passo di ognuna è il doppio di quello della precedente. In questo modo

si riesce a disegnare la isosuperficie prendendo i ‘pezzi’ approssimati nei vari livelli di

griglia. Si usano multipli di 2 del passo di griglia in modo che sia facile gestire la

giunzione tra ‘pezzi’ appartenenti a livelli diversi.

Figura 39. Rappresentazione di pezzi di una forma ai diversi livelli LOD


______________________________________________________________________________________________

______________________________________________________________________________________________


3.2.6 Requisiti di memoria I requisiti di memoria per la versione correntemente implementata nell’applicativo a

corredo della tesi sono piuttosto elevati. Non poteva essere altrimenti dato che si

salvano quasi tutti i risultati del MC per ogni cella, in modo da poter ‘recuperare’ al

frame successivo i risultati precedenti. Calcoliamo i vari contributi e sommiamo:

v1 = vettore della energia di griglia;

v2 = vettore dei ‘possibili’ vertici e normali per un dato lato della griglia

v1 = (n+1)^3 elementi di tipo float

v2= (n*(n+1)^2)*6 elementi di tipo float (3 scalari per la posizione + 3 scalari per la

normale)

Ovviamente v2 va moltiplicato per 3 poiché ogni triangolo è composto da 3 vertici.

Facendo un veloce calcolo considerando una griglia fatta di 100 posizioni e

considerando un float fatto da 32 bit (4 byte) come è comune sui moderni PC:

Memoria = 4*[(100+1)^3+3*(100*(100+1)^2)*6] = 77568404 bytes = 74 Mb

Una dimensione poco più grande farebbe esplodere i requisiti di memoria necessari

all’applicazione.

3.2.7 Considerazioni Finali I Marching cubes sono sicuramente la soluzione giusta per generare un rendering

‘ottimale’ che calcoli con precisione tutti i poligoni della isosuperficie. Purtroppo la

natura stessa dell’algoritmo lo vede relegato sempre dal lato CPU. Tra l’altro le quantità

enormi di triangoli passati via AGP dalla CPU causano degli stalli incredibili, e

situazioni in cui la CPU è ferma poiché sta aspettando che la GPU si liberi per poterle

mandare i nuovi triangoli. Il campo del rendering volumetrico va sicuramente

approfondito con l’utilizzo delle nuove tecniche dal lato GPU. A tale proposito si vuole

citare la nascita di siti dedicati all’utilizzo delle GPU consumer a basso costo,

provenienti dal mercato dei videogames di ultima generazione, per scopi di

visualizzazione puramente scientifici [33].


______________________________________________________________________________________________

______________________________________________________________________________________________


Capitolo 4: Rendering in GPU Nell’ultimo quinquennio, l’avvento degli acceleratori 3D nel campo dell’entertainment

ha avuto una notevole risonanza anche nel campo della visualizzazione scientifica. E’

possibile infatti oggi avere sul proprio PC a prezzi assolutamente abbordabili,

tecnologie che fino a poco tempo fa erano relegate nelle mani dei pochi fortunati che

avevano accesso alle costose workstation SGI o simili.

E’ nata, ed è in costante evoluzione, una generazione di acceleratori hardware video che

ha delle capacità di programmabilità incredibili. Non si parla più di pura potenza di

calcolo, di ‘numero di triangoli per secondo’, o altre cifre prestazionali poco

significative, quanto della possibilità di creare dei microprogrammi che vengono

eseguiti in un ambiente dedicato, detto GPU.

Questo termine sta per Graphical Processing Unit, ed è sostanzialmente una CPU

specificamente studiata per risolvere problemi visualizzativi in tempo reale. Già da

tempo le GPU hanno superato la complessità delle moderne CPU, tanto da avere un

numero di transistor superiore a queste ultime.

Figura 40. Curva del numero di transistor installati sulle GPU recenti.

Inoltre l’evoluzione con la quale si stanno succedendo versioni sempre più evolute di

questo hardware a basso costo permette ai ricercatori di non fossilizzarsi sull’utilizzo di

una macchina costosa tipo workstation come in passato. Non è quindi più necessario

spendere ingenti somme in questo ambito specifico della visualization technology, le

moderne GPU assolvono egregiamente il questo compito. Peraltro in ambito


______________________________________________________________________________________________

______________________________________________________________________________________________


prettamente scientifico ci si sta orientando quasi ovunque a sfruttarne la potenza non più

semplicemente aumentandone le frequenze di lavoro o il numero di transistor, ma

facendole lavorare in parallelo attraverso cluster di PC.

Figura 41. Confronto qualitativo tra performance GPU/CPU

Le aziende che portano avanti la ricerca su questo tipo di hardware (al momento

principalmente NVidia e ATI con PowerVR e Matrox con minore importanza) hanno a

disposizione know-how e fondi enormi, derivanti dal mercato dell’entertainment per il

quale producono la loro tecnologia. Fino a poco tempo fa queste aziende hanno avuto il

ruolo di portare a basso costo tecnologie già ampiamente disponibili e documentate in

ambito scientifico. Ultimamente però si sta avendo un’inversione di tendenza, dove è la

ricerca nel campo dell’entertainment a ‘dare qualcosa’ al campo scientifico. Non a caso

poi molti gruppi di ricerca universitari sono ‘gemellati’ con qualcuna in particolare delle

aziende citate. Entrambe le parti hanno visto gli innumerevoli vantaggi di una tale

collaborazione. Non sfruttare queste opportunità significa non essere al passo con i

tempi.

Prima di analizzare i dettagli tecnici, cerchiamo di capire quali sono e come sono

classificate le molteplici tecniche di visualizzazione disponibili


______________________________________________________________________________________________

______________________________________________________________________________________________


4.1 Classificazione delle moderne tecniche di visualizzazione Le classiche tecniche di visualizzazione non sono appropriate per visualizzare grandi

dataset, dove per dataset si intendono moli enormi di dati provenienti da vari ambiti

scientifici.

Alcuni esempi [34] potrebbero essere le simulazioni al calcolatore nell’ambito della

fisica, fluidodinamica, chimica etc., oppure da strumentazioni di misura. Queste ultime

sono particolarmente importanti nel campo medico, sismico, satellitare ma anche

prettamente meccanico.

Ogni struttura implica un diverso approccio per la progettazione di algoritmi di

visualizzazione.

Figura 42. Il tipo di dati implica un determinato algoritmo [34]

Le tecniche di visualizzazione avanzata possono essere suddivise in base al dominio di

definizione ed al relativo condominio.

Figura 43. Classificazione delle varie tipologie di rendering [34]


______________________________________________________________________________________________

______________________________________________________________________________________________


4.1.1 Campi vettoriali I campi vettoriali invece possono rappresentarsi ‘sezionando’ il volume contenente il

campo stesso, per poter visualizzare vari tipi di linee ‘interessanti’ per il ricercatore

fluidodinamico, quali ad esempio le strikelines [34] o i flussi di vario tipo, sicuramente

più interessanti nel campo dell’ingegneria meccanica. Lo scopo di queste tecniche però

è puramente visualizzativo, non si pensa assolutamente di ricavare risultati numerici

approssimati da questo genere di algoritmi, per i quali bisogna ripiegare ai sistemi

classici basati su CPU e programmi simulativi in linguaggi tipo C o Fortran.

Figura 44. Visualizzazione di Campi vettoriali fluidodinamici e streaklines [34]


______________________________________________________________________________________________

______________________________________________________________________________________________


4.1.2 Simboli matematici Un’altra categoria di oggetti visualizzabili è quella di strumenti matematici quali

tensori, quaternioni o altri simboli più o meno astratti.

Figura 45. Visualizzazione di un simbolo matematico [34]

4.1.3 Campi scalari 2D E’ possibile disegnare campi scalari bidimensionali, f(x,y) = k, detti anche heightfields,

colorando diversamente zone separate magari da curve di iso-livello.

Figura 46. Visualizzazione di un HeightField [34]


______________________________________________________________________________________________

______________________________________________________________________________________________


4.1.4 Campi scalari 3D

Il caso proprio di questo tipo di simulazione sono i campo scalari tridimensionali,

definiti, come già descritto nel capitolo di rendering poligonale, come f(x,y,z) oppure

acquisiti da strumenti di misura sensibili a determinate grandezze fisiche rappresentabili

come scalari definiti in ogni punto dello spazio, ottenuti per sampling o per

ricostruzione. La terminologia utilizzata in questo caso è quella di ‘voxel’ [38] per

indicare una singola cella della griglia tridimensionale. La griglia invece viene chiamata

Texture 3D, che sarà formalizzata più avanti.

Figura 47. Campo scalare 3D ‘Voxel’ [38]


______________________________________________________________________________________________

______________________________________________________________________________________________


4.2 Teoria del Rendering Volumetrico Lo scopo del rendering volumetrico è quello di ricostruire graficamente un segnale

proveniente da scansioni effettuate ad esempio con delle apparecchiature mediche

attraverso convoluzione con filtri appropriati, oppure di rappresentare funzioni

analitiche, come nel caso del virtual sculpturing, oggetto di questa tesi.

Figura 48. Filtri comunemente usati nella ricostruzione [38]

A differenza degli approcci poligonali, queste tecniche di rendering operano tutte in

DVR (Direct Volume rendering), ovverosia senza estrarre superfici di alcun tipo [38]

Quello che invece si fa è appunto mappare direttamente range di valori scalari con

proprietà ottiche quali colore ed opacità, attraverso una funzione chiamata Transfer

Function. [40]. In [37] possiamo vedere una trattazione delle transfer function anche

multidimensionali, capaci di mostrare i risultati di diversi dataset acquisiti

contemporaneamente

Figura 49. Applicazione di una TF per evidenziare l’interno di una scansione[37]


______________________________________________________________________________________________

______________________________________________________________________________________________


Per ottenere il colore di un singolo pixel sullo schermo, si immagina di creare un raggio

che parte dall’osservatore e attraversa tutto il volume del dataset in considerazione.

Questo raggio prende i valori scalari che incontra sul suo percorso li trasforma con la

funzione di trasferimento e li ‘compone’. La composizione, come vedremo potrà essere

additiva ma anche di altri tipi. Se la composizione è additiva, questa operazione

equivale a risolvere un integrale che viene appunto chiamato Volume Rendering

Integral (equazione )(28 ) [38]

Figura 50. Integrale di Volume Rendering (parzialmente preso da [38])

Questo integrale è composto da due contributi. Il primo è quello che tiene in

considerazione l’emissione iniziale, dipendente sostanzialmente dalle condizioni di

illuminamento dell’ambiente virtuale nel quale si immagina di visualizzare il modello.

Questa emissione perde intensità lungo il percorso del raggio a causa dello scattering

ambientale. A questo contributo si ‘sovrappone’ quello specifico della tranfer function

in questione. La I(s) rappresenta quindi la parte di ‘opacità’ della transfer function

integrata lungo il percorso del raggio, alla quale va accoppiata ovviamente l’altra parte

della transfer function che si occupa del ‘colore’.

L’ultima cosa da notare è che per una data transfer function, questo integrale può essere

precalcolato dal lato CPU ed inserito in una tabella di look-up, ovverosia una texture di

dimensione pari a quella della transfer function.


______________________________________________________________________________________________

______________________________________________________________________________________________


4.3 Tecniche di Rendering Volumetrico Si cercherà ora di spiegare come queste tecniche di rendering possano essere

praticamente implementate in hardware [43]. E’ da dire che molte di queste tecniche

sono nate nelle loro forme primitive già alla fine degli anni ’80 in campo medico, ma le

implementazioni attuali differiscono lievemente dai modelli ‘base’ ideati circa 25 anni

fa. Si introdurrà prima il concetto di texture tridimensionale, e successivamente le due

famiglie di rendering volumetrico principali: le texture slices ed il raytracing.

4.3.1 Texture 3D Nonostante le texture tridimensionali [42] siano state implementate da circa un decennio

nelle workstation High End, e da circa 3 anni nelle schede ‘consumer’, solo poche

applicazioni sono state create per sfruttarle. La realtà è che il rendering volumetrico in

tempo reale è ancora per larga parte una scienza che tenta di ‘incollare’ diverse

immagini 2D su superfici 3D. I pixel program hanno finalmente dato la possibilità di

implementare algoritmi grafici avanzati e le texture 3D sono un modo conciso di

implementare tecniche che sarebbero assolutamente scomode o impossibili da

programmare con texture 1D o 2D.

Definizione Il modo più semplice per immaginare una texture 3D è una pila di texture 2D. Noi

indicheremo con le lettere s, t, r e q le coordinate in larghezza,altezza,profondità e

coordinate nello spazio di proiezione.

Figura 51. Comparazione tra una Texture 2D e una Texture 3D [42]


______________________________________________________________________________________________

______________________________________________________________________________________________


Filtering Con l’aggiunta di una dimensione in più, la nomenclatura del filtering cambia

leggermente. Per una texture 3D il semplice filtraggio lineare di un singolo divello

della catena di mipmap è riferito come trilinear filtering, poiché esegue un filtraggio

lineare su tutti i tre assi della texture, leggendo colori da otto diversi textel. In genere si

è abituati a definire il trilinear filtering come diverse mipmap fuse insieme, ma non è il

caso riguardante le texture 3D. Con le texture 3D ci si riferisce a questo caso con la

terminologia di quadrilineaar filtering in quanto aggiunge un filtering lineare sulla

quarta dimensione delle mipmap volumetriche, leggendo 16 diversi texel. Come

suggerisce questo cambio di nomenclatura, il filtering diventa due volte più pesante

quando si passa da texture 2D a texture 3D. Con le texture 2D, l’anisotropic filtering è

pensato come un filtro su una linea lungo la texture con una larghezza di 2 e di

lunghezza limitata solo dal massimo grado di anisotropia. Questo significa che il filtro

anisotropico con il massimo livello (8) richiede 16 texel. In tre dimensioni, la linea

diventa una ‘fetta’ con laghezza 2, che si estende in due dimensioni fino al massimo

grado di anisotropia. Ora il caso peggiore per il livello massimo di anisotropia diventa

8*8*2=128 texel. Questo è pari a 32 volte il costo della semplice lettura bilineare di una

texture 2D che la maggior parte dei programmatori considera come unità base delle

performance delle texture.

Figura 52. Comparazione Filtering Anisotropico su Texture 2D e Texture 3D [42]


______________________________________________________________________________________________

______________________________________________________________________________________________


Memoria Le texture tridimensionali crescono in modo estremamente rapido con l’aumento della

risoluzione. Per esempio una texture 32bit di 256x256x256 occupa circa 64 mega di

meoria video. Questo problema viene in genere aggirato tenendo il numero di canali per

pixel al minimo. Se i dati da salvare in una texture 3D rappresentano quantità scalri

coem densità o intesità, si provvederà a salvarla come canale di luminanza o alfa e non

come una texture in RGBA. In secondo luogo le imperfezioni degli stili di

compressione DXTC non sono un problema, le texture 3D possono essere compresse

per salvare memoria.

Applicazioni Nonostante i suoi costi di memoria, le texture 3D rendono semplici moltissimi compiti e

possono frequentemente essere più efficienti dell’utilizzo di insiemi di texture 2D per

implementare lo stesso algoritmo, come vedremo successivamente a riguardo delle

texture slices. Tra l’altro i problemi citati sono destinati a diventare meno rilevanti col

passare del tempo, come è già successo. La memoria sugli acceleratori grafici è

raddoppiata di anno in anno, e il numero di letture da texture (texture fetch) per ciclo

aumenta allo stesso modo sugli acceleratori di moderna generazione.

La prima applicazione è quella che vede la texture 3D come una funzione di tre variabili

indipendenti. Possiamo quindi rappresentare una funzione siffatta in una texture 3D,

ottenendo automaticamente approssimazione lineare sui singoli punti di discretizzazione

della griglia.

La seconda applicazione, che poi è quella che interessa questo lavoro di tesi più da

vicino, è quella di rappresentare direttamente campi scalari o vettoriali. Se il campo è

scalare può bastare una 3D Texture con un solo canale, altrimenti si possono affiancare

fino a 4 altri canali, inserendoli nelle componenti RGBA della texture stessa. Possiamo

quindi rappresentare anche campi vettoriali a valori in R4 senza alcun problema. Si può

dire che fu [44] a ipotizzare l’uso delle 3D texture nell’ambito della visualizzazione

scientifica al posto della solita tecnica delle texture slice 2D utilizzata da tempo. Nel

seguito vedremo due tecniche che usano le texture 3D nelle modalità ipotizzate da

Drebin, ma con qualche piccola modifica dovuta all’architettura reale delle GPU

correnti, diversa di quella ipotizzata nell’88 da Drebin stesso.


______________________________________________________________________________________________

______________________________________________________________________________________________


4.3.2 Texture slices La tecnica delle texture slices [38] consiste nel disegnare delle ‘sezioni’ del dataset

volumetrico utilizzando una geometria intermedia o ‘proxy’. Queste sezioni vengono

combinate attraverso operazioni di alpha blending classico per poter ottenere le

trasparenze volute dalla transfer function.

Figura 53. Composizione di ‘fette’ 2D (slice) per il rendering 3D [38]

Texture Slices 2D Nella tecnica delle texture slices 2D si prende il campo scalare e si creano 3 pile di

‘sezioni’ lungo gli assi coordinati. Ogni sezione viene rappresentata con una texture

bidimensionale contenente i valori del campo. A seconda del punto di vista poi viene

scelta la pila che ‘si allinea meglio’ con l’osservatore, ovverosia quella che espone più

superficie verso la camera rispetto alle altre

Figura 54. Composizione di slice allineate all’oggetto [38]

Questa tecnica è semplice e veloce ma consuma molta memoria, in quanto bisogna

tenere in memoria le varie sezioni. Tra l’altro essendo queste sezioni orientate con

l’oggetto, è facile per l’osservatore incontrare bruschi cambiamenti di qualità del


______________________________________________________________________________________________

______________________________________________________________________________________________


rendering. Questi ‘scatti’ derivano dal fatto che in prossimità dei 45 gradi su uno dei 3

assi coordinati, si ha il ‘cambio’ di una pila con un’altra meglio allineate all’osservatore.

Figura 55. Problema del ‘cambio’ del set di slices [38]

Si può ovviare parzialmente a questo problema con l’utilizzo di schemi di cambio ad

isteresi ma il problema fondamentale permane egualmente.

La qualità del rendering non è eccezionale poiché la distanza di sampling non è costante

per i raggi idealmente mandati dal punto di vista dell’osservatore. Un aumento della

qualità si può avere ovviamente con l’aumento del numero di sezioni disegnate, ed è

possibile farlo anche interpolando due slices vicine in alphablending.

Figura 56. Blending di una slice intermedia interpolata tra due diverse slices [38]

Texture Slices 3D Se si ha a disposizione dell’hardware capace di gestire le Texture 3D si può fare

l’inverso di quanto si faceva nello slicing bidimensionale [39]. Si genera una sola pila di

sezioni allineate con la camera, e la si ‘taglia’ (clipping) lungo la forma geometrica di

contenimento scelta (in genere un cubo).


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 57. Slices Orientate alla camera [39].

Ogni sezione ha i vertici correttamente mappati direttamente sulla texture 3D e quindi

provvederà automaticamente a visualizzare sulla sua superficie i giusti pixel presi dal

dataset. Il compositing dei layer si fa esattamente come nello slicing 2D, anche qui per

aumentare la qualità di rendering basta aumentare il numero di sezioni.

Figura 58. A destra una singola slice. A sinistra una composizione di slice

Se proviamo a giocare con la transfer function possiamo avere diversi output di

rendering:


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 59. Due Transfer Function diverse sulo stesso oggetto

Lo slicing 3D è più lento dello slicing bidimensionale, ma con l’hardware di ultimo

grido questa differenza si va assottigliando. In compenso però si risparmia un notevole

quantitativo di memoria, che può essere dedicata ad altri fini, principalmente quello di

aumentare la risoluzione massima dei dataset gestibili.

4.3.3 Pre integrated volume rendering In [45] è presentate una interessante tecnica di accelerazione del rendering utilizzando le

transfer function. Nel contesto delle textured slice questa tecnica richiede i pixel

program, ma supporta a pieno anche transfer function di qualità nettamente superiore a

quella derivante da una semplice composizione di sezioni, in particolare per transfer

function con salti e cuspidi o picchi improvvisi.

Il ‘rendering volumetrico pre-integrato’ riduce il numero di sezioni disegnando le ‘fette’

(slab) di volume tra di esse al posto delle slice stesse. Uno slab è quindi uno spazio tra

due slice consecutive. Per un singolo pixel il problema si riduce al calcolo del

contribuito di colore di un raggio tra due slices.


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 60. Un raggio che passa attraverso due slice

Se si assume di usare l’interpolazione lineare dei dati all’interno del volume lungo

questo raggio, il colore dipenderà solamente da tre numeri: dal valore sulla slice frontale

(punto F della figura), dal valore della slice posteriore (al punto B) e dalla distanza tra le

slice stesse. Se la distanza tra di esse p anche assunta costante, il colore del raggio

dipenderà esclusivamente dai due valori del volume, e quindi una tabella di lookup

bidimensionale è tutto il necessario per descrivere questa dipendenza.

Nel contesto del GPU rendering, questa tabella è ovviamente implementata come una

texture bidimensionale, o meglio una dependant texture per chi è a conoscenza dei

dettagli tecnici del GPU rendering. Questa texture è generata come preprocesso

calcolando l’integrale di rendering volumetrico della transfer function lungo un raggio

che ha come punti inziali e finali (F e B) appartenenti a tutte le possibili coppie di slab.

Questo precalcolo è appunto la ragione del nome pre-integrated volume rendering.

Il corrispondente vertex program dovrà prendere il valore della volume texture sulla

slice frontale e quello sulla slice successiva. Poi questi due valori sono usati come

coordinate texture per indicizzare nella texture che funge da lookup dell’integrale

precalcolato come appena descritto. Nell’appendice sulle implementazioni al paragrafo

0 sarà mostrato il codice pixel shader e quello necessario alla generazione della texture

pre-integrata. Come si vedrà Esso è una semplice applicazione del VRI.


______________________________________________________________________________________________

______________________________________________________________________________________________


4.3.4 GPU Raytraced volumetric Rendering Da quanto esposto sappiamo che è possibile rappresentare la nostra funzione f(x,y,z)

con una texture 3D. Nell’ambito del volumetric rendering [39] però, bisogna cambiare il

modo di considerare la funzione f. Ora infatti non si parla più di f(x,y,z) = k, ma di una

f(x,y,z)<k che definisce lo spazio ‘racchiuso’ dalla nostra superficie. Questo spazio

viene detto spazio pieno mentre il restante è chiamato spazio vuoto. L’algoritmo che

verrà presentata fa uso di loop all’interno dei pixel shader, di istruzioni di controllo

dinamico del flusso di esecuzione (IF,BREAK) e di numero illimitato di letture da

texture, anche con coordinate generate nel pixelshader stesso (dependant read). Infine si

utilizzeranno anche le chiamate di funzione. Prima di descrive l’algoritmo però si

provvede a dare un brevissimo rinfresco sulla nozione di raytracing.

Introduzione al raytracing Il ray-tracing [28] è una tecnica di rendering nella quale i raggi sono usati per

determinare la visibilità dei vari elementi. Il meccanismo base è molto semplice, e

infatti come vedremo è possibile implementarlo in hardware. Nei raytracer ‘classici’ i

raggi sono ‘sparati’ dall’osservatore per ogni pixel del rettangolo che contiene la scena.

Per ogni raggio, si trova il punto di intersezione più vicino all’osservatore stesso. Da qui

si possono mandare altri raggi per ottenere ulteriori informazioni di shading su quel

punto, applicando ad esempio le leggi della riflessione o della rifrazione in base alle

caratteristiche del ‘materiale’ che compone la superficie alla quale esso appartiene.

Figura 61. Funzionamento basilare di un raytracer classico


______________________________________________________________________________________________

______________________________________________________________________________________________


Come si può notare questa nozione di raytracing è applicata più a delle primitive

poligonali o definire da equazioni parametriche. Nel nostro caso è la transfer function a

prendere il sopravvento sulle leggi di riflessione o rifrazione, si potrebbero addirittura

programmare due transfer function diverse per i raggi sparati in questi due casi.

L’algoritmo L’algoritmo di base è molto semplice:

1. Si disegni un cubo che abbia gli otto vertici con coordinate texture 3D scelte in

modo da mappare l’intero volume della texture stessa, cioè con un vertice a

coordinate [0 0 0] e l’angolo opposto a coordinate [1 1 1]. Se il piano di clipping

frontale della camera taglia il cubo, bisogna aggiungere dei poligoni di chiusura

sul ‘buco’ creato.

Figura 62. Cubo con vertici mappati sulla texture 3D

2. Nel vertex shader si emetta la coordinata texture 3D ed un vettore indicante la

direzione del raggio che va dall’osservatore al vertice (step vector).

3. Nel pixel shader, cominciando alla coordinata texture data, si proceda per piccoli

passi lungo questo raggio, entrando in sempre più in profondità nel volume. Ad

ogni passo si legga il contenuto del volume.

La lunghezza del vettore dalla camera al vertice, calcolata nel vertex shader, controlla

direttamente la distanza camminata attraverso il volume ad ogni loop. Di conseguenza

le volume texture contenenti piccoli oggetti avranno bisogno di essere disegnate con


______________________________________________________________________________________________

______________________________________________________________________________________________


lunghezze di questo vettore inferiori a quella di un voxel. Per esempio per una texture

100x100x100 è bene usare una lunghezza dello step vector pari a 100-1.

Problemi dello Step Vector E’ da sottolineare che l’interpolazione fatta di questi vettori lungo il triangolo è lineare,

di conseguenza i pixel che si trovano verso la zona centrale del poligono riceveranno

degli step vector più corti del previsto.

Figura 63. Lo step vector interpolato dal PS si ‘accorcia’ nei pixel centrali

La soluzione ovviamente è quella di imporre la lunghezza di questo vettore nel pixel

shader, prima normalizzandolo e poi scalandolo per la lunghezza effettivamente scelta.

Raytracing con Transfer Function Questa tecnica è anche chiamata dei ‘voxel accumulati’ (accumulated voxels)

riferendosi all’operazione di somma cumulata in cui si trasforma l’integrale di volume

rendering successivamente alla sua discretizzazione.Il funzionamento di base è il

seguente:

1. Man mano che il raggio attraversa il volume, si calcola un altro ‘pezzo’

dell’integrale di VR.

2. Se si è deciso di utilizzare un numero fissato di intervallini di integrazione da

effettuare attraverso il volume, il loop non è più dinamico ma diventa statico, in

quanto ha un numero fisso di iterazioni da compiere. Di conseguenza questo

algoritmo potrebbe essere implementato anche con Shader Model 2.0, stando

attenti però a fare in modo che i pixel del bordo texture siano ‘vuoti’, altrimenti

si vedrà il material ‘stirato’ verso l’infinito.

3. Se il raggio viene fermato quando lascia il range texture [0..1], allora saranno

effettuati meno letture da texture, ma non si potrà ovviamente utilizzare lo


______________________________________________________________________________________________

______________________________________________________________________________________________


shader mode 2.0 per questo scopo, a causa della mancanza di strutture di

controllo dinamico.

4. Infine, se necessario, si può ridurre l’intensità del colore calcolato per un fattore

costante su tutti i pixel. Questo fattore è costante e non calcolato sul singolo

pixel per evitare la formazione di ‘bande’ di colore, tipiche delle procedure di

discretizzazione.

Uno pseudocodice che implemnta la procedura descritta sopre è riportato a seguire,

dove la funzione VRI(s1,s2), che sta per Volume Rendering Integral, è stata già

precalcolata e messa in una texture con la tecnica del pre-integrated volume rendering

descritta precedentemente:

p0 = v; //Posizione corrente sulla superficie

finalColor = 0

while in [0…1] range

p1 = p0+stepVector

finalColor += VRI(p0,p1)*f(p0);

p0=p1

Output finalColor

Raytracing ‘First Hit’ Se nel rendering non si ha la necessità di dover visualizzare il dataset utilizzando le

transfer function, l’algoritmo può essere velocizzato enormemente. Immaginando quindi

a avere a che fare con delle celle ‘solide’ (solid voxel) della griglia 3D, il problema si

riduce al solo avanzare attraverso il volume fino a quando non si trovi una cella piena.

Questo quindi implica che il raggio termini quando lascia il volume o trova materia.

Figura 64. Step successivi nel pixel shader per trovare la superficie


______________________________________________________________________________________________

______________________________________________________________________________________________


Abbiamo due possibilità ora per rappresentare la texture 3D nel caso specifico della

simulazione presentata.

La prima è quella di calcolare tutte le zone dentro/fuori in cpu e questo si fa facilmente

senza costi addizionali rispetto a quello che già si fa. Per tenere in memoria una mappa

del genere basterebbe 1 bit per pixel, che ci dica se il voxel in questione è pieno o vuoto,

ma il minimo che possiamo avere sono gli 8 bit del formato LUMINANCE8. E’

possibile se lo si vuole ‘comprimere’ questa mappa di ‘opacità’ in una texture3D 8 bit,

ma a costo di istruzioni aggiuntive per la decompressione nel pixel shader.

La seconda possibilità è quella di passare semplicemente il campo scalare riportato al

range 8 bit [0 – 255] e lasciar fare al pixel shader la comparazione f(x,y,z) < k per

identificare i punti ‘solidi’ che interrompano il cammino del raggio.

La differenza tra queste due soluzioni è il solito compromesso tra velocità di esecuzione

e memoria utilizzata, dove il vantaggio dell’una implica uno svantaggio dell’altra. Sono

in corso delle sperimentazioni per valutare l’alternativa più vantaggiosa per quanto

riguarda il rendering specifico del sistema particellare descritto che si basa su questa

tecnica.

Calcolo della normale Finora forse non ci siamo resi conto di aver disegnato con gli algoritmi descritti, oggetti

completamente opachi. Mancando informazioni sulla normale alla superficie nel punto

disegnato ci si trova nell’impossibilità di eseguire un qualunque algoritmo di

illuminazione, anche il più banale.

Una soluzione potrebbe essere creare un’altra texture 3D contenente le normali per tutti

i punti della griglia, calcolate in modi diversi a seconda di come è stata generata la

texture 3D principale. Nel nostro caso dato che deriva da una funzione analitica nota, è

possibile calcolare la normale differenziando la funzione stessa nel punto, come è stato

ampiamente descritto nella sezione del rendering poligonale.

Se invece non è possibile calcolare queste normali in modo ‘esatto’, oppure non si è

disposti a spendere la memoria texture per la mappa delle normali, allora bisogna

calcolare una normale per ogni pixel nel pixel shader. Questo è un problema non da

poco e richiederebbe una descrizione approfondita, ma di seguito sarà descritta una

tecnica che possa dare una buona approssimazione del calcolo della normale. Ad ogni

modo è bene sapere che esistono metodi migliori per svolgere questo compito, con

maggiore qualità visuale di rendering. Ad esempio si possono tracciare tre raggi per


______________________________________________________________________________________________

______________________________________________________________________________________________


pixel e calcolare il prodotto vettoriale dei punti di intersezione con il solido, ma è una

procedura molto costosa in termini di calcolo e ricca di casi singolari da gestire.

Una soluzione migliore è quella illustrata in figura, e di seguito descritta:

Figura 65. Algoritmo approssimato per il calcolo della normale nel PS

1. Il punto di intersezione viene trovato con l’algoritmo descritto precedentemente

2. Si prende un numero di punti ‘campione’ su una sfera centrata attorno al punto.

Per un risultato soddisfacente può bastare prendere otto punti di sampling ben

disposti.

3. Si sommano gli offset vector, ovverosia i vettori che vanno dal punto di

intersezione al punto campione (i vettori a,b,c, e di nella figura) che non toccano

il solido e si normalizza il vettore risultante (N nel diagramma).

I punti a,b,c etc., verranno di seguito chiamati punti di sampling per la ricerca della

normale. Con la lettera N, ovviamente, si indica la normale approssimata cercata

Migliorie I punti di sampling scelti per la generazione della normale possono essere ottimizzati

per ottenere risultati migliori.

1. Gli offset vector utilizzati per scegliere i punti di sampling sulla sfera, possono

essere salvati in costanti del pixel shader. Queste posizioni verranno salvate per

un solo ottante della sfera, i restanti sette ottanti possono essere calcolati

invertendo i segni degli offset vector su ogni asse coordinato. Di conseguenza

tre offset vector, con segni cambiati applicati al punto di intersezione del raggio


______________________________________________________________________________________________

______________________________________________________________________________________________


con la superficie, danno coordinate per 24 possibili sample. In figura si vede

l’analogo caso 2D, dove salvando due vettori di offset per quadrante avremmo

un totale di 8 punti di sampling.

2. Si possono raggiungere risultati migliori effettuando il sampling da diverse sfere

con raggi variabili. Ovviamente questo si fa facilmente variando la lunghezza

degli offset vector.

Bisogna però stare attenti ai possibili errori che si possono incontrare nell’utilizzo di

questa tecnica. Si immagini ad esempio di avere una risoluzione come in figura:

Figura 66. Caso ‘difficile’ da trattare nel calcolo della normale

I punti di sampling 2,3,4,6,7 e 8 ci diranno di trovarsi nello spazio vuoto. La risultante

di questi vettori ovviamente ci darà una normale di lunghezza zero. Per ovviare a questo

problema,bisogna introdurre un ulteriore controllo: quando si calcola la normale

bisogna tenere in considerazione solo i punti dove il prodotto scalare tra gli offset vector

e lo step vector è negativo ( in figura ad esempio passerebbero questo test solo i punti

5,6,7 e 8).

Questo controllo porta con se anche il positivo effetto di dimezzare il numero di letture

da texture, che è un grandissimo vantaggio.

Ad ogni modo si sta valutando se sia più conveniente precalcolare le normali in CPU e

passarle alla GPU frame per frame poiché questa tecnica per il calcolo ‘al volo’ della

normale è risultata abbastanza pesante dai test.


______________________________________________________________________________________________

______________________________________________________________________________________________


Z-Buffering Nel caso si utilizzino transfer function, e quindi si disegnino volumi ‘traslucenti’ nella

scena, bisogna stare attenti a disegnarli dopo ogni altro oggetto. Nel caso specifico di

questa tesi ad esempio bisogna effettuare il disegno sicuramente dopo i tool utilizzati

dal simulatore.

Figura 67. Un oggetto esterno dovrebbe ‘bloccare’ i raggi

In figura si può notare come un raggio idealmente tracciato dalla faccia del volume di

fronte all’osservatore, non debba terminare alla fine del volume,ma anche quando

incontra un oggetto. Quindi se l’oggetto si muove gradualmente attraverso il volume, da

dietro verso l’osservatore, l’intensità visibile di fronte all’oggetto non diminuirà

dolcemente, bensì passera da una intensità piena a nessuna intensità. Per risolvere

questo problema bisognerebbe avere conoscenza dei contenuti dello Z-Buffer, cosa

possibile solamente tenendone una copia in una texture,aggiornata ogni frame.

Quand’anche non si usassero transfer function, e quindi si abbia a che fare con voxel

solidi provenienti da raytracing del tipo ‘first hit’, bisognerebbe che il pixel shader

calcolasse anche il valore di profondità corretto (nel registro oDepth) all’interno del

pixel shader stesso. Se è possibile evitare questa cosa, magari per gli oggetti solidi per i

quali si è sicuri di non avere mai la possibilità di avere intersezioni con altri oggetti,

questo non è necessario e va assolutamente evitato per questioni di performance.

L’hardware 3D corrente infatti non gradisce molto che i pixel shader scrivano anche

nello Z-Buffer, perché impediscono una serie di importantissime ottimizzazioni tra le

quali citiamo solo lo Hierarchical Z-Buffer [46]. Ad ogni modo se c’è assolutamente

bisogno di avere il valore corretto nello Z-Buffer non si può fare molto,bisogna scrivere

il registro oDepth dal pixel shader, scrivendo la Z del punto di intersezione.

E’ da notare che a differenza del rendering con transfer function, non è necessario avere

il depth buffer anche in input al pixel shader per il rendering first hit.


______________________________________________________________________________________________

______________________________________________________________________________________________


Capitolo 5: Il Software Implementato Tutto il sistema qui presentato è stato implementato e testato in linguaggio C++ su un

normale PC (AMD Athlon XP 2000+, ATI Radeon 9600 Pro, 512 Mb di RAM). Il

compilatore utilizzato è il Visual C++ 6 della Microsoft insieme al kit di sviluppo delle

DirectX 9, in modo da assicurare le migliori performance per il rendering sul target

scelto. E’ stato programmato anche il manipolatore 3D visibile in quasi tutte le

immagini, poiché DirectX, in quanto API a basso livello, ne è sprovvista.

Per il futuro sono già previsti porting in ambienti VR Immersivi come quelli descritti in

[51].

Le soluzioni presentate in questa tesi sono il risultato di un processo iterativo nel quale

sono stati provate diverse strutture dati e schemi di aggiornamento. E’ stato speso

anche del tempo con l’interfaccia visuale. E’ stato usato il framework wxWidgets per la

creazione delle GUI (una libreria Open Source Multi-platform). Tra le caratteristiche

interessanti dell’applicazione c’è quella di poter prototipare la funzione di risposta

salvandola e leggendola da disco, molto utile in fase di sperimentazione. E’ possibile

selezionare le varie modalità di calcolo della normale, decidere la costante k di

isosuperficie, il numero di aggiornamenti per frame e creare tools di forme predefinite.

Ecco uno screen -shot dell’applicazione:

Figura 68. L’applicazione realizzata per questa tesi


______________________________________________________________________________________________

______________________________________________________________________________________________


Ecco visualizzata la rappresentazione interna del blocco così come essa è vista dal

simulatore:

Figura 69. La rappresentazione interna del blocco di materiale

C’è la possibilità di caricare modelli provenienti da scansioni di vario genere. Di

seguito ad esempio c’è la rappresentazione particellare di una scansione MRI di una

testa umana:

Figura 70. Un teschio proveniente da una scansione MRI


______________________________________________________________________________________________

______________________________________________________________________________________________


Esiste anche una modalità che riporta tutta la simulazione su un piano, modalità

utilizzata in fase di prototipazione del sistema:

Figura 71. Modalità di deformazione 2D

E’ possibile ‘staccare’ un pezzo di materiale ed usarlo come tool:

Figura 72. Materiale originario usato come tool

Di seguito è possibile vedere le particelle dopo una deformazione applicata dal tool.

C’e’ da dire che vedere in azione il programma è tutt’altra cosa, le immagini statiche


______________________________________________________________________________________________

______________________________________________________________________________________________


non rendono assolutamente l’idea di quello che sta succedendo. Nella seguente figura si

può vedere la funzione di risposta utilizzata:

Figura 73. Risultati della deformazione

Ora per provare il funzionamento della parte ‘attrattiva’ della funzione LJ, si

selezionino alcune particelle definendole come ‘tool’ e si provi a staccarle.

Figura 74. Selezione delle particelle del tool

Allontanando il tool lentamente, le particelle rimangono attaccate ad esso per effetto

della funzione LJ:


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 75. Le particelle vicine al tool rimangono ‘attaccate’ grazie ad LJ

Si crei un blocco di materiale per provare ad ‘incidere’ qualcosa su di esso. Da notare

che questo è un tool fatto di sfere, come descritto nell’apposita sezione dedicata al

paragrafo 2.7.1 .

Figura 76. Il blocco di materiale e la punta pronta per ‘incidere’ qualcosa

Si usi il tool trascinandolo per incidere il materiale:


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 77. Si trascini la punta e si incida il materiale

Figura 78. Si continui il percorso


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 79. Si chiuda il percorso

Per vedere bene come l’effetto di tensione superficiale sia presente grazie alla funzione

LJ, creiamo un blocco compsto da relativamente poche particelle:

Figura 80. Un blocco composto da poche particelle


______________________________________________________________________________________________

______________________________________________________________________________________________


Provando a spingere con il tool sul materiale si possono vedere gli effetti della ‘tensione

superficiale’ imposta dalla funzione di LJ. La superficie non si ‘rompe’ subito, ma

bisogna applicare una certa forza affinché questo avvenga:

Figura 81. Iniziamo a premere sul materiale. La superficie non si rompe.

Dopo aver premuto un po’, la superficie si rompe creando un buco nel materiale

Figura 82. La superficie si è rotta.


______________________________________________________________________________________________

______________________________________________________________________________________________


Ecco il ‘foro’ creato:

Figura 83. Il foro creato

Figura 84. Il foro visto da lontano


______________________________________________________________________________________________

______________________________________________________________________________________________


Si provi a tirare una sfera sul materiale:

Figura 85. Ecco la ‘vittima’ da rompere

Figura 86. La pallina tocca la superficie


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 87. La pallina entra nella superficie

Figura 88. Oramai il danno è fatto


______________________________________________________________________________________________

______________________________________________________________________________________________


Figura 89. Alcune immagini random


______________________________________________________________________________________________

______________________________________________________________________________________________


Capitolo 6: Appendice A In questo capitolo saranno dati dettagli tecnici implementativi delle tecnologie citate o

realizzate effettivamente all’interno del programma a corredo della tesi. E’ stato scelto

di fare un capitolo a parte su questo aspetto per evitare di spezzare i discorsi portati

avanti nei precedenti capitoli con dettagli tecnici. Questo non vuol dire che i dettagli

tecnici non siano importanti, anzi è fondamentale per un buon ricercatore avere una

solida conoscenza di come funzionano ‘dietro le quinte’ determinati algoritmi o test

geometrici e conoscerne la complessità. In questo modo ci si può creare un bagaglio di

conoscenze e di algoritmi che saranno la base per nuove idee e nuovi algoritmi

nell’inarrestabile mondo della ricerca.

Questo capitolo è una collezione quindi di spezzoni di codice citati negli altri capitoli e

spiegati in maggior dettaglio in questa sede.

6.1 Test di Intersezione Sfera-Cubo Un algoritmo per testare una sfera e un AABB è stato presentato da Arvo [47] ed è

sorprendentemente semplice.

In questo codice si usa una versione ulteriormente semplificata: //Rappresentazione intera di un valore in floating point #define MUSE_IR(x) ((unsigned int&)(x)) //Rappresentazione in floating point di un valore intero #define MUSE_FR(x) ((float&)(x)) //Procedura ottimizzata per 'ripulire' un numero floating point //dal bit di segno. Può essere più o meno veloce a seconda del contesto static __forceinline float FastFabs(float x) { unsigned int mfloatBits = MUSE_IR(x)&0x7fffffff; return MUSE_FR(mfloatBits); } bool LooseCell::RecursiveContactingObjects(SimulationSphere* spheres, int nSpheres,LooseSphereQuery* query,LooseOctree* octree,LooseQueryCallback callback,void* userData) {


______________________________________________________________________________________________

______________________________________________________________________________________________


//Controlliamo se la sfera e l'AABB non si toccano const float test = octree->HalfSizeTable[depth]+query->radius; if( FastFabs(cx - query->x)>test|| FastFabs(cy - query->y)>test|| FastFabs(cz - query->z)>test) { return false; } //...missing }

Come si può vedere dal codice, al di là di ottimizzazioni specifiche che sfruttano la

particolarità dei numeri floating point di mantenere il segno in un bit ben preciso, quello

che si fa è trasformare lo spazio nel quale si effettua la comparazione. In due dimensioni

il discorso è quello in figura.

Figura 90. Trasformazioni geometriche nel Test AABB-Sfera

Come prima operazione si trasporta la circonferenza al centro degli assi. La seconda

cosa che si fa è quella di aumentarne il raggio di una quantità pari a metà diagonale del

quadrato. A questo punto basterà controllare che il centro del quadrato si trovi

all’interno della ‘nuova’ circonferenza creata. In 3D è tutto analogo,con la sola aggiunta

della terza coordinata.


______________________________________________________________________________________________

______________________________________________________________________________________________


6.2 Funzione di risposta void CMetaballs::RecalcGridInfluencedRadius() { int counter=1; float fDist = 0; const float thresold=m_fLevel/1000.0f; while(EnergyFunction(fDist*fDist)>thresold) { fDist+=m_fVoxelSize; counter++; } m_nGridInfluencedRadius = (counter+1)/2; m_fInfluenceRadius = m_nGridInfluencedRadius*m_fVoxelSize; }

6.3 Calcolo raggio funzione di influenza del MC Il calcolo del raggio di influenza della funzione di blending delle particelle è molto

semplice, semplicemente ci si allontana dalla funzione fino a quando essa non void CMetaballs::RecalcGridInfluencedRadius() { int counter=1; float fDist = 0; const float thresold=m_fLevel/1000.0f; while(EnergyFunction(fDist*fDist)>thresold) { fDist+=m_fVoxelSize; counter++; } m_nGridInfluencedRadius = (counter+1)/2; m_fInfluenceRadius = m_nGridInfluencedRadius*m_fVoxelSize; }


______________________________________________________________________________________________

______________________________________________________________________________________________


6.4 Ricerca dell’ottante giusto del looseoctree LooseCell* LooseCell::RecursiveInsert(SimulationSphere* spheres, int nSpheres,int sphereIndex,LooseOctree* octree) { SimulationSphere* sphere=spheres+sphereIndex; //...Missing // Selezioniamo l'ottante in base al segno del vettore che // va dal centro della sfera al centro del cubo dell'octree const int i = (sphere->mCenter.x <= this->cx) ? 0 : 1; const int j = (sphere->mCenter.y <= this->cy) ? 0 : 1; const int k = (sphere->mCenter.z <= this->cz) ? 0 : 1; //...Missing }

Questo banalissimo pezzo di codice riesce a dire con pochissime istruzioni in quale

ottante di un cubo di octree suddiviso, cioè in quale nodo figlio deve venire posizionata.

La cosa importante da tenere a mente però è la convenzione di indici usata quando si ha

la necessità di indicare un determinato nodo figlio.

Figura 91. Numerazione dei nodi figli dell’octree


______________________________________________________________________________________________

______________________________________________________________________________________________


6.5 Codice di crazione texture pre-integrata Il seguente codice precalcola per tutte le possibili coppie di punti iniziali e finali di

integrazione la pre-integration map definita dalla seguente formula già illustrata: void createLookupTable(float thickness, float tfunc[TABLE_SIZE][4], float lutable[TABLE_SIZE][TABLE_SIZE][4]) { for (int x = 0; x < TABLE_SIZE; x++) { for (int y = 0; y < TABLE_SIZE; y++) { int n = 10 + 2 * abs(x-y); double step = thickness / n; double dr = 0.0, dg = 0.0, db = 0.0,dtau = 0.0; for (int i = 0; i < n; i++) { double w = x + (y-x)*(double)i/n; if ((int)(w + 1) >= TABLE_SIZE) w = (double)(TABLE_SIZE - 1) - 0.5 / n; double tau = step * (tfunc[(int)w][3]*(w-floor(w))+ tfunc[(int)(w+1)][3]*(1.0-w+floor(w))); double r = exp(-dtau)*step*( tfunc[(int)w][0]*(s-floor(w))+tfunc[(int)(w+1)][0]*(1.0-w+floor(w))); double g = exp(-dtau)*step*( tfunc[(int)w][1]*(w-floor(w))+ tfunc[(int)(w+1)][1]*(1.0-w+floor(w))); double b = exp(-dtau)*step*( tfunc[(int)w][2]*(w-floor(w))+ tfunc[(int)(w+1)][2]*(1.0-w+floor(w))); dr += r; dg += g; db += b; dtau += tau; } } lutable[x][y][0]=(dr > 1.0 ? 1.0f : (float)dr); lutable[x][y][1]=(dg > 1.0 ? 1.0f : (float)dg); lutable[x][y][2]=(db > 1.0 ? 1.0f : (float)db); lutable[x][y][3]=(1.0 - exp(-dtau)); } }


______________________________________________________________________________________________

______________________________________________________________________________________________


6.6 Dettagli tecnici sul GPU rendering attuale Come una qualunque macchina, prima di poter sfruttare una GPU bisogna capire come

essa funziona e cosa è capace di fare. Il secondo punto è strettamente legato ai modelli

particolare di schede che si considera, ma gran parte delle caratteristiche base sono

condivise dai modelli della stessa ‘fascia’ promossa dai singoli Hardware Vendors

(HV). In questa sede verranno esposte potenzialità e tecniche che si applicano ad

hardware programmabile di seconda e terza generazione [48] (per intenderci con

supporto pixelshader 2.0 e 3.0), che sono le prime veramente ‘utilizzabili’ in ambito

scientifico. Si cercherà ora di spiegare in breve in cosa consistono le tecnologie citate,

che hanno dei nomi sicuramente molto d’effetto ma che sicuramente non ne esplicano la

funzionalità.

La terminologia ‘vertex program’ e ‘pixel program’ sono sicuramente più corrette dei

corrispondenti commerciali vertex shader e pixel shader. Nel seguito si farà riferimento

a queste due tecnologie indifferentemente con le abbreviazioni VP e PP o VS e PS.

6.6.1 La pipeline

Figura 92. Diagramma schematico di una Pipeline di GPU moderna

Per pipeline si intende la sequenza di operazioni che porta una primitiva ad essere

disegnata a schermo. In figura è riportato uno schema generico di pipline. Quello che ci

interessa sapere è solo quanto descritto di seguito. Alla base di tutto il processo ci sono i

vertici, che formano le primitive, ad esempio un triangolo o una superficie parametrica,


______________________________________________________________________________________________

______________________________________________________________________________________________


e delle mappe di pixel N-dimensionali (texture 2D, 3D etc). La pipeline prende

singolarmente i vertici delle primitive ed esegue su ognuno di loro un vertex shader

attraverso il Vertex Processor.Il VS prenderà un vertice in input e lo trasformerà in un

certo output, con la possibilità di aggiungere delle quantità calcolate ‘al volo’ su di esso.

Dopo alcuni processi tecnici necessari (culling, clipping etc) la primitiva approda al

pixel shader, che provvede a ‘colorare’ un pixel alla volta, eseguendo tramite il

fragment processor un microprogramma su ognuno di essi . Di seguito si parlerà di

primitiva riferendosi principalmente a dei triangoli.

Probabilmente con le poche parole che è possibile dedicare in questa introduzione

all’argomento non si riesce a capire il perché della necessità di ‘due livelli’ di

programmabilità, uno per vertice ed uno per pixel,ma gli tecniche che verranno esposte

attireranno sicuramente anche l’attenzione del più scettico.

6.6.2 Vertex Program

Figura 93. Il Vertex Program all’internod della pipeline

Come abbiamo appena detto, un Vertex Program è un piccolo programmino che viene

nativamente eseguito sulla GPU, compilato in un suo linguaggio macchina.

L’applicazione che sembra più naturale per un vertex program, è quello di modificare le

proprietà geometriche di un vertice, sostanzialmente di ‘trasformarlo’. Si potrebbe

obiettare che anche prima dell’avvento della programmabilità della GPU, era possibile

impostare delle matrici di trasformazione per vertice, ma già ci si rende conto di essere

su due piani diversi. Un vertex program ha la possibilità di eseguire una qualunque

operazione di trasformazione sul vertice, non solo una trasformazione lineare come può

essere quella effettuata tramite una matrice. Ad ogni modo vediamo come è strutturato

un semplice vertex program:


______________________________________________________________________________________________

______________________________________________________________________________________________


float4x4 view_proj_matrix; float4x4 texture_matrix0; struct VS_OUTPUT { float4 Pos : POSITION; float3 Pshade : TEXCOORD0; }; VS_OUTPUT main (float4 vPosition : POSITION) { VS_OUTPUT Out = (VS_OUTPUT) 0; // Si trasforma il vertice in projspace Out.Pos = mul (view_proj_matrix, vPosition); // si crea una coordinata texture a partire dalla posizione // da una matrice Out.Pshade = mul (texture_matrix0, vPosition); return Out; }

Questo piccolo vertex shader trasforma il vertice in questione proiettandolo a schermo e

generando dinamicamente una coordinata texture al suo interno. Quest’ultima viene

ottenuta moltiplicando la posizione in object space dell’oggetto per una matrice di

trasformazione texture.

Quello che è importante capire di questo esempio è che i nuovi ‘campi’ che vengono

creati per vertice possono avere dei significati arbitrari. In questo caso specifico i valori

inseriti nel campo Pshade potranno in seguito essere utilizzati come tali nei conti oppure

per indicizzare una texture. Astraendosi un attimo dall’ambito specifico dello shading a

fini visualizzativi, si possono osservare alcune cose:

1. Disponibilità di un linguaggio di programmazione evoluto

2. Programmi eseguiti su una architettura dedicata

3. Possibilità di indicizzare tabelle, eseguire operazioni SIMD (single instruction

multiple data) per processare parallelamente diversi dati.

Tutta questa potenza è dedicata ai soli fini visualizzativi. La CPU potrà in questo modo

esser scaricata di una mole di lavoro enorme ed occuparsi di compiti per i quali è

certamente più portata, quali ad esempio gli aspetti simulativi e di gestione. Il segreto di

questa architettura, ma contemporaneamente la sua debolezza sta nell’unidirezionalità

delle comunicazioni tra CPU e GPU. La CPU manda continuamente dati alla GPU, che


______________________________________________________________________________________________

______________________________________________________________________________________________


provvede a processarli e a disegnarli a video. Non è assolutamente possibile però per la

CPU riutilizzare l’output della GPU per un ulteriore processing, o meglio, è possibile

ma a costo di una irreparabile rottura del parallelismo del sistema. Questa rottura si

manifesta infatti con un calo vertiginoso delle performance. Non avrebbe senso tra

l’altro ipotizzare un sistema del genere poiché a conti fatti equivarrebbe ad un sistema

con due CPU, con tutti i limiti che una scelta del genere comporta. La difficoltà del

progettare algoritmi per la GPU sta proprio in questo: ricordarsi costantemente che si sta

lavorando su uno stream processor. Il termine stream (flusso) indica che è in grado di

gestire un enorme flusso di dati unidirezionale, come un fiume in piena.

Figura 94. Comunicazione unidirezionale da CPU a GPU

Questo concetto è perfettamente esemplificato nella lettura [49]

6.6.3 Pixel (o fragment) program

Figura 95. Posizionamento del Pixel Program nella Pipeline


______________________________________________________________________________________________

______________________________________________________________________________________________


Una volta usciti dal vertex program, i triangoli sono composte da un insieme di vertici

con delle quantità numeriche ben definite al loro interno, che possono venire

interpretate nei modi più disparati per riempire l’area coperta a schermo dal triangolo

stesso.

Figura 96. Rasterizing di un triangolo

Nella fase di disegno vera e propria a schermo, queste quantità vengono interpolate

lungo tutto il triangolo. Su ogni pixel ‘interpolato’ viene eseguito un microprogramma

che permette di effettuare ulteriori calcoli prima di plottarlo a schermo.

Figura 97. Il Pixel Program viene applicato ad ogni pixel disegnato

In figura ad esempio si prende il colore calcolato per vertice, ed interpolato sull’i-esimo

pixel e si scrive a schermo l’esponenziale. Quella appena descritta potrebbe benissimo

essere una procedura per riportare in un range visualizzabile su un normale monitor PC

(LDR o Low Dynamic Range), un colore calcolato in un altro spazio, quale ad esempio


______________________________________________________________________________________________

______________________________________________________________________________________________


l’HDR (High Dynamic Range). Per gli interessati questa procedura si chiama Tone

Mapping [50].

6.6.4 Novità dei modelli VS e PS 3.0 I Modelli VS e PS 3.0 sono quelli implementati dalle nuovissime schede in commercio

al tempo della scrittura di questa tesi. I prezzi di questo tipo di hardware non sono

ancora a livello di consumer, ma sicuramente alla portata di chi vuole utilizzarle per

scopi scientifici o di ricerca. Il pregio di questi nuovi modelli è quello della

‘convergenza’ di features delle diverse proposte da parte dei vari Hardware Vendors. Si

è ragionevolmente certi di poter sviluppare algoritmi che funzionino indifferentemente

sui diversi modelli da essi proposti.

Dichiarazioni Input/Output Nel modello shader 3.0 è obbligatorio dichiarare i registri di input (sia per VS che PS) e

quelli di output (solo per i VS) prima che essi vengano utilizzati. In compenso però i

registri dichiarati hanno la caratteristica di essere in grado di contenere più di un tipo di

dichiarazione. Questo significa che anche se il numero di input di uno shader program

eccede il limite di registri di input a disposizione, gli input possono essere impacchettati

insieme nei registri nella fase di dichiarazione dello shader. Per esempio una

dichiarazione di VS 3.0 potrebbe accettare posizione,normale,tangente, binormale,due

color, otto pesi di blending, otti indici di blending e 14 coodinate texture bidimensionali

attraverso un uso ottimale dell’impacchettamente dei registri di input. Con lo shader

model 2.0 una tale dichiarazione sarebbe stata possibile solo pre-impaccando tutti questi

dati in una struttura vertice utilizzando 16 input totali. Gli outputs sono dichiarati allo

stesso modo degli input. Ad esempio: vs_3_0 ; Declare inputs dcl_position0 v0.xyzw dcl_normal0 v1.xyz dcl_tangent0 v2.xyz dcl_binormal0 v3.xyz dcl_blendweight0 v1.w dcl_blendweight1 v2.w dcl_blendweight2 v3.w dcl_texcoord0 v4.xy dcl_texcoord1 v4.wz


______________________________________________________________________________________________

______________________________________________________________________________________________


; Declare outputs dcl_position o0.xyzw dcl_texcoord0 o1.xy dcl_texcoord1 o1.zw dcl_texcoord2 o2.xyz dcl_fog o2.w dcl_psize o3

Poichhe il numero di VS output e di pixel shader input è lo stesso (12) questa

caratteristica non è vantaggiosa nel pixel shader come lo è nel vertex shader. Il metodo

preferito per selezionare gli input nei PS è lo swizzling sorgente arbitrario, di cui si

parlerà in seguito.

Queste dichiarazioni di input ed output così flessibile permettono di accoppiare diversi

VS e PS senza assicurarsi di utilizzare esattamente lo stesso assegnamento dei registri.

Questa può essere una caratteristica utile quando si ha a che fare con un grande numero

di vertex e pixel program.

Predicazione Il registro dei predicati (p0) è un insieme di 4 flag booleani (uno per x,y,z, e la w) che è

sostanzialmente una ‘maschera di scrittura dinamca’. Esso permette alle istruzioni

shader di essere eseguite per ogni singolo canale, basandosi sui risultati dei calcoli

precedenti. I flag nel registro dei predicati sono impostati con il comando setp_comp p0,

src1,src2, dove comp è la modalità di comparazione (maggiore di, minore di, etc.), p0 è

il registro dei predicati e src1 e src2 sono due registri di input. La comparazione viene

eseguita quatrro volte sulle corrispondenti component dei registri sorgenti, e i risultati

sono salvati come flag booleani nel registro dei predicati.

Per esempio il seguente codice imposta i registri dei predicati a (falso,vero,falso falso): def c0, 0.0f, 2.0f, -4.0f, 1.0f def c1, 4.0f, -8.0f, -2.0f, 1.0f mov r0, c0 setp_gt p0, r0, c1 Una volta che il registro dei predicati è stato impostato, il suo contenuto può essere

utilizzato per permettere o impedire che vengano portate a termine operazioni ‘per

canale’. Per abilitare la predicazione, (p0) è aggiunto di fronte alla corrispondente

istruzione aritmetica o di texture. Per esempio, basandosi sul contenuto del registro dei


______________________________________________________________________________________________

______________________________________________________________________________________________


predicati definito precedentemente, solo la compnente y del registro di destinazione r0 è

presa in considerazione dalla seguente istruzione: (p0) mul r0, r1.x, c1 Usare la predicazione come una maschera di scrittura dinamica ha vari campi di

applicazione. Per sequenze di branching corte, bisognerebbe preferire istruzioni del tipo

if_comp. Il registro dei predicati usa anche meno registri temporanei in comparazione

alla versione equivalente che non usa i predicati, che può aiutare il compilatore ad

eseguire alcune ottimizzazioni. Istruzioni di controllo condizionale come loop, if_comp,

etc.,non possono essere usate insieme al predicate mode, anche se i registri di predicate

possono essere usati come condizionale usando le istruzioni dedicate if_pred,

callnz_pred, e break_pred.

Controllo di flusso statico e dinamico Il modello vs_2_0 supporta il controllo statico del flusso del programma (cioè istruzioni

di branching come loop e sottofunzioni o subroutine che sono chiamate in base a

determinati valori presenti nei registri costanti). Il controllo di dlusso statico permette di

combinare differenti shaders in shaders più lunghi, riducendo il numero di cambi di

shader durante il processo. I loop statici posso anche essere utili con un numero fisso di

iterazioni (ad esempio effettuare un ciclo su un certo numero di luci). L’unica differenza

nel modello 3.0 sul controllo statico di flusso, è il numero di loop che possono essere

innestati l’uno dentro l’altro (nesting). Mentre il l vs_2_0 non supporta il neesting (cioè

avere loop dentro loop), sia vs_3.0 e ps_3_0 supportano il controllo di flusso statico con

una profondità di nesting di 4. Comunque la vera potenza del modello 3.0 deriva

dall’abilita di supportare il controllo di flusso dinamico.

Il controllo di flusso dinamica è un modo di specificare diversi percorsi di codice in

base alla comparazione di registri il cui contenuto è modificato dinamicamente dentro lo

shader program. Ci sono due vantaggi principali con questa nuova carrteristica:

flessibilità e performance. La flessibilità deriva dal fatto che diverse istruzioni di

branching possono ora essere eseguite a livello di vertice o di pixel, permettendo

l’implementazione di percorsi di codice molto complessi.

Le performance derivano invece dal fatto che il codice può venire eseguito solo sui

vertici o pixel che lo richiedono (anche se la performance guadagnata da codice non

eseguito può variare a seconda dell’implementazione hardware). Le istruzioni di branch

dinamico possono essere innestate fino a 24 livelli di ricorsione. Una tipica applicazione


______________________________________________________________________________________________

______________________________________________________________________________________________


del dynamic branching sta nel classico calcolo dell’illuminazione di una superficie. In

base al risultato del prodotto scalare N.L, dove N è la normale alla superficie e L è il

vettore che collega il punto sulla superficie stessa alla luce, possiamo sapere se il punto

è rivolto dalla parte opposta rispetto alla luce. In questo modo con un condizionale

possiamo scartare la pesante equazione per il calcolo dell’illuminazione che tiene conto

dell’attenuazione, del colore speculare, etc etc. a seconda del modello usato. Il seguente

esempio mostra questa cosa ps_3_0 ; costanti def c0, 0.0f, 0.0f, 0.0f, 1.0f ; Dichiariamo I vari registri di sampling dcl_2d s1 ; Normal map ; Dichiariamo l’input dcl_texcoord0 v0.xy ; Coordinate texture dcl_texcoord1 v1.xyz ; Vettore Punto->Luce texld r2, v0, s1 ; Leggiamo la normale alla superficie nel punto nrm r1, v1 ; Normalizziamo il vettore punto->Luce dp3 r0.w, r2, r1 ; Calcoliamo il fattore(N.L) if_gt r0.w, c0.x ; Se (N.L)>0 ; Calcoliamo il resto dell’equazione di illuminazione: colore ; speculare, attenuazione, lightmap, etc. ; r3 conterrà il valore finale di tutta questo calcolo else ; Altrimenti scriviamo nero (o un colore ambient) ; saltando tutta il complesso calcolo di illuminazione! mov r3, c0 endif mov oC0, r3 ; Scriviamo a schermo il colore Il dynamic branching può quindi servire ad eseguire forti ottimizzazioni. Bisogna solo

stare attenti a non usarlo con numeri di istruzioni molto piccoli, poiché risulterebbe più

lento che eseguendo tutte le istruzioni.

Le istruzioni di break all’interno dei loop hanno la stessa funzione, quando si è trovato

il giusto risultato si può stoppare l’iterazione. Come vedremo si farà uso di questa

istruzione esattamente nel nostro caso di volumetric rendering, per effettuare il ray-

tracing del campo scalare generato dalla simulazione per estrarne una iso-superficie.


______________________________________________________________________________________________

______________________________________________________________________________________________


Altre novità Esistono altre caratteristiche nuove dello shader model 3.0 ma sono state esposte solo

quelle fondamentali per la scrittura degli algoritmi necessari al tipo di rendering

utilizzato in questa sede. Il modello shader 3.0 è un grande passo in avanti rispetto al

2.0. Nuove caratteristiche sono state introdotte mentre i limiti di registri e numero di

istruzioni sono stati aumentati enormemente. La semplicità aggiunta dall’unificazione

dei modelli vertex e pixel shader permettono una maggiore flessibilità sulle istruzioni e

sui registri.


______________________________________________________________________________________________

______________________________________________________________________________________________


Conclusioni In questa tesi è stata presentata una nuova tecnica per la simulazione e rendering in

tempo reale di corpi deformabili utilizzando una ridotta potenza di calcolo. La

semplicità del modello, le ottimizzazioni nello schema di partizionamento e le strategie

di aggiornamento hanno portato la simulazione ai livelli del real-time su numero di

particelle molto alti per un normale PC.

Il collo di bottiglia attuale è la fase di rendering attraverso Marching Cubes, ma i

risultati dei primi test sul GPU rendering fanno ben sperare per soluzioni alternative

valide.

Il metodo qui presentato può anche essere esteso per utilizzare periferiche di force

feedback. In questo modo gli utenti possono sentire realisticamente la presenza fisica

del materiale virtuale nel corso del processo di modellazione.


______________________________________________________________________________________________

______________________________________________________________________________________________


Bibliografia

[1] Dewaele G., Cani M.-P.,, “Interactive Global and Local Deformations for Virtual

Clay”, in Pacific Graphics, Canmore, Canada, October 2003

[2] A. Raviv and G. Elber, “Three-dimensional freeform sculpting via zero sets of

scalar trivariate functions”, Computer-Aided Design, 32(8-9):513–526, August

2000.

[3] Arata H., Takai Y., Takai N. K., Yamamoto T.: “Free-form shape modeling by 3d

cellular automata”. In Shape Modeling International (1999), pp. 242–247.

[4] Druon S., A.Crosnier, Brigandat L.: “Efficient cellular automata for 2d / 3d free-

form modeling”. In WSCG (Feb. 2003).

[5] Ferley E., Cani M.-P., Gascuel J.-D.: “Practical volumetric sculpting”. the Visual

Computer 16, 8 (Dec. 2000), 469–480.

[6] IX F. D., El-Sana J., Qin H., Kaufman A.: “Haptic sculpting of dynamic

surfaces”. 1999 ACM Symposium on Interactive 3D Graphics (April 1999), 103–

110.

[7] McDonnell K. T., Qin H., Wlodarczyk R. A.: “Virtual clay: A real-time sculpting

system with haptic toolkits”. 2001 ACM Symposium on Interactive 3D Graphics

(2001).

[8] Onoue K., Nishita T.: “Virtual sandbox”. In Pacific Graphics 2003 (2003), pp.

252–259.

[9] Rossignac J.: “Finger sculpting with digital clay: 3d shape input and output

through a computer-controlled real surface”. In Shape Modeling International

(2003), pp. 229–231.

[10] T. W. Sederberg and S. R. Parry. “Free-form deformation of solid geometric

models”. In Computer Graphics (Proceedings of SIGGRAPH 86), volume 20,

pages 151–160, Aug. 1986.

[11] C. Shaw and M. Green., “THRED: A two-handed design system. Multimedia

Systems”, 5(2):126–139, 1997.

[12] J. R. Bill and S. Lodha. “Computer sculpting of polygonal models using virtual

tools”. In Tech. Rep. UCSC-CRL-94- 27, Baskin Center for Computer


______________________________________________________________________________________________

______________________________________________________________________________________________


Engineering and Information Sciences, University of California, Santa Cruz, US,

July 1994.

[13] Kuester, F., Duchaineau, M.A., Hamann, B., Joy, K.I. and Uva, A.E., (1999),

“3DIVS: 3-dimensional immersive virtual sculpting”, in Workshop on New

Paradigms in Information Visualization and Manipulation (NPIV '99), ACM

Press, New York, New York, pp. 92-95

[14] Thatcher U., “Notes on spatial partitioning”

http://tulrich.com/geekstuff/partitioning.html

[15] Fiorentino M., De Amicis R., Stork A., Monno G., “Spacedesign: A Mixed

Reality Workspace for Aesthetic Industrial Design”, Proc. of ISMAR 2002 IEEE

and ACM International Symposium on Mixed and Augmented Reality, pp. 86-94,

Sept. 30 - Oct. 1, 2002, Darmstadt, Germany.

[16] J. M. Haile, “Molecular Dynamics Simulation: Elementary Methods”, Wiley

Professional.

[17] N. J. A. Sloane, R. H. Hardin, T. S. Duff and J. H. Conway, “Minimal-Energy

Clusters of Hard Spheres” Discrete Computational Geom., 14 (1995), pp. 237-

259.

[18] D. L. Tonnesen. “Dynamically Coupled Particle Systems for Geometric Modelling, Reconstruction, and Animation”

[19] M. Matyka, M. Ollila "Pressure Model of Soft Body Simulation" [20] J. Krivanek “Representing and Rendering Surfaces with Points”

[21] D. B. Murray "Molecular Dynamics Simulation of an Elastic Material" [22] J. M. Haile, “Molecular Dynamics Simulation: Elementary Methods”, Wiley

Professional [23] Press, William H. et al, Numerical Recipes, Cambridge University Press, 1993.

http://www.nr.com/nronline_switcher.html [24] T. Jakobsen "Advanced Character Physics". GDC

[25] Verlet, L. "Computer experiments on classical fluids. I. Thermodynamical

properties of Lennard-Jones molecules", Phys. Rev., 159, 98-103 (1967).

[26] Stewart, D. E., and J. C. Trinkle, “An Implicit Time-Stepping Scheme for Rigid

Body Dynamics with Inelastic Collisions and Coulomb Friction”, International

Journal of Numerical Methods in Engineering

http://tulrich.com/geekstuff/partitioning.html

http://www.nr.com/nronline_switcher.html


______________________________________________________________________________________________

______________________________________________________________________________________________


[27] Eric W. Weisstein. "Cubic Close Packing." From MathWorld--A Wolfram Web

Resource. http://mathworld.wolfram.com/CubicClosePacking.html

[28] Tomas Akenine-Moller; Eric Haines. “Real-Time Rendering”

[29] T.H. Cormen, C.E. Leiserson; R.L. Rivest, C. Stain, "Introduction to

Algorithms, Second Edition". The MIT Press

[30] J. W. Ratcliff “Sphere Trees for Fast Visibility Culling, Ray Tracing, And

Range Searching” from Game Programming Gems

[31] R.Green "Spherical Harmonic Lighting: The Gritty Details"

[32] Andreas Jonsson. http://www.angelcode.com/

[33] http://www.gpgpu.org

[34] K. Engel, M. Hadwiger, J. Kniss, A. Lefohn, D. Weiskopf "Interactive

Visualization of Volumetric Data On Consumer PC Hardware.". Slides. IEEE

Visualization 2003 Tutorial

[35] D. Weiskopf “Flow Visualization”. Slides. IEEE Visualization 2003 Tutorial:

Interactive Visualization of Volumetric Data On Consumer PC Hardware

[36] J. Kniss “Transfer Functions: Optical Properties”. Slides. IEEE Visualization

2003 Tutorial: Interactive Visualization of Volumetric Data On Consumer PC

Hardware

[37] J. Kniss “Transfer Functions: Classification” Slides. IEEE Visualization 2003

Tutorial: Interactive Visualization of Volumetric Data On Consumer PC

Hardware

[38] M. Hadwiger “Introduction to GPU Volume Rendering: Part 1” Slides. IEEE

Visualization 2003 Tutorial: Interactive Visualization of Volumetric Data On

Consumer PC Hardware

[39] M. Hadwiger “Introduction to GPU Volume Rendering: Part 2” Slides. IEEE



[40] K. Engel “Basic Illumination Techniques”. Slides. IEEE Visualization 2003

Tutorial: Interactive Visualization of Volumetric Data On Consumer PC

Hardware

http://mathworld.wolfram.com/CubicClosePacking.html

http://www.angelcode.com/

http://www.gpgpu.org


______________________________________________________________________________________________

______________________________________________________________________________________________


[41] K. Engel “Advanced Volume Rendering Techniques”. Slides. IEEE



[42] E. Hart “3D Textures and Pixel Shaders” on “ShaderX. Vertex And Pixel Shader

Tips and Tricks” edited By W. F. Engel.

[43] M. Kraus “Truly Volumetric Effects” on “ShaderX. Vertex And Pixel Shader

Tips and Tricks” edited By W. F. Engel.”

[44] Robert A. Drebin, Loren Carpenter, and Pat Hanrahan, “Volume

Rendering,”SIGGRAPH Proceedings, Vol. 22, No. 4, August 1988, pp. 65-74.

[45] Klaus Engel, Martin Kraus, and Thomas Ertl, “High-Quality Pre-Integrated

Volume Rendering Using Hardware-Accelerated Pixel Shading,” SIGGRAPH

Proceedings/Eurographics Workshop on Graphics Hardware, August 2001, pp. 9-

16 (http://wwwvis.informatik.uni-stuttgart.de/~engel/pre-integrated/)

[46] N. Greene, M. Kass, G. Miller “Hierarchical Z-Buffer Visibility”. Computer

Graphics (SIGGRAPH 93 Proceedings), pp. 231-238, August 1993

[47] Arvo, James, “A simple Method for Box-Sphere Intersection Testing” in Andrei

S. Glassner, ed., Graphics Gems, Academic Press, pp. 335-339 1990.

http://www.graphicsgems.org/

[48] N. Thibieroz, K. Beets and A. Burton. “Introduction to the vs_3_0 and ps_3_0

Shader Models” on “ShaderX2 Introduction&Tutorials with DirectX9” edited By

W. Engel.

[49] T.J. Purcell “Ray Tracing on a Stream Processor” 2004

[50] Reinhard, Erik, Mike Stark, Peter Shirley, and Jim Ferwerda. "Photographic

Tone Reproduction for Digital Images" .

[51] Fiorentino M., De Amicis R., Stork A., Monno G., “Spacedesign: A Mixed

Reality Workspace for Aesthetic Industrial Design”, Proc. of ISMAR 2002 IEEE

and ACM International Symposium on Mixed and Augmented Reality, pp. 86-94,

Sept. 30 - Oct. 1, 2002, Darmstadt, Germany.

http://wwwvis.informatik.uni-stuttgart.de/~engel/pre-integrated/

http://www.graphicsgems.org/

Tesi - Simulazione e Rendering di materiali deformabili in ... · Questa caratteristica rende...

Documents

Transcript of Tesi - Simulazione e Rendering di materiali deformabili in ... · Questa caratteristica rende...