Коррекция движения при томографических и планарных радионуклидных исследованиях Текст научной статьи по специальности «Математика»

Сер. 10. 2011. Вып. 1

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

УДК 517.9+519.6

Е. Д. Котина, К. М. Максимов

КОРРЕКЦИЯ ДВИЖЕНИЯ ПРИ ТОМОГРАФИЧЕСКИХ И ПЛАНАРНЫХ РАДИОНУКЛИДНЫХ ИССЛЕДОВАНИЯХ

Введение. Важным этапом при обработке данных радионуклидных исследований является обнаружение и коррекция движения пациента во время обследования, поскольку даже небольшое смещение пациента или его внутреннего органа во время сбора проекционных данных может повлиять на достоверность результатов диагностики, а полностью избежать перемещений различного характера при этом практически невозможно. Были рассмотрены радионуклидные изображения, полученные методами однофотонной эмиссионной томографии (ОФЭКТ) и динамической сцинтиграфии при помощи программного комплекса «Диагностика». В настоящей работе предлагаются алгоритмы коррекции движения пациента и его отдельных внутренних органов. Впервые предложено использовать понятие оптического потока для данных радионуклидных исследований.

Коррекция движения при ОФЭКТ. Сбор диагностических данных здесь, как правило, занимает более 30 мин. Для достижения наилучшего качества получаемого изображения требуется полная неподвижность пациента, поскольку даже небольшое смещение исследуемого органа во время сбора проекционных данных может стать причиной искажения изображения или появления «артефактов» [1, 2]. В случае ОФЭКТ важную для коррекции движения роль играет сам характер сбора данных, так как заранее подразумевается присутствие движения, обусловленного вращением камеры вокруг пациента. Потому для математической постановки и дальнейшего рассмотрения задачи для этого случая введем две системы координат (рис. 1): систему координат (x,y), связанную с детектором, вращающимся по круговой орбите вокруг центра неподвижной системы координат (x',y',z'), связанной с гантри гамма-камеры. В таком случае для точки 3-мерного пространства проекционная координата будет определяться следующим образом:

{x = A sin ($ — ш),

У = z',

где со = arctan (у1 / х')] $ - угол наблюдения; А = \]х/2 + у/2. Рассмотрим относительное перемещение точки в проекционных координатах между двумя последовательными

Котина Елена Дмитриевна — кандидат физико-математических наук, доцент кафедры теории управления факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 63. Научные направления: математическое моделирование, численные методы, методы оптимизации. E-mail: ekotina123@mail.ru.

Максимов Константин Михайлович — аспирант кафедры теории систем управления электрофизической аппаратурой факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 3. Научный руководитель: кандидат физико-математических наук, доц. Е. Д. Котина. Научные направления: математическое моделирование, численные методы и комплексы программ. E-mail: maximovkm@gmail.com.

© Е. Д. Котина, К. М. Максимов, 2011

кадрами:

{ Ах = А сов ($ — ш)Ад, \Ду = 0.

Таким образом, в отсутствие движения траектория проекции точки системы (хХ, у', г') должна представлять собой синусоиду относительно оси х и прямую относительно оси у.

Далее под поперечным движением будем подразумевать смещения положения исследуемого органа, параллельные плоскости (х',у'), под продольным - параллельные оси г'. К заметным искажениям реконструированного изображения приводит только движение, амплитуда которого превосходит приблизительно 3.5 мм в продольном направлении и 13 мм в поперечном [3].

Предлагается разбить процесс коррекции на следующие этапы: 1) определение наличия движения и момента (группы кадров), когда оно произошло; 2) произведение коррекции при необходимости. Такое разбиение позволяет экономить время и вычислительные ресурсы, необходимые для коррекции движения.

Для решения задачи определения наличия и момента движения будем использовать метод синограмм и линограмм [4]. При достаточной малости отверстий коллиматора гамма-камеры можно предположить, что один ряд проекционного изображения представляет собой проекцию распределения радиофармацевтического препарата (РФП) в слое на той же высоте. Если фиксировать эту высоту и последовательно расположить друг под другом соответствующие ряды проекционных кадров, то получившееся изображение, синограмма, будет отображать синусоидальные траектории движения точек относительно оси х проекционной системы координат. При поперечном движении изучаемого органа происходит отклонение траекторий от синусоид и, как следствие, нарушается плавность переходов между рядами синограммы, что становится заметно визуально. Таким образом, на основании построенной синограммы можно определить, между сбором данных с каких углов наблюдения произошло поперечное движение. Для обнаружения продольного движения, в свою очередь, строится линограм-ма - изображение, столбцы матрицы которого есть последовательно расположенные проекционные кадры, для которых было произведено суммирование значений счета для каждого ряда. Таким образом, полученное изображение, в силу сохранения продольной проекционной координаты, должно представлять собой набор прямых линий. Следовательно, наличие продольного движения можно определить как визуально

заметный переход от изображенных на линограмме прямых линий к ломаным. В итоге на основании полученных синограммы и линограммы (рис. 2) можно сделать вывод

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

Рис. 2. Синограмма (а) и линограмма (б)

Существует ряд методов, используемых для количественной оценки движения и автоматизированной коррекции данных ОФЭКТ. Наиболее широко применяются методы расходящихся квадратов [5], двумерного соответствия [6] и методы, основанные на использовании функции взаимной корреляции [7, 8]. Так как первые два способа коррекции могут быть применены только для коррекции данных исследований перфузии миокарда, а также учитывая результаты сравнения [9, 10] методов между собой, предлагается использовать метод функции взаимной корреляции, как наиболее достоверный и независимый от типа исследования. Дискретная функция взаимной корреляции между двумя одномерными последовательностями данных Ря и РN-1 может быть записана в виде

м

(в) = ^2 Рн(з)Рм-1 (з + в),

3 = 1

где М - число точек в последовательности; —К ^ 5 ^ К и Ря-1 (з + в) =0, когда (з + в) < 1 или (3 + в) > М. Пусть исходные данные томографического исследования представляют собой набор из q проекционных изображений размером I х I пикселей. Таким образом, будем иметь набор матриц А^\ к = 1,д; г, 3 = 1,/, элементами которых являются значения плотности распределения препарата в точках с координатами (хъ ,уг). Для анализа проекционных кадров в качестве наборов данных используются суммарные профили, полученные из планарных изображений для каждого из углов наблюдения:

I

Сгк = ^2л(£)^=~1’к=Т^’ (!)

3=1

^ = Еа(?^' = м,^ = м. (2)

ъ=1

Формулы (1), (2) представляют суммарные профили вдоль осей х и у соответственно.

Рассмотрим корреляцию между заданным и предшествующим ему кадрами. Параметр К обычно выбирается равным 20 из предположения, что два последовательных кадра будут коррелированны в пределах 20 пикселей в каждом направлении. Таким образом, функции взаимной корреляции для двух последовательных кадров с номерами к и (к — 1) будут иметь вид

I

йк)(з) = '52сзксз+8,к-1, к = 1^, -К^в^К, Ке Z>0, (3)

3=1

где Оз+3,к-1 = 0, если з + в > I или З + в < 1,

I

= к=Т^Ч, -К ^8 ^К, К ех>0, (4)

1=1

где Въ+3,к-1 = 0, если г + в > I или г + в < 1. Формулы (3), (4) представляют собой вид функции взаимной корреляции относительно профилей х и у соответственно.

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

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

Реализация [11] решения задачи коррекции движения при ОФЭКТ была произведена в виде отдельного модуля (клинической программы), совместимой с программным комплексом «Диагностика» [12]. В качестве входных данных использовались отдельные файлы в формате DICOM или осуществлялась связь с базой данных томографических исследований. Корректировалось движение более 1 пикселя в поперечном и более 0.5 пикселя в продольном направлениях. Алгоритм подтвердил теоретические рассуждения, показав наличие движения для некоторых групп кадров, и смог его удачно скорректировать. Кроме того, для некоторых исследований алгоритм определил небольшой сдвиг между кадрами в ^-направлении, обусловленный неправильным объединением данных с двух детекторов. Основываясь на полученных с помощью программы результатах, можно сделать вывод, что приведенный алгоритм позволяет с высокой точностью определить наличие движения пациента в ходе ОФЭКТ обследования и произвести необходимый сдвиг проекционных данных.

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

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

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

В настоящее время широкое распространение получили методы оптического потока [13-15], которые широко применяются в телевизионных системах, системах видеонаблюдения, анализа оптических изображений поверхности материалов и т. д. В статье предлагается использовать понятие оптического потока при обработке радионуклидных исследований. В рассматриваемом случае под оптическим потоком будем понимать двумерное поле скоростей (поле векторов перемещения), описывающее наблюдаемое в изображении смещение точек, происходящее при движении изображаемых объектов относительно детектора гамма-камеры.

Будем рассматривать модель, предложенную в статье [16]:

x(k +1) = x(k) + u(k, x(k)), k = 0, 1, 2, ..., где x(k) - 2-мерный вектор пространственных координат;

u(k,x(k))=<ф)) = (uitlk)))

