CC BY
40
Оптическая томография преломляющих
объектов
Афанасьев В.В., Лебедев A.C.
Институт прикладной математики им. М.В. Келдыша РАН, Лаборатория компьютерной графики и мультимедиа ВМК МГУ {vafanasjev I alebedev} @ graphics, es. msu. ru
Игнатенко A.B.
Лаборатория компьютерной графики и мультимедиа ВМК МГУ ignatenko @ graphics, es. msu. ru
Аннотация. В данной статье рассмотрены основные существующие методы вычислительной томографии. Описываются особенности оптической томографии, и модифицированные для этого случая методы.
Ключевые слова: томография, оптическая томография, ART.
1 Введение
Вычислительная томография - класс задач восстановления внутренней структуры некоторого объёма по набору его проекций. Самой распространённой областью применения методов вычислительной томографии является рентгеновская томография.
Установка для рентгеновской томографии состоит из излучателя, приёмника и места для исследуемого объекта между ними. Излучатель испускает рентгеновские лучи с фиксированной интенсивностью, внутри объекта происходит их поглощение, и приёмник регистрирует остаточную интенсивность лучей. Таким образом получается одна проекция объекта. Излучатель и приёмник вращаются вокруг объекта и снимают набор проекций. Как правило, из-за конструктивных особенностей установок вращение происходит вокруг одной оси на 360°.
Рис. 1. Схема рентгеновского томографа
Возможные конфигурации пучка лучей:
■ Плоский параллельный
■ Плоский расходящийся
■ Объёмный параллельный
■ Объёмный расходящийся
В зависимости от типа пучка лучей, могут применяться различные типы методов вычислительной томографии:
■ Послойная томография. Структура трёхмерного объёма восстанавливается независимо в каждом срезе, перпендикулярном оси вращения приёмника.
■ Объёмная томография. Восстановление структуры объекта происходит сразу во всём его объёме.
Поглощение излучения в материале подчиняется закону Бугера-Ламберта-Бера:
/ =
где I - интенсивность излучения на приёмнике, /0 ~~ интенсивность излучения источника в данном направлении, к (х) - распределение показателя поглощения среды вдоль луча.
2 Классическая вычислительная томография
Пусть на плоскости ХУ задана двумерная функция ^х, у). Рассмотрим всевозможные прямые, лежащие в плоскости ХУ.
Рис. 2. Преобразование Радона для одной прямой
Преобразованием Радона от двумерной функции называется несобственный интеграл от этой функции вдоль прямой:
Новые информационные технологии в автоматизированных системах 2014
fco
f(P + fit) dt
-со
Предполагается, что несобственные интегралы от этой функции вдоль всех прямых сходятся, т.е. функция достаточно близка к 0 на бесконечности. На практике рассматривают функции, заданные в некоторой области. Аналогичным образом вводится преобразование Радона от функции трёх переменных. Физический смысл преобразования Радона в том, что оно отражает показатель поглощения некоторого вещества, распределённого в объёме, вдоль одного луча.
Рассмотрим все направления проекции, для каждого направления проекции - все прямые с этим направлением. Получим двумерную функцию R(k, а), где к - смещение прямой относительно центра координат, а - угол наклона прямой. Такая функция называется синограммой для параллельной проекции.
Рис. 3. Преобразование Радона
Обратное преобразование Радона - восстановление функции f(xly) по функции й(к, се). Оно выражается точной формулой:
^ г 2'я. (" Ш
(2л)2, где
ÎCO
R(s, a)e~™sds
-со
Ha практике все прямые, вдоль которых берётся интеграл, проходят через площадку фиксированного конечного размера (приёмник). Поэтому синограмма для такого случая задаётся в некотором прямоугольнике на плоскости (к, а).
3 Дискретизация
Реальные данные томографии изначально дискретны. Дискретизуются также остальные величины для записи в память и расчётов.
■ Дискретизация по углу поворота объекта относительно приёмника. Имеется конечный набор ракурсов. Как правило, используется выборка ракурсов с равномерным шагом.
■ Дискретизация смещения луча относительно центра координат. Приёмник излучения всегда имеет фиксированное число чувствительных элементов.
■ Дискретизация значений остаточной интенсивности излучения.
■ Дискретизация восстанавливаемой функции / по аргументам. Функция задаётся на сетке с равномерным шагом.
■ Дискретизация значений функции /.
4 Существующие методы
Методы восстановления изображения по проекциям делятся на прямые и итерационные. [Губарени, 1997] К прямым методам относятся методы свёртки и обратной проекции, методы Фурье-синтеза. Итерационных методов существует множество, наиболее известный из них -ART.
4.1 Метод свёртки и обратной проекции
Метод основан на утверждении теоремы о центральном сечении: Фурье-образ каждой проекции является сечением двумерного Фурье-образа исходной функции, проходящим через центр координат. Поэтому задача решается путём применения преобразования Фурье к исходным проекциям, инетрполяции между получившимися сечениями и взятия обратного двумерного преобразования Фурье от получившегося образа. Также к одномерным Фурье-образам проекций могут применяться фильтры.
4.2 Методы ART [Gordon et al, 1970]
Идея этой группы методов состоит в том, чтобы путём минимизации отклонения наблюдаемого значения текущей проекции от истинного значения проекции на каждом отдельном ракурсе и последовательным учётом разных ракурсов добиться минимизации отклонения неизвестной изначально функции / от её истинного значения, которое даёт наблюдаемые значения проекций.
Рассмотрим наиболее простой, аддитивный метод ART в двумерном варианте для параллельной проекции. На вход алгоритму подаётся синограмма (набор из N проекций по М значений в каждой) и геометрические характеристики каждой проекции (направление хода лучей и размер приёмника). На выходе должна получиться искомая функция
[(х,у), заданная на сетке фиксированного размера W х Н. Функция считается кусочно-постоянной, т.е. внутри одной ячейки сетки она равна константе.
Изначально искомая функция инициализируется произвольными значениями (например, 0) на всей сетке. Затем итеративно производятся следующие действия либо заранее заданное число раз, либо до выполнения некоторого условия:
16) Выбор очередного значения проекции из некоторого места некоторого ракурса. Как именно они выбираются, может быть определено специальным порядком обхода проекций и ракурсов, самый простой вариант - последовательная выборка.
17) Генерация прямой, соответствующей данному значению проекции. Параметры прямой вычисляются из геометрических параметров ракурса.
18) Расчёт интеграла Я0 от текущей функции вдоль прямой
19) Вычисление знаковой разности Д = Я — Е0 между известным наблюдаемым значением и текущим интегралом.
20) Проход по сетке и коррекция существующих значений функции вдоль прямой. Коррекция значения функции в одной ячейке сетки производится следующим образом:
Лк* к* I ¡к+1 = /к ^ £
где Лк - коэффициент регуляризации, I - длина пересечения прямой с текущей ячейкой сетки, Ь - общая длина пересечения прямой со всей сеткой.
5 Недостатки классической модели
Классическая модель имеет ряд принципиальных ограничений:
■ Единственный фактор, изменяющий мощность излучения, - это поглощение среды
■ Восстанавливаемая функция должна быть конечна в любой точке, значения проекций также должны быть конечны
6 Проблемы классического метода ART
Классический метод ART основан на последовательной "подгонке" искомой функции к наблюдаемым значениям интеграла. Но, несмотря на то, что на каждой отдельной итерации возможно так изменить функцию, что интеграл от неё вдоль одной прямой в точности совпадёт с наблюдаемым, результат работы алгоритма зависит от последовательности выбора ракурсов и пикселов проекции для выполнения очередного шага. Возможны ситуации, когда суммарная коррекция, например, вдвое превысит необходимую даже при последовательной выборке пикселов с одной проекции. Это связано с тем, что при пересечении соседними лучами одной ячейки сетки на первом шаге коррекция вычисляется точно, а на следующем шаге в эту же ячейку может добавиться ещё одна коррекция от следующей прямой. Ниже приведён пример некорректной работы алгоритма ART для случая восстановления функции, равной 100, на сетке 2x2.
Также классический ART не подходит для реализации на параллельных архитектурах в силу необходимости последовательной выборки отдельных пикселов. Наличие условных операций в алгоритме прохода сетки с пересечением границ, необходимого для расчёта длины пересечения с каждой ячейкой, приводит к тому, что реализация на SIMD архитектуре (например, на видеокарте) была бы неэффективна.
Источник света
л * Ракурсы
/ ^—я X 1 Камера
Рис. 5. Схема установки для оптической томографии
7 Особенности оптической томографии
Оптическая томография принципиально похожа на рентреновскую
томографию, но имеет свою специфику.
■ Каждое значение проекции - не значение интеграла от функции, а остаточная яркость источника света с учётом затухания из-за поглощения среды и пропускания границы между средами по Френелю. Также вклад в яркость вносит отражение от всех поверхностей.
■ На проекции возможно нулевое значение яркости, что соответствует бесконечному поглощению. Бесконечное поглощение согласованно выражается нулевыми значениями яркости на всех ракурсах в соответствующей области. В классической же модели берётся интеграл от показателя поглощения вдоль луча. Вследствие этого, любое значение проекции конечно. Поэтому в линейной модели не существует
такого набора проекций, которые согласованно отображали бы бесконечное поглощение в какой-то области объёма.
■ Существует эффект преломления света на поверхностях раздела разных сред, часть направлений обзора внутри преломляющей среды отсутствует. Внутри преломляющей среды ракурсы распределены неравномерно, шаг по углу переменный.
■ В общем случае на одном ракурсе лучи, которые в воздухе шли параллельно, преломляются на границе среды и проходят через целевой объём сразу с нескольких направлений.
■ Шаг между лучами внутри преломляющей среды переменный и зависит от угла падения лучей на преломляющую поверхность.
8 Источник освещения
В качестве источника освещения лучше всего подходит некоторая светящаяся площадка, расположенная за исследуемым объектом. Таким образом, минимизируется количество света, отражённого от поверхности объекта. Распределение освещённости этой площадки должно быть известно. Существуют следующие подходы к моделированию освещения:
■ Настройка параметрической модели освещения на основе фотографий фона [Афанасьев и др., 2013]. Такой способ эффективен в случае существенного вклада в яркость объекта областей фона, лежащих за пределами области видимости камеры. Однако сама параметрическая модель может вносить некоторую погрешность в освещение видимой области.
■ Использование фотографии фона напрямую. Строится некоторая геометрическая поверхность, которая является источником света, яркость каждой её точки определяется фотографией фона с камеры сканера, сделанной до установки сканируемого объекта. Такая модель точно описывает видимую часть источника света.
■ На практике используется комбинация первых двух подходов.
9 Предложенный метод
Предлагается модифицированный алгоритм ART на основе алгоритма трассировки лучей. Процедуры, связанные с томографией, встраиваются в процесс трассировки одного луча, для этого подходит любая конфигурируемая система трассировки лучей. В качестве такой взята система трассировки лучей на графических процессорах NVidia Optix. [Parker et al, 2010].
Входные данные:
■ Набор фотографий объекта, снятый с разных ракурсов при повороте объекта на 360° вокруг вертикальной оси
■ Геометрические характеристики кадров (положение и направление ракурсов, размеры плоскости проекции)
■ Фотография источника света (фона) с той же камеры и того же ракурса
■ Геометрические модели объектов сцены
■ Показатели преломления прозрачных объектов
■ Положение и размеры исследуемого объёма, заданного параллелепипедом
Результат:
■ Функция распределения показателя поглощения среды в целевом объёме, заданная на сетке
Изначально вся сетка инициализируется значением 0. Алгоритм заключается в последовательном учёте кадров. Одна итерация состоит из следующих этапов:
21) Фотореалистичная визуализация сцены с учётом текущего распределения показателя поглощения. Первый шаг трассировки лучей. Для каждого луча производится:
a. Пересечение луча с геометрией сцены
b. Расщепление на отражённый/преломлённый, вычисление коэффициентов Френеля для преломляющих поверхностей
c. Вычисление интеграла от функции на сетке вдоль луча
d. Пересечение с источником света, вычисление его яркости
e. Вычисление итоговой яркости, запись в пиксел на плоскости проекции
22) Вычисление на основе визуализированного кадра и фотографии карты коррекции - изображения размером с фотографию, содержащего значения необходимой коррекции для каждого пиксела
23)Коррекция значений показателя поглощения на основе имеющейся карты коррекции. Второй шаг трассировки лучей. Для каждого луча производятся те действия из первого шага трассировки, которые необходимы для коррекции искомой функции: пересечение с геометрией и проход сетки.
24)Применение глобальной коррекции к функции (сложение).
Основные отличия алгоритма от классического ART заключаются в следующем:
■ Коррекция функции накапливается в отдельной сетке за один кадр, нормируется на число модификаций каждой ячейки сетки. Это делается для того, чтобы избежать проблемы ART с многократной коррекцией одной ячейки сетки, и для работы в параллельном режиме.
■ Коррекция вычисляется на основе значений яркости фотографии и визуализации, а не интегралов от функции поглощения. Это приводит к большей устойчивости алгоритма при нулевых значениях яркости (что соответствует "бесконечному" значению функции в сетке)
■ При проходе лучом сетки не используется честное пересечение луча с
клетками сетки. Вместо этого равномерно выбираются точки вдоль луча, в них происходит взятие значения функции, или коррекция соответствующей ячейки сетки. Это позволяет отказаться от условных переходов в основном цикле работы, что увеличивает скорость на 81МБ архитектуре.
10 Результаты
Предложенный алгоритм был опробован на реальных данных для построения моделей дефектов в драгоценных камнях. Алгоритм является частью более сложной системы, включающей также предобработку входных данных, сегментацию результата оптической томографии и полигонизацию моделей дефектов. Во многих случаях результат алгоритма имеет такое же качество, как результат ручного построения моделей по тем же входным данным.
Рис. 6. Результат алгоритма оптической томографии в одном срезе
Также проведено тестирование алгоритма на синтетических данных,
которое выявило некоторые проблемы:
■ В случае полного поглощения света областью дефекта на восстановленной функции поглощение распределяется неравномерно в этой области: меньше по центру и больше по краям.
■ Дискретность сетки вызывает ложные коррекции для лучей, близких к резким краям тёмных областей на фотографиях. Это приводит к образованию «вееров» рядом с некоторыми дефектами на реконструированной функции.
■ Низкая скорость сходимости в случае наличия тонких продолговатых дефектов.
11 Список литературы
[Афанасьев и др., 2013] Афанасьев В.В., Игнатенко А.В., Тисевич И.О. Калибровка источника света в сцене с прозрачным объектом по фотографии // 23-я Международная конференция по компьютерной графике и зрению ГрафиКон'2013. Владивосток., 2013. С. 215-218.
[Губарени, 1997] Губарени Н.М. Вычислительные методы и алгоритмы малоракурсной компьютерной томографии. Киев: Наукова думка, 1997.
[Gordon et al, 1970] Gordon, R; Bender, R; Herman, GT. Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography // Journal of Theoretical Biology 29 (3), 1970. pp. 471-81
[Parker et al, 2010] Parker S., Bigler J., Dietrich A. OptiX: a general purpose ray tracing engine //ACM SIGGRAPH 2010 papers, Article No. 66
