CC BY
58
УДК 620.179.15
С. А. ЗОЛОТАРЕВ1, М. А. МИРЗАВАНД2
БЫСТРАЯ ИТЕРАЦИОННАЯ КИЛОВОЛЬТНАЯ ТОМОГРАФИЯ В КОНИЧЕСКОМ ПУЧКЕ ЛУЧЕЙ
1ГНУ «Институт прикладной физики НАН Беларуси» 2Белорусский национальный технический университет
Создание быстрых параллельных итерационных томографических алгоритмов, основанных на использовании графических ускорителей, которые одновременно обеспечивают минимизацию невязки и тотальной вариации реконструированного изображения является важной и актуальной задачей, имеющей большое научное и прикладное значение. Такие алгоритмы могут использоваться, например, при осуществлении радиационной терапии пациентов, потому что предварительно всегда производится компьютерная томография пациентов с целью более точного определения областей, которые потом будут подвергнуты лучевому воздействию.
Введение
Как в промышленности, так и в медицине широкое распространение получила компьютерная томография с использованием конического пучка рентгеновских лучей (Cone Beam Computed Tomography (CBCT)). Самым распространенными на сегодняшний день методами реконструкции по коническим проекци ям являются методы, основанные на использовании алгоритма Фельдкампа (FDK). Вследствие наметившейся в последние годы тенденции к значительному уменьшению рентгеновской дозы для пациентов, подвергаемых томографическому обследованию, возникает необходимость сокращения, как числа снимаемых проекций, так и уменьшения анодного тока рентгеновской трубки, при котором производится съемка. Например, перед осуществлением радиационной терапии производится CBCT пациентов, для того, чтобы более точно определить области, которые потом будут подвергнуты мощному лучевому воздействию. При компьютерной томографии пациента, как правило, используется порядка 200-300 рентгеновских проекций с лучевой нагрузкой 0.4 mAs на одну проекцию (0.4 mAs/projection). При уменьшении рентгеновской дозы до 0.1 mAs/projection и числа снимаемых проекций до 30-90 алгоритм Фель-дкампа реконструирует сильно деградированные изображения пациентов, которые непригодны для дальнейшего медицинского использования.
Опыт показывает, что исправить эту ситуацию можно путем применения итерационных томографических алгоритмов с использованием в качестве регуляризирующего функционала тотальной вариации. Использование минимизации тотальной вариации реконструируемого изображения позволяет устранить артефакты, вызванные зашумлением рентгеновских проекций, возникающего вследствие малого тока рентгеновской трубки. Кроме того, использование тотальной вариации обеспечивает практически полное устранение полосовых артефактов, вызванных уменьшением числа проекций. Задача состоит в том, чтобы разработать итерационные алгоритмы томографической реконструкции с использованием минимизации тотальной вариации, позволяющие по ограниченному числу проекций (30-90) и рентгеновской дозе 0.1-0.2 mAs/projection обеспечить практически такое же качество реконструируемого изображения, как и метод FDK для 200300 проекций с лучевой нагрузкой на одну проекцию 0.4 mAs.
Однако, итерационные методы томографии затрачивают на реконструкцию трехмерного изображения пациента неприемлемо большое время. Использование графических ускорителей позволяет разработать сверхбыстрые параллельные алгоритмы для трехмерной томографической реконструкции, которые решают эту задачу за временной промежуток от нескольких секунд до одной-двух минут.
Основная часть
1. Эвристический метод итерационной реконструкции
Рентгеновские проекции объекта «голова-шея» для нескольких пациентов были получены на многофункциональном ускорителе для лучевой терапии с модуляцией интенсивности (ЛТМИ) и расширенным визуальным контролем (ЛТВК) Электа Синержи с использованием киловольтного источника. Помимо мега-вольтного рентгеновского источника в этой рентгеновской системе имеется киловольтный источник рентгеновского излучения, предназначенный для осуществления конической томографической реконструкции изображений пациентов с целью их правильного позиционирования перед осуществлением лучевой терапии. В данной работе исследована возможность улучшения качества стандартной томографической реконструкции (FDK) путем использования вместо нее параллельной итерационной реконструкции с использованием библиотеки OpenGL и регуляризацией. Кроме того, было оценено качество реконструкции по малому числу проекций - 90 проекционных видов. Итерационные способы томографии основаны на решении системы линейных алгебраических уравнений (СЛАУ)
Ах = р, (1.1)
где А = (Лу) - проекционная матрица, x = (xi, x2, ..., Xj) -вектор изображения,р = (pi, ..., pI) -вектор проекций. Таким образом, решение основной задачи восстановления изображения по заданному набору проекций р, сводится к решению СЛАУ вида (1), при этом вектор р задан заведомо с некоторой погрешностью. Давайте определим
A1+=Tj A для i = 1, .., I, (1.2) A+j=Y.i A для j = i,..., J. (1.3)
Тогда формула для реконструкции методом SART записывается следующим образом
х(к+1)= х(к)+ Х(к)/Л+] А /Al+((pl - (A1, х(к))),
(1.4)
где параметры релаксации Х(к) представляют собой последовательность вещественных чисел (0 < Х(к) < 2). Мы испытаем в этой работе новый эвристический итерационный алгоритм по теоретическим соображениям близкий к ал-
горитму, опубликованному в работе [1], но отличающийся прежде всего тем, что мы используем коррекцию текущего приближения с помощью минимизации невязки для лучевых интегралов, а не для интенсивностей. Предлагаемый алгоритм реконструкции можно записать так
х(к+1) = х(к)*[1 +Х(к) X, Aj/ (15)
/Ai+((pi - (Al, х(к)))/ (Al, х(к))], .
где 0 < Х(к) < 1. Приведем ниже фрагмент кода на языке СИ, для расчета отношения разности между экспериментальной и текущей модельной проекций к модельной проекции:
for(i= 0; i<pos; i++)
{if(proj_mod[i]> = 1.0) {proj_mult[l]=1.0 +X *(proJ_exp[i]-proJ_mod[i])/proj_mod[i]; if(proJ_mult[i]>2.0)proJ_mult[i]= 2.0; if(proJ_ mult[i]<0.0)proj_mult[i]=fabs(proj_mult [i]); }else proj_mult[i]= 1.0;}.
Алгоритм коррекции, представленный формулой (1.5) является практически таким же удобным для распараллеливания с помощью графической библиотеки OpenGL, как и известный алгоритм SART [2], описываемый формулой (1.4). Прямое проецирование трехмерной текстуры, содержащей изображение текущего приближения осуществляется так же, как и для алгоритма SART, которое подробно описано в монографии [3]. Отличие заключается в том, что для его реализации приходится ввести две дополнительных двухмерных текстуры. В одной из них содержится корректирующее изображение, полученное согласно приведенному выше фрагменту программного кода, а во второй текстуре содержится изображение корректируемого слоя трехмерной текстуры. Корректирующее изображение сначала перспективно проецируется [4] на место, где находится корректируемый слой, а потом перемножается с изображением корректируемого слоя и записывается на его место. Для перемножения двухмерных текстур следует использовать наложение текстуры на объект с режимом GL_MODULATE.
2. Минимизация тотальной вариации
Модель анизотропной тотальной вариации для удаления шума в исходном изображении f размерностью NxN может быть представлена как задача минимизации выражения [5]:
Рис. 2. Ортогональные сечения изображения хоу (слева), yoz (посередине) и xoz (справа). Верхний ряд - реконструкция по 373 проекциям без регуляризации. Нижний ряд - реконструкция с ^ регуляризацией
тт
и
II Итт II -II /-1 ГЩ. + V,,и\\ +— и - / х 111 11 у 11 2"
(2.1)
где - соответственно выбранный положительный параметр, здесь
\\и\\ =
II ||р
1
( \2 I |и(/, ] )| Р
^ 1</, ]<м
для 1 < р < да . (2.2)
Градиенты по х:
Vхи^,у) = и(/, у) - и(/ -1,у), i = 2,...,N, у = 1,...,N.
Ухи(1,у) = 0, у = 1,...,N.
Градиенты по у:
Vуи(/, у) = и(/, у) - и(/, у -1), /' = 1,...,N, у = 2,...,N.
Vуи(/,1) = 0, i = 1,...,N Итерационная схема:
Ькх = си1(У хик + Ькх-1,}{),
Ьку = си1(У уих + ЬХ-1,1
у
Х+1
у
X,
и"" = /--^Х +ууЬу^)
-
где для параметра X > 0
си«с,-) =
1/
Л
с для -1/
(2.3)
для с > у^ _ 1/<с< 1/
для с <
-1/
Рис. 3. Ортогональные сечения изображения хоу (слева), yoz (посередине) и xoz (справа). Реконструкция по 90 проекциям с ^ регуляризацией
объекта «голова и шея»: напряжение на трубке -100 кВ, анодный ток - 10 мА, время экспозиции - 10 мс. На рис. 1 показаны проекции одного из пациентов для углов обзора 0° и 90°. Для реконструкций использовались 373 проекции и 90 проекций, расположенных в угловом диапазоне 201°. Результаты томографических реконструкций показаны ниже на рис. 2-3.
Заключение
В статье предложен разработанный авторами новый параллельный эвристический алгоритм реконструкции для круговой схемы сканирования в коническом пучке рентгеновских лучей. С помощью реконструкций по рентгеновским проекциям реального пациента показано, что он может быть использован в качестве альтернативы алгоритму Фельдкампа, который используется на линейном ускорителе Электа Синержи. Показано, что использование регуляризации методом тотальной вариации уменьшает уровень шума в реконструируемом трехмерном изображении и улучшает его качество, одновременно сохраняя четкие границы.
Литература
1. Lange K., Fessler J. A. Globally Convergent Algorithms for Maximum a Posteriori Transmission Tomography / K. Lange, J. A. Fessler // IEEE Trans. on Image Processing - 1995. - Vol. 8. - No 10. - P. 1430-1438.
2. Andersen A. H., Algebraic reconstruction in CT from limited views / A. H. Andersen // IEEE Trans. Medical Imaging -1989. - Vol. 8. - P. 50-55.
3. Венгринович В. Л. Итерационные методы томографии / В. Л. Венгринович, С. А. Золотарев // Минск: «Белорусская наука», - 2009. - 227 с.
4. Segal M., Fast shadows and lighting effects using texture mapping / M. Segal, C. Korobkin, R. van Widenfelt, J. Foran, and P. E. Haeberli // SIGGRAPH'92. - 1992. - Vol. 26. - P. 249-252.
5. Jia R. Q., A fast algorithm for the total variation model of image denoising / H. Q. Zhao // Adv. Comput. Math. -2010. - Vol. 33. - P. 231-241.
Поступила 10.07.15 После доработки 22.09.15
Сопряженные градиенты для Vx, Vy: -w(2, j) если i -1,
w(i,j) - w(i +1, j) если i = 2,...,N -1, w(N, j) если i = N, -w(2',2)e&riH j -1,
w(hj) ~ MjJ +1) если 7 = 2 ,...,N-1, w(i, N) если j = N,
Итерационная схема обеспечивает сходимость последовательности uk к искомому изображению и*, представляющему собой «обес-шумленное» изображение f [5].
limk^ uk = u* если 0 < < ^ .
Минимизация тотальной вариации для трехмерного изображения осуществлялась на обычном процессоре между корректирующими итерациями, которые были реализованы на графической видеокарте NVIDIA GeForce GTX470 c 1280 Mb текстурной памяти.
3. Результаты
Все проекции были получены на линейном ускорителе последнего поколения Электа Синержи. Параметры для рентгеновской съемки
S. A. ZOLOTAREV, M. A. MIRZAVAND FAST ITERATIVE KILOVOLTAGE CONE BEAM TOMOGRAPHY
Creating a fast parallel iterative tomographic algorithms based on the use ofgraphics accelerators, which simultaneously provide the minimization of residual and total variation of the reconstructed image is an important and urgent task, which is of great scientific and practical importance. Such algorithms can be used, for example, in the implementation of radiation therapy patients, because it is always done pre-computed tomography of patients in order to better identify areas which can then be subjected to radiation exposure.