- 2-мерный вектор перемещения, который соответствует номеру кадра. Считаем заданными x(0) - координаты точек контура на начальном кадре. Надо найти координаты точек контура на следующем кадре или результирующий вектор перемещения u(x(k)). В этой связи рассмотрим распределение плотности (концентрации) РФП по кадрам, которое обозначим р = p(k, x(k)).

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

p(k +1,x(k + 1)) = p(k, x(k)), k = 0,1, 2,....

Далее рассмотрим один шаг процесса, т. е. два соседних кадра, и обозначим x(0) = x, u(0, x(0)) = u(x) = u, x(1) = x + u(x) = x + u, р(0, x(0)) = p(x), р(1, x(1)) = p(x + u(x)) = p(x + u).

В статье [16] для нахождения вектора перемещения используется метод вложенных итераций

uk+1 = uk + duk,

где для определения duk получена система, которая в случае постоянства плотности РФП вдоль траектории принимает вид

(рХс1)2duk — a2 Aduk + pXipx2dv‘k = 0/2 Auk + (р — рк )рхЪ (5)

(р12)2duk - a2Aduk + pkXipkX2duk = a2 Auk + (р - рк)р’ Здесь р = p(x), pk = p(x + uk), a2 - весовой коэффициент.

Заменяя в системе (5) производные соответствующими конечными разностями, получаем линейную систему, которую решаем итерационным блочным методом Гаусса-Зейделя [17].

