Visualizzazione dei dati 3D con strumenti open source: un tutorial che utilizza VTK

Nel suo recente articolo sul blog di Toptal, l'esperto data scientist Charles Cook ha scritto di calcolo scientifico con strumenti open source. Il suo tutorial fa un punto importante sugli strumenti open source e sul ruolo che possono svolgere nell'elaborazione semplice dei dati e nell'acquisizione dei risultati.

Ma non appena abbiamo risolto tutte queste complesse equazioni differenziali si presenta un altro problema. Come comprendiamo e interpretiamo le enormi quantità di dati che escono da queste simulazioni? Come si visualizzano potenziali gigabyte di dati, ad esempio dati con milioni di punti della griglia all'interno di una simulazione di grandi dimensioni?

Durante il mio lavoro su problemi simili per la mia tesi di laurea, sono entrato in contatto con il Visualization Toolkit, o VTK, una potente libreria grafica specializzata per la visualizzazione dei dati.

In questo tutorial darò una rapida introduzione a VTK e alla sua architettura della pipeline e proseguirò discutendo un esempio di visualizzazione 3D reale utilizzando i dati di un fluido simulato in una pompa a girante. Infine elencherò i punti di forza della libreria, così come i punti deboli che ho riscontrato.

Visualizzazione dei dati e pipeline VTK

La libreria open source VTK contiene una solida pipeline di elaborazione e rendering con molti sofisticati algoritmi di visualizzazione. Le sue capacità, tuttavia, non si fermano qui, poiché nel tempo sono stati aggiunti anche algoritmi di elaborazione di immagini e mesh. Nel mio attuale progetto con una società di ricerca dentale, sto utilizzando VTK per attività di elaborazione basate su mesh all'interno di un'applicazione simile a CAD basata su Qt. I casi di studio VTK mostrano l'ampia gamma di applicazioni adatte.

L'architettura di VTK ruota attorno a un potente concetto di pipeline. Lo schema di base di questo concetto è mostrato qui:

• Le fonti sono proprio all'inizio della pipeline e creano "qualcosa dal nulla". Ad esempio, un vtkConeSource crea un cono 3D e un vtkSTLReader legge i file di geometria 3D *.stl .
• I filtri trasformano l'output delle sorgenti o di altri filtri in qualcosa di nuovo. Ad esempio un vtkCutter taglia l'output dell'oggetto precedente negli algoritmi utilizzando una funzione implicita, ad esempio un piano. Tutti gli algoritmi di elaborazione forniti con VTK sono implementati come filtri e possono essere liberamente concatenati.
• I mappatori trasformano i dati in primitive grafiche. Ad esempio, possono essere utilizzati per specificare una tabella di ricerca per colorare i dati scientifici. Sono un modo astratto per specificare cosa visualizzare.
• Gli attori rappresentano un oggetto (geometria più proprietà di visualizzazione) all'interno della scena. Cose come colore, opacità, ombreggiatura o orientamento sono specificate qui.
• Renderer e Windows descrivono infine il rendering sullo schermo in modo indipendente dalla piattaforma.

Una tipica pipeline di rendering VTK inizia con una o più sorgenti, le elabora utilizzando vari filtri in diversi oggetti di output, che vengono quindi renderizzati separatamente utilizzando mappatori e attori. Il potere dietro questo concetto è il meccanismo di aggiornamento. Se le impostazioni dei filtri o delle sorgenti vengono modificate, tutti i filtri dipendenti, i mappatori, gli attori e le finestre di rendering vengono aggiornati automaticamente. Se, invece, un oggetto più in basso nella pipeline ha bisogno di informazioni per svolgere i suoi compiti, può facilmente ottenerle.

Inoltre, non è necessario occuparsi direttamente di sistemi di rendering come OpenGL. VTK incapsula tutte le attività di basso livello in modo indipendente dalla piattaforma e (parzialmente) dal sistema; lo sviluppatore lavora a un livello molto più alto.

Esempio di codice con un set di dati di una pompa a rotore

Diamo un'occhiata a un esempio di visualizzazione dei dati utilizzando un set di dati del flusso di fluido in una pompa a girante rotante dell'IEEE Visualization Contest 2011. I dati stessi sono il risultato di una simulazione di fluidodinamica computazionale, proprio come quella descritta nell'articolo di Charles Cook.

