Визуализация 3D-данных с помощью инструментов с открытым исходным кодом: учебное пособие с использованием VTK

В своей недавней статье в блоге Toptal опытный специалист по данным Чарльз Кук написал о научных вычислениях с помощью инструментов с открытым исходным кодом. В его учебном пособии делается важный вывод об инструментах с открытым исходным кодом и о той роли, которую они могут играть в простой обработке данных и получении результатов.

Но как только мы решили все эти сложные дифференциальные уравнения, возникают другие проблемы. Как мы понимаем и интерпретируем огромные объемы данных, получаемых в результате этих симуляций? Как визуализировать потенциальные гигабайты данных, например, данные с миллионами точек сетки в рамках большой симуляции?

Во время работы над аналогичными задачами для магистерской диссертации я столкнулся с Visualization Toolkit, или VTK — мощной графической библиотекой, специализирующейся на визуализации данных.

В этом руководстве я кратко расскажу о VTK и его архитектуре трубопровода, а затем перейду к обсуждению реального примера 3D-визуализации с использованием данных смоделированной жидкости в крыльчатом насосе. Наконец, я перечислю сильные стороны библиотеки, а также слабые места, с которыми я столкнулся.

Визуализация данных и конвейер VTK

Библиотека с открытым исходным кодом VTK содержит надежный конвейер обработки и рендеринга со множеством сложных алгоритмов визуализации. Однако его возможности на этом не заканчиваются, так как со временем были добавлены алгоритмы обработки изображений и сетки. В моем текущем проекте с стоматологической исследовательской компанией я использую VTK для задач обработки на основе сетки в приложении на основе Qt, похожем на САПР. Тематические исследования VTK показывают широкий спектр подходящих приложений.

Архитектура VTK основана на концепции мощного конвейера. Основная схема этой концепции показана здесь:

• Исходники находятся в самом начале пайплайна и создают «что-то из ничего». Например, vtkConeSource создает трехмерный конус, а vtkSTLReader читает файлы трехмерной геометрии *.stl .
• Фильтры преобразуют вывод источников или других фильтров во что-то новое. Например, vtkCutter вырезает вывод предыдущего объекта в алгоритмах, используя неявную функцию, например плоскость. Все алгоритмы обработки, поставляемые с VTK, реализованы в виде фильтров и могут свободно объединяться в цепочку.
• Картографы преобразуют данные в графические примитивы. Например, их можно использовать для указания справочной таблицы для раскрашивания научных данных. Это абстрактный способ указать, что отображать.
• Актеры представляют объект (геометрию плюс свойства отображения) в сцене. Здесь задаются такие вещи, как цвет, непрозрачность, затенение или ориентация.
• Рендереры и Windows , наконец, описывают рендеринг на экране независимым от платформы способом.

Типичный конвейер рендеринга VTK начинается с одного или нескольких источников, обрабатывает их с помощью различных фильтров в несколько выходных объектов, которые затем рендерятся отдельно с помощью картографов и акторов. Сила этой концепции заключается в механизме обновления. При изменении настроек фильтров или источников автоматически обновляются все зависимые фильтры, мапперы, акторы и окна рендеринга. Если, с другой стороны, объекту, находящемуся дальше по конвейеру, нужна информация для выполнения своих задач, он может легко ее получить.

Кроме того, нет необходимости напрямую иметь дело с системами рендеринга, такими как OpenGL. VTK инкапсулирует все низкоуровневые задачи платформо- и (частично) независимым от системы рендеринга способом; разработчик работает на гораздо более высоком уровне.

Пример кода с набором данных роторного насоса

Давайте рассмотрим пример визуализации данных с использованием набора данных потока жидкости в насосе с вращающейся крыльчаткой из конкурса визуализации IEEE 2011. Сами данные являются результатом вычислительного гидродинамического моделирования, очень похожего на то, что описано в статье Чарльза Кука.

Заархивированные данные моделирования показанного насоса имеют размер более 30 ГБ. Он содержит несколько частей и несколько временных шагов, отсюда и большой размер. В этом руководстве мы поиграем с роторной частью одного из этих временных шагов, размер которого в сжатом виде составляет около 150 МБ.

