Vizualizarea datelor 3D cu instrumente open source: un tutorial folosind VTK

În articolul său recent de pe blogul lui Toptal, expertul în știință de date Charles Cook a scris despre calculul științific cu instrumente open source. Tutorialul său face un punct important despre instrumentele open source și rolul pe care acestea îl pot juca în procesarea cu ușurință a datelor și obținerea de rezultate.

Dar, de îndată ce am rezolvat toate aceste ecuații diferențiale complexe, apare o altă problemă. Cum înțelegem și interpretăm cantitățile uriașe de date care ies din aceste simulări? Cum vizualizăm potențiali gigaocteți de date, cum ar fi datele cu milioane de puncte de grilă într-o simulare mare?

În timpul lucrului meu pe probleme similare pentru lucrarea mea de master, am intrat în contact cu Visualization Toolkit, sau VTK - o bibliotecă grafică puternică specializată pentru vizualizarea datelor.

În acest tutorial, voi oferi o introducere rapidă în VTK și arhitectura conductei sale și voi continua să discut despre un exemplu de vizualizare 3D real, folosind date dintr-un fluid simulat într-o pompă cu rotor. În cele din urmă, voi enumera punctele tari ale bibliotecii, precum și punctele slabe pe care le-am întâlnit.

Vizualizarea datelor și conducta VTK

Biblioteca open source VTK conține o conductă solidă de procesare și randare cu mulți algoritmi de vizualizare sofisticați. Cu toate acestea, capacitățile sale nu se opresc aici, deoarece de-a lungul timpului au fost adăugați și algoritmi de procesare a imaginilor și a rețelei. În proiectul meu actual cu o companie de cercetare stomatologică, folosesc VTK pentru sarcini de procesare bazate pe ochiuri într-o aplicație CAD bazată pe Qt. Studiile de caz VTK arată o gamă largă de aplicații adecvate.

Arhitectura VTK se învârte în jurul unui concept puternic de conductă. Schița de bază a acestui concept este prezentată aici:

• Sursele sunt chiar la începutul conductei și creează „ceva din nimic”. De exemplu, un vtkConeSource creează un con 3D, iar un vtkSTLReader citește fișierele de geometrie 3D *.stl .
• Filtrele transformă ieșirea fie a surselor, fie a altor filtre în ceva nou. De exemplu, un vtkCutter taie ieșirea obiectului anterior în algoritmi folosind o funcție implicită, de exemplu, un plan. Toți algoritmii de procesare care vin cu VTK sunt implementați ca filtre și pot fi înlănțuiți liber împreună.
• Cartografele transformă datele în primitive grafice. De exemplu, ele pot fi folosite pentru a specifica un tabel de căutare pentru colorarea datelor științifice. Ele sunt o modalitate abstractă de a specifica ceea ce trebuie afișat.
• Actorii reprezintă un obiect (geometrie plus proprietăți de afișare) în cadrul scenei. Lucruri precum culoarea, opacitatea, umbrirea sau orientarea sunt specificate aici.
• Renderers și Windows descriu în cele din urmă redarea pe ecran într-un mod independent de platformă.

O conductă tipică de randare VTK începe cu una sau mai multe surse, le procesează folosind diferite filtre în mai multe obiecte de ieșire, care sunt apoi randate separat folosind mapper și actori. Puterea din spatele acestui concept este mecanismul de actualizare. Dacă setările filtrelor sau surselor sunt modificate, toate filtrele, mapperele, actorii și ferestrele de randare dependente sunt actualizate automat. Dacă, pe de altă parte, un obiect mai jos în conductă are nevoie de informații pentru a-și îndeplini sarcinile, le poate obține cu ușurință.

În plus, nu este nevoie să te ocupi direct de sisteme de randare precum OpenGL. VTK încapsulează toate sarcinile de nivel scăzut într-un mod independent de platformă și (parțial) de sistem; dezvoltatorul lucrează la un nivel mult mai înalt.