I dati di simulazione compressi della pompa in primo piano hanno una dimensione di oltre 30 GB. Contiene più parti e più passaggi temporali, da qui le grandi dimensioni. In questa guida, giocheremo con la parte del rotore di uno di questi timestep, che ha una dimensione compressa di circa 150 MB.

La mia lingua preferita per l'utilizzo di VTK è C++, ma ci sono mappature per molti altri linguaggi come Tcl/Tk, Java e Python. Se l'obiettivo è solo la visualizzazione di un singolo set di dati, non è necessario scrivere codice e si può invece utilizzare Paraview, un front-end grafico per la maggior parte delle funzionalità di VTK.

Il set di dati e perché 64 bit sono necessari

Ho estratto il set di dati del rotore dal set di dati da 30 GB fornito sopra, aprendo un timestep in Paraview ed estraendo la parte del rotore in un file separato. È un file griglia non strutturato, ovvero un volume 3D costituito da punti e celle 3D, come esaedri, tetraedri e così via. Ciascuno dei punti 3D ha valori associati. A volte le celle hanno anche valori associati, ma non in questo caso. Questo allenamento si concentrerà su pressione e velocità nei punti e cercherà di visualizzarli nel loro contesto 3D.

La dimensione del file compresso è di circa 150 MB e la dimensione in memoria è di circa 280 MB quando caricato con VTK. Tuttavia, elaborandolo in VTK, il set di dati viene memorizzato nella cache più volte all'interno della pipeline VTK e raggiungiamo rapidamente il limite di memoria di 2 GB per i programmi a 32 bit. Esistono modi per risparmiare memoria quando si utilizza VTK, ma per semplificare, compileremo ed eseguiremo l'esempio a 64 bit.

Ringraziamenti : Il set di dati è reso disponibile per gentile concessione dell'Istituto di Meccanica Applicata, Università di Clausthal, Germania (Dipl. Wirtsch.-Ing. Andreas Lucius).

Il bersaglio

Quello che otterremo usando VTK come strumento è la visualizzazione mostrata nell'immagine qui sotto. Come contesto 3D, il profilo del set di dati viene mostrato utilizzando un rendering wireframe parzialmente trasparente. La parte sinistra del set di dati viene quindi utilizzata per visualizzare la pressione utilizzando una semplice codifica a colori delle superfici. (Per questo esempio salteremo il rendering del volume più complesso). Per visualizzare il campo della velocità, la parte destra del set di dati è riempita di linee di flusso, codificate a colori in base all'entità della loro velocità. Questa scelta di visualizzazione non è tecnicamente l'ideale, ma volevo mantenere il codice VTK il più semplice possibile. Inoltre, c'è una ragione per cui questo esempio fa parte di una sfida di visualizzazione, cioè molta turbolenza nel flusso.

Passo dopo passo

Discuterò passo dopo passo il codice VTK, mostrando come apparirà l'output di rendering in ogni fase. Il codice sorgente completo può essere scaricato al termine del corso.

Iniziamo includendo tutto ciò di cui abbiamo bisogno da VTK e apriamo la funzione principale.

 #include <vtkActor.h> #include <vtkArrayCalculator.h> #include <vtkCamera.h> #include <vtkClipDataSet.h> #include <vtkCutter.h> #include <vtkDataSetMapper.h> #include <vtkInteractorStyleTrackballCamera.h> #include <vtkLookupTable.h> #include <vtkNew.h> #include <vtkPlane.h> #include <vtkPointData.h> #include <vtkPointSource.h> #include <vtkPolyDataMapper.h> #include <vtkProperty.h> #include <vtkRenderer.h> #include <vtkRenderWindow.h> #include <vtkRenderWindowInteractor.h> #include <vtkRibbonFilter.h> #include <vtkStreamTracer.h> #include <vtkSmartPointer.h> #include <vtkUnstructuredGrid.h> #include <vtkXMLUnstructuredGridReader.h> int main(int argc, char** argv) {

