Visualização de dados 3D com ferramentas de código aberto: um tutorial usando VTK

Em seu artigo recente no blog da Toptal, o habilidoso cientista de dados Charles Cook escreveu sobre computação científica com ferramentas de código aberto. Seu tutorial faz uma observação importante sobre ferramentas de código aberto e o papel que elas podem desempenhar no processamento fácil de dados e na obtenção de resultados.

Mas assim que resolvemos todas essas equações diferenciais complexas, surgem outros problemas. Como entendemos e interpretamos as enormes quantidades de dados que saem dessas simulações? Como visualizamos possíveis gigabytes de dados, como dados com milhões de pontos de grade em uma grande simulação?

Durante meu trabalho em problemas semelhantes para minha Dissertação de Mestrado, entrei em contato com o Visualization Toolkit, ou VTK - uma poderosa biblioteca gráfica especializada em visualização de dados.

Neste tutorial, darei uma rápida introdução ao VTK e sua arquitetura de pipeline, e continuarei discutindo um exemplo de visualização 3D da vida real usando dados de um fluido simulado em uma bomba de rotor. Finalmente vou listar os pontos fortes da biblioteca, assim como os pontos fracos que encontrei.

Visualização de dados e o pipeline VTK

A biblioteca de código aberto VTK contém um sólido pipeline de processamento e renderização com muitos algoritmos de visualização sofisticados. Seus recursos, no entanto, não param por aí, pois ao longo do tempo os algoritmos de processamento de imagem e malha também foram adicionados. Em meu projeto atual com uma empresa de pesquisa odontológica, estou utilizando o VTK para tarefas de processamento baseadas em malha em um aplicativo semelhante a CAD baseado em Qt. Os estudos de caso da VTK mostram a ampla gama de aplicações adequadas.

A arquitetura do VTK gira em torno de um poderoso conceito de pipeline. O esboço básico deste conceito é mostrado aqui:

• As fontes estão no início do pipeline e criam “algo do nada”. Por exemplo, um vtkConeSource cria um cone 3D e um vtkSTLReader lê arquivos de geometria 3D *.stl .
• Os filtros transformam a saída de fontes ou outros filtros em algo novo. Por exemplo, um vtkCutter corta a saída do objeto anterior nos algoritmos usando uma função implícita, por exemplo, um plano. Todos os algoritmos de processamento que acompanham o VTK são implementados como filtros e podem ser encadeados livremente.
• Os mapeadores transformam dados em primitivos gráficos. Por exemplo, eles podem ser usados ​​para especificar uma tabela de consulta para colorir dados científicos. Eles são uma maneira abstrata de especificar o que exibir.
• Os atores representam um objeto (geometria mais propriedades de exibição) dentro da cena. Coisas como cor, opacidade, sombreamento ou orientação são especificadas aqui.
• Renderizadores e Windows finalmente descrevem a renderização na tela de maneira independente de plataforma.

Um pipeline de renderização VTK típico começa com uma ou mais fontes, processa-as usando vários filtros em vários objetos de saída, que são renderizados separadamente usando mapeadores e atores. O poder por trás desse conceito é o mecanismo de atualização. Se as configurações de filtros ou fontes forem alteradas, todos os filtros, mapeadores, atores e janelas de renderização dependentes serão atualizados automaticamente. Se, por outro lado, um objeto mais abaixo no pipeline precisar de informações para realizar suas tarefas, ele poderá obtê-las facilmente.

Além disso, não há necessidade de lidar diretamente com sistemas de renderização como o OpenGL. O VTK encapsula toda a tarefa de baixo nível de uma maneira independente de plataforma e (parcialmente) de renderização do sistema; o desenvolvedor trabalha em um nível muito mais alto.

Exemplo de código com um conjunto de dados de bomba de rotor

Vejamos um exemplo de visualização de dados usando um conjunto de dados de fluxo de fluido em uma bomba de rotor rotativo do IEEE Visualization Contest 2011. Os dados em si são o resultado de uma simulação computacional de dinâmica de fluidos, muito parecida com a descrita no artigo de Charles Cook.

Os dados de simulação compactados da bomba apresentada têm mais de 30 GB de tamanho. Ele contém várias partes e várias etapas de tempo, daí o tamanho grande. Neste guia, vamos brincar com a parte do rotor de um desses timesteps, que tem um tamanho compactado de cerca de 150 MB.

