CC BY
44
УДК 681.3.06
Р. Ф. Гибадуллин, А. А. Максимов, А. А. Новиков, М. Ю. Перухин
РЕКОНСТРУКЦИЯ ТОМОГРАФИЧЕСКИХ СНИМКОВ С ПРИМЕНЕНИЕМ МНОГОПРОЦЕССОРНЫХ СИСТЕМ
Ключевые слова: мультипроцессоры, параллельные алгоритмы, компьютерная томография, реконструкция.
Реконструкция томографического изображения имеет широкий спектр применений в медицине, технике, промышленности. Алгебраические методы реконструкции, использующиеся для решения реконструкции томографических изображений, довольно медленные. Статья охватывает один из видов мультипликативного алгебраического метода реконструкции. Было показано, что при применении параллельной реализации этого алгоритма увеличивается скорость вычислений без ущерба качества реконструкции.
Keywords: Multiprocessors, parallel algorithms, computed tomography, reconstruction.
Reconstruction of a tomographic image has wide range of applications in medicine, engineering, industry. Algebraic reconstruction methods used to solve the reconstruction of tomographic images are rather slow. The article is covered one of the types of the multiplicative algebraic reconstruction method. The parallel implementation of the algorithm is increased the speed of computation without compromising of the reconstruction quality.
Введение
В широком классе задач медицины, физики и техники исследователей интересует пространственная структура неоднородных объектов, недоступных для локальных измерений. Информацию о локальных характеристиках таких объектов приходится получать по результатам измерений их косвенных проявлений, решая обратную задачу. Исследуемый объект облучается извне или сам является излучателем [1-3].
Реконструкция томографического изображения один из известных и доступных методов неразру-шающего контроля. Этот метод включает в себя восстановление поперечного сечения исходного объекта по полученным проекционным данным. Проблема реконструкции образа имеет разные решения. Одним из них является метод мультипликативной алгебраической реконструкции (MART -multiplicative algebraic reconstruction technique). [4-5]
Принцип компьютерной томографии
Если М общее количество лучей от всех проекций, полученных от сканеров, N - общее количество вокселей в сетке, если N >> M, то уравнение для реконструкции изображения проекции можно записать в следующем виде:
где ^ является проекцией данных, а W представляет собой матрицу порядка MxN, содержащая вес Wij для каждого вокселя. Вектор f (вектор неизвестных) это столбец сохранения перестраиваемого изображения и e - критерий останова (требуемая точность) [6].
Вышеуказанное выражение можно записать в развернутом виде, как систему линейных уравнений:
WMlfl+WMlh + WM?.f?. Н-----^WM*tfN = Фм
Вес wij представляет собой вклад i-го луча из точки зрения проекции на j-м вокселе креста объек-
та поперечного сечения. Понятно, что веса несут решающую роль в решение этой системы уравнений. Вес wij точно представляет влияние вокселя £ на луче г- проходящего через пиксель р (рис. 1).
Рис. 1 - Луч Г/ проходящий через пиксель
Значение wij геометрии вычисляется с помощью сканирования геометрических параметров. Используем для расчета метод, основанный на следующей идее. Если луч г- проходит через пиксель рь то соответствующий wij равен 1, в противном случае равняется 0.
Это позволяет использовать методы хранения матрицы веса, которые уменьшают объем памяти требуемый для хранения веса матрицы. В полученной системе число неизвестных не совпадает с числом уравнения. Следует отметить, что полученная система линейных уравнений, некорректна. Эта система неопределённостей и, следовательно, имеет бесконечное множество решений. Некорректные условия подразумевают, что мы не можем использовать прямые методы решения системы линейных уравнений [1].
Последовательный алгоритм мультипликативный алгебраической реконструкции (SvMART)
Пусть S некоторое множество лучей: 5 = 4,4,4, — . Множество S имеет р-подмножеств одинакового размера где, =
{4,4,-,4), ^2 = (4+1,4+2, -,4Л и так далее, где
t _ ^/р . p является суммарным числом видов проекции. В MART только одно подмножество следует принимать во внимание при реконструкции [6-8].
Выходной сигнал этого подмножества работает в качестве входных данных для следующего подмножества. Основные шаги MART конкретной итерации приведены ниже в виде псевдоалгоритма:
MART&U),f0,Npa,k,fj)
Входные данные-^,/0,Npa, Я; Выходные данные-fj <pig, - проекционные данные полученные из эксперимента. /0 - Начальное приближение значений пикселей. Я - параметр релаксации. Npa- Общее количество просмотров проекционных. Sg, - набор уравнений. For j = 1 to Npa:
Для подмножества (Se) Шаг А:
Вычислить W в соответствии с уравнением
W:
1, если луч проходит через 0 , в противном случае
= Г1, есл = 1 0 ,
Шаг В:
Для каждого луча (¿б)подмножества 5е Подсчитать числовую проекцию
Рассчитать коррекцию как Аф^д = —
Шаг С:
Для каждой ячейки (])
Если в Wie j ф 0;
^new _ ^old * (i
-Я*(1- Афш)) fj_f]n
Параллельный алгоритм мультипликативный алгебраической реконструкции (P-SvMART)
В этом разделе рассматривается реализация параллельного алгоритма P-SvMART для конусного сканирования объекта по круговой траектории сканирования.
Алгоритм, P-SvMART, можно разделить на три части:
а) генерация весовой матрицы,
б) вычисление коррекции элементаОегт)
в) обновление значения поля.
Основные этапы параллельного P-SvMART могут быть описаны следующим образом.
1. Общее число лучей Предположим, что каждый вид проекции имеет Ыгау число лучей и существуют процессоры N доступные для параллельного применения. Предположим, что £ представляет собой набор общего числа лучей, присутствующих в свете проекции так, что:5 = {51,52,53,^,5и} где
= {1,2,3.....р}
$2 = {р + 1, Р + 2, р + 3,... ,2 * р},
SN _ {(N - 1) * р + 1,(N - 1) * р + 2,(N - 1) * р + 3, ...,N * р]
Nrav
где р _-.
r N
2. Формирования веса матриц.
Веса можно рассчитать с помощью параметров сканирования,
где Npa - общее число просмотров проекции, Nray -общее количество лучей в каждом ракурсе проекции, Npix - общее количество вокселей на каждой стороне сетки объекта, Scob -расстояние между источником и центром объекта, Scd - расстояние между источником и детектором, d - расстояние между двумя детекторными пикселями, и поле, над которым были собраны проекционные данные.
Используем распределение один-к-одному для подмножеств S устанавливаемых на процессоры. Каждый процессор вычисляет вес каждого луча, принадлежащего к подгруппе. Эти веса будут в дальнейшем использоваться для расчета терминов коррекции. Алгоритмический шаг (а) можно записать в виде:
P-SvMART (Npa, NMy, Npix, Scob, Scd, d, span, W) Входные данные: Npa, Nray, Npix, Scob, Scd, d, span Выходные данные: W
Npa - общее количество просмотров проекционных окон. Nray - общее количество лучей в каждом ракурсе проекции окон. Npix - общее количество вок-селов на каждой стороне сетки. Scob - расстояние от источника до объекта. Scd - расстояние от источника до детектора. d - расстояние между двумя детекторными пикселями. W - вес матрицы. Span - расстояние над спроецированными данными. Шаг A: Для каждого луча (id)
Вычислить уравнение луча Вычислить точку пересечения луча с объектом сетки веса.
w,
= Г1, если луч проходит через у = 1 0 , в противном случае
3. Вычисление коррекция Веса полученные из шага (А) используются для расчета поправочного члена.
Каждый луч имеет свой собственный поправочный член и этот расчет независим от других лучей. Мы будем использовать распределение один-к-одному для назначения процессора на этот шаг. Алгоритмический шаг (В) можно записать в виде: Р-МАЯТ (фш,Г0^ра,Л,Г1). Входные данные: ф1д,и,№ра,Х. Выходные- fj Шаг Б: Для каждого луча (¿б)подмножества Бд Рассчитать числовую проекцию по
ФГв =т>1= гЩЬ
Рассчитать коррекцию по Аф^д = —
Рассчитанная на шаге 2 коррекция - это эффект луча на воксели, через которые он проходит. Следующий шаг (С) включает в себя операции записи значения поля объекта, следовательно, он будет выполняться последовательно. Покажем это в следующей форме: Шаг 3:
Для каждой ячейки (]) Если ф 0
fnew _ f
Jj _ Jj fj_fj
j
new
* (1 — Я *
(wie,j)max
* (1 - Афш))
Результаты
Алгоритм был протестирован на численном фантоме размером 512х512х512. Исходные изображения числового фантома, показано на рисунке 2а. Численный фантом протестирован с различными разрешениями для проверки точности, качества обработки и скорости алгоритма.
а б
Рис. 2 - Изображения числового фантома: a) исходное изображение числового фантома; б) реконструированное изображение числового фантома
Алгоритм был реализован на мультипроцессорной системе. Качество реконструкции остается таким же, при обработке одним процессором, или при использовании нескольких процессоров. Полученное реконструированное изображение показано на рисунке 2б. Используя для реконструкции изображения параллельный алгоритм мультипликативный алгебраической реконструкции (P-MART), мы получили увеличение скорости обработки изображения более чем в 3 раза при использовании 8 процессоров (типа) CPU. При реализации этого алгоритма даже на одном процессоре GPU (использовались Tesla M2070) ускорение составит более 16,5 раз. Более подробные результаты показаны в таблице 1 и изображены на рисунках.
Таблица
Количество использованных CPU Время выполнения реконструкции с использованием CPU ^е^иды) Количество ис-пользованных GPU Время выполнения реконструкции с использованием GPU (cекунды)
1 4796 1 290
2 2915 2 242
4 1989 4 191
8 1585 8 169
с -!-1-1-1-'-'-1-
о 1 г 3 4 & 6 т 8
Нолличество GPU
Рис. 3 - Время выполнения от количества используемых GPU для P-SvMART
5СО
О--J-1-L- '
о 2 ¿t a а
Колличестзо CPU
Рис. 4 - Время выполнения от количества используемых CPU для P-SvMART
Литература
1. Терещенко С.А. Методы вычислительной томографии. -М.: ФИЗМАТЛИТ, 2004. - 320 с.
2. Хермен Г.Т. Восстановление изображений по проекциям. Основы реконструктивной томографии. - М.: Мир, 1983 (перевод с английского). - 352 с.
3. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. - М.: Мир, 1991. - 368 с.
4. Mishra D., Muralidhar K. and Munshi P. A robust MART algorithm for tomographic applications. - Numerical Heat Transfer Part B 35, 1999: 485-506.
5. Бакушинский А.Б., Гончарский А.В. Итеративные методы решения некорректных задач. - М.: Наука, 1989. - 130 с.
6. Gordon R., Bender R., Herman G. Algebraic reconstruction technique (ART) for three-dimensional electron microscopy and X-ray photography. - Journal of Theoretical Biology 29, 1970: pp. 471-481.
7. Гибадуллин Р.Ф., Новиков А.А., Хевронин Н.В., Перухин М.Ю. Разработка параллельного модуля генерации защищенной картографической базы данных // Вестник технол. ун-та. - 2016. - № 10. - С.102-105.
8. Гибадуллин Р.Ф., Савельев А.Г., Перухин М.Ю. Ускорение обработки SQL-запросов к базам данных на GPU посредством аппаратно-программной платформы NVIDIA CUDA // Вестник технол. ун-та. - 2016. - № 20. - С.110-116.
© Р. Ф. Гибадуллин - к.т.н., доц. каф. компьютерных систем КНИТУ им. А.Н. Туполева - КАИ (КНИТУ-КАИ), landwatersun@mail.ru; А. А. Максимов - магистрант, асс. той же кафедры; А. А. Новиков - асс. той же кафедры; М. Ю. Перухин - к.т.н., доц. каф. автоматизированных систем сбора и обработки информации КНИТУ, perukhin@inbox.ru.
© R. Gibadullin - PhD associate professor of computer system department of KNRTU named after A.N.Tupolev - KAI (KNRTU-KAI), leonov.artem.d@gmail.com; А. Maksimov - master student; assistant of the same Department; А. Novikov- assistant of the same Department; M. Perukhin - PhD, associate professor of automated systems for the collection and processing of information department of KNRTU, perukhin@inbox.ru.