Successivamente, impostiamo il renderer e la finestra di rendering per visualizzare i nostri risultati. Impostiamo il colore di sfondo e la dimensione della finestra di rendering.

 // Setup the renderer vtkNew<vtkRenderer> renderer; renderer->SetBackground(0.9, 0.9, 0.9); // Setup the render window vtkNew<vtkRenderWindow> renWin; renWin->AddRenderer(renderer.Get()); renWin->SetSize(500, 500);

Con questo codice potremmo già visualizzare una finestra di rendering statica. Invece, scegliamo di aggiungere un vtkRenderWindowInteractor per ruotare, ingrandire e fare una panoramica della scena in modo interattivo.

 // Setup the render window interactor vtkNew<vtkRenderWindowInteractor> interact; vtkNew<vtkInteractorStyleTrackballCamera> style; interact->SetRenderWindow(renWin.Get()); interact->SetInteractorStyle(style.Get());

Ora abbiamo un esempio in esecuzione che mostra una finestra di rendering grigia e vuota.

Successivamente, carichiamo il set di dati utilizzando uno dei tanti lettori forniti con VTK.

 // Read the file vtkSmartPointer<vtkXMLUnstructuredGridReader> pumpReader = vtkSmartPointer<vtkXMLUnstructuredGridReader>::New(); pumpReader->SetFileName("rotor.vtu");

Breve escursione nella gestione della memoria VTK : VTK utilizza un comodo concetto di gestione automatica della memoria che ruota attorno al conteggio dei riferimenti. Diversamente dalla maggior parte delle altre implementazioni, tuttavia, il conteggio dei riferimenti viene mantenuto all'interno degli oggetti VTK stessi, anziché nella classe del puntatore intelligente. Questo ha il vantaggio che il conteggio dei riferimenti può essere aumentato, anche se l'oggetto VTK viene passato come un puntatore non elaborato. Esistono due modi principali per creare oggetti VTK gestiti. vtkNew<T> e vtkSmartPointer<T>::New() , con la differenza principale che un vtkSmartPointer<T> è in grado di eseguire il cast implicito per il puntatore non elaborato T* e può essere restituito da una funzione. Per le istanze di vtkNew<T> dovremo chiamare .Get() per ottenere un puntatore non elaborato e possiamo restituirlo solo avvolgendolo in un vtkSmartPointer . All'interno del nostro esempio, non torniamo mai dalle funzioni e tutti gli oggetti risiedono per tutto il tempo, quindi useremo l'abbreviazione vtkNew , con solo l'eccezione sopra a scopo dimostrativo.

A questo punto, non è stato ancora letto nulla dal file. Noi o un filtro più in basso nella catena dovremmo chiamare Update() affinché la lettura del file avvenga effettivamente. Di solito è l'approccio migliore lasciare che le classi VTK gestiscano gli aggiornamenti da sole. Tuttavia, a volte vogliamo accedere direttamente al risultato di un filtro, ad esempio per ottenere l'intervallo di pressioni in questo set di dati. Quindi dobbiamo chiamare Update() manualmente. (Non perdiamo le prestazioni chiamando Update() più volte, poiché i risultati vengono memorizzati nella cache.)

 // Get the pressure range pumpReader->Update(); double pressureRange[2]; pumpReader->GetOutput()->GetPointData()->GetArray("Pressure")->GetRange(pressureRange);

Successivamente, è necessario estrarre la metà sinistra del set di dati, utilizzando vtkClipDataSet . Per ottenere questo, definiamo prima un vtkPlane che definisce lo split. Quindi, vedremo per la prima volta come la pipeline VTK è collegata insieme: successor->SetInputConnection(predecessor->GetOutputPort()) . Ogni volta che richiediamo un aggiornamento da clipperLeft , questa connessione ora garantirà che anche tutti i filtri precedenti siano aggiornati.

 // Clip the left part from the input vtkNew<vtkPlane> planeLeft; planeLeft->SetOrigin(0.0, 0.0, 0.0); planeLeft->SetNormal(-1.0, 0.0, 0.0); vtkNew<vtkClipDataSet> clipperLeft; clipperLeft->SetInputConnection(pumpReader->GetOutputPort()); clipperLeft->SetClipFunction(planeLeft.Get());

Infine, creiamo i nostri primi attori e mappatori per visualizzare il rendering wireframe della metà sinistra. Si noti che il mapper è collegato al suo filtro esattamente nello stesso modo dei filtri tra loro. Il più delle volte, il renderer stesso attiva gli aggiornamenti di tutti gli attori, i mappatori e le catene di filtri sottostanti!

