Visualización de datos 3D con herramientas de código abierto: un tutorial usando VTK

En su artículo reciente en el blog de Toptal, el experto científico de datos Charles Cook escribió sobre computación científica con herramientas de código abierto. Su tutorial destaca un punto importante sobre las herramientas de código abierto y el papel que pueden desempeñar en el procesamiento de datos y la adquisición de resultados con facilidad.

Pero tan pronto como hemos resuelto todas estas ecuaciones diferenciales complejas, surge otro problema. ¿Cómo entendemos e interpretamos las enormes cantidades de datos que surgen de estas simulaciones? ¿Cómo visualizamos gigabytes potenciales de datos, como datos con millones de puntos de cuadrícula dentro de una gran simulación?

Durante mi trabajo en problemas similares para mi tesis de maestría, entré en contacto con Visualization Toolkit o VTK, una poderosa biblioteca de gráficos especializada para la visualización de datos.

En este tutorial, daré una introducción rápida a VTK y su arquitectura de tubería, y continuaré discutiendo un ejemplo de visualización 3D de la vida real utilizando datos de un fluido simulado en una bomba impulsora. Finalmente, enumeraré los puntos fuertes de la biblioteca, así como los puntos débiles que encontré.

Visualización de datos y el canal VTK

La biblioteca de código abierto VTK contiene una sólida canalización de procesamiento y representación con muchos algoritmos de visualización sofisticados. Sin embargo, sus capacidades no terminan ahí, ya que con el tiempo también se han agregado algoritmos de procesamiento de imágenes y mallas. En mi proyecto actual con una empresa de investigación dental, estoy utilizando VTK para tareas de procesamiento basadas en malla dentro de una aplicación similar a CAD basada en Qt. Los estudios de casos de VTK muestran la amplia gama de aplicaciones adecuadas.

La arquitectura de VTK gira en torno a un poderoso concepto de tubería. El esquema básico de este concepto se muestra aquí:

• Las fuentes están al principio de la canalización y crean "algo de la nada". Por ejemplo, un vtkConeSource crea un cono 3D y un vtkSTLReader lee archivos de geometría 3D *.stl .
• Los filtros transforman la salida de fuentes u otros filtros en algo nuevo. Por ejemplo, un vtkCutter corta la salida del objeto anterior en los algoritmos utilizando una función implícita, por ejemplo, un plano. Todos los algoritmos de procesamiento que vienen con VTK se implementan como filtros y se pueden encadenar libremente.
• Los mapeadores transforman los datos en primitivas gráficas. Por ejemplo, se pueden usar para especificar una tabla de consulta para colorear datos científicos. Son una forma abstracta de especificar qué mostrar.
• Los actores representan un objeto (geometría más propiedades de visualización) dentro de la escena. Aquí se especifican cosas como el color, la opacidad, el sombreado o la orientación.
• Renderers & Windows finalmente describen el renderizado en la pantalla de una manera independiente de la plataforma.

Una canalización típica de renderizado de VTK comienza con una o más fuentes, las procesa usando varios filtros en varios objetos de salida, que luego se renderizan por separado usando mapeadores y actores. El poder detrás de este concepto es el mecanismo de actualización. Si se cambia la configuración de los filtros o las fuentes, todos los filtros, mapeadores, actores y ventanas de renderizado dependientes se actualizan automáticamente. Si, por el contrario, un objeto más abajo en la canalización necesita información para realizar sus tareas, puede obtenerla fácilmente.

Además, no hay necesidad de tratar directamente con sistemas de renderizado como OpenGL. VTK encapsula todas las tareas de bajo nivel de una manera independiente de la plataforma y (parcialmente) del sistema de renderizado; el desarrollador trabaja en un nivel mucho más alto.

Ejemplo de código con un conjunto de datos de bomba de rotor

Veamos un ejemplo de visualización de datos utilizando un conjunto de datos de flujo de fluido en una bomba de impulsor giratorio del Concurso de visualización IEEE 2011. Los datos en sí son el resultado de una simulación de dinámica de fluidos computacional, muy parecida a la descrita en el artículo de Charles Cook.

Los datos de simulación comprimidos de la bomba presentada tienen un tamaño de más de 30 GB. Contiene múltiples partes y múltiples pasos de tiempo, de ahí su gran tamaño. En esta guía, jugaremos con la parte del rotor de uno de estos intervalos de tiempo, que tiene un tamaño comprimido de unos 150 MB.

Mi lenguaje de elección para usar VTK es C++, pero hay asignaciones para varios otros lenguajes como Tcl/Tk, Java y Python. Si el objetivo es solo la visualización de un solo conjunto de datos, no es necesario escribir ningún código y, en su lugar, puede utilizar Paraview, una interfaz gráfica para la mayoría de las funciones de VTK.