Exemplu de cod cu un set de date pompe cu rotor

Să ne uităm la un exemplu de vizualizare a datelor folosind un set de date privind fluxul de fluid într-o pompă cu rotor rotativ de la IEEE Visualization Contest 2011. Datele în sine sunt rezultatul unei simulări computaționale de dinamică a fluidelor, la fel ca cea descrisă în articolul lui Charles Cook.

Datele de simulare arhivate ale pompei prezentate au o dimensiune de peste 30 GB. Conține mai multe părți și mai mulți pași de timp, de unde și dimensiunea mare. În acest ghid, ne vom juca cu partea rotorului a unuia dintre acești pași de timp, care are o dimensiune comprimată de aproximativ 150 MB.

Limbajul meu ales pentru utilizarea VTK este C++, dar există mapări pentru mai multe alte limbi precum Tcl/Tk, Java și Python. Dacă ținta este doar vizualizarea unui singur set de date, nu este nevoie să scrieți cod deloc și, în schimb, puteți utiliza Paraview, un front-end grafic pentru majoritatea funcționalității VTK.

Setul de date și de ce este necesar 64 de biți

Am extras setul de date rotor din setul de date de 30 GB furnizat mai sus, deschizând un interval de timp în Paraview și extragând partea rotorului într-un fișier separat. Este un fișier grilă nestructurat, adică un volum 3D format din puncte și celule 3D, cum ar fi hexaedre, tetraedre și așa mai departe. Fiecare dintre punctele 3D are valori asociate. Uneori celulele au și valori asociate, dar nu în acest caz. Acest antrenament se va concentra pe presiunea și viteza în puncte și va încerca să le vizualizeze în contextul lor 3D.

Dimensiunea fișierului comprimat este de aproximativ 150 MB și dimensiunea în memorie este de aproximativ 280 MB când este încărcat cu VTK. Cu toate acestea, procesându-l în VTK, setul de date este stocat în cache de mai multe ori în conducta VTK și atingem rapid limita de memorie de 2 GB pentru programele pe 32 de biți. Există modalități de a economisi memorie când folosiți VTK, dar pentru a rămâne simplu, vom compila și rula exemplul pe 64 de biți.

Mulțumiri : Setul de date este pus la dispoziție prin amabilitatea Institutului de Mecanică Aplicată, Universitatea Clausthal, Germania (Dipl. Wirtsch.-Ing. Andreas Lucius).

Ținta

Ceea ce vom realiza folosind VTK ca instrument este vizualizarea prezentată în imaginea de mai jos. Ca context 3D, conturul setului de date este afișat folosind o redare wireframe parțial transparentă. Partea din stânga a setului de date este apoi utilizată pentru a afișa presiunea utilizând codificarea simplă prin culori a suprafețelor. (Vom sări peste randarea volumului mai complexă pentru acest exemplu). Pentru a vizualiza câmpul de viteză, partea dreaptă a setului de date este umplută cu linii fluide, care sunt codificate în culori după mărimea vitezei lor. Această alegere de vizualizare nu este ideală din punct de vedere tehnic, dar am vrut să păstrez codul VTK cât mai simplu posibil. În plus, există un motiv pentru ca acest exemplu să facă parte dintr-o provocare de vizualizare, adică multă turbulență în flux.

Pas cu pas

Voi discuta codul VTK pas cu pas, arătând cum ar arăta rezultatul de randare în fiecare etapă. Codul sursă complet poate fi descărcat la sfârșitul cursului.

Să începem prin a include tot ce avem nevoie de la VTK și să deschidem funcția principală.

 #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) {

Apoi, configuram redarea și fereastra de randare pentru a afișa rezultatele noastre. Setăm culoarea de fundal și dimensiunea ferestrei de randare.

 // 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);