На рис. 3 приведен пример построения оптического потока для проекционного кадра в случае движения различного рода: смещения всех точек изображения вдоль оси х или у, в случае реального движения органа во время исследования.

Г

I \Л 7Л

чч'.ч^п1-1! ’’■л"!.

■’-V

. ... ^■|А-іЧК\Ч,лЛ ■ ■ ■ ■

ь

■л.чЧкуЛ'л;,1;

—ЛЗй

-л1л\, ■. -ч V- ■

КВДгВДНВДВДЗД*'

\+М■ ч>,Л^Л Л1.

■л -ь-н ->-4. т. \ ^ \ уч % \ч 4% V

Vі 'л '-1 V 'Л V \

? ччл. чч\.г\\

\\кЛтлч-ч. ч,чч. к ч-уч^ч^-лл

Ч. \ % V-1-

*л-ччч*а\

1|Н^,-.НН’.Н 4.-% ■ Ч1 г ■ -і-иЧ %V -Ч ' ь

■. ■.,Ч Ч. Ч Ч| Ч \Ч. ■ ■*■ і Ч.-ЧИ.-Ц.-Ч Ч'-.ч-Ч*

"ь ^ 'і N "і Ч %

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

Для повышения чувствительности предложенного метода, а именно для возможности правильно определить движения большой (более 2 пикселей) амплитуды, будем использовать пирамиду изображений [18]. Она представляет собой набор изображений в уменьшающемся масштабе, организованный в виде пирамиды. Основу пирамиды составляет исходное изображение высокого разрешения; вершина пирамиды состоит из приближения низкого разрешения (рис. 4). По мере движения вверх по пирамиде масштабы (размер и разрешение) уменьшаются. Если нижний уровень Ь имеет размеры 2Ь х 2Ь или N х N, где N = 1с^2 Ь, то размеры промежуточного уровня I будут

21 х 21, где 0 ^ I ^ Ь.

Рис. 4- Структура пирамиды изображений

Рассмотрим алгоритм получения различных уровней пирамиды изображений: для построения следующего уровня пирамиды, имеющего меньший масштаб,

необходимо произвести фильтрацию исходного изображения и применить прореживающую выборку с фактором 2 (под такой выборкой будем понимать операцию, заключающуюся в отбрасывании каждого второго отсчета). Для этого могут быть использованы различные виды фильтрации: усреднение по окрестности (пирамида средних значений), гауссова низкочастотная фильтрация (гауссова пирамида) или фильтрация не проводится вовсе (пирамида прореженных значений). Для изучаемого случая проекционных изображений предлагается использовать пространственный усредняющий фильтр. Поскольку, как правило, в качестве проекционных кадров выступают изображения размером 64 х 64 или 128 х 128, а амплитуда движения между кадрами в большинстве случаев не превышает 5 пикселей, то на данном этапе ограничимся рассмотрением двух уровней пирамиды. С помощью описанного метода коррекции движения для уменьшенного изображения можно определить, существует ли необходимость применения пирамид изображений к данной группе кадров или же можно работать сразу с полномасштабными изображениями.