El conjunto de datos y por qué es necesario 64 bits

Extraje el conjunto de datos del rotor del conjunto de datos de 30 GB proporcionado anteriormente, abriendo un paso de tiempo en Paraview y extrayendo la parte del rotor en un archivo separado. Es un archivo de cuadrícula no estructurado, es decir, un volumen 3D que consta de puntos y celdas 3D, como hexaedros, tetraedros, etc. Cada uno de los puntos 3D tiene valores asociados. A veces, las celdas también tienen valores asociados, pero no en este caso. Este entrenamiento se concentrará en la presión y la velocidad en los puntos y tratará de visualizarlos en su contexto 3D.

El tamaño del archivo comprimido es de aproximadamente 150 MB y el tamaño en memoria es de aproximadamente 280 MB cuando se carga con VTK. Sin embargo, al procesarlo en VTK, el conjunto de datos se almacena en caché varias veces dentro de la canalización de VTK y alcanzamos rápidamente el límite de memoria de 2 GB para programas de 32 bits. Hay formas de ahorrar memoria cuando se usa VTK, pero para simplificar, solo compilaremos y ejecutaremos el ejemplo en 64 bits.

Agradecimientos : El conjunto de datos está disponible por cortesía del Instituto de Mecánica Aplicada, Universidad de Clausthal, Alemania (Dipl. Wirtsch.-Ing. Andreas Lucius).

El objetivo

Lo que lograremos usando VTK como herramienta es la visualización que se muestra en la imagen a continuación. Como contexto 3D, el contorno del conjunto de datos se muestra mediante una representación de estructura alámbrica parcialmente transparente. Luego, la parte izquierda del conjunto de datos se usa para mostrar la presión usando una codificación de color simple de las superficies. (Omitiremos la representación de volumen más compleja para este ejemplo). Para visualizar el campo de velocidad, la parte derecha del conjunto de datos se llena con líneas de corriente, que están codificadas por colores según la magnitud de su velocidad. Esta opción de visualización técnicamente no es ideal, pero quería mantener el código VTK lo más simple posible. Además, hay una razón para que este ejemplo sea parte de un desafío de visualización, es decir, muchas turbulencias en el flujo.

Paso a paso

Discutiré el código VTK paso a paso, mostrando cómo se vería la salida de renderizado en cada etapa. El código fuente completo se puede descargar al final de la capacitación.