Мой предпочтительный язык для использования VTK — C++, но есть сопоставления для нескольких других языков, таких как Tcl/Tk, Java и Python. Если целью является просто визуализация одного набора данных, вам вообще не нужно писать код, и вместо этого можно использовать Paraview, графический интерфейс для большей части функций VTK.

Набор данных и зачем нужны 64-битные версии

Я извлек набор данных ротора из предоставленного выше набора данных объемом 30 ГБ, открыв один временной шаг в Paraview и извлек часть ротора в отдельный файл. Это файл неструктурированной сетки, т. е. трехмерный объем, состоящий из точек и трехмерных ячеек, таких как шестигранники, тетраэдры и т. д. Каждая из трехмерных точек имеет связанные значения. Иногда ячейки также имеют связанные значения, но не в этом случае. Это обучение будет сосредоточено на давлении и скорости в точках и попытается визуализировать их в трехмерном контексте.

Размер сжатого файла составляет около 150 МБ, а размер в памяти — около 280 МБ при загрузке с помощью VTK. Однако, обрабатывая его в VTK, набор данных многократно кэшируется в конвейере VTK, и мы быстро достигаем предела памяти в 2 ГБ для 32-битных программ. Есть способы сэкономить память при использовании VTK, но для простоты мы просто скомпилируем и запустим пример в 64-битной версии.

Благодарности : набор данных предоставлен Институтом прикладной механики Клаустальского университета, Германия (дипл. Wirtsch.-Ing. Andreas Lucius).

Цель

Чего мы добьемся, используя VTK в качестве инструмента, так это визуализацию, показанную на изображении ниже. В качестве 3D-контекста схема набора данных отображается с использованием частично прозрачного каркасного рендеринга. Затем левая часть набора данных используется для отображения давления с помощью простого цветового кодирования поверхностей. (Для этого примера мы пропустим более сложный объемный рендеринг). Чтобы визуализировать поле скорости, правая часть набора данных заполнена линиями тока, которые имеют цветовую кодировку по величине их скорости. Такой вариант визуализации технически не идеален, но я хотел, чтобы код VTK был максимально простым. Кроме того, есть причина, по которой этот пример является частью задачи визуализации, т. е. много турбулентности в потоке.

Шаг за шагом

Я буду обсуждать код VTK шаг за шагом, показывая, как результат рендеринга будет выглядеть на каждом этапе. Полный исходный код можно скачать в конце обучения.

Начнем с того, что включим все необходимое из VTK и откроем основную функцию.

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

Затем мы настраиваем рендерер и окно рендеринга, чтобы отобразить наши результаты. Мы устанавливаем цвет фона и размер окна рендера.

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

С помощью этого кода мы уже могли отображать статическое окно рендеринга. Вместо этого мы решили добавить vtkRenderWindowInteractor для интерактивного вращения, масштабирования и панорамирования сцены.

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

Теперь у нас есть работающий пример, показывающий серое пустое окно рендеринга.

Затем мы загружаем набор данных с помощью одного из множества ридеров, поставляемых с VTK.

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

Краткий экскурс в управление памятью VTK : VTK использует удобную концепцию автоматического управления памятью, основанную на подсчете ссылок. Однако, в отличие от большинства других реализаций, счетчик ссылок хранится в самих объектах VTK, а не в классе интеллектуальных указателей. Это имеет то преимущество, что счетчик ссылок может быть увеличен, даже если объект VTK передается как необработанный указатель. Существует два основных способа создания управляемых объектов VTK. vtkNew<T> и vtkSmartPointer<T>::New() , с основным отличием в том, что vtkSmartPointer<T> неявно может быть приведен к необработанному указателю T* и может быть возвращен из функции. Для экземпляров vtkNew<T> нам нужно вызвать .Get() , чтобы получить необработанный указатель, и мы можем вернуть его, только завернув его в vtkSmartPointer . В нашем примере мы никогда не возвращаемся из функций, и все объекты живут все время, поэтому мы будем использовать короткий vtkNew , за исключением только приведенного выше исключения в демонстрационных целях.

На данный момент из файла еще ничего не было прочитано. Нам или фильтру, находящемуся дальше по цепочке, пришлось бы вызывать Update() , чтобы чтение файла действительно произошло. Обычно лучше всего позволить классам VTK самим обрабатывать обновления. Однако иногда мы хотим получить прямой доступ к результату фильтра, например, чтобы получить диапазон давлений в этом наборе данных. Затем нам нужно вызвать Update() вручную. (Мы не теряем производительность, вызывая Update() несколько раз, так как результаты кэшируются.)

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