Cu acest cod am putea deja afișa o fereastră de randare statică. În schimb, alegem să adăugăm un vtkRenderWindowInteractor pentru a roti, mări și deplasăm scena în mod interactiv.

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

Acum avem un exemplu de rulare care arată o fereastră de randare gri, goală.

În continuare, încărcăm setul de date folosind unul dintre numeroasele cititoare care vin cu VTK.

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

Scurtă excursie în managementul memoriei VTK : VTK utilizează un concept convenabil de gestionare automată a memoriei care se învârte în jurul numărării referințelor. Totuși, diferit de majoritatea celorlalte implementări, numărul de referințe este păstrat în obiectele VTK înseși, în locul clasei pointer inteligent. Acest lucru are avantajul că numărul de referințe poate fi crescut, chiar dacă obiectul VTK este transmis ca un pointer brut. Există două moduri majore de a crea obiecte VTK gestionate. vtkNew<T> și vtkSmartPointer<T>::New() , cu principala diferență fiind că un vtkSmartPointer<T> poate fi transmis implicit către pointerul brut T* și poate fi returnat de la o funcție. Pentru instanțe de vtkNew<T> , va trebui să apelăm .Get() pentru a obține un pointer brut și îl putem returna doar prin împachetarea lui într-un vtkSmartPointer . În exemplul nostru, nu ne întoarcem niciodată de la funcții și toate obiectele trăiesc tot timpul, prin urmare vom folosi scurtul vtkNew , cu doar excepția de mai sus în scopuri demonstrative.

În acest moment, încă nu a fost citit nimic din fișier. Noi sau un filtru mai jos în lanț ar trebui să apelăm Update() pentru ca citirea fișierului să se întâmple cu adevărat. De obicei, este cea mai bună abordare de a lăsa clasele VTK să gestioneze actualizările. Cu toate acestea, uneori dorim să accesăm direct rezultatul unui filtru, de exemplu pentru a obține intervalul de presiuni din acest set de date. Apoi trebuie să apelăm manual Update() . (Nu pierdem performanța apelând Update() de mai multe ori, deoarece rezultatele sunt stocate în cache.)

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

Apoi, trebuie să extragem jumătatea stângă a setului de date, folosind vtkClipDataSet . Pentru a realiza acest lucru, definim mai întâi un vtkPlane care definește diviziunea. Apoi, vom vedea pentru prima dată cum este conectată conducta VTK: successor->SetInputConnection(predecessor->GetOutputPort()) . Ori de câte ori solicităm o actualizare de la clipperLeft , această conexiune se va asigura acum că toate filtrele precedente sunt, de asemenea, actualizate.

 // 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());

În cele din urmă, creăm primii noștri actori și mapper pentru a afișa redarea wireframe a jumătății din stânga. Observați că mapper-ul este conectat la filtrul său exact în același mod ca filtrele unul la celălalt. De cele mai multe ori, redarea în sine declanșează actualizările tuturor actorilor, cartografilor și lanțurilor de filtre care stau la baza!

Singura linie care nu se explică de la sine este probabil leftWireMapper->ScalarVisibilityOff(); - interzice colorarea wireframe-ului prin valori de presiune, care sunt setate ca matrice activă în prezent.

 // 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());

În acest moment, fereastra de randare arată în sfârșit ceva, adică, wireframe pentru partea din stânga.

Redarea wireframe pentru partea dreaptă este creată într-un mod similar, prin comutarea normală a planului unui vtkClipDataSet (nou creat) în direcția opusă și modificarea ușor a culorii și opacității mapperului și actorului (nou creat). Observați că aici conducta noastră VTK se împarte în două direcții (dreapta și stânga) din același set de date de intrare.

 // 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());

Fereastra de ieșire arată acum ambele părți wireframe, așa cum era de așteptat.