L'unica riga che non è autoesplicativa è probabilmente leftWireMapper->ScalarVisibilityOff(); - vieta la colorazione del wireframe per valori di pressione, che sono impostati come array attualmente attivo.

 // Create the wireframe representation for the left part vtkNew<vtkDataSetMapper> leftWireMapper; leftWireMapper->SetInputConnection(clipperLeft->GetOutputPort()); leftWireMapper->ScalarVisibilityOff(); vtkNew<vtkActor> leftWireActor; leftWireActor->SetMapper(leftWireMapper.Get()); leftWireActor->GetProperty()->SetRepresentationToWireframe(); leftWireActor->GetProperty()->SetColor(0.8, 0.8, 0.8); leftWireActor->GetProperty()->SetLineWidth(0.5); leftWireActor->GetProperty()->SetOpacity(0.8); renderer->AddActor(leftWireActor.Get());

A questo punto, la finestra di rendering mostra finalmente qualcosa, ovvero il wireframe per la parte sinistra.

Il rendering wireframe per la parte destra viene creato in modo simile, spostando la normale del piano di un vtkClipDataSet (appena creato) nella direzione opposta e modificando leggermente il colore e l'opacità del mapper e dell'attore (appena creati). Si noti che qui la nostra pipeline VTK si divide in due direzioni (destra e sinistra) dallo stesso set di dati di input.

 // Clip the right part from the input vtkNew<vtkPlane> planeRight; planeRight->SetOrigin(0.0, 0.0, 0.0); planeRight->SetNormal(1.0, 0.0, 0.0); vtkNew<vtkClipDataSet> clipperRight; clipperRight->SetInputConnection(pumpReader->GetOutputPort()); clipperRight->SetClipFunction(planeRight.Get()); // Create the wireframe representation for the right part vtkNew<vtkDataSetMapper> rightWireMapper; rightWireMapper->SetInputConnection(clipperRight->GetOutputPort()); rightWireMapper->ScalarVisibilityOff(); vtkNew<vtkActor> rightWireActor; rightWireActor->SetMapper(rightWireMapper.Get()); rightWireActor->GetProperty()->SetRepresentationToWireframe(); rightWireActor->GetProperty()->SetColor(0.2, 0.2, 0.2); rightWireActor->GetProperty()->SetLineWidth(0.5); rightWireActor->GetProperty()->SetOpacity(0.1); renderer->AddActor(rightWireActor.Get());

La finestra di output ora mostra entrambe le parti wireframe, come previsto.

Ora siamo pronti per visualizzare alcuni dati utili! Per aggiungere la visualizzazione della pressione alla parte sinistra, non è necessario fare molto. Creiamo un nuovo mappatore e lo colleghiamo anche a clipperLeft , ma questa volta coloriamo in base all'array di pressione. È anche qui che finalmente utilizziamo il pressureRange che abbiamo derivato sopra.

 // Create the pressure representation for the left part vtkNew<vtkDataSetMapper> pressureColorMapper; pressureColorMapper->SetInputConnection(clipperLeft->GetOutputPort()); pressureColorMapper->SelectColorArray("Pressure"); pressureColorMapper->SetScalarRange(pressureRange); vtkNew<vtkActor> pressureColorActor; pressureColorActor->SetMapper(pressureColorMapper.Get()); pressureColorActor->GetProperty()->SetOpacity(0.5); renderer->AddActor(pressureColorActor.Get());

L'output ora appare come l'immagine mostrata di seguito. La pressione al centro è molto bassa, aspirando materiale nella pompa. Quindi, questo materiale viene trasportato all'esterno, guadagnando rapidamente pressione. (Ovviamente dovrebbe esserci una legenda della mappa dei colori con i valori effettivi, ma l'ho omessa per mantenere l'esempio più breve.)