Minha linguagem de escolha para usar VTK é C++, mas existem mapeamentos para várias outras linguagens como Tcl/Tk, Java e Python. Se o objetivo é apenas a visualização de um único conjunto de dados, não é necessário escrever código e, em vez disso, pode-se utilizar o Paraview, um front-end gráfico para a maioria das funcionalidades do VTK.

O conjunto de dados e por que 64 bits é necessário

Extraí o conjunto de dados do rotor do conjunto de dados de 30 GB fornecido acima, abrindo um timestep no Paraview e extraindo a parte do rotor em um arquivo separado. É um arquivo de grade não estruturado, ou seja, um volume 3D composto por pontos e células 3D, como hexaedros, tetraedros e assim por diante. Cada um dos pontos 3D tem valores associados. Às vezes, as células também têm valores associados, mas não neste caso. Este treinamento se concentrará na pressão e velocidade nos pontos e tentará visualizá-los em seu contexto 3D.

O tamanho do arquivo compactado é de cerca de 150 MB e o tamanho da memória é de cerca de 280 MB quando carregado com VTK. No entanto, ao processá-lo no VTK, o conjunto de dados é armazenado em cache várias vezes no pipeline do VTK e rapidamente alcançamos o limite de memória de 2 GB para programas de 32 bits. Existem maneiras de economizar memória ao usar o VTK, mas para simplificar vamos apenas compilar e executar o exemplo em 64 bits.

Agradecimentos : O conjunto de dados é disponibilizado por cortesia do Institute of Applied Mechanics, Clausthal University, Alemanha (Dipl. Wirtsch.-Ing. Andreas Lucius).

O alvo

O que vamos conseguir usando o VTK como ferramenta é a visualização mostrada na imagem abaixo. Como um contexto 3D, o contorno do conjunto de dados é mostrado usando uma renderização de estrutura de arame parcialmente transparente. A parte esquerda do conjunto de dados é então usada para exibir a pressão usando um código de cores simples das superfícies. (Vamos pular a renderização de volume mais complexa para este exemplo). Para visualizar o campo de velocidade, a parte direita do conjunto de dados é preenchida com linhas de corrente, que são codificadas por cores pela magnitude de sua velocidade. Essa opção de visualização não é tecnicamente ideal, mas eu queria manter o código VTK o mais simples possível. Além disso, há uma razão para este exemplo fazer parte de um desafio de visualização, ou seja, muita turbulência no fluxo.

Passo a passo

Discutirei o código VTK passo a passo, mostrando como a saída de renderização ficaria em cada estágio. O código-fonte completo pode ser baixado no final do treinamento.

Vamos começar incluindo tudo o que precisamos do VTK e abrir a função main.

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

Em seguida, configuramos o renderizador e a janela de renderização para exibir nossos resultados. Definimos a cor de fundo e o tamanho da janela de renderização.

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

Com este código já poderíamos exibir uma janela de renderização estática. Em vez disso, optamos por adicionar um vtkRenderWindowInteractor para girar, ampliar e deslocar interativamente a cena.

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

Agora temos um exemplo em execução mostrando uma janela de renderização cinza e vazia.

Em seguida, carregamos o conjunto de dados usando um dos muitos leitores que acompanham o VTK.

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

Curta excursão ao gerenciamento de memória VTK : VTK usa um conceito conveniente de gerenciamento automático de memória que gira em torno da contagem de referências. No entanto, diferente da maioria das outras implementações, a contagem de referência é mantida dentro dos próprios objetos VTK, em vez da classe de ponteiro inteligente. Isso tem a vantagem de que a contagem de referência pode ser aumentada, mesmo se o objeto VTK for passado como um ponteiro bruto. Há duas maneiras principais de criar objetos VTK gerenciados. vtkNew<T> e vtkSmartPointer<T>::New() , com a principal diferença sendo que um vtkSmartPointer<T> pode ser convertido implícito para o ponteiro bruto T* e pode ser retornado de uma função. Para instâncias de vtkNew<T> teremos que chamar .Get() para obter um ponteiro bruto e só podemos devolvê-lo envolvendo-o em um vtkSmartPointer . Em nosso exemplo, nunca retornamos de funções e todos os objetos permanecem ativos o tempo todo, portanto, usaremos o short vtkNew , com apenas a exceção acima para fins de demonstração.

Neste ponto, nada foi lido do arquivo ainda. Nós ou um filtro mais abaixo na cadeia teríamos que chamar Update() para que a leitura do arquivo realmente acontecesse. Geralmente é a melhor abordagem deixar as classes VTK lidarem com as atualizações por conta própria. No entanto, às vezes queremos acessar o resultado de um filtro diretamente, por exemplo, para obter a faixa de pressões neste conjunto de dados. Então precisamos chamar Update() manualmente. (Não perdemos desempenho chamando Update() várias vezes, pois os resultados são armazenados em cache.)

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