Acum suntem gata să vizualizăm câteva date utile! Pentru a adăuga vizualizarea presiunii în partea stângă, nu trebuie să facem mare lucru. Creăm un nou mapper și îl conectăm și la clipperLeft , dar de data aceasta colorăm după matricea de presiune. Tot aici folosim în sfârșit intervalul de pressureRange pe care l-am derivat mai sus.

 // 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());

Ieșirea arată acum ca imaginea de mai jos. Presiunea la mijloc este foarte scăzută, aspirând material în pompă. Apoi, acest material este transportat spre exterior, câștigând rapid presiune. (Desigur, ar trebui să existe o legendă a hărții de culori cu valorile reale, dar am omis-o pentru a menține exemplul mai scurt.)

Acum începe partea mai dificilă. Vrem să desenăm linii de viteză în partea dreaptă. Liniile de fluidizare sunt generate prin integrarea într-un câmp vectorial din punctele sursă. Câmpul vectorial este deja parte a setului de date sub forma matricei vectoriale „Velocities”. Deci trebuie doar să generăm punctele sursă. vtkPointSource generează o sferă de puncte aleatorii. Vom genera 1500 de puncte sursă, deoarece cele mai multe dintre ele nu se vor afla oricum în setul de date și vor fi ignorate de urmăritorul fluxului.

 // 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);

Apoi, creăm streamtracer-ul și îi setăm conexiunile de intrare. „Stai, conexiuni multiple ?”, ai putea spune. Da - acesta este primul filtru VTK cu intrări multiple pe care îl întâlnim. Conexiunea normală de intrare este utilizată pentru câmpul vectorial, iar conexiunea sursă este folosită pentru punctele de inițiere. Deoarece „Velocities” este matricea vectorială „activă” în clipperRight , nu trebuie să-l specificăm aici în mod explicit. În cele din urmă precizăm că integrarea trebuie efectuată în ambele direcții de la punctele semințe și setăm metoda de integrare la Runge-Kutta-4.5.

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

Următoarea noastră problemă este colorarea liniilor de curgere după mărimea vitezei. Deoarece nu există o matrice pentru mărimile vectorilor, vom calcula pur și simplu mărimile într-o nouă matrice scalară. După cum ați ghicit, există și un filtru VTK pentru această sarcină: vtkArrayCalculator . Preia un set de date și îl scoate neschimbat, dar adaugă exact o matrice care este calculată din unul sau mai multe dintre cele existente. Configuram acest calculator matrice pentru a lua magnitudinea vectorului „Velocity” și a-l scoate ca „MagVelocity”. În cele din urmă, apelăm manual Update() din nou, pentru a deriva intervalul noii matrice.

 // 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 scoate direct polilinii și vtkArrayCalculator le transmite neschimbate. Prin urmare, am putea afișa direct rezultatul magCalc folosind un nou mapper și actor.

În schimb, în ​​acest antrenament optăm să facem rezultatul puțin mai plăcut, afișând în schimb panglici. vtkRibbonFilter generează celule 2D pentru a afișa panglici pentru toate poliliniile de intrare.

 // 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());

Ceea ce acum lipsește și este de fapt necesar pentru a produce și randările intermediare, sunt ultimele cinci linii pentru a reda scena și a inițializa interactorul.

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

În sfârșit, ajungem la vizualizarea finalizată, pe care o voi prezenta încă o dată aici:

Codul sursă complet pentru vizualizarea de mai sus poate fi găsit aici.

Cel bun cel rau si cel urat

Voi încheia acest articol cu ​​o listă cu avantajele și dezavantajele mele personale ale cadrului VTK.

• Pro : Dezvoltare activă : VTK este în curs de dezvoltare activă din partea mai multor colaboratori, în principal din cadrul comunității de cercetare. Aceasta înseamnă că sunt disponibili unii algoritmi de ultimă generație, multe formate 3D pot fi importate și exportate, erorile sunt remediate în mod activ și problemele au de obicei o soluție gata făcută în forumurile de discuții.

