CC BY
24
Серия «Математика» 2013. Т. 6, № 3. С. 48—59
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 519.6
Определение поля скоростей в задачах обработки изображений *
Е. Д. Котина
Санкт-Петербургский государственный университет
Г. А. Пасечная
Санкт-Петербургский государственный университет
Аннотация. В работе рассматриваются два подхода к определению поля скоростей в задачах обработки радионуклидных изображений. В обоих случаях задача сведена к решению линейных систем специального вида блочными итерационными методами, показана сходимость данных методов к единственному решению систем.
Ключевые слова: оптический поток; поле скоростей; блочные итерационные методы; сходимость блочных итерационных методов, радионуклидные изображения.
В статье рассматриваются математические модели для определения поля скоростей при обработке радионуклидных изображений. Ради-онуклидная диагностика — это активно развивающаяся область медицинской функциональной диагностики заболеваний человека с использованием радиоактивных изотопов [3, 8]. Радионуклидные исследования проводятся с помощью гамма-камер или гамма-томографов [1]. Динамический режим сбора данных позволяет наблюдать распределение радиофармпрепарата (РФП) в исследуемой системе организма в зависимости от времени. Таким образом, в результате исследования можно получить плотность распределения РФП, зависящую от времени и пространственных координат р = р(Ь,х,у),Ь € [0,Т], (х,у) € Б или с учетом дискретизации по времени и пространственным координатам имеем последовательность матриц р1 (г, у), р2(г,з), рт (г,3), г= 0,... ,п + 1.
1. Введение
* Работа выполнена при финансовой поддержке СПбГУ: тема № 9.39.1065.2012, тема № 9.38.673.2013.
Задача построения поля скоростей или определения оптического потока рассматривалась в работах многих авторов [10-12, 16, 17, 21, 22]. Оптический поток между парой изображений есть векторное поле, задающее естественную трансформацию первого изображения во второе. В данном случае под оптическим потоком будем понимать двумерное поле скоростей, описывающее наблюдаемое в изображении смещение точек, происходящее при движении изображаемых объектов относительно детектора гамма-камеры.
В работе рассматриваются два подхода к определению поля скоростей, а именно подход, основанный на предположении о постоянстве плотности распределения РФП вдоль траекторий движения и на предположении о постоянстве ее градиента. Задача нахождения поля скоростей сводится к решению линейной системы специального вида. Линейные системы решаются блочными итерационными методами, показывается сходимость этих методов к единственному решению системы.
2. Постановка задачи
Рассмотрим систему дифференциальных уравнений:
x = u(t, x,y), ,
. (/ ) (2.1) y = v(t,x,y).
Предполагается, что транспорт РФП (индикатора) описывается системой 2.1. Здесь x и y пространственные координаты, t - время, функции u и v определяют поле скоростей.
Рассмотрим два случая. В первом будем предполагать, что плотность распределения РФП вдоль траекторий системы 2.1 остается постоянной:
Pt + pxU + pyV = 0. (2.2)
Здесь pt, px, py — обозначения для частных производных первого порядка по t, x,y соответственно.
Во втором случае, что градиент плотности распределения вдоль траекторий системы 2.1 остается постоянным:
I pxx U + pxyV = -pxt, ,0 ол
\ pyxU + pyyV = -pyt, ( . )
где pxx, pxt, pxy, pyy, pyt - обозначения для частных производных второго порядка, (u,v)T = f - поле скоростей системы 2.1.
Исходя из описанного выше модельного представления транспорта индикатора, рассмотрим обратную задачу определения поля скоростей
в системе 2.1, т. е. нахождение функций и,у по известной функции р(Ь,х,у). В общем случае - это некорректная задача. В связи с этим будем использовать метод регуляризации, описанный А. Н. Тихоновым [9]. Рассмотрим задачу построения поля скоростей в области М в некоторый момент времени Ь. Составим интегральный функционал:
.] (и, у) = ! (ч>2 + а2 ф2^хйу, (2.4)
м
где
V2 = р + рхи + руь)2в случае предположения 2.2, V2 = (рхл+рххи+рхуУ)2 + (руг+рухи+рууу)2в случае предположения 2.3,
ф2 = иХ + и^ + ьХ + у^,
а2 - параметр регуляризации, М — область ненулевой меры из Я2.
Таким образом, будем рассматривать задачу нахождения поля скоростей как задачу минимизации функционала 2.4.
3. Уравнения Эйлера-Лагранжа
Выпишем для интегрального функционала 2.4 уравнения Эйлера -Лагранжа. В первом случае, когда рассматривается уравнение 2.2, они будут иметь вид:
-а2Аи + рХи + рхру у = -ргрх,
2 . 2 (3.1) -а Ау + ру у + рхру и = -ргру.
Во втором случае при рассмотрении системы 2.3 имеем:
-(Х2Аи + (р2хх + р2ух)и + рху(рхх + руу)у = РххРхЬ - рyxрyt,
2 2 2 (3.2)
-а Ау + (рху + руу)и + рху(рхх + руу)и = -рхурxt - рууру^'
Здесь А — оператор Лапласа, А и = + и А у = + тг?.
" ^ ^ ' ах2 ау2 ах2 ау2
В результате данного подхода задача нахождения поля скоростей системы 2.1 сводится к решению системы 3.1 или 3.2 при соответствующих граничных условиях.
4. Разреженные системы специального вида
Будем рассматривать системы 3.1, 3.2, где (х,у) € М и функции и и V заданы на границе области М, которую обозначим - О. Правая часть представляет собой известную функцию от х, у .
Учитывая дискретный характер измерений, как уже говорилось во введении, обозначаем плотность распределения РФП в точке, находящейся на пересечении г-й строки, j-го столбца и к-го момента времени за рк (г, у), г, у = 1,...,п. Под единичным изменением расстояния будем понимать приращение вдоль какой-либо оси величиной в один пиксель изображения. Будем искать решения систем 3.1, 3.2 в узлах квадратной сетки с шагом в один пиксель.
Если обозначить через и(г,.])^(г,у) приближение к решению задач 3.1, 3.2 в точке (г,]) , принадлежащей сетке, то, заменяя оператор Лапласа конечными разностями, получим для системы 3.1 систему линейных уравнений:
-а2(и(г - 1,э) + и(г + 1,у) + и(г- 1) + и(г+ 1))+
+(4а2 + рХ(г,з))и(г,з) + Рх(г,3)Ру(г,3Мьз) = —рх(г,3)рь(г,3), .лл.
о (4.1)
-а2^(г - 1,у)+ v(г + 1,у)+ v(г,j - 1) + v(г,j + 1))+
^ +(4а2 + ру(г,j))v(г,j) + рх(г. )ру(г,j)u(г,j) = -ру(г,j )рь(г,j).
г. = 1,...,п.
Для системы 3.2, соответственно, получаем:
' -а2(и(г - 1^)+ и(г + 1^) + и(г^ - 1) + и(г^ + 1)) +
+ (4а2 + р2хх(г,Л + рух^^ ЪФ. ) + (рхх(Ь. ) + руу(г,j))Х Хрху(г^^(г^) = -рхх(г^)рхг(г^) - рух(г^)руЛг^),
-а2^(г - 1^) + v(г + 1^) + v(г,j - 1) + v(г,j + 1)) +
+(4а2 + р2ху(г,j) + ру^^. + (рхх(Ь. ) + руу(г,j ))х
Хрху (г,j )и(г,Л = -рху (г,j)рxt(г,j) - руу (г,j )рyt(г,j)■
г^ = 1,...,п.
Принимая во внимание введенную дискретизацию, можем лить приближенные значения частных производных р^ рх, ру, ных производных второго порядка рхх, р^, рху, руу, р^, используя значения плотности в соседних точках сетки по выбранной схеме.
Согласно нашему предположению о том, что функции и и v заданы на границе области, только 2п2 переменных u(г,j), v(г,j), г,j = 1,...,п, во внутренних точках сетки являются неизвестными в уравнениях 4.1,
(4.2)
вычис-и част-
или, соответственно, 4.2. Таким образом, получается линейная система из 2п2 уравнений, решение которой дает приближение к решению системы 3.1(или 3.2) в точках сетки.
В результате дискретизации систем дифференциальных уравнений в частных производных 3.1 и 3.2 получили линейные системы разностных уравнений 4.1 и 4.2. Рассмотрим далее их решение.
Системы 4.1, 4.2 можно записать в виде:
A B B C
(4.3)
где
u = (ui, . . .,un2 ) = (uii, . . .,uin ,U21, . . .,U2n, ■ ■ -,Unl, . . . ,Unn)
V = (vi, . . . ,Vn2 )T = (vii, . . . ,Vin,V2i,. . . ,V2n, . . . ,Vni, . . . ,Vnn) значения искомых функций;
т
т
d = ^i^.., dn^)T = (da,..., dln, d2i,..., d2^ ..., dni,..., dnn)
T
е — (e1,■■ ■,еп2) — (e11, ■ ■ ■ , е1п1 e21,■ ■ ■ , e2n,■ ■ ■ , en1,■ ■ ■ , епп) ■
Запишем подробнее для системы 4.1:
díj = -px(i,j)pt(i,j),i,j = 1 ,п.
ehi = ~py(hj)pt(hj),i,j = T~ñ.
Матрица В — диагональная, т.е. Ьзг=0, в,т=1,п2, Ьзз = рх(г)ру(г, у)■ Матрицы А, С - разреженные и отличаются диагональными элементами, то есть азг = езг и ненулевые из них азг = езг = -а2■ Диагональные элементы матрицы А: а33 = 4а2 + рх(г^), матрицы : с33 = 4а2 + =Т~п,,в = 1 ,п2. Теперь запишем подробно для системы 4.2:
= -Рхх{г,з)рхг{г,з) - РУх(г, .])руг(г,.]), г,] = 1 ,п. = -рхУ(г,Лрхл(г,Л ~ РУУ(г, з)рУЛг, 3 = 1 ,п. Матрица В — диагональная, т.е.
Ьзг = 0,з,г = 1,п2,Ьзз = рху (г,3 )(рхх (г,з) + руу(г,3)) ■
Матрицы А, С - разреженные и отличаются диагональными элементами, то есть азг = сзг и ненулевые из них азг = сзг = -а2■ Диагональные элементы матрицы А: азз = 4а2 + рхх(г,3) + рух(г,3), матрицы : сзз =
4 а2 + р2уу{1,]) + р1х(ьз),ЬЗ =Т~п,,в = 1 ,п2.
Таким образом, матрица системы 4.3 имеет вид: рис. 1.
Пусть
■•••V
■- ••■л,
Рис. 1. Вид матрицы системы 4.3 для сетки 6х6.
г = (гг,...,гп2 )т ,гг =
Vl
п
—
й1
д = (ц1,...,Чп2) ,91 = ^ ¡,...,Яп2 = Тогда система 4.3 примет вид:
Нг = д,
где компоненты матрицы Н:
(1п2
(4.4)
Ня
, Нвт
,в = Г, в, Г = 1, п2.
Представим матрицу Н следующим образом:
Н = Б - Е - Г, где Б - матрица, состоящая из диагональных блоков:
Б
(Ни 0 0 Н22
0 0
Е + Г = -
V 0 ■■■ 0 Нп2п2 /
( 0 Н12 ••• Н21 0 •••
\Нп21 ■■■ Н,
п2 ,п2 — 1
Н1п2 Н2п2
0 У
(4.5)
Будем решать систему 4.4 с матрицей вида 4.5 блочными итерационными методами Гаусса - Зейделя и последовательной верхней релаксации.
20 30 40 Э0 00 70
ип2
^2
е п2
Е. Д. КОТИНА, Г. А. ПАСЕЧНАЯ 5. Блочные методы. Сходимость
Блочный метод Гаусса-Зейделя решения системы 4.4, отвечающий введенному разбиению, имеет вид:
n2 n2
Hsszks+1 = -Y, HsrZkr+1 Hsrzk + qs, (5.1)
r<s r>s
s = l,...,n2, k = 0,l, 2,...,
или с учетом 4.5, используя матричные обозначения
(D - E)zs+1 = Fzs + q, k = 0, l, 2,.... (5.2)
Для его сходимости в случае, когда матрица B — диагональная, матрицы A и C отличаются диагональными элементами, достаточно, чтобы выполнялись условия, полученные в статье [4]: 1) asscss - b2ss > 0, ass > 0, css > 0, s = 1, п2,
2 Г ~
Ю £ ||iw|| + \ ( _££___££ I + b2ss,s = l,n2 (5.3)
r=1,r=s
и для некоторого s
т=1,т=8 '
3) условие «блочной» непроводимости.
Блочный итерационный метод последовательной верхней релаксации для системы 4.4 имеет вид:
п2 п2
Н^1 = Нвг гкг+1 - Е Н33гк + д8) + (1 - ш)Н33гк3, (5.4)
Г<8 Г>8
в = 1,...,п2, к = 0,1,2,.... В матричном виде решение можно представить следующим образом:
Бгк+1 = ш(Егк+1 + Ггк + д) + (1 - ш)Бгк, к = 0,1,2,.... (5.5)
В случае, когда матрица В — диагональная, матрицы А и С симметричны и отличаются диагональными элементами, и выполняются перечисленные условия 5.3, метод последовательной верхней релаксации сходится при любом ш € (0, 2) и при любом начальном приближении.
Нетрудно проверить, что для рассматриваемых систем условия 5.3 выполняются, и таким образом метод Гаусса - Зейделя и метод последовательной верхней релаксации сходятся к единственному решению при любом начальном приближении.
6. Построение поля скоростей для радионуклидных
изображений
Рассмотрим примеры с использованием радионуклидных изображений, полученных в результате исследований проводимых на двухде-текторном гамма-томографе «ЭФАТОМ» [1] в «Федеральном научно-клиническом центре специализированных видов медицинской помощи и медицинских технологий Федерального медико-биологического агентства» (ФНКЦ ФМБА России, Москва). На данном томографе обработка данных производится с помощью программного комплекса «Диагностика» [5, 6].
Важными этапом при обработке данных радионуклидных исследований является обнаружение и коррекция движения пациента во время исследования, поскольку даже небольшое смещение пациента или исследуемого органа во время сбора проекционных данных может повлиять на достоверность результатов диагностики [13, 15, 18, 22], а полностью избежать изменения положения пациента или его отдельных внутренних органов во время исследования практически невозможно.
В результате проведенных опытов с использованием двух построенных алгоритмов получены значения поля скоростей в узлах сетки. В некоторых случаях метод, основанный на предположении о постоянстве градиента плотности распределения, дает лучшие результаты, в частности, в случае изменения общего уровня яркости изображения.
Рассмотрим два изображения, полученных при исследовании гепа-тобилиарной системы человека: рк(грк+1(г^),г,j = 0,...,п + 1, на которых имеет место сдвиг исследуемого органа (рис. 2).
В области интереса, выделенной пунктиром на рис. 2 а), поле скоростей имеет вид: рис. 3.
На данном примере мы видим, что предложенный метод может использоваться для обнаружения движения исследуемого органа и его последующей коррекции [2, 7].
Рассмотрим теперь возможность использования метода при построении контуров исследуемых объектов. Проиллюстрируем это на примере кардиологических изображений (рис. 4).
Рис. 2. а) Изображение рк(г,]) в момент времени к, г,] = 0,...,п + 1, Ь) Изображение рк+г(г,з) в момент времени к + 1,г,] =0,...,п + 1.
Рис. 3. Поле скоростей
На рис. 4 а) контуром выделен левый желудочек сердца. Поле скоростей в точках данного контура представлено на рис. 5.
7. Заключение
Приведенные в статье методы могут быть использованы для обработки изображений, полученных при радионуклидных исследованиях для обнаружения движения и его коррекции, а также при построении контуров исследуемых объектов. Следует отметить, что эти методы могут найти применение и при обработке других данных. В частности, при построении поля скоростей в задачах анализа динамики заряженных частиц [19, 20].
Рис. 4. а) Изображение сердца рк(%, ]) в момент времени к, %,] = 0,...,п + 1, Ь) Изображение сердца рк+1(г,]) в момент времени к + 1,%,] = 0,...,п + 1.
Рис. 5. Поле скоростей
Список литературы
1. Двухдетекторный однофотонный эмиссионный гамма-томограф ЭФАТОМ / М. А. Арлычев, В. Л. Новиков, А. В. Сидоров, А. М. Фиалковский, Е. Д. Котина, Д. А.Овсянников, В. А. Плоских // Журн. техн. физики. - 2009. - Т. 79, вып. 10. - С. 138.
2. Котина Е. Д. К теории определения поля вектора перемещения на основе уравнения переноса для дискретного случая / Е. Д. Котина // Вестн. С.-Петерб. ун-та. Сер. 10, Прикл. математика, информатика, процессы управления. -2010. - № 3. - С. 38-43.
3. Котина Е. Д. Обработка данных радионуклидных исследований / Е. Д. Котина // Вопр. атомной науки и техники. Сер. Ядер.-физ. исслед. - 2012. - № 3(79). - С. 195-198.
58
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15
16
17
18
19
20
21
Котина Е. Д. О сходимости блочных итерационных методов / Е. Д. Котина
// Изв. Иркут. гос. ун-та. - 2012. - Т. 5, вып. 3. - С. 41-55.
Котина Е. Д. Программный комплекс «Диагностика» для обработки радио-
нуклидных исследований / Е. Д. Котина // Вестн. С.-Петерб. ун-та. Сер. 10,
Прикл. математика, информатика, процессы управления. - 2010. - № 2. - С.
100-113.
Котина Е. Д. Программный комплекс обработки радионуклидных исследований / Е. Д. Котина // Вестн. С.-Петерб. гос. ун-та технологии и дизайна. Сер. 1, Естеств. и техн. науки. - 2010. - № 1. - С. 43-51.
Котина Е. Д. Коррекция движения при томографических и планарных ра-дионуклидных исследованиях / Е. Д. Котина, К. М. Максимов // Вестн. С.-Петерб. Ун-та. Сер. 10, Прикл. математика, информатика, процессы управления. - 2011. - С. 29-36.
Радионуклидные методы в кардиологической клинике / Е. Н. Остроумов, Е. Д. Котина, О. Р. Сенченко, А. Б. Миронков // Сердце : журн. для практ. врачей. - 2010. - Т. 9, вып. 3. - С. 190-195.
Тихонов А. Н. Методы решения некорректных задач / А. Н. Тихонов, В. Я. Арсенин. - М. : Наука, 1974. - 285 с.
Anandan P. A computational framework and an algorithm for the measurement of visual motion / P. Anandan // International Journal of Computer Vision. - 1989. - Vol. 2. - P. 283-310.
Barron J. Performance of optical flow techniques / J. Barron, D. Fleet // International Journal of Computer Vision. - 1994. - Vol. 12. - P. 43-77. Bergholm F. A theory of optical flow / F. Bergholm, S. Carlsson // Computer vision. Graphics and Image Processing: Image Understanding. - 1991. - Vol. 53, N 2. - P. 171-188.
Cooper J. A. Detection of patient motion during tomographic myocardial perfusion imaging / J. A. Cooper, P. H. Neumann, B. K. McCandless // Journal of Nuclear Medicine. - 1993. - Vol. 34. - P. 1341-1348.
Fleet D. J. Optical flow estimation / D. J. Fleet, Y. Weiss // Mathematical models for Computer Vision: The Handbook. - Springer, 2005. - 24 p. Detection and correction of patient motion in dynamic and static myocardial SPECT using a multi-detector camera / G. Germano [et al.] // Journal of Nuclear Medicine. - 1993. - Vol. 34. - P. 1394-1395.
Horn B. K. P. Determining optical flow / B. K. P. Horn, B.G. Schunck // Artificial intelligence. - 1981. - Vol. 17. - P. 185-203.
Kotina E. D. Data Processing and Quantitation in Nuclear Medicine / E. D. Kotina, V. A. Ploskikh // Proceedings of RuPAC 2012, http://accelconf.web.cern.ch/AccelConf/rupac2012/. - 2012. - P. 526-528. Quantitative assessment of motion artifacts and validation of a new motion-correction program for myocardial perfusion SPECT / N. Matsumoto [et al.] // Journal of Nuclear Medicine. - 2001. - Vol. 42. - P. 687-694.
Ovsyannikov D. A. Determination of velocity field by given density distribution of charged particles / D. A. Ovsyannikov, E. D. Kotina // Problems of Atomic Science and Technology. - 2012. - Vol. 3(79). - P. 122-125.
Ovsyannikov D. A. Reconstruction of velocity field / D.A. Ovsyannikov, E. D. Kotina // Proceedings of ICAP2012, Rostock-Warnemunde, Germany. - 2012. -P. 256—258. - URL: http://accelconf.web.cern.ch/AccelConf/ICAP2012. Highly Accurate Optic Flow Computation with Theoretically Justified Warping / N. Papenberg [et al.] // International Journal of Computer Vision. - 2006. - Vol. 67, N 2. - P. 141-158.
22. Effect of motion on Thallium-201 SPECT / M. F. Prigent, M. Hyun, D. S. Berman, A. Rozanski // Journal of Nuclear Medicine. - 1993. - Vol. 34. - P. 1845-1850.
E. D. Kotina, G. A. Pasechnaya Determining of velocity field for image processing problems
Abstract. In this work two methods on determining of velocity field for radionuclide image processing problems are considered. In both cases the problem add up to solving linear systems of special sort by block iterative methods, the convergence of this methods is investigated.
Keywords: optical flow; velocity field; block iterative methods; convergence of block iterative methods; radionuclide images.
Котина Елена Дмитриевна, доктор физико-математических наук, профессор, Санкт-Петербургский государственный университет, 198504, Санкт-Петербург, Университетский пр., 35, тел.: (812)4284868 (ekotina123@mail.ru)
Пасечная Галина Александровна, студентка, Санкт-Петербургский государственный университет, 198504, Санкт-Петербург, Университетский пр., 35, тел.: (812)4284868 (pasechnayagalina@gmail.com)
Kotina Elena, Saint-Petersburg State University, 35, Universitetskij pr., Saint-Petersburg, 198504, professor, Phone: (812)4284868 (ekotina123@mail.ru)
Pasechnaya Galina, Saint-Petersburg State University, 35, Universitetskij pr., Saint-Petersburg, 198504, student, Phone: (812)4284868 (pasechnayagalina@gmail.com)
