Современные томографы позволяют получать высокоточные трехмерные изображения внутренней структуры большого размера, предоставляя ценную информацию о геометрии, составе и дефектах исследуемых образцов. Размеры одной реконструкции обычно колеблются от 512 мегабайт до 1 терабайта. Для анализа таких данных используются специализированные инструменты, но традиционная визуализация трёхмерного объема реконструкции до сих пор является важным этапам оценки качества реконструкции и её интерпретации специалистом.

В сегодняшней статье мы расскажем, как нам удалось обойти все трудности и сделать визуализатор, работающий в режиме реального времени. Одна из важных задач визуализатора, созданного командой Smart Engines, - отображение процесса реконструкции, которая тоже потребовала креативного подхода, так как все вычислители нагружены выполнением другой вычислительно затратной задачи - выполнением реконструкции, а самих данных для отображения еще частично нет.

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

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

2. Реконструированные изображения содержат различные искажения (например, кольцевые артефакты, полосы) или шумы. Наличие артефактов и шумов затрудняют визуальную интерпретацию и требуют предварительной обработки объема для визуализации.

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

4. Трехмерный массив данных трудно целиком воспринимать человеку. Ограничения человеческого восприятия для интуитивно понятной визуализации требуют комбинирования различных методов - сечений, проекций, трехмерного рендеринга и т.д.

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

При разработке собственного 3D визуализатора мы решили остановиться на методе визуализации трехмерных изображений (другими словами, raycast), как наиболее универсального метода. Метод можно применять для визуализации как в непрозрачном, так и прозрачном виде. В его основе лежит расчет прохождения множества виртуальных лучей света сквозь объем, который надо визуализировать. Результат выглядит как полупрозрачный трехмерный объект. Можно интерактивно менять параметры визуализации, подбирая нужный контраст и прозрачность для выделения нужных деталей в реальном времени. В данной статье мы описали, как мы смогли реализовать такой метод, и с какими ограничениями используемых подходов мы столкнулись.

Сначала опишем пошагово сам метод raycast.

1. Для каждого пикселя результирующего 2D изображения испускается виртуальный луч, который проходит через объемные данные (3D массив вокселей).

2. Вдоль траектории луча с определенным шагом выбираются точки сэмплирования. В этих точках интерполируются значения соседних вокселей для получения локальной плотности/интенсивности.

3. Полученные значения плотности преобразуются с помощью передаточной функции (transfer function) в оптические свойства - цвет и прозрачность.

4. Финальный цвет пикселя определяется путем композитинга (смешивания) оптических свойств вдоль луча согласно модели объемного рендеринга, учитывающей ослабление и излучение света в полупрозрачной среде.

5. Дополнительно могут учитываться эффекты освещения и тени для повышения реалистичности и выразительности визуализации.

Ключевыми параметрами метода raycast являются передаточная функция (transfer function), определяющая цветовое отображение значений плотности, а также шаг сэмплирования (volume sampling), влияющий на качество и скорость рендеринга, как изображено на Рис. 3.

При выборе инструментов для создания визуализатора мы отталкивались от необходимости пользователям с разными вычислительными ресурсами запускать нашу программу для томографической реконструкции STE. Вычислительные ресурсы пользователей могут варьироваться от самых слабых компьютеров без видеокарт до серверных станций с несколькими мощнейшими видеокартами или станций с российскими вычислительными платформами “Эльбрус” или КОМДИВ. Для поддержки такого разнообразия вычислителей мы реализовали метод визуализации raycast как на CPU, так и на GPU с помощью OpenGL. В дальнейшем будем говорить о нашей реализации raycast на GPU.

Подготовка к визуализации

Перед выполнением каких-либо вычислений  на видеокарте необходимо подготовить и заполнить структуры с данными. OpenGL предоставляет широкий набор структур данных на GPU, из которых для хранения трехмерного массива томографической реконструкции лучше всего подходит 3D текстура (GL_TEXTURE_3D в терминах OpenGL).

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

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

Стандартные объемы оперативной памяти видеокарт потребительского сегмента равны 4Gb, 6Gb, 8Gb, 12Gb, 24Gb. Для хранения широкодиапазонных значений реконструкции нам требуется использовать значения с плавающей точкой. OpenGL позволяет использовать 16-битные, 32-битные и 64-битные значения. На основе размера значения одного вокселя и длины ребра визуализируемого куба можно оценить потребляемую память и, как следствие, максимальный размер визуализируемого куба реконструкции относительно размера видеопамяти.