Рис. 5. Исходный кадр с выделенным контуром (а), последующий кадр, координаты контура не изменены (б), вид оптического потока для точек контура (в) и последующий кадр с откорректированным положением контура (г)

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

Заключение. Рассмотрена задача коррекции движения при различных радионуклидных исследованиях. Для исследований ОФЭКТ решена задача коррекции движения

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

Литература

1. Cooper J. A., Neumann P. H., McCandless D. K. Effect of patient motion on Thallium-201 SPECT / / J. of Nucl. Medicine. 1992. Vol. 33, N 8. P. 1566-1571.

2. Botvinivk E., Zhu Y., O’Connell W., Dae M. A quantitative assessment of patient motion and its effect on myocardial perfusion SPECT image SPECT // J. of Nucl. Medicine. 1993. Vol. 34, N 2. P. 303-310.

3. Eisner R., Churchwell A., Nowak D., Cloninger K. Quantitative analysis of the tomographic Thallium-201 myocardial bullseyedisplay: critical role of correcting for patient motion // J. of Nucl. Medicine. 1989. Vol. 29. P. 91-97.

4. Mommennezhad M., Zakavi S. R., Sadeghi R., Kakhki V. Review of the linogram and sinogram: an easy way to detect off-peak artifacts in myocardial perfusion SPECT // J. Nucl. Med. Technology. 2009. Vol. 37, N 3. P. 188-190.

5. Geckle W. J., Frank T. L., Links J. M. Correction for patient and organ movement in SPECT: Application to exercise Thallium-201 cardiac imaging // J. of Nucl. Medicine. 1988. Vol. 29, N 4. P. 441450.

6. Cooper J. A., Neumann P. H., McCandless B. K. Detection of patient motion during tomographic myocardial perfusion imaging // J. of Nucl. Medicine. 1993. Vol. 34, N 8. P. 1341-1348.

7. Noever T., Nowak D., Carlson W. et al. Use of cross-correlation function to detect patient motion during SPECT imaging // J. of Nucl. Medicine. 1987. Vol. 28. P. 97-101.

8. Appledorn C. R., Oppenheimand B. E., Wellman H. An automated method for the alignment of image pairs // J. of Nucl. Medicine. 1980. Vol. 21. P. 165-167.

9. O’Connor M. K., Kanal K. M., Gebhard M. W., Rossman P. J. Comparison offour motion correction

techniques in SPECT imaging of the heart: A cardiac phantom study // J. of Nucl. Medicine. 1998. Vol. 39, N 12. P. 2027-2034.

10. Tabrizi M. H. N., Dong K., Movahed A. Performance comparison of different motion correction

algorithms used in cardiac SPECT imaging // ITCC. 2001. Vol. 1. P. 488.

11. Максимов К. М. Программная реализация алгоритмов коррекции движения при однофотонной эмиссионной томографии // Процессы управления и устойчивость: Труды 40-й Междунар. науч. конференции студентов и аспирантов / под ред. Н. В. Смирнова, Г. Ш. Тамасяна. СПб.: Издат. дом С.-Петерб. гос. ун-та, 2009. С. 347-352.

12. Котина Е. Д. Программный комплекс «Диагностика» для обработки радионуклидных исследований // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2010. Вып. 2. С. 100-113.

13. Anandan P. A computational framework and an algorithm for the measurement of visual motion // Intern. Journal of Computer Vision. 1989. Vol. 2. P. 283-310.

14. Horn B. K. P., Schunck B. G. Determining of Optical Flow // Artificial Intelligence. 1981. Vol. 17. P. 185-203.

15. Ibrahim N., Riyadi S., Zakaria N. et al. Implementation of differential optical flow algorithms in natural rigid video motion // Proc. of IMECS. 2009. Vol. 1. P. 795-798.

16. Котина Е. Д. К теории определения поля вектора перемещения на основе уравнения переноса для дискретного случая // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2010. Вып. 3. C. 38-43.

17. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем / пер. с англ. Х. Д. Икрамова, И. Е. Капорина; под ред. Х. Д. Икрамова. М.: Мир, 1991. 367 с.

18. Гонсалес Р., Вудс Р. Цифровая обработка изображений / пер. с англ.; под ред. П. А. Чочиа. М.: Техносфера, 2005. 1072 с.

Статья рекомендована к печати проф. Д. А. Овсянниковым.

Статья принята к печати 14 октября 2010 г.