Сенсорные системы, 2019, T. 33, № 2, стр. 166-172

ВВЕДЕНИЕ

Метод компьютерной томографии позволяет восстанавливать внутреннюю структуру объекта без его физического разрушения, т.е. описывать локальную способность объекта ослаблять рентгеновское излучение. Если зондирующее рентгеновское излучение монохроматично, то решением задачи компьютерной томографии является пространственное распределение линейного коэффициента ослабления рентгеновского излучения для используемой длины волны. При полихроматическом зондировании, если не применяются специализированные методы реконструкции (Zhao et al., 2018), то восстанавливается распределение “усредненного по спектру линейного коэффициента”. В медицинских томографах применяется полихроматическое зондирование, что в некоторых случаях затрудняет интерпретацию получаемых результатов (Birnbacher et al., 2018).

Задача томографической реконструкции в методе рентгеновской томографии формулируется следующим образом: по зарегистрированным 2D-изображениям восстановить пространственное (3D) распределение коэффициента ослабления, описывающего способность объекта ослаблять рентгеновское излучение. Зарегистрированные позиционно-чувствительным оборудованием пиксельные изображения называют томографическими проекциями. Значение в пикселе – это число рентгеновских квантов, зарегистрированных ячейкой детектора, установленного за объектом. Излучение, сформированное рентгеновским источником, распространяясь прямолинейно, проходит через объект, на своем пути ослабляется материалом объекта, и ослабленное излучение регистрируется позиционно-чувствительным детектором, т.е. матрицей регистрирующих элементов. Каждый из элементов будем описывать пикселем изображения. Говоря о размере элемента матрицы (она же ячейка детектора), мы опустим из рассмотрения функцию размытия точки. Будем считать, что размер ячейки определяет пространственное разрешение метода. Тогда для повышения пространственного разрешения метода томографии следует увеличивать число ячеек регистрирующего оборудования, уменьшая их геометрические размеры. Сегодня в лабораторных приборах используются детекторы с числом пикселей порядка четырех тысяч в одном направлении, т.е. число регистрируемых данных на одну проекцию – 1.6 × 10⁷, а число вокселов восстанавливаемого объема – 6.4 × 10¹⁰. Это и есть число неизвестных в задаче реконструкции. Чтобы уменьшить ошибку при решении обратной некорректно поставленной задачи, потребуем, чтобы число неизвестных в задаче было приблизительно равно числу наблюдаемых значений, т.е. число проекционных углов приблизительно равно числу ячеек детектора в одном из направлений.

Разнотипные подходы применяются для решения задачи томографической реконструкции. Рассмотрим подробнее два подхода – интегральный и алгебраический (Наттерер, 1990). Их принципиальное различие заключается в том, что в интегральных методах реконструкции, использующих аналитическое обращение, задача дискретизуется только на последнем шаге при вычислении интегралов. Аналитическое решение построено в предположении, что функции пространства Радона (пространства проекций) и функция пространства объекта (восстанавливаемое распределение) – непрерывны. Алгебраические подходы к решению задачи реконструкции используют сразу дискретные данные. Измерения – это дискретные изображения. Восстанавливаемый объем покрывается сеткой вокселей, на которой восстанавливается кусочно-постоянная функция. Каждый из этих типов реконструкции обладает своим набором достоинств и недостатков. Интегральные методы выполняются за приемлемое время, сравнимое со временем проведения измерений, и не налагают строгих требований к ресурсам памяти, алгебраические – являются сильно ресурсно-затратными. Однако именно алгебраические методы сегодня привлекают внимание многих групп исследователей. Это обусловлено стремлением снизить дозовые нагрузки при проведении измерений. Есть два способа это сделать – уменьшить число проекционных углов или уменьшить экспозицию. Первый ведет к потере пространственного разрешения, второй ухудшает отношение сигнал-шум. Сильная сторона алгебраических подходов в том, что основанные на них методы реконструкции с регуляризацией решения устойчивы к шуму входных данных (томографических проекций). Поскольку эти методы реконструкции итеративные, то время вычислений значительно больше по сравнению с интегральными методами. Создаются оптимальные версии итерационных алгоритмов (Прун и др., 2013), но существенно ускорить томографическую реконструкцию позволяет использование графических ускорителей (GPU). Их применение позволяет параллельно обрабатывать тысячи пикселей и вокселей (Xu, Mueller, 2007; Buzmakov et al., 2011). Появляются программные пакеты с открытым исходным кодом для выполнения операций прямого и обратного проецирования (основные операции в алгоритмах реконструкции) на GPU, эффективно реализованные и поддерживающиеся сообществом ученых, работающих в области решения задач реконструкции (Van Aarle, 2016). Однако, если все измеренные томографические данные не помещаются в память GPU, то необходимо изменять реализацию существующих алгоритмов реконструкции (Andersson et al., 2016).

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

