Санкт-Петербургский государственный университет
Факультет прикладной математики ‒ процессов управления
Кафедра теории управления
Пасечная Галина Александровна
Выпускная квалификационная работа
Построение поля скоростей для обработки данных радионуклидных исследований
Специальность 05.13.01 ‒ Системный анализ, управление и
обработка информации
(по прикладной математике и процессам управления)
Заведующий кафедрой,
доктор физ.-мат. наук,
профессор
Жабко А. П.
Научный руководитель,
доктор физ.-мат. наук,
профессор
Котина Е. Д.
Рецензент,
Гордеев Д.Ф.
кандидат физ.-мат. наук
г. Санкт-Петербург
2016
2
Оглавление
Введение……………………………………………………………………………3
Постановка задачи…………………………………………………………………5
Обзор литературы………………………………………………………………….7
Глава 1. Определение поля скоростей в задачах обработки радионуклидных
изображений………………………………………………………………………12
§1.1. Уравнения Эйлера-Лагранжа………………………………………………12
§1.2. Разреженные системы специального вида………………………………..12
§1.3. Блочные методы. Сходимость……………………………………………..16
Глава 2. Практическое применение методов определения поля скоростей для
обработки радионуклидных исследований……………………………………..19
§2.1. Построение поля скоростей для радионуклидных изображений………..19
§2.2. Построение поля скоростей для последовательностей радионуклидных
изображений……………………………………………………………………….21
Заключение………………………………………………………………………...25
Список литературы………………………………………………………………..27
3
Введение
Радионуклидная диагностика (ядерная медицина) — современный метод
лучевой диагностики для оценки функционального состояния различных органов
и систем организма с помощью меченных радионуклидами веществ —
радиофармпрепаратов.
Наиболее широкое распространение получил метод сцинтиграфии — метод
функциональной визуализации, заключающийся во введении в организм
радиоактивных изотопов и получении изображений путём детектирования
испускаемого ими излучения.
Существует несколько режимов сбора данных радионуклидного
исследования в зависимости от его целей: планарное статическое или
динамическое сканирование, сцинтиграфия всего тела, томографическое
сканирование, «Синхронизация» и «Томография с синхронизацией». Они
отличаются положением пациента относительно детектора гамма-камеры,
количеством формируемых изображений, способностью наблюдать распределение
индикатора в организме в зависимости от времени и получать картины объемного
распределения радиофармпрепарата.
Важную роль в процессе радионуклидной диагностики играет аппаратное
средство, с помощью которого она производится, выбор подходящего в
конкретном случае радиофармпрепарата и, конечно же,
обработка данных с
использованием полученных изображений.
В данной работе рассматривается построение поля скоростей для обработки
данных радионуклидных исследований. Предложенные методы основаны на
понятии оптического потока, который представляет собой распределение видимых
скоростей движения объектов, получаемое на основе их изображений в разные
моменты времени. Здесь под оптическим потоком будем понимать двумерное или
трехмерное поле векторов перемещения, описывающее наблюдаемое в
4
изображении смещение точек, происходящее при движении изображаемых
объектов относительно детектора гамма-камеры.
В первой главе рассматриваются методы определения поля скоростей для
двумерных изображений и их последовательностей в предположении о
постоянстве плотности распределения радиофармпрепарата вдоль траекторий
движения и о постоянстве ее градиента. Для решения используется метод
регуляризации по А.Н. Тихонову и исследуется вариационная задача.
Составляется интегральный функционал и рассматривается задача его
минимизации. Выписываются уравнения Эйлера-Лагранжа, которые
представляют собой дифференциальные уравнения в частных производных
второго порядка с заданными граничными условиями. Полученная система
сводится заменой частных производных конечными разностями к системе
линейных уравнений, выполняются некоторые преобразования, в результате
которых получается линейная система специального вида с разреженными
матрицами. Далее она решается известными итерационными методами. Матрицы
системы преобразуются в блочные, с блоками второго порядка. Полученная
система решается блочными итерационными методами Гаусса-Зейделя и
последовательной верхней релаксации (SOR), после предварительного
доказательства их сходимости к единственному решению системы.
Во второй главе представлены результаты экспериментов, проведенных с
использованием радионуклидных исследований, в ходе которых были
апробированы разработанные алгоритмы. Предложенные в работе методы
используются для определения движения областей интереса на изображениях, а
также для их оконтуривания.
5
Постановка задачи
Целью данной работы является построение двумерного поля скоростей в
задачах цифровой обработки радионуклидных изображений. Будем рассматривать
алгоритм для анализа перемещений объектов на основе информации о плотности
распределения радиофармпрепарата в разные моменты времени. Динамический
режим сбора данных позволяет наблюдать распределение показателя (РФП) при
изучении организма в зависимости от времени и пространственных координат
(Рис. 1).
Обозначим функцию плотности распределения
ρ=ρ( t , x , y ),
t ∈[0 ,T ],
( x , y )∈M ,
M ∈R 2 .
Рис.1. Распределение радиофармпрепарата в гепатобилиарной системе
Рассмотрим систему дифференциальных уравнений:
Ẋ=U (t , X ),
T
г д е X=( x , y )
— вектор пространственных координат, t
(1)
— время,
U =(u , v )T — определяет искомое поле скоростей.
Далее будем использовать два предположения: о том, что плотность
распределения РФП вдоль траекторий системы (1) остается постоянной и что
6
градиент плотности распределения вдоль траекторий системы (1) остается
постоянным.
Рассмотрим первое из них:
ρt + ρ x u+ ρ y v=0 .
Здесь ρt , ρ x , ρ y
— частные производные по t ,
x
(2)
, y соответственно.
Запишем второе:
{
ρ xx u+ ρ xy v=−ρ xt ,
ρ yx u+ ρ yy v=−ρ yt ,
(3)
г д е ρ xx , ρ xy , ρ xt , ρ yt , ρ yy
— обозначения для частных производных второго
порядка.
Такие модели позволяют рассматривать задачу определения поля скоростей
как некорректную задачу нахождения функций u,v
в системе (1) по известной
функции ρ(t , x , y ) . Наиболее распространенным методом решения такой задачи
является метод регуляризации по А. Н. Тихонову [24].
Согласно предложенному методу, введем интегральный функционал и
рассмотрим задачу построения поля скоростей в некоторый момент времени t .
J ( u , v )=∫ ( ϕ2 +α 2 ψ 2 )dxdy ,
M
(4)
где
ϕ2 =( ρt + ρ x u+ ρ y v )2 — в первом случае,
ϕ2 =( ρ xt + ρ xx u+ ρ xy v )2 +( ρ yt + ρ yx u+ ρ yy v )2 — во втором,
ψ 2=u2x +u2y + v 2x +v 2y ,
α 2 — параметр регуляризации, M — область ненулевой меры из R 2 .
Таким образом, для решения поставленной задачи о нахождении двумерного
поля скоростей необходимо минимизировать полученный функционал (4).
7
Обзор литературы
Радионуклидная диагностика — активно развивающаяся область
современной диагностической и функциональной медицины [6, 13, 15, 21, 22, 37].
Исследованиям по данной теме посвящено множество трудов, как в России, так и
за рубежом. Наряду с другими методами диагностики, ядерная медицина
получила широкое распространение и применяется для обнаружения заболеваний
различных органов, для наблюдения за их развитием, а также для уточнения
поставленных диагнозов. Особенно же важна функциональность методов
радионуклидной диагностики. Полученные изображения способны отображать
различные изменения, происходящие в организме человека, и динамику
п р о ц е с с о в , ч то д о с т и г а е т с я з а сч е т в в ед е н и я с о от в е т с т ву ю щ и х
радиофармпрепаратов, способных определенным образом накапливаться в
морфологических структурах. Возможность диагностировать функциональные
дефекты работы органов на самых ранних стадиях позволяет эффективно
вылечивать болезни, экономя при этом средства, требуемые для лечения.
Основными компонентами радионуклидного исследования являются:
прибор для регистрации излучения, соответствующий радиофармпрепарат и
обработка результатов с использованием математического моделирования и
компьютерных технологий.
В ядерной медицине используются такие основные аппаратные средства как
гамма-камера, гамма-томограф и позитронно-эмиссионный томограф [3].
Современная гамма-камера состоит из одного или нескольких блоков
детектирования, предназначенных для регистрации гамма излучения, штативноповоротного устройства, обеспечивающего крепление блоков детектирования,
ложа пациента, биосинхронизатора (например, кардиосинхронизатора),
8
электронной системы позиционирования, монитора укладки, компьютера сбора
данных и компьютера обработки изображения (Рис. 2). Если гамма-камера
снабжена более чем одним блоком детектирования и соответствующей системой
многоракурсового сбора информации, то такую гамма-камеру называют гамматомографом. Она позволяет с помощью компьютерной реконструкции по набору
проекций, сделанных под различными углами, осуществить получение
трехмерного распределения радионуклида в исследуемом органе.
Рис. 2. Функциональная блок-схема гамма-камеры (гамма-томографа)
О с н о в н ы м и ком п о н е н т а м и бл о ко в д е т е кт и р о ва н и я я вл я ют с я
многоканальный коллиматор, сцинтилляционный кристалл NaI(Tl) с большой
площадью поверхности, набор фотоумножителей и электронное устройство,
обеспечивающее определение координат и амплитуд сигналов.
Для получения радиоизотопного изображения необходимо применять
систему коллимации, которая способна выделять направление -квантов,
падающих на камеру. Наиболее распространены коллиматоры с параллельными
отверстиями, пропускающие только лучи, которые движутся перпендикулярно его
поверхности. Коллиматор определяет также геометрическое поле зрения камеры и
обусловливает пространственное разрешение и чувствительность всей системы.
Помимо коллиматоров с параллельными отверстиями существуют и коллиматоры
с единственным отверстием малого размера, предназначенные для визуализации
малых, приповерхностных органов, а так же коллиматоры со сходящимися или
9
расходящимися отверстиями для получения изображений всего тела и органов
средних размеров.
В сцинтилляторе при полном или частичном поглощении энергии
падающих на него -квантов возникают световые вспышки (сцинтилляции) очень
низкой интенсивности. Чтобы зарегистрировать такие вспышки, необходимо
специальное вакуумное устройство фотоэлектронный умножитель (ФЭУ).
Далее импульсы преобразуются, формируются координаты сигналов,
которые затем обрабатываются компьютером.
Процесс визуализации распределения гамма-излучения от РФП с помощью
гамма-камеры или гамма-томографа называется однофотонной эмиссионной
компьютерной томографией (ОФЭКТ).
Позитронно-эмиссионный томограф используется для регистрации
пространствнно-временного распределения позитронно-излучающего
радиофармпрепарата в теле пациента по аннигиляционному излучению. Важной
характеристикой любого ПЭТ сканера является тип материала, используемого для
сцинтиллятора. Наиболее часто применяемые в сцинтилляторах материалы:
детекторы на основе монокристаллов оксиортосиликата лютеция (Lu2SiO5, LSO).
В последнее время появляются приборы, соединяющие в себе гамматомограф с КТ и МРТ, позитронно-эмиссионный томограф с КТ. Таким образом,
объединяются лучшие стороны разных методов, используется преимущество
одних при проведении структурных обследований и других в получении
функциональных изображений.
Радионуклидные методы диагностики могут использоваться для
исследований с различными целями. В зависимости от этого выбирается
соответствующий режим сбора данных [5, 13]:
Планарное статическое сканирование. В процессе сбора информации
перемещения пациента и блоков детектирования не происходит.
Позволяет оценивать статическое распределение РФП в исследуемом
объекте.
10
Планарное динамическое сканирование. Данный режим аналогичен
предыдущему, но в процессе сбора для каждого детектора
формируется последовательность изображений с фиксированной
выдержкой. Динамический режим позволяет наблюдать кинетику
РФП в исследуемой системе организма.
Сканирование в режиме «все тело». Для сбора информации пациент
перемещается в горизонтальном направлении, но угловое положение
блоков не изменяется. Данный режим позволяет оценивать
статическое распределение индикатора во всём теле.
Планарное статическое сканирование с синхронизацией. Сбор
информации ведётся с использованием сигнала внешнего устройства
(кардиосинхронизатора). Такие исследования позволяют получать
изображения сердца в различных фазах сердечного цикла.
Томографическое сканирование. Формируется последовательность
изображений, соответствующих различным угловым ракурсам, и с
помощью томографической реконструкции получают объемное
распределение препарата. Перемещение пациента в процессе
сканирования не происходит.
Томографическое сканирование с синхронизацией. Перемещения
блоков осуществляются в томографическом режиме, но в каждой
позиции формируется серия изображений.
Полученные в результате изображения обрабатываются специальными
компьютерными комплексами программ, которые вычисляют различные
диагностические параметры, строят графики [8, 10, 39, 41].
Важными задачами в обработке информации, полученной в ходе
диагностики, являются обнаружение и коррекция движения пациента и органов во
время движения, а также построение контуров исследуемых объектов. Для их
решения используются различные методы, один из них – определение поля
скоростей. При таком подходе используется понятие оптического потока,
11
которому посвящено множество работ современных российских и зарубежных
авторов [25, 27, 28, 29, 32, 36, 43, 44, 45, 47, 49].
Основоположниками теории оптического потока считаются Horn B. K. P.,
Schunck B. G., которые в своей статье [36] дали определение этому термину и
предложили метод вычисления оптического потока. Представленный ими способ
основан на предположении о постоянстве яркости изображений и использует
алгоритм минимизации интегрального функционала с помощью уравнений
Эйлера-Лагранжа.
Другой базовый алгоритм определения оптического потока — алгоритм
Лукаса-Канаде, представленный в статье [45]. В данном случае используется то же
предположение, записывается основное уравнение оптического потока для всех
пикселей окрестности и полученная система уравнений решается методом
наименьших квадратов.
Для записи основного уравнения оптического потока могут использоваться
различные предположения о постоянстве некоторой величины, например о
постоянстве Лапласиана, Гессиана, нормы Гессиана и т.д. Большинство из них
описано в статье Papenberg N. [47].
Теория оптического потока активно развивалась вместе с развитием теле- и
видеосистем. Его методы используются при обработке различных изображений
[17]. Одной из важнейших сфер применения оптического потока является
радионуклидная диагностика. С появлением задач, связанных с обработкой
изображений, полученных в ходе радионуклидных исследований, методы
оптического потока стали применяться для обнаружения движения наблюдаемых
органов и его коррекции [34, 46].
Предпринимаются попытки построения трехмерного поля скоростей для
кардиологических исследований [12, 30, 35].
12
Глава 1. Определение поля скоростей в задачах обработки
радионуклидных изображений
§1.1. Уравнения Эйлера-Лагранжа
Выпишем для интегрального функционала (4) уравнения Эйлера-Лагранжа
в предположении (2):
−α 2 Δu+ ρ 2x u+ ρ x ρ y v=−ρ t ρ x ,
−α 2 Δv+ ρ 2y v+ ρ x ρ y u=−ρt ρ y .
(5)
Теперь с учетом предположения (3) для функционала (4) уравнения ЭйлераЛагранжа будут иметь вид:
−α 2 Δu+( ρ 2xx + ρ 2yx )u+ ρ xy ( ρ xx + ρ yy ) v=−ρ xx ρ xt −ρ yx ρ yt ,
−α 2 Δv+( ρ 2xy + ρ2yy ) v+ ρ xy ( ρ xx + ρ yy )u=−ρ xy ρ xt −ρ yy ρ yt .
(6)
∂2 u ∂2 u
∂2 v ∂2 v
Δu= 2 + 2
Δu= 2 + 2
∂
x
∂
y
∂x ∂ y .
Δ
Здесь
— оператор Лапласа,
и
В результате данного подхода, задача нахождения поля скоростей системы
(1) сводится к решению систем (5) и (6) при соответствующих граничных
условиях.
§1.2. Разреженные системы специального вида
Будем рассматривать системы (5) и (6), где ( x , y )∈M
и функции
u
и
v
заданы на границе области M , которую обозначим — Ω . Правые части
систем представляют собой известную функцию от x , y .
Учитывая дискретный характер измерений, обозначаем плотность
распределения РФП в точке, находящейся на пересечении i -ой строки, j -го
13
столбца и k -го момента времени за ρ k (i , j),
i , j=1,…,n . Под единичным
изменением расстояния будем понимать приращение вдоль какой-либо оси
величиной в один пиксель изображения. Будем искать решение систем в узлах
квадратной сетки с шагом в один пиксель (Рис. 3).
(i+1,j)
(i+1,j+1)
(i,j)
(i+1,j-1)
(i,j+1)
(i-1,j)
(i,j-1)
(i-1,j+1)
(i-1,j-1)
Рис. 3. Квадратная сетка для дискретизации
Если обозначить через u(i , j),
v(i, j)
приближение к решению задач в
точке (i , j) , принадлежащей сетке, то, заменяя оператор Лапласа конечными
разностями [14, 19], получим для системы (5) систему линейных уравнений:
{
−α 2 ( u(i−1 , j )+u(i +1 , j )+u(i , j−1)+u (i , j +1))+
+( 4 α 2 + ρ 2x (i , j ))u(i , j )+ ρ x (i , j ) ρ y (i , j )v (i , j )=−ρ x (i , j ) ρ t (i , j ),
−α 2 ( v (i−1 , j )+v (i +1 , j )+v (i , j−1 )+v (i , j+1))+
+( 4 α 2 + ρ 2y (i , j )) v( i , j )+ ρ x (i , j ) ρ y (i , j )u(i , j )=−ρ y (i , j ) ρ t (i , j ),
i , j=1,…, n.
А для системы (6)
(7)
14
{
−α 2 (u(i−1 , j )+u (i+1 , j )+u (i , j−1 )+u(i , j +1 ))+(4 α 2 + ρ2xx (i , j )+
ρ 2yx (i , j ))u(i , j )+ ρ xy (i , j )( ρ xx (i , j )+ ρ yy (i , j ))v (i , j )=
¿−ρ xx (i , j ) ρ xt ( i , j )− ρ yx (i , j ) ρ yt (i , j ),
−α 2 ( v (i−1 , j )+ v( i+1 , j )+v (i , j−1)+v (i , j +1))+( 4 α 2 + ρ 2xy ( i , j )+
+ ρ 2yy ( i , j ))v (i , j )+ ρ xy (i , j )( ρ xx (i , j )+ ρ yy (i , j ))u(i , j )=
¿− ρ xy (i , j ) ρ xt (i , j )− ρ yy ( i , j ) ρ yt ( i , j ),
i , j=1,…, n.
(8)
Принимая во внимание введенную дискретизацию, можем вычислить
приближенные значения частных производных ρ x ,
ρy ,
ρt
и частных
производных второго порядка ρ xx , ρ yy , ρ xt , ρ yt , ρ xy , используя значения
плотности в соседних точках сетки по выбранной схеме.
Согласно нашему предположению о том, что функции
2
границе области, только 2n
переменных u(i , j),
u
v(i, j) ,
и
v
заданы на
i , j=1,…,n , во
внутренних точках сетки являются неизвестными в уравнениях (7) и (8). Таким
2
образом, получается линейная система из 2n уравнений, решение которой дает
приближение к решению систем (5) и (6) в точках сетки.
В результате дискретизации систем дифференциальных уравнений в
частных производных (5) и (6) получили линейные системы разностных
уравнений (7) и (8), рассмотрим далее их решение.
Эти системы можно записать в виде
(
A
B
)( ) ( )
B u
d
=
C v
e
(9
,
)
где
T
T
u=(u1 ,…,un2 ) =( u11 ,…, u1n , u21 ,…,u2 n ,…un1 ,…,unn ) ,
T
T
v=( v1 ,…, v n2 ) =( v 11 ,…, v 1 n , v21 ,…, v2 n ,…v n 1 ,…, v nn )
—
значения искомых функций;
T
T
d =( d1 ,…, d n2 ) =( d11 ,…, d 1n , d 21 ,…, d 2n ,…dn 1 ,…, d nn ) ,
15
T
T
e=( e1 ,…,e n2 ) =(e11 ,…, e1 n ,e21 ,…, e2n ,…en1 ,…,enn ) .
Запишем подробнее для системы (7)
d i , j =−ρ x (i , j ) ρ t (i , j ),i , j=1,n.
ei , j =−ρ y (i , j) ρ t (i , j),i , j=1,n.
— д и а г о н а л ь н а я , т. е . b sr =0 ,
Матрица B
s≠r ,
s,r=1,n2 ,
b ss=ρ x (i , j) ρ y (i , j ).
Матрицы A ,C — разреженные и отличаются диагональными элементами, то
есть a sr =csr
2
и ненулевые из них a sr =c sr =−α .
2
2
a
=4
α
+
ρ
A
ss
x (i , j ),
Диагональные элементы матрицы
:
2
2
матрицы C : c ss =4 α + ρ y (i , j),
s=1,n2 .
i, j=1,n,
Для системы (8):
d i , j =−ρ xx (i , j) ρ xt (i , j)−ρ yx (i, j) ρ yt (i , j ),i , j=1,n.
ei , j =−ρ xy ( i, j) ρ xt (i, j)−ρ yy (i , j) ρ yt (i , j),i , j=1,n.
Матрица B
— д и а г о н а л ь н а я , т. е . b sr =0 ,
s≠r ,
s,r=1,n2 ,
b ss=ρ xy (i , j)( ρ xx (i , j )+ ρ yy (i , j)).
Матрицы A ,C — разреженные и отличаются диагональными элементами, то
есть a sr =csr
2
и ненулевые из них a sr =c sr =−α .
2
2
2
Диагональные элементы матрицы A : a ss=4 α + ρ xx (i , j )+ ρ yx (i , j),
2
2
2
c
=4
α
+
ρ
(i
,
j)+
ρ
C
ss
yy
yx (i , j ),
матрицы
:
Матрица системы (13) имеет вид: рис. 4.
i, j=1,n,
s=1,n2 .
16
0
10
20
30
40
50
60
70
0
01
02
03
40
05
60
07
384
=nz
Рис. 4. Вид матрицы системы (9) для сетки 6х6.
Пусть
( )
( )
( )
z =( z 1 ,…, z n 2 )T , z 1 =
u2
u1
,… , z n2 = n ,
v1
v n2
q=( q 1 ,… ,q n2 )T ,q1 =
( )
d 2
d1
,… ,q n2 = n .
e1
en2
Тогда система (9) примет вид [18, 48]:
Hz=q,
(10)
где компоненты матрицы H :
(
H ss =
ass
bss
)
bss
,
c ss
H sr =
(
a sr
b sr
Представим матрицу H
)
bsr
,
c sr
s≠r ,
s,r=1,n2 .
следующим образом:
H =D−E−F ,
где D — матрица, состоящая из диагональных блоков:
(11)
17
(
⋯
⋱
⋯
H 11
D= ⋮
0
(
E + F=−
0
⋮
H n2 n2
)
0
H 12
H 21
0
,
⋯
H 2 n2
⋮
⋮
⋱
H n2 n2 −1
⋯
H n2 1
H 1 n2
0
)
.
Будем решать систему (10) с матрицей H
вида (11) блочными
итерационными методами Гаусса-Зейделя и последовательной верхней
релаксации [1].
§1.3. Блочные методы. Сходимость.
Блочный метод Гаусса-Зейделя решения системы (10), отвечающий
введенному разбиению, имеет вид
H ss z ks +1=−
n2
∑
r <s
H sr z kr +1 −
n2
∑ H sr zkr +q s , s=1,…,n2 , k=0,1, …
r> s
(12)
или, с учетом (11), используя матричные обозначения
( D−E ) zk +1=Fz k +q , k=0,1,…
(13)
Для его сходимости в случае, когда матрица B — диагональная, матрицы
и C
A
отличаются диагональными элементами, достаточно, чтобы
выполнялись условия, полученные в статье [7]:
2
2
a
c
−b
>0
,
a
>0
,c
>0
,
s=1
,
n
,
ss
ss
ss
ss
ss
a)
ass +c ss
б)
2
2
n
≥
∑
r=1 ,r≠s
и для некоторого
|asr|+
s
√(
)
2
a ss −c ss
+b2ss , s=1 , n2 ,
2
(14)
18
ass +c ss
2
n2
>
∑
r=1 , r≠s
|asr|+
√(
)
2
a ss −c ss
+b2ss , s=1 , n2 ,
2
в) условие блочной неприводимости.
Блочный итерационный метод последовательной верхней релаксации для
системы (10) имеет вид:
(
n2
n
2
r <s
r> s
)
H ss z ks +1=ω −∑ H sr z kr +1 −∑ H sr z r +q s +(1−ω) H ss z s ,
k
k
(15)
s=1,…, n2 , k=0,1,…
В матричном виде решение можно представить следующим образом:
Dz k +1=ω( Ez k +1 + Fzk +q )+(1−ω) Dzk , k=0,1, …
В случае, когда матрица B
— диагональная, матрицы A
(16)
и C
симметричны и отличаются диагональными элементами, и выполняются
перечисленные условия (14),
метод последовательной верхней релаксации
сходится при любом ω∈( 0,2) и при любом начальном приближении.
Нетрудно проверить, что для рассматриваемых систем условия (14)
выполняются, и таким образом, метод Гаусса-Зейделя и метод последовательной
верхней релаксации сходятся к единственному решению при любом начальном
приближении.
19
Глава 2. Практическое применение методов определения поля
скоростей для обработки радионуклидных исследований
§2.1. Построение поля скоростей для радионуклидных изображений
Рассмотрим использование представленного метода для обнаружения
движения на радионуклидных изображениях. Изображения, полученные в ходе
радионуклидной диагностики, обладают определенной спецификой, являются
зашумленными и негладкими. Представленный в данной работе алгоритм
включает в себя предварительную обработку данных с помощью низкочастотных
пространственных фильтров с маской 3х3 [2, 20].
20
Будем обрабатывать два радионуклидных изображения, полученных в
разные моменты времени, на которых происходит перемещение и сжатие
исследуемого объекта, выделенного пунктиром.
Рис. 5. Радионуклидные изображения с выделенным объектом
исследования
В результате применения предложенного метода поле скоростей в
выделенной области интереса имеет вид, представленный на рис. 6.
Рис. 6. Поле скоростей в области интереса
Рассмотрим изображения, полученные при исследовании гепатобилиарной
системы человека [5, 22] . Печень — самый большой орган человеческого тела,
гепатобилиарная система включает в себя гепатоциты, желчный пузырь и
желчевыводящие пути. Динамическая сцинтиграфия гипатобилиарной системы
представляет собой комплексное исследование, включающее в себя оценку
функционального состояния печени, концентрационной и двигательной функции
желчного пузыря, проходимости желчных путей и наличия дисфункции
сфинктера Одди. Метод позволяет также исследовать и измерить рефлюкс желчи
из двенадцатиперстной кишки в желудок. Показаниями к проведению
исследования являются воспалительные и обменные заболевания печени,
21
желчного пузыря, дискинезии билиарного тракта, аномалии и пороки развития
желчевыделительной системы, заболевания желудочно-кишечного тракта и т.д.
Рис. 7. Радионуклидные изображения гепатобилиарной системы с выделенным
объектом исследования
Возьмем для обработки два радионуклидных изображения динамической
сцинтиграфии гипатобилиарной системы, на которых имеет место сдвиг желчного
пузыря (Рис. 7) с течением времени.
После применения представленного метода получим, что в области
интереса, выделенной пунктиром на Рис. 7, поле скоростей имеет вид,
изображенный на Рис. 8.
Рис. 8. Поле скоростей
Таким образом, метод определения поля скоростей может быть использован
для обнаружения движения органа или пациента в процессе проведения
динамических и томографических исследований и для дальнейшей его коррекции.
Решение данной задачи является важным этапом при радионуклидной
диагностике, так как даже небольшое смещение во время сбора данных может
повлиять на достоверность результатов.
22
§2.2. Построение поля скоростей для последовательностей
радионуклидных изображений
Предложенный метод определения поля скоростей может использоваться
для обработки последовательностей радионуклидных изображений, в частности,
для построения контуров исследуемых объектов.
Рассмотрим в качестве примера изображения, полученные в процессе
равновесной вентрикулографии сердца [5, 22]. Такое исследование используется
для качественной и количественной оценки состояния левого желудочка сердца,
определения подвижности его стенок, определения синхронности вступления
отдельных участков миокарда в сокращение [16, 40]. При его проведении
изучается влияние физической нагрузки и лекарственных препаратов на
сократимость миокарда. Проводится исследование, синхронизированное с
сигналом ЭКГ, в результате которого мы получаем серию сцинтиграмм,
соответствующую «представительному» сердечному циклу, таким образом, мы
получаем изображения сердца в различные временные интервалы сердечного
цикла.
Такое обследование подразумевает вычисление целого набора параметров
по кривой кровенаполнения левого желудочка сердца [33, 39]. А это значит, что от
того на сколько точно и правильно будет определен его контур, зависит
достоверность результатов диагностики.
Пусть на первом изображении произведено оконтуривание левого желудочка
вручную или автоматически (Рис. 9.).
Рис. 9. Первое изображение
23
последовательности с оконтуренной
областью интереса
Рис. 10. Поле скоростей в точках
контура
Далее, используя предложенный в работе способ, определим поле скоростей
в точках полученного контура (Рис. 10). Это поле дает точки нового контура —
контура области интереса на втором изображении.
Произведенные действия можно представить в виде алгоритма
оконтуривания областей интереса на последовательности изображений (Рис. 11).
24
Рис. 11. Алгоритм построения контуров области интереса на последовательности
изображений.
Применяя последовательно к каждой паре изображений представленный
способ вычисления оптического потока, получаем его контуры на всех
изображениях (Рис. 12).
Рис. 12. Последовательность изображений равновесной вентрикулографии сердца
с оконтуренным левым желудочком
Как правило, при проведении томографических исследований важно
обрабатывать большое количество изображений, что делает автоматизацию
оконтуривания их последовательности актуальной.
25
Заключение
В данной работе рассматривалась задача определения поля скоростей в
задачах ц и фровой обработки радионуклидных изображений и их
последовательностей. Актуальность проведенного исследования обусловлена
активным развитием ядерной медицины. Сегодня практически ни одна
диагностика заболеваний различного типа не обходится без проведения ОФЭКТ
или ПЭТ. Применение радионуклидных методов исследования помогает
обнаруживать болезни на ранних стадиях, локализовывать их очаг или уточнять
диагноз. Средства ядерной медицины и ее программное обеспечение постоянно
совершенствуются, что приводит к появлению все новых математических задач, в
том числе связанных с обработкой полученных изображений.
В ходе данного исследования была изучена литература, связанная с его
тематикой, достижения российских и зарубежных ученых в данной области. В
основной части работы были предложены два метода определения поля скоростей.
Один из них основан на предположении о постоянстве плотности распределения
радиофармпрепарата вдоль траектории движения объектов, другой предполагает
постоянство градиента плотности. Проблема нахождения поля скоростей сведена
к решению больших разреженных систем линейных уравнений. Полученные
системы линейных уравнений решены итерационными методами Гаусса-Зейделя
и последовательной верхней релаксации, показана их сходимость.
Полученные в процессе исследования методы использовались в качестве
алгоритмов, применяемых при цифровой обработке различных изображений, в
26
том числе и радионуклидных, и были реализованы в среде MATLAB. Таким
образом, в работе проведен системный анализ информации на основе
компьютерных методов обработки данных.
Результаты данного исследования были апробированы с использованием
радионуклидных изображений, полученных в результате исследований
проводимых на двухдетекторном гамма-томографе «ЭФАТОМ» [23, 26] в
«Федеральном научно-клиническом центре специализированных видов
медицинской помощии медицинских технологий Федерального медикобиологического агентства» (ФНКЦ ФМБА России, Москва). На данном томографе
обработка данных производится с помощью программного комплекса
«Диагностика» [8].
Рассмотренные методы могут быть использованы при обработке данных,
полученных в процессе радионуклидной диагностики, а именно для обнаружения
движения во время исследования и его коррекции, а также для построения
контуров областей интереса [4, 9, 11, 40].
27
Список литературы
1.
Голуб Дж., Ван Лоун Ч. Матричные вычисления. М.: Мир, 1999. 548 с.
2. Гонсалес Р., Вудс Р. Цифровая обработка изображений.
М.: Техносфера,
2006. 1072 с.
3. Гребенщиков В.В., Котина Е.Д. Физико-технические основы ядерной
медицины. СПб.: СПбГУ, 2007. 172 с.
4. Котина Е.Д. К теории определения поля перемещений на основе уравнения
переноса в дискретном случае // Вестник Санкт-Петербургского
университета. Серия 10: Прикладная математика. Информатика. Процессы
управления. 2010. № 3. С. 38–43.
5. Котина Е.Д. Некоторые вопросы моделирования динамических процессов в
радионуклидных исследованиях. СПб.: ВВМ, 2013. 150 c.
6. Котина Е.Д. Обработка данных радионуклидных исследований // Вопросы
атомной науки и техники. 2012. № 3. С. 195‒198.
7. Котина Е.Д. О сходимости блочных итерационных методов // Известия
Иркутского государственного университета. 2012. Т. 5. № 3. С. 41–55.
8. Котина Е.Д. Программный комплекс «Диагностика» для обработки
радионуклидных исследований // Вестник Санкт-Петербургского
университета. Серия 10: Прикладная математика. Информатика. Процессы
управления. 2010. № 2. С. 100‒113.
9. Котина Е.Д., Максимов К.М. Коррекция движения при томографических и
планарных радионуклидных исследованиях // Ве стник СанктПетербургского университета. Серия 10: Прикладная математика.
Информатика. Процессы управления. 2011. № 1. С. 29‒36.
10. Котина Е.Д., Овсянников Д.А., Плоских В.А., Бабин А.В., Тузикова О.Ф.
Обработка данных в радионуклидной диагностике // Ульяновский медикобиологический журнал. 2014. № 1. С. 174‒175.
28
11. Котина Е. Д., Пасечная Г. А. Определение поля скоростей в задачах
обработки изображений // Известия Иркутского государственного
университета. 2013. Т. 6. № 1. С. 48–59.
12. Котина Е.Д., Пасечная Г.А. Построение трёхмерного поля скоростей для
томографических исследований сердца // Устойчивость и процессы
управления. Материалы III международной конференции. 2015. С. 589–590.
13. Марусина М.Я., Казначеева А.О. Современные виды томографии. СПб.:
СПбГУ ИТМО, 2006. 132 с.
14. Марчук Г. И. Методы вычислительной математики. М.: Наука, 1980. 536 c.
15. Остроумов Е.Н., Котина Е.Д., Сенченко О.Р., Миронков А.Б.
Радионуклидные методы в кардиологической клинике // Сердце: журнал для
практикующих врачей. 2010. Т. 9. № 3. С. 190‒195.
16. Остроумов Е.Н., Котина Е.Д., Шмыров В.А., Слободяник В.В.,
Тонкошкурова В.В., Можейко Н.П., Ильинский И.М., Шумаков Д.В.
Кардиоресинхронизирующая терапия и перфузия миокарда левого и правого
желудочков // Вестник трансплантологии и искусственных органов. 2012. Т.
XIV. № 3. С. 60‒68.
17. Овсянников Д.А., Котина Е.Д. Определение поля скоростей по заданной
плотности заряженных частиц // Вопросы атомной науки и техники. 2012. №
3. С. 122‒125.
18. Ортега Дж. Введение в параллельные и векторные методы решения
линейных систем. М.: Мир, 1991. 367 c.
19. Панов Д.Ю. Справочник по численному решению дифференциальных
уравнений в частных производных. М.: Государственное издательство
технико-теоретической литературы, 1943. 128 c.
20. Прэтт У. Цифровая обработка изображений: в 2-х томах; Пер. с англ. — М.:
Мир, 1982. Т. 1: 312с., Т. 2: 480 с.
21. Радионуклидная диагностика. Под редакцией Лясса Ф.Н. — М.: Медицина,
1983. 304 с.
29
22. Радионуклидная диагностика для практических врачей. Под редакцией
Лишманова Ю.Б., Чернова В.И. — Томск: STT, 2004. 394 с.
23. Сидоров А. В., Новиков В. Л., Гребенщиков В. В. и др. Опытный образец
двухдетекторной томографической гамма-камеры // Вопросы атомной науки
и техники. Сер.: Электрофизическая аппаратура. 2006. № 4. С. 30–34.
24. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.:
Наука, 1974. 285 c.
25. Anandan P.A. A computational framework and an algorithm for the measurement
of visual motion // International Journal of Computer Vision. 1989. V. 2. P. 283–
310.
26. Arlychev M. A., Novikov V. L, Sidorov A. V., Fialkovskii A. M., Kotina E.D.,
Ovsyannikov D.A., Ploskikh V. A. EFATOM Two-Detector One-Photon Emission
Gamma Tomograph // Technical Physics. 2009. V. 54, No 10. P. 1539–1547.
27.Barron J., Fleet D. Performance of optical flow techniques // International Journal
of Computer Vision. 1994. V. 12. P. 43–77.
28. Bergholm F., Carlsson S. A theory of optical flow // Computer vision. Graphics
and Image Processing: Image Understanding. 1991. V. 53, No 2. P. 171–188.
29. Fleet D., Weiss J. Optical Flow Estimation // Mathematical Models in Computer
Vision: The Handbook. Chapter 15. Springer. 2005. P. 239–258.
30. Becciu A., van Assen H., Florack L., Kozerke S., Roode V., ter Haar Romeny B.
M. A Multi-scale Feature Based Optic Flow Method for 3D Cardiac Motion
Estimation // Proceedings of SSVM 2009, LNCS. 2009. V. 5567. P. 588–599.
31. Canclini S., Terzi A., Rossini P., Vignati A., La Canna G., Magri J.C., Pizzocaro
C., Giubbini R. Gated blood pool tomography for the evaluation of global and
regional left ventricular function in comparison to planar techniques and
echocardiography // Ital Heart. 2001. V. 2. P. 42–48.
32. Cremers D. A Multiphase Levelset Framework for Variational Motion
Segmentation // Lecture Notes in Computer Science, 2003. V. 2695. P. 599‒614.
30
33. DePuey E. G., Garcia E. Optimal specificity of thallium-201 SPECT through
recognition of imaging artifacts // Journal of Nuclear Medicine. 1989. Vol. 30. Pp.
441–449.
34. Germano G., Chua T., Kavanagh P. et al. Detection and correction of patient
motion in dynamic and static myocardial SPECT using a multi-detector camera //
Journal of Nuclear Medicine. 1993. V. 34. P. 1394–1395.
35. Germano G., Kavanagh P. et al. Quantitation in gated perfusion SPECT imaging:
The Cedars-Sinai approach // Journal of Nuclear Cardiology. 2007. V. 14, No 4. P.
433–454.
36. Horn B.K.P., Schunck B.G. Determining optical flow // Artificial intelligence.
1981. V. 17. No 11. P. 185–203.
37. Kotina E.D., Latypov V.N., Ploskikh V.A. Universal system for tomographic
reconstruction on GPUS // Вопросы атомной науки и техники. 2013. Т. 88. №
6. С. 175‒178.
38. Kotina E.D., Ostroumov E.N., Ploskih V.A. Left and right ventricular phase
analysis of gated SPECT myocardial perfusion // European Journal of Nuclear
Medicine and Molecular Imaging. 2012. V. 39. No 2. P. 213.
39. Kotina E., Ovsyannikov D., Ploskikh V., Latypov V., Babin A., Shirokolobov A.
Data processing in nuclear medicine // Cybernetics and Physics. 2014. Т. 3. № 2.
С. 55‒61.
40. Kotina E.D., Pasechnaya G.A. Optical flow-based approach for the contour
detection in radionuclide images processing // Cybernetics and physics. 2014. V.
3. No 2. P. 62–65.
41. Kotina E. D., Ploskih V. A. Data Processing and Quantitation in Nuclear
Medicine // Proceedings of RuPAC 2012. P. 526–528.
42. Kotina E.D., Ploskih V.A., Babin A.V. Radionuclide methods application in
cardiac studies // Problems of atomic science and technology. 2013. V. 88. No 6.
P. 179–182.
31
43. Kumar A., Tannenbaum A. R., Balas G. J. Optic Flow: A Curve Evolution
Approach // IEEE Transactions on Image Processing, 1996. V. 5. No 4. P. 598‒
610.
44. Lefebure M., Cohen L. D. Image Registration, Optical Flow and Local
Rigidity // Journal of Mathematical Imaging and Vision, 2001. V. 14. No 2. P.
131‒147.
45. Lucas B.D., Kanade T. An Iterative Image Registration Technique with an
Application to Stereo Vision // Proceedings of Imaging Understanding Workshop.
1981. P. 121‒130.
46. Ovsyannikov D. A., Kotina E. D., Shirokolobov A. Yu. Mathematical methods of
motion correction in radionuclide studies // Problems of atomic science and
technology. 2013. V. 88. No 6. P. 137–140.
47. Papenberg N. Highly accurate optic flow computation with theoretically justified
warping // International Journal of Computer Vision. 2006. V. 67. No 2. P. 141–
158.
48. Saad Y. Iterative Methods for Sparse Linear Systems. Philadelphia: Society for
Industrial and Applied Mathematics, 2003. 547 p.
49. Uras S., Girosi F., Verri A., Torre V. A Computational Approach to Motion
Perception // Biological Cybernetics, 1988. V. 60. P. 79‒87.
Отзывы:
Авторизуйтесь, чтобы оставить отзыв