Em seguida, precisamos extrair a metade esquerda do conjunto de dados, usando vtkClipDataSet . Para conseguir isso, primeiro definimos um vtkPlane que define a divisão. Então, veremos pela primeira vez como o pipeline VTK está conectado: successor->SetInputConnection(predecessor->GetOutputPort()) . Sempre que solicitarmos uma atualização do clipperLeft , essa conexão agora garantirá que todos os filtros anteriores também estejam atualizados.

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

Por fim, criamos nossos primeiros atores e mapeadores para exibir a renderização de wireframe da metade esquerda. Observe que o mapeador está conectado ao seu filtro exatamente da mesma maneira que os filtros entre si. Na maioria das vezes, o próprio renderizador está acionando as atualizações de todos os atores, mapeadores e as cadeias de filtros subjacentes!

A única linha que não é autoexplicativa é provavelmente leftWireMapper->ScalarVisibilityOff(); - proíbe a coloração do wireframe por valores de pressão, que são definidos como o array atualmente ativo.

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

Neste ponto, a janela de renderização finalmente mostra algo, ou seja, o wireframe da parte esquerda.

A renderização de wireframe para a parte direita é criada de maneira semelhante, alternando o plano normal de um vtkClipDataSet (recém-criado) para a direção oposta e alterando ligeiramente a cor e a opacidade do mapeador e ator (recém-criado). Observe que aqui nosso pipeline VTK se divide em duas direções (direita e esquerda) do mesmo conjunto de dados 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());

A janela de saída agora mostra ambas as partes do wireframe, como esperado.

Agora estamos prontos para visualizar alguns dados úteis! Para adicionar a visualização de pressão à parte esquerda, não precisamos fazer muito. Criamos um novo mapeador e o conectamos ao clipperLeft também, mas desta vez colorimos pela matriz de pressão. É também aqui que finalmente utilizamos o pressureRange que derivamos acima.

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

A saída agora se parece com a imagem mostrada abaixo. A pressão no meio é muito baixa, sugando material para dentro da bomba. Em seguida, esse material é transportado para o exterior, ganhando pressão rapidamente. (É claro que deve haver uma legenda do mapa de cores com os valores reais, mas deixei de fora para manter o exemplo mais curto.)

Agora começa a parte mais complicada. Queremos desenhar linhas de corrente de velocidade na parte direita. As linhas de fluxo são geradas pela integração dentro de um campo vetorial a partir de pontos de origem. O campo vetorial já faz parte do conjunto de dados na forma do vetor-array “Velocities”. Portanto, precisamos apenas gerar os pontos de origem. vtkPointSource gera uma esfera de pontos aleatórios. Geraremos 1500 pontos de origem, porque a maioria deles não estará dentro do conjunto de dados e será ignorada pelo rastreador de fluxo.

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

Em seguida, criamos o streamtracer e definimos suas conexões de entrada. “Espere, várias conexões?”, você pode dizer. Sim - este é o primeiro filtro VTK com várias entradas que encontramos. A conexão de entrada normal é usada para o campo vetorial e a conexão de origem é usada para os pontos de semente. Já que “Velocities” é o array de vetores “ativo” em clipperRight , não precisamos especificá-lo aqui explicitamente. Finalmente, especificamos que a integração deve ser realizada em ambas as direções a partir dos pontos de semente e configuramos o método de integração para Runge-Kutta-4.5.

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

Nosso próximo problema é colorir as linhas de corrente pela magnitude da velocidade. Como não há matriz para as magnitudes dos vetores, simplesmente calcularemos as magnitudes em uma nova matriz escalar. Como você adivinhou, também existe um filtro VTK para esta tarefa: vtkArrayCalculator . Ele pega um conjunto de dados e o produz inalterado, mas adiciona exatamente uma matriz que é calculada a partir de uma ou mais das existentes. Configuramos esta calculadora de matriz para obter a magnitude do vetor “Velocity” e produzi-lo como “MagVelocity”. Finalmente, chamamos Update() manualmente novamente, para derivar o intervalo do novo 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 gera polilinhas diretamente e vtkArrayCalculator as transmite inalteradas. Portanto, poderíamos apenas exibir a saída do magCalc diretamente usando um novo mapeador e ator.