Мы составили таблицу максимальной длины ребра визуализируемого куба реконструкции для различных объемов видеопамяти и размеров значения вокселя. Значения в ячейках таблицы представляют собой максимальную длину ребра куба, которую можно визуализировать при заданных параметрах.

Весь столбец "1 байт" был написан курсивом, так как по нашим экспериментам его диапазона значений от 0 до 255 не хватало для качественной визуализации. Разницу же между использованием 2 байт и 8 байт на значение вокселя визуально сложно отличить. Таким образом, формат float16 (2 байта на значение вокселя) для реализации был выбран по следующим причинам:

1. Обеспечивает достаточный диапазон значений для качественной полутоновой визуализации.

2. Позволяет визуализировать кубы большего размера по сравнению с форматами float и double при том же объеме видеопамяти.

3. Визуально разница между float16 и float/double незаметна.

При этом, обычно мы работаем с детекторами размером от 512 на 512 пикселей до 3000 на 3000. Так что мы почти полностью покрываем одной потребительской видеокартой потребности в размере визуализируемого объема.  Однако просто взять и потребовать у видеокарты отдать всю свою оперативную память нельзя, так как всегда есть некоторые фоновые процессы, потребляющие часть памяти. Для динамического определения максимального размера возможного визуализируемого куба мы используем GL_PROXY_TEXTURE_3D, который позволяет опросить видеокарту на возможность выделения 3D текстуры.

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

Для реализации метода объемного рендеринга с помощью рейкастинга мы использовали шейдеры OpenGL - специальные программы, исполняющиеся на графическом процессоре (GPU), что обеспечивает значительно более высокую производительность по сравнению с рендерингом на центральном процессоре. В нашей реализации используется фрагментный шейдер, вызываемый для каждого пикселя выходного изображения. Он выполняет трассировку соответствующего луча через трехмерный объем реконструированных данных.

Ключевыми этапами этого шейдера являются:

1. Вычисление направления луча на основе координат пикселя и параметров виртуальной камеры.

2. Определение точек входа и выхода луча в объеме данных на основе расчета пересечения луча с ограничивающим параллелепипедом (bounding box) визуализируемого объема.

3. Трассировка луча внутри объема с фиксированным шагом. На каждом шаге:

   • Выборка значения плотности из 3D текстуры реконструкции в текущей точке.

   • Применение передаточной функции для получения цвета и прозрачности.

   • Смешивание полученного цвета с уже накопленным цветом луча согласно уравнениям модели объемного рендеринга.

4. Запись финального цвета и прозрачности луча в выходной пиксель.

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

Первый режим визуализации моделирует построение линеаризованных проекций, и вид передаточной функции выражается формулой

где C(r) - значение цветового вектора в точке пространства выражаемого вектором r, R(r) - значение реконструкции. И функция смешения вдоль семплируемого луча

Изображение с первым режимом нашей визуализации ниже.

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

а функция смешивания имеет вид

где αi - виртуальный коэффициент прозрачности начинающийся от 0, когда этот коэффициент достигает близкое к 1 значение, то семплирование луча останавливается. Изображение со вторым режимом нашей визуализации ниже.

Третий наш режим продолжает идею alpha наложения цветных вокселей друг на друга, лишь изменяя передаточную функцию на спектральную функцию, где 0 - это синий цвет, 1 - красный, 0.5 примерно желтый

где R_max - максимальное значение во всем визуализируемом объеме. Изображение с третьим режимом нашей визуализации ниже.

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

Создание эффективных шейдеров для визуализации объемных данных требует тщательного баланса между качеством изображения и производительностью рендеринга. Увеличение числа шагов трассировки лучей, сложности передаточной функции повышает (но может и не повысить) визуальное качество, но может снизить интерактивность из-за большего числа выполняемых инструкций на GPU. 

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

Такой шаг по нашим экспериментам позволяет ускорять визуализацию без значимого ухудшения визуализации. 

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

Заключение

В итоге нам удалось создать эффективный инструмент интерактивной визуализации томографических данных на основе метода рейкастинга с использованием возможностей OpenGL. Визуализатор показал высокую производительность и гибкость в настройке на различных GPU, обеспечивая качественное отображение сложных внутренних структур исследуемых образцов. Разработанные подходы, такие как динамическое обновление текстур и потоковая визуализация, позволили органично интегрировать инструмент в нашу программу для обработки, визуализации и анализа томографических данных STE. Видео с результатами томографической реконструкцией можно посмотреть в нашем ютуб-канале.