ОПТИМИЗАЦИОННАЯ ЗАДАЧА ДЛЯ РЕКОНСТРУКЦИИ ПОЛНОГО ОБЪЕМА

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

Математическая модель формирования проекции при монохроматическом зондировании, применяемая для решения обратной задачи, одинакова для обеих, описанных выше, схем сбора томографических проекций и включает следующие предположения. Источник рентгеновского излучения является точечным. Каждой ячейкой детектора регистрируется ровно один бесконечно тонкий луч. Запишем выражение, которое связывает величину сигнала I(r, φ), регистрируемого ячейкой позиционно-чувствительного детектора, с восстанавливаемым распределением величины коэффициента линейного ослабления m(x, y, z) при монохроматическом зондировании:

Здесь r – радиус-вектор, определяющий положение наблюдаемой ячейки детектора или пикселя на изображении 2D-томографической проекции, φ – проекционный угол, под которым зондируется объект, ${{I}_{0}}$ – начальная интенсивность бесконечно тонкого зонда, ${{L}_{{r,\varphi }}}$ – траектория луча, формирующего сигнал в ячейке детектора, и $\mu \left( l \right)$ – коэффициент линейного ослабления текущей точки объема, через которую проходит траектория. Интегрирование ведется вдоль направления распространения луча. Мы ввели декартову систему координат, связанную с неподвижным объектом, описав объект функцией $\mu \left( {x,y,z} \right)\quad$, и описали процедуру получения набора проекций движением системы источник-детектор по круговой траектории. После нормирования зарегистрированного сигнала на начальную интенсивность и взятия отрицательного логарифма приходим к выражению, известному в непрерывном случае, как преобразование Радона:

Важнейшим свойством этого преобразования является его обратимость при выполнении некоторых условий. Это свойство использовано в интегральных методах реконструкции, о преимуществах и недостатках которых говорилось выше. Чтобы перейти к рассмотрению алгебраического подхода, вспомним, что томографические проекции представляют собой изображения с пикселями конечных размеров. Введем вектор измерений P размера ${{N}^{{dx}}} \times {{N}^{{dz}}} \times {{N}^{\varphi }}$, координатами которого являются левые части выражения (2) для всех пикселей детектора и всех проекционных углов. Здесь ${{N}^{{dx}}}$ – число пикселей в линейке детектора, ${{N}^{{dz}}}$ – число строк детектора и ${{N}^{\varphi }}$ – число проекционных углов. Будем использовать воксельное представление реконструируемого объекта. Покроем куб сеткой вокселов. Физический размер воксела свяжем геометрическим соотношением с размером пикселя изображения, формируемого детектором

где SO – расстояние источник – ось вращения, SD – расстояние источник–плоскость окна детектора. Количество вокселов в восстанавливаемом объеме ${{N}^{x}} \times {{N}^{y}} \times {{N}^{z}}$, где ${{N}^{x}} = {{N}^{y}} = {{N}^{{dx}}}$, а ${{N}^{z}} = {{N}^{{dz}}}$. Введем вектор решения $\vec {x}$ размерности ${{N}^{x}} \times {{N}^{y}} \times {{N}^{z}}$, координаты которого есть значения восстанавливаемой кусочно-постоянной функции. Перепишем выражение (2) где – проекционная матрица. Каждая строка матрицы содержит ${{N}^{x}} \times {{N}^{y}} \times {{N}^{z}}$ элементов. Число элементов строки равно числу вокселов восстанавливаемого объема. Число строк равно полному числу измерений ${{N}^{{dx}}} \times {{N}^{{dz}}} \times {{N}^{\varphi }}$. В самом простом случае элементами этой матрицы являются 0 и 1. Если луч, формирующий сигнал в текущей ячейке детектора, пересекает воксел, то элемент матрицы – “1”, в противном случае – “0”. Модели отрезков и весовых фракций также используются при составлении проекционных матриц (Kak, Slaney,1988). Размер системы линейных алгебраических уравнений (СЛАУ) (3) при рассматриваемых нами размерах детекторов составляет порядка нескольких миллиардов уравнений, поэтому прямые методы не используются для их решения. Переформулируем задачу поиска решения СЛАУ как задачу оптимизации. Будем искать такой вектор $\vec {x}$, который минимизирует функционал

Здесь α – параметр, определяющий вклад каждого из слагаемых в оптимизируемое выражение, $R\left( x \right)\quad$ – линейная комбинация регуляризирующих операторов, областью определения которых является векторное пространство решений. Поскольку искажения восстанавливаемого изображения могут иметь разную природу (Nikolaev et al., 2016), параметр α введен таким образом, чтобы регулировать вклад каждого из слагаемых в решаемую оптимизационную задачу.