Em vez disso, neste treinamento, optamos por tornar a saída um pouco melhor, exibindo fitas. vtkRibbonFilter gera células 2D para exibir fitas para todas as polilinhas de sua 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());

O que ainda está faltando, e é realmente necessário para produzir as renderizações intermediárias também, são as últimas cinco linhas para realmente renderizar a cena e inicializar o interator.

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

Por fim, chegamos à visualização finalizada, que apresentarei mais uma vez aqui:

O código-fonte completo para a visualização acima pode ser encontrado aqui.

O bom, o Mau e o Feio

Vou fechar este artigo com uma lista dos meus prós e contras pessoais do framework VTK.

• Pro : Desenvolvimento ativo : VTK está em desenvolvimento ativo de vários colaboradores, principalmente da comunidade de pesquisa. Isso significa que alguns algoritmos de ponta estão disponíveis, muitos formatos 3D podem ser importados e exportados, os bugs são corrigidos ativamente e os problemas geralmente têm uma solução pronta nos fóruns de discussão.

• Contras : Confiabilidade : Acoplar muitos algoritmos de diferentes contribuidores com o design de pipeline aberto do VTK, no entanto, pode levar a problemas com combinações de filtros incomuns. Eu tive que entrar no código-fonte do VTK algumas vezes para descobrir por que minha complexa cadeia de filtros não está produzindo os resultados desejados. Eu recomendaria fortemente configurar o VTK de uma maneira que permita a depuração.

• Prós : Arquitetura de software : O design do pipeline e a arquitetura geral do VTK parecem bem pensados ​​e é um prazer trabalhar com eles. Algumas linhas de código podem produzir resultados surpreendentes. As estruturas de dados integradas são fáceis de entender e usar.

• Contras : Microarquitetura : Algumas decisões de design de microarquitetura escapam ao meu entendimento. A correção de const é quase inexistente, as matrizes são passadas como entradas e saídas sem distinção clara. Eu aliviei isso para meus próprios algoritmos, desistindo de algum desempenho e usando meu próprio wrapper para vtkMath que utiliza tipos 3D personalizados como typedef std::array<double, 3> Pnt3d; .

• Pro : Micro Documentação : A documentação do Doxygen de todas as classes e filtros é extensa e utilizável, os exemplos e casos de teste no wiki também são uma grande ajuda para entender como os filtros são usados.

• Contras : Documentação Macro : Existem vários bons tutoriais e introduções ao VTK na web. No entanto, até onde eu sei, não há uma grande documentação de referência que explique como as coisas específicas são feitas. Se você quer fazer algo novo, espere pesquisar como fazê-lo por algum tempo. Além disso, é difícil encontrar o filtro específico para uma tarefa. No entanto, uma vez que você o encontrou, a documentação do Doxygen geralmente será suficiente. Uma boa maneira de explorar o framework VTK é baixar e experimentar o Paraview.

• Pró : Suporte à Paralelização Implícita : Se suas fontes podem ser divididas em várias partes que podem ser processadas independentemente, a paralelização é tão simples quanto criar uma cadeia de filtros separada dentro de cada thread que processa uma única parte. A maioria dos grandes problemas de visualização geralmente se enquadram nesta categoria.

• Contras : Não Suporte explícito de paralelização : Se você não é abençoado com problemas grandes e divisíveis, mas deseja utilizar vários núcleos, está por conta própria. Você terá que descobrir quais classes são thread-safe, ou mesmo reentrantes por tentativa e erro ou lendo a fonte. Certa vez, rastreei um problema de paralelização em um filtro VTK que usava uma variável global estática para chamar alguma biblioteca C.

• Pro : Buildsystem CMake : O meta-build-system multiplataforma CMake também é desenvolvido pela Kitware (os criadores do VTK) e usado em muitos projetos fora do Kitware. Ele se integra muito bem ao VTK e torna a configuração de um sistema de compilação para várias plataformas muito menos dolorosa.

• Pro : Independência de plataforma, licença e longevidade : VTK é independente de plataforma e está licenciado sob uma licença de estilo BSD muito permissiva. Além disso, o suporte profissional está disponível para os projetos importantes que o exigem. Kitware é apoiado por muitas entidades de pesquisa e outras empresas e estará disponível por algum tempo.

Última palavra

No geral, o VTK é a melhor ferramenta de visualização de dados para os tipos de problemas que eu amo. Se você encontrar um projeto que exija visualização, processamento de malha, processamento de imagem ou tarefas semelhantes, tente iniciar o Paraview com um exemplo de entrada e avalie se o VTK pode ser a ferramenta para você.