• Contra : Fiabilitate : Cuplarea multor algoritmi de la diferiți contribuitori cu designul conductei deschise a VTK poate duce totuși la probleme cu combinațiile neobișnuite de filtre. A trebuit să intru de câteva ori în codul sursă VTK pentru a-mi da seama de ce lanțul meu complex de filtre nu produce rezultatele dorite. Aș recomanda cu tărie configurarea VTK într-un mod care să permită depanarea.

• Pro : Arhitectura software : Designul conductei și arhitectura generală a VTK par bine gândite și este o plăcere să lucrezi. Câteva linii de cod pot produce rezultate uimitoare. Structurile de date încorporate sunt ușor de înțeles și de utilizat.

• Contra : Micro-Arhitectură : Unele decizii de proiectare micro-arhitecturale îmi scapă înțelegerii. Corectitudinea constantă este aproape inexistentă, tablourile sunt transmise ca intrări și ieșiri fără nicio distincție clară. Am atenuat acest lucru pentru proprii algoritmi, renunțând la anumite performanțe și folosind propriul meu wrapper pentru vtkMath , care utilizează tipuri 3D personalizate, cum ar fi typedef std::array<double, 3> Pnt3d; .

• Pro : Micro Documentație : Documentația Doxygen pentru toate clasele și filtrele este extinsă și utilizabilă, exemplele și cazurile de testare de pe wiki sunt, de asemenea, de mare ajutor pentru a înțelege cum sunt utilizate filtrele.

• Contra : Documentație macro : Există câteva tutoriale bune pentru și introduceri în VTK pe web. Cu toate acestea, din câte știu, nu există o documentație mare de referință care să explice cum se fac anumite lucruri. Dacă doriți să faceți ceva nou, așteptați-vă să căutați cum să o faceți pentru ceva timp. În plus, este greu să găsești filtrul specific pentru o sarcină. Cu toate acestea, odată ce l-ați găsit, documentația Doxygen va fi de obicei suficientă. O modalitate bună de a explora cadrul VTK este să descărcați și să experimentați cu Paraview.

• Pro : Suport pentru paralelizare implicită : Dacă sursele dvs. pot fi împărțite în mai multe părți care pot fi procesate independent, paralelizarea este la fel de simplă ca și crearea unui lanț de filtre separat în fiecare fir care procesează o singură parte. Majoritatea problemelor mari de vizualizare se încadrează de obicei în această categorie.

• Contra : Fără suport pentru paralelizare explicită : Dacă nu sunteți binecuvântat cu probleme mari, divizibile, dar doriți să utilizați mai multe nuclee, sunteți pe cont propriu. Va trebui să vă dați seama ce clase sunt sigure pentru fire sau chiar reintrante prin încercare și eroare sau citind sursa. Odată am depistat o problemă de paralelizare la un filtru VTK care a folosit o variabilă globală statică pentru a apela o bibliotecă C.

• Pro : Buildsystem CMake : Sistemul de meta-build multi-platformă CMake este dezvoltat și de Kitware (producătorii VTK) și folosit în multe proiecte în afara Kitware. Se integrează foarte bine cu VTK și face mult mai puțin dureroasă configurarea unui sistem de compilare pentru mai multe platforme.

• Pro : Independența platformei, licența și longevitatea : VTK este independent de platformă și este licențiat sub o licență în stil BSD foarte permisivă. În plus, sprijinul profesional este disponibil pentru acele proiecte importante care o necesită. Kitware este susținut de multe entități de cercetare și alte companii și va exista de ceva timp.

Ultimul cuvant

În general, VTK este cel mai bun instrument de vizualizare a datelor pentru tipurile de probleme pe care le iubesc. Dacă întâlniți vreodată un proiect care necesită vizualizare, procesare mesh, procesare imagini sau sarcini similare, încercați să porniți Paraview cu un exemplu de intrare și evaluați dacă VTK ar putea fi instrumentul pentru dvs.