Comencemos por incluir todo lo que necesitamos de VTK y abramos la función 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) {

A continuación, configuramos el renderizador y la ventana de renderizado para mostrar nuestros resultados. Establecemos el color de fondo y el tamaño de la ventana de renderizado.

 // 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 este código ya podríamos mostrar una ventana de renderizado estática. En su lugar, optamos por agregar un vtkRenderWindowInteractor para rotar, ampliar y desplazar la escena de forma interactiva.

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

Ahora tenemos un ejemplo en ejecución que muestra una ventana de representación gris y vacía.

Luego, cargamos el conjunto de datos usando uno de los muchos lectores que vienen con VTK.

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

Breve excursión a la gestión de memoria VTK : VTK utiliza un conveniente concepto de gestión de memoria automática que gira en torno al recuento de referencias. Sin embargo, a diferencia de la mayoría de las otras implementaciones, el recuento de referencias se mantiene dentro de los propios objetos VTK, en lugar de la clase de puntero inteligente. Esto tiene la ventaja de que se puede aumentar el recuento de referencias, incluso si el objeto VTK se pasa como un puntero sin procesar. Hay dos formas principales de crear objetos VTK administrados. vtkNew<T> y vtkSmartPointer<T>::New() , con la principal diferencia de que un vtkSmartPointer<T> se puede convertir de forma implícita al puntero sin procesar T* y se puede devolver desde una función. Para instancias de vtkNew<T> tendremos que llamar a .Get() para obtener un puntero sin formato y solo podemos devolverlo envolviéndolo en un vtkSmartPointer . Dentro de nuestro ejemplo, nunca regresamos de las funciones y todos los objetos viven todo el tiempo, por lo tanto, usaremos el vtkNew corto, con solo la excepción anterior para fines de demostración.

En este punto, todavía no se ha leído nada del archivo. Nosotros o un filtro más abajo en la cadena tendríamos que llamar a Update() para que la lectura del archivo realmente suceda. Por lo general, el mejor enfoque es dejar que las clases VTK manejen las actualizaciones por sí mismas. Sin embargo, a veces queremos acceder directamente al resultado de un filtro, por ejemplo, para obtener el rango de presiones en este conjunto de datos. Luego, debemos llamar a Update() manualmente. (No perdemos rendimiento llamando a Update() varias veces, ya que los resultados se almacenan en caché).

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

A continuación, necesitamos extraer la mitad izquierda del conjunto de datos, usando vtkClipDataSet . Para lograr esto, primero definimos un vtkPlane que define la división. Luego, veremos por primera vez cómo se conecta la canalización VTK: successor->SetInputConnection(predecessor->GetOutputPort()) . Siempre que solicitemos una actualización de clipperLeft , esta conexión ahora garantizará que todos los filtros anteriores también estén actualizados.

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

Finalmente, creamos nuestros primeros actores y mapeadores para mostrar la representación alámbrica de la mitad izquierda. Tenga en cuenta que el mapeador está conectado a su filtro exactamente de la misma manera que los filtros entre sí. La mayoría de las veces, el propio renderizador activa las actualizaciones de todos los actores, mapeadores y las cadenas de filtros subyacentes.

La única línea que no se explica por sí misma es probablemente leftWireMapper->ScalarVisibilityOff(); - prohíbe la coloración de la estructura alámbrica por valores de presión, que se establecen como la matriz actualmente activa.

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

En este punto, la ventana de renderizado finalmente muestra algo, es decir, la estructura alámbrica de la parte izquierda.

La representación alámbrica para la parte derecha se crea de manera similar, cambiando el plano normal de un vtkClipDataSet (recién creado) a la dirección opuesta y cambiando ligeramente el color y la opacidad del asignador y actor (recién creado). Tenga en cuenta que aquí nuestra tubería VTK se divide en dos direcciones (derecha e izquierda) desde el mismo conjunto de datos de entrada.

 // 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 ventana de salida ahora muestra ambas partes de estructura alámbrica, como se esperaba.

¡Ahora estamos listos para visualizar algunos datos útiles! Para agregar la visualización de presión a la parte izquierda, no necesitamos hacer mucho. Creamos un nuevo mapeador y lo conectamos a clipperLeft también, pero esta vez coloreamos según la matriz de presión. También es aquí donde finalmente utilizamos el rango de pressureRange que hemos obtenido anteriormente.

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

La salida ahora se parece a la imagen que se muestra a continuación. La presión en el centro es muy baja y succiona material hacia la bomba. Luego, este material es transportado al exterior, ganando rápidamente presión. (Por supuesto, debería haber una leyenda de mapa de colores con los valores reales, pero la omití para acortar el ejemplo).

Ahora empieza la parte más complicada. Queremos dibujar líneas de corriente de velocidad en la parte derecha. Las líneas de corriente se generan mediante la integración dentro de un campo vectorial desde puntos de origen. El campo vectorial ya forma parte del conjunto de datos en forma de matriz vectorial "Velocidades". Entonces solo necesitamos generar los puntos de origen. vtkPointSource genera una esfera de puntos aleatorios. Generaremos 1500 puntos de origen, porque la mayoría de ellos no se encontrarán dentro del conjunto de datos de todos modos y serán ignorados por el rastreador de flujo.

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

A continuación, creamos el trazador de flujo y establecemos sus conexiones de entrada. “Espera, ¿conexiones múltiples ?”, podrías decir. Sí, este es el primer filtro VTK con múltiples entradas que encontramos. La conexión de entrada normal se usa para el campo vectorial y la conexión de origen se usa para los puntos semilla. Dado que "Velocidades" es la matriz vectorial "activa" en clipperRight , no necesitamos especificarlo aquí explícitamente. Finalmente, especificamos que la integración debe realizarse en ambas direcciones desde los puntos semilla y establecemos el método de integración en Runge-Kutta-4.5.

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

Nuestro siguiente problema es colorear las líneas de corriente según la magnitud de la velocidad. Dado que no existe una matriz para las magnitudes de los vectores, simplemente calcularemos las magnitudes en una nueva matriz escalar. Como habrá adivinado, también hay un filtro VTK para esta tarea: vtkArrayCalculator . Toma un conjunto de datos y lo genera sin cambios, pero agrega exactamente una matriz que se calcula a partir de una o más de las existentes. Configuramos esta calculadora de matriz para tomar la magnitud del vector "Velocity" y generarlo como "MagVelocity". Finalmente, volvemos a llamar a Update() manualmente para derivar el rango de la nueva matriz.

 // 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 genera directamente polilíneas y vtkArrayCalculator las pasa sin cambios. Por lo tanto, podríamos simplemente mostrar la salida de magCalc directamente usando un nuevo mapeador y actor.

En cambio, en esta capacitación optamos por hacer que la salida sea un poco más agradable, mostrando cintas en su lugar. vtkRibbonFilter genera celdas 2D para mostrar cintas para todas las polilíneas de su entrada.

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

Lo que ahora todavía falta, y en realidad se necesita para producir las renderizaciones intermedias también, son las últimas cinco líneas para renderizar la escena e inicializar el interactor.

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

Finalmente, llegamos a la visualización terminada, que presentaré una vez más aquí:

El código fuente completo para la visualización anterior se puede encontrar aquí.

Lo bueno, lo malo y lo feo

Cerraré este artículo con una lista de mis ventajas y desventajas personales del marco VTK.

• Pro : Desarrollo activo : VTK está en desarrollo activo por parte de varios colaboradores, principalmente dentro de la comunidad de investigación. Esto significa que hay algunos algoritmos de vanguardia disponibles, se pueden importar y exportar muchos formatos 3D, los errores se corrigen activamente y los problemas generalmente tienen una solución preparada en los foros de discusión.

• Con : Confiabilidad : Sin embargo, combinar muchos algoritmos de diferentes colaboradores con el diseño de canalización abierta de VTK puede generar problemas con combinaciones de filtros inusuales. Tuve que ingresar al código fuente de VTK varias veces para descubrir por qué mi compleja cadena de filtros no produce los resultados deseados. Recomiendo encarecidamente configurar VTK de una manera que permita la depuración.

• Pro : Arquitectura de software : El diseño de tubería y la arquitectura general de VTK parecen bien pensados ​​y es un placer trabajar con ellos. Unas pocas líneas de código pueden producir resultados sorprendentes. Las estructuras de datos integradas son fáciles de entender y usar.

• Con : Micro Arquitectura : Algunas decisiones de diseño de micro-arquitectura escapan a mi comprensión. La corrección constante es casi inexistente, las matrices se transmiten como entradas y salidas sin una distinción clara. Alivié esto para mis propios algoritmos renunciando a algo de rendimiento y usando mi propio contenedor para vtkMath que utiliza tipos 3D personalizados como typedef std::array<double, 3> Pnt3d; .

• Pro : Micro Documentación : La documentación de Doxygen de todas las clases y filtros es extensa y utilizable, los ejemplos y casos de prueba en la wiki también son de gran ayuda para entender cómo se usan los filtros.

• Con : Documentación de macros : Hay varios buenos tutoriales e introducciones a VTK en la web. Sin embargo, que yo sepa, no hay una gran documentación de referencia que explique cómo se hacen las cosas específicas. Si quieres hacer algo nuevo, espera buscar cómo hacerlo durante algún tiempo. Además, es difícil encontrar el filtro específico para una tarea. Sin embargo, una vez que lo haya encontrado, la documentación de Doxygen normalmente será suficiente. Una buena manera de explorar el marco VTK es descargar y experimentar con Paraview.

• Pro : soporte de paralelización implícita : si sus fuentes se pueden dividir en varias partes que se pueden procesar de forma independiente, la paralelización es tan simple como crear una cadena de filtro separada dentro de cada subproceso que procesa una sola parte. La mayoría de los grandes problemas de visualización suelen caer en esta categoría.

• Desventaja : No hay soporte de paralelización explícito : si no está bendecido con problemas grandes y divisibles, pero quiere utilizar múltiples núcleos, está solo. Tendrá que averiguar qué clases son seguras para subprocesos, o incluso reingresar por prueba y error o leyendo la fuente. Una vez localicé un problema de paralelización en un filtro VTK que usaba una variable global estática para llamar a alguna biblioteca C.

• Pro : Buildsystem CMake : Kitware (los creadores de VTK) también desarrolla el meta-sistema de compilación multiplataforma CMake y se usa en muchos proyectos fuera de Kitware. Se integra muy bien con VTK y hace que la configuración de un sistema de compilación para múltiples plataformas sea mucho menos dolorosa.

• Pro : independencia de la plataforma, licencia y longevidad : VTK es una plataforma independiente lista para usar, y tiene una licencia de estilo BSD muy permisiva. Además, se dispone de apoyo profesional para aquellos proyectos importantes que lo requieran. Kitware está respaldado por muchas entidades de investigación y otras empresas y estará disponible por algún tiempo.

Ultima palabra

En general, VTK es la mejor herramienta de visualización de datos para los tipos de problemas que amo. Si alguna vez se encuentra con un proyecto que requiere visualización, procesamiento de mallas, procesamiento de imágenes o tareas similares, intente iniciar Paraview con un ejemplo de entrada y evalúe si VTK podría ser la herramienta para usted.