Далее нам нужно извлечь левую половину набора данных, используя vtkClipDataSet . Для этого мы сначала определяем vtkPlane , который определяет разделение. Затем мы впервые увидим, как конвейер VTK связан вместе: successor->SetInputConnection(predecessor->GetOutputPort()) . Всякий раз, когда мы запрашиваем обновление от clipperLeft это соединение теперь гарантирует, что все предыдущие фильтры также обновлены.

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

Наконец, мы создаем наших первых акторов и мапперов для отображения каркасного рендеринга левой половины. Обратите внимание, что маппер связан со своим фильтром точно так же, как фильтры друг с другом. Большую часть времени рендерер сам запускает обновления всех акторов, картографов и базовых цепочек фильтров!

Единственная строка, которая не говорит сама за себя, вероятно, leftWireMapper->ScalarVisibilityOff(); - запрещает окрашивание каркаса по значениям давления, которые заданы как текущий активный массив.

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

В этот момент окно рендеринга наконец-то что-то показывает, т. е. каркас для левой части.

Рендеринг каркаса для правой части создается аналогичным образом, путем переключения нормали плоскости (вновь созданного) vtkClipDataSet в противоположном направлении и небольшого изменения цвета и непрозрачности (вновь созданного) маппера и актера. Обратите внимание, что здесь наш конвейер VTK разделяется на два направления (правое и левое) из одного и того же набора входных данных.

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

Окно вывода теперь показывает обе части каркаса, как и ожидалось.

Теперь мы готовы визуализировать некоторые полезные данные! Чтобы добавить визуализацию давления в левую часть, нам не нужно много делать. Мы создаем новый маппер и также подключаем его к clipperLeft , но на этот раз мы раскрашиваем массив давления. Также здесь мы, наконец, используем диапазон pressureRange , который мы получили выше.

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

Вывод теперь выглядит как изображение, показанное ниже. Давление в середине очень низкое, материал всасывается насосом. Затем этот материал транспортируется наружу, быстро набирая давление. (Конечно, должна быть легенда цветовой карты с фактическими значениями, но я не учел ее, чтобы сделать пример короче.)

Теперь начинается более сложная часть. Мы хотим нарисовать линии тока скорости в правой части. Линии тока генерируются путем интегрирования в векторном поле из исходных точек. Векторное поле уже является частью набора данных в виде векторного массива «Скорости». Итак, нам нужно только сгенерировать исходные точки. vtkPointSource генерирует сферу случайных точек. Мы сгенерируем 1500 исходных точек, потому что большинство из них все равно не будут лежать в наборе данных и будут проигнорированы трассировщиком потока.

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

Затем мы создаем трассировщик потока и устанавливаем его входные соединения. «Подождите, несколько подключений?», — можете сказать вы. Да, это первый фильтр VTK с несколькими входами, с которым мы сталкиваемся. Обычное входное соединение используется для векторного поля, а исходное соединение используется для исходных точек. Поскольку «Velocities» является «активным» векторным массивом в clipperRight , нам не нужно указывать его здесь явно. Наконец, мы указываем, что интегрирование должно выполняться в обоих направлениях от исходных точек, и устанавливаем метод интегрирования Рунге-Кутта-4.5.

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

Наша следующая задача — раскрасить линии тока по величине скорости. Поскольку нет массива для значений векторов, мы просто вычислим значения в новом скалярном массиве. Как вы уже догадались, для этой задачи тоже есть фильтр VTK: vtkArrayCalculator . Он берет набор данных и выводит его без изменений, но добавляет ровно один массив, который вычисляется из одного или нескольких существующих. Мы настраиваем этот калькулятор массива так, чтобы он брал величину вектора «Скорость» и выводил ее как «MagVelocity». Наконец, мы снова вызываем Update() вручную, чтобы получить диапазон нового массива.

 // 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 напрямую выводит полилинии, а vtkArrayCalculator передает их без изменений. Поэтому мы могли бы просто отображать вывод magCalc напрямую, используя новый преобразователь и актор.

Вместо этого в этом обучении мы решили сделать вывод немного лучше, вместо этого отображая ленты. vtkRibbonFilter генерирует 2D-ячейки для отображения лент для всех входных полилиний.

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