Будем решать задачу алгебраическими методами на графическом процессоре. Для метода наискорейшего спуска решение на i-м итерационном шаге будет иметь вид

Здесь β – шаг в направлении наискорейшего спуска. Чтобы реализовать итерацию алгоритма на GPU потребуется память под несколько областей размера вектора измерений и несколько областей размера вектора решений. Точное количество областей зависит от метода, используемого для решения оптимизационной задачи. Таблица 1 позволяет провести сравнение требуемого количества памяти для метода градиентного спуска с постоянным шагом, метода наискорейшего спуска и метода сопряженных градиентов. Для эффективных расчетов на GPU для хранения каждой переменной использовалось четыре байта (т.е. переменные хранились в формате FLOAT). Объем памяти, доступный графическому процессору, на сегодняшний день не превышает нескольких гигабайт, т.е. позволяет одновременно загрузить приблизительно миллиард переменных. Как упоминалось выше, количество уравнений в системе имеет примерно такой же порядок (размер вектора измерений), плюс такое же количество неизвестных (размер вектора решений). С учетом значений, указанных в таблице, приходим к выводу о необходимости разбивать итерируемый объем на субобъемы и работать с каждым из них по отдельности.

ПРОЦЕДУРА ВЫДЕЛЕНИЯ СУБОБЪЕМА

На рис. 2 проиллюстрирован предложенный нами ранее метод разбиения объема на субобъемы. Чтобы рассчитать количество слоев субобъема, запишем выражение, которое связывает его $N_{S}^{z}$ и требуемый для вычислений объем памяти ${\Pi }$ в байтах

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

ОПТИМИЗАЦИОННАЯ ЗАДАЧА ДЛЯ СУБОБЪЕМА

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

Поскольку на GPU доступны только данные субобъема, рассчитаем для них лучевые суммы

и будем решать следующую систему уравнений:

Тогда оптимизационная задача для субобъема формулируется так

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

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Для расчета модельных томографических проекций в качестве фантома был использован цилиндр, все поперечные сечения которого одинаковы и представляют собой фантом Шеппа-Логана (рис. 3). Для расчета модельных проекций использовались следующие параметры: число строк детектора – 100, в каждой строке 525 ячеек; физический размер ячейки детектора – 11 мкм; расстояние источник-объект – 55 мм; расстояние источник-детектор – 225 мм; схема измерения – конусная. Проекционные данные рассчитаны для интервала углов 0–180°, величина проекционного угла – 0.2°. Чтобы протестировать места сшивки субобъемов на результате восстановления, мы провели три эксперимента. Восстановили объем целиком, провели реконструкцию, разделив объем на два и четыре субобъема. Расчеты проведены для десяти итераций. Следов сшивки субобъемов на идеальных данных визуально обнаружить не удалось. Чтобы визуально оценить качество реконструкции, мы представляем центральное горизонтальное сечение объекта, восстановленного с использованием четырех субобъемов (рис. 4, а) и двух профилей (б), соответствующих 250 линиям этого сечения (серая линия) и центральному сечению фантома. Чтобы продемонстрировать наличие слабого эффекта сшивки, мы приводим участки вертикальных центральных сечений трех восстановленных объемов, наложив на изображения белые линии равного уровня (рис. 5). На а представлен участок сечения восстановленного целиком объема, на б – результат восстановления с использованием двух субобъемов, место сшивки в районе 50-го слоя. На в представлен результат восстановления с использованием четырех субобъемов, места сшивки в районе 25-, 50- и 75-го слоев. Полученные результаты подтверждают приемлемое качество реконструкции при разбиении описанным методом восстанавливаемого объема на субобъемы.

ЗАКЛЮЧЕНИЕ

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

В заключение хотелось бы отметить, что в июле 2018 г. CERN поместил на своей странице новость о создании регистрирующего устройства нового поколения для использования при решении задач компьютерной томографии в медицинских приложениях (Greenwood, 2018). Новый прибор позволяет регистрировать в каждом пикселе детектора рентгеновский спектр. Это означает, что число измерений, с которыми придется работать при решении задачи томографической реконструкции, увеличится в число раз, равное числу регистрируемых спектральных линий. В связи с этим мы уверены, что создаваемые для работы с субобъемами данных алгоритмы в дальнейшем найдут свое применение при решении задач спектральной компьютерной томографии.

Исследования, связанные с анализом требуемых вычислительных ресурсов для решения задачи томографии, проводились при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ “Кристаллография и фотоника” РАН. Работы по оптимизации необходимых расчетов выполнялись при частичной финансовой поддержке РФФИ (проекты № 17-29-03492 и № 18-52-7819).