Ora inizia la parte più complicata. Vogliamo disegnare linee di flusso di velocità nella parte giusta. Le linee di flusso sono generate dall'integrazione all'interno di un campo vettoriale da punti di origine. Il campo vettoriale fa già parte del set di dati sotto forma di vettore-array "Velocities". Quindi abbiamo solo bisogno di generare i punti sorgente. vtkPointSource genera una sfera di punti casuali. Genereremo 1500 punti sorgente, perché la maggior parte di essi non si troverà comunque all'interno del set di dati e verrà ignorata dal tracciante del flusso.

 // Create the source points for the streamlines vtkNew<vtkPointSource> pointSource; pointSource->SetCenter(0.0, 0.0, 0.015); pointSource->SetRadius(0.2); pointSource->SetDistributionToUniform(); pointSource->SetNumberOfPoints(1500);

Quindi creiamo lo streamtracer e impostiamo le sue connessioni di input. "Aspetta, più connessioni?", potresti dire. Sì, questo è il primo filtro VTK con ingressi multipli che incontriamo. La normale connessione di ingresso viene utilizzata per il campo vettoriale e la connessione di origine viene utilizzata per i punti seme. Poiché "Velocities" è l'array vettoriale "attivo" in clipperRight , non è necessario specificarlo qui in modo esplicito. Infine specifichiamo che l'integrazione deve essere eseguita in entrambe le direzioni dai seed point e impostiamo il metodo di integrazione su Runge-Kutta-4.5.

 vtkNew<vtkStreamTracer> tracer; tracer->SetInputConnection(clipperRight->GetOutputPort()); tracer->SetSourceConnection(pointSource->GetOutputPort()); tracer->SetIntegrationDirectionToBoth(); tracer->SetIntegratorTypeToRungeKutta45();

Il nostro prossimo problema è colorare le linee di flusso in base all'intensità della velocità. Poiché non esiste un array per le grandezze dei vettori, calcoleremo semplicemente le grandezze in un nuovo array scalare. Come avrai intuito, esiste anche un filtro VTK per questa attività: vtkArrayCalculator . Prende un set di dati e lo restituisce invariato, ma aggiunge esattamente un array calcolato da uno o più di quelli esistenti. Configuriamo questo calcolatore di array in modo che prenda la grandezza del vettore "Velocity" e lo emetta come "MagVelocity". Infine, chiamiamo di nuovo Update() manualmente, per ricavare l'intervallo del nuovo array.

 // Compute the velocity magnitudes and create the ribbons vtkNew<vtkArrayCalculator> magCalc; magCalc->SetInputConnection(tracer->GetOutputPort()); magCalc->AddVectorArrayName("Velocity"); magCalc->SetResultArrayName("MagVelocity"); magCalc->SetFunction("mag(Velocity)"); magCalc->Update(); double magVelocityRange[2]; magCalc->GetOutput()->GetPointData()->GetArray("MagVelocity")->GetRange(magVelocityRange);

vtkStreamTracer emette direttamente le polilinee e vtkArrayCalculator le passa invariate. Pertanto potremmo semplicemente visualizzare l'output di magCalc direttamente utilizzando un nuovo mappatore e attore.

Invece, in questa formazione scegliamo di rendere l'output un po' più gradevole, visualizzando invece dei nastri. vtkRibbonFilter genera celle 2D per visualizzare i nastri per tutte le polilinee del suo input.

 // Create and render the ribbons vtkNew<vtkRibbonFilter> ribbonFilter; ribbonFilter->SetInputConnection(magCalc->GetOutputPort()); ribbonFilter->SetWidth(0.0005); vtkNew<vtkPolyDataMapper> streamlineMapper; streamlineMapper->SetInputConnection(ribbonFilter->GetOutputPort()); streamlineMapper->SelectColorArray("MagVelocity"); streamlineMapper->SetScalarRange(magVelocityRange); vtkNew<vtkActor> streamlineActor; streamlineActor->SetMapper(streamlineMapper.Get()); renderer->AddActor(streamlineActor.Get());

Ciò che ora manca ancora, ed è effettivamente necessario per produrre anche i rendering intermedi, sono le ultime cinque righe per rendere effettivamente la scena e inizializzare l'interattore.

 // Render and show interactive window renWin->Render(); interact->Initialize(); interact->Start(); return 0; }

Infine, arriviamo alla visualizzazione finita, che presenterò ancora una volta qui:

Il codice sorgente completo per la visualizzazione di cui sopra può essere trovato qui.

Il buono il brutto e il cattivo

Chiuderò questo articolo con un elenco dei miei pro e contro personali del framework VTK.

• Pro : Sviluppo attivo : VTK è in fase di sviluppo attivo da parte di diversi contributori, principalmente all'interno della comunità di ricerca. Ciò significa che sono disponibili alcuni algoritmi all'avanguardia, molti formati 3D possono essere importati ed esportati, i bug vengono risolti attivamente e i problemi di solito hanno una soluzione già pronta nei forum di discussione.

• Contro : Affidabilità : L'accoppiamento di molti algoritmi di diversi contributori con il design della pipeline aperta di VTK, tuttavia, può portare a problemi con combinazioni di filtri insolite. Ho dovuto entrare nel codice sorgente VTK alcune volte per capire perché la mia complessa catena di filtri non sta producendo i risultati desiderati. Consiglio vivamente di configurare VTK in un modo che consenta il debug.

• Pro : Architettura del software : il design della pipeline e l'architettura generale di VTK sembrano ben congegnati ed è un piacere lavorare con loro. Poche righe di codice possono produrre risultati sorprendenti. Le strutture dati integrate sono facili da capire e utilizzare.

• Contro : Microarchitettura : Alcune decisioni di progettazione micro-architettura sfuggono alla mia comprensione. La correttezza di cost è quasi inesistente, gli array vengono passati come input e output senza una chiara distinzione. Ho alleviato questo problema per i miei algoritmi rinunciando ad alcune prestazioni e usando il mio wrapper per vtkMath che utilizza tipi 3D personalizzati come typedef std::array<double, 3> Pnt3d; .

• Pro : Micro Documentazione : La documentazione Doxygen di tutte le classi e filtri è ampia e utilizzabile, gli esempi e i casi di test sul wiki sono anche un grande aiuto per capire come vengono utilizzati i filtri.

• Contro : Documentazione Macro : Ci sono molti buoni tutorial e introduzioni a VTK sul web. Tuttavia, per quanto ne so, non esiste una grande documentazione di riferimento che spieghi come vengono fatte le cose specifiche. Se vuoi fare qualcosa di nuovo, aspettati di cercare come farlo per un po' di tempo. Inoltre è difficile trovare il filtro specifico per un'attività. Una volta trovato, tuttavia, la documentazione di Doxygen di solito sarà sufficiente. Un buon modo per esplorare il framework VTK è scaricare e sperimentare con Paraview.

• Pro : supporto per la parallelizzazione implicita : se i tuoi sorgenti possono essere suddivisi in più parti che possono essere elaborate indipendentemente, la parallelizzazione è semplice come creare una catena di filtri separata all'interno di ogni thread che elabora una singola parte. La maggior parte dei grandi problemi di visualizzazione di solito rientrano in questa categoria.

• Contro : Nessun supporto per la parallelizzazione esplicita : se non sei benedetto da problemi grandi e divisibili, ma desideri utilizzare più core, sei da solo. Dovrai capire quali classi sono thread-safe o addirittura rientranti per tentativi ed errori o leggendo la fonte. Una volta ho rintracciato un problema di parallelizzazione su un filtro VTK che utilizzava una variabile globale statica per chiamare alcune librerie C.

• Pro : Buildsystem CMake : anche il meta-sistema di build multipiattaforma CMake è sviluppato da Kitware (i creatori di VTK) e utilizzato in molti progetti al di fuori di Kitware. Si integra molto bene con VTK e rende molto meno doloroso l'impostazione di un sistema di compilazione per più piattaforme.

• Pro : Indipendenza dalla piattaforma, licenza e longevità : VTK è indipendente dalla piattaforma e viene concesso in licenza con una licenza in stile BSD molto permissiva. Inoltre, è disponibile un supporto professionale per quei progetti importanti che lo richiedono. Kitware è supportato da molti enti di ricerca e altre società e sarà in circolazione per un po' di tempo.

Ultima parola

Nel complesso, VTK è il miglior strumento di visualizzazione dei dati per i tipi di problemi che amo. Se ti imbatti in un progetto che richiede visualizzazione, elaborazione mesh, elaborazione di immagini o attività simili, prova ad avviare Paraview con un esempio di input e valuta se VTK potrebbe essere lo strumento per te.