Чего сейчас все еще не хватает, и что действительно необходимо для создания промежуточного рендеринга, так это последних пяти строк для фактического рендеринга сцены и инициализации интерактора.

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

Наконец, мы подошли к готовой визуализации, которую я еще раз представлю здесь:

Полный исходный код для приведенной выше визуализации можно найти здесь.

Хороший, плохой, злой

Я закончу эту статью списком своих личных плюсов и минусов фреймворка VTK.

• Плюсы : Активная разработка : VTK находится в стадии активной разработки несколькими участниками, в основном из исследовательского сообщества. Это означает, что доступны некоторые передовые алгоритмы, многие 3D-форматы могут быть импортированы и экспортированы, баги активно исправляются, а проблемы обычно имеют готовое решение в форумах.

• Минусы : Надежность . Объединение множества алгоритмов от разных участников с открытым конвейером VTK, однако, может привести к проблемам с необычными комбинациями фильтров. Мне пришлось несколько раз зайти в исходный код VTK, чтобы понять, почему моя сложная цепочка фильтров не дает желаемых результатов. Я настоятельно рекомендую настроить VTK таким образом, чтобы разрешить отладку.

• Pro : Архитектура программного обеспечения : Дизайн конвейера и общая архитектура VTK кажутся хорошо продуманными, и с ними приятно работать. Несколько строк кода могут дать потрясающие результаты. Встроенные структуры данных просты для понимания и использования.

• Минусы : микроархитектура : некоторые микроархитектурные дизайнерские решения ускользают от моего понимания. Const-корректность почти отсутствует, массивы передаются как входы и выходы без четкого различия. Я облегчил это для своих собственных алгоритмов, отказавшись от некоторой производительности и используя собственную оболочку для vtkMath , которая использует пользовательские 3D-типы, такие как typedef std::array<double, 3> Pnt3d; .

• Pro : Micro Документация : документация Doxygen по всем классам и фильтрам обширна и удобна для использования, примеры и тестовые случаи на вики также очень помогают понять, как используются фильтры.

• Минусы : документация по макросам : в Интернете есть несколько хороших руководств и вводных материалов по VTK. Однако, насколько я знаю, нет большой справочной документации, объясняющей, как выполняются конкретные действия. Если вы хотите сделать что-то новое, ожидайте, что какое-то время будете искать, как это сделать. Кроме того, сложно найти конкретный фильтр для задачи. Однако, как только вы его найдете, документации Doxygen обычно будет достаточно. Хороший способ изучить структуру VTK — загрузить и поэкспериментировать с Paraview.

• Pro : неявная поддержка параллелизации : если ваши источники можно разделить на несколько частей, которые можно обрабатывать независимо, распараллеливание так же просто, как создание отдельной цепочки фильтров в каждом потоке, обрабатывающем одну часть. Большинство крупных задач визуализации обычно попадают в эту категорию.

• Минусы : нет явной поддержки параллелизации : если у вас нет больших делимых проблем, но вы хотите использовать несколько ядер, вы сами по себе. Вам придется выяснить, какие классы являются потокобезопасными или даже допускают повторный вход методом проб и ошибок или путем чтения исходного кода. Однажды я обнаружил проблему распараллеливания в фильтре VTK, который использовал статическую глобальную переменную для вызова некоторой библиотеки C.

• Pro : система сборки CMake : многоплатформенная мета-система сборки CMake также разработана Kitware (создатели VTK) и используется во многих проектах за пределами Kitware. Он очень хорошо интегрируется с VTK и значительно упрощает настройку системы сборки для нескольких платформ.

• Плюсы : Независимость от платформы, лицензия и долговечность : VTK изначально не зависит от платформы и распространяется под очень разрешающей лицензией в стиле BSD. Кроме того, профессиональная поддержка доступна для тех важных проектов, которые в ней нуждаются. Kitware поддерживается многими исследовательскими организациями и другими компаниями и будет существовать еще какое-то время.

Последнее слово

В целом, VTK — лучший инструмент визуализации данных для задач, которые мне нравятся. Если вы когда-нибудь столкнетесь с проектом, который требует визуализации, обработки сетки, обработки изображений или подобных задач, попробуйте запустить Paraview с входным примером и оцените, может ли VTK быть инструментом для вас.