Реконструкция МР-изображений с учетом неоднородностей полей Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

ISSN ÜS6S-5SS6

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2ÜÜ3, том ІЗ, № 2, с. І7-2І

300 лет Санкт-Петербургу —

Материалы XXXII конференции СПбГИТМО(ТУ)

УДК 538.69.083.2: 621.318.4

© В. А. Иванов, М. Я. Марусина, Н. Г. Рущенко, В. С. Сизиков

РЕКОНСТРУКЦИЯ МР-ИЗОБРАЖЕНИЙ С УЧЕТОМ НЕОДНОРОДНОСТЕЙ ПОЛЕЙ

Рассматривается задача реконструкции МР-изображений в условиях неоднородностей основного магнитного поля В и градиентных полей Ох, Оу, Ог. Неоднородности ДВ(х, у), ДGx (х, у), ДGy (х, у),

ДGz (х, у) должны быть известными детерминированными монотонными функциями. Задача решается модифицированным методом Кумара—Велти—Эрнста и методом регуляризации Тихонова.

ВВЕДЕНИЕ

В настоящее время существует несколько методов реконструкции МР-изображений. Реконструкция МР-изображения — это определение распределения плотности протонов (или спиновой плотности) по сечению с(х, у) или по объему с(х, у, г) на основе измеренного эхо-сигнала 8 (^) [1-5]. Одним из наиболее эффективных является метод Кумара—Велти—Эрнста [2-5], в котором плотность с определяется посредством двумерного или трехмерного преобразования Фурье над эхо-сигналом 8(^). Для реализации этого метода существует несколько практических схем [4, 5]. Однако в этом случае требуется высокая степень однородности основного магнитного поля

(относительная неоднородность ~ 10 6 +10 7) [6]

и градиентных полей (~10 3) [7], что приводит к большому количеству технических проблем.

Между тем в работах [8-10] предложен метод, математически учитывающий неоднородности полей и позволяющий решать задачу реконструкции МР-изображений в условиях значительной неоднородности магнитных полей (~ 10 1).

В данной статье проводится дальнейшее развитие метода, предложенного в работах [8-10] (см. также [5, с. 53-55]).

ЭКСПЕРИМЕНТ В СЛУЧАЕ ОДНОРОДНОСТИ ПОЛЕЙ И ЕГО МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ

Рассмотрим так называемую 1-ю практическую схему эксперимента [4, с. 272], [5, с. 48-49] (см. рисунок). Согласно этой схеме, вдоль г включается градиентное поле Ог = gzz . Вдоль х включается одно градиентное поле Ох = gxx (в прин-

ципе бесконечной длительности), т. е. создается частотное кодирование, поскольку этому соответствует постоянство лишь частоты Y gxx , где Y — гиромагнитное отношение. Вдоль у включается поочередно несколько градиентных полей Gy = gyy с различными gy продолжительностью Ту, т.е. при каждом значении gy создается

фазовое кодирование, поскольку каждому эхо-сигналу соответствует постоянство фазы Y gyyTy . При этом для селекции полагается, как

обЫчно что gz >> gx и gz >> gy [5, с. 46]. Помимо трех градиентных полей включается также статическое поляризующее однородное магнитное поле В0 = const.

Для получения эхо-сигнала от ансамбля протонов (или сигнала спада свободной индукции (ССИ)) создаются п/2- и п-импульсы Карра—Пар-селла (согласно методу эха Ганна). При этом, чтобы устранить расфазирование магнитных моментов протонов (спинов) вдоль z-координаты в пределах эффективной толщины слоя, используется методика Хоулта — изменение направления градиентного поля Gz на противоположное [5, с. 48].

Данный эксперимент требует идеальной однородности полей. Под этим будем подразумевать не только однородность основного поля (В0 = = const), но и изменение градиентных полей по точным законам: Gz = gzz , Gx = gxx ,

Gy = gy y .

В условиях такого эксперимента от z-слоя протонов будет исходить ряд эхо-сигналов S (столько сигналов в виде временных процессов S (t), сколько задано значений g y ). Двумерный эхо-сигнал, идущий от всего z-слоя, с учетом его волновой

Практическая схема эксперимента

природы можно записать в виде (ср. [5, с. 49])

S(X gy) = A J J c(x y) x

x exp[ iY(gxxt + Bot + gvyTv )] ] (1)

где 8 ^, gy ) = (?, gy ) — измеренный сильный

эхо-сигнал, с(х, у) = сг (х, у) — искомая плотность в г-слое, А — некоторая константа. Запишем (1) иначе:

8^, gy) = А | | с(х у) х

х ехр{/ {( В0 + gxx) t + У Tyygy )]} }у. (2)

Введя новые переменные £ = у(В0 + gxx), П = 7 Туу, получим (2) в виде двумерного преобразования Фурье:

S(t, gy) = D j j

П

Y Tv

х

где

Ygy і

V J

х exp[z[t + ngy)] ] dn,

D = A /(Y 2Tygx), ao = Y Bo-

(3)

Используя одно из свойств преобразования Фурье — масштабирование (scaling) [4, с. 24], полу-

чим:

F j j S(t, gy) х

X exp[- І [t + ngy )] dt dgy ,

Y gv Y Ty

(4)

где Р = 1/(4п2 Ц).

Формула (4) позволяет в принципе получить распределение плотности с = сг по измеренному эхо-сигналу 8 = 8Ю . При этом двумерное непрерывное преобразование Фурье (НПФ) (4) должно реализовываться через дискретное преобразование Фурье (ДПФ) или быстрое преобразование Фурье (БПФ). В этом случае задание равномерных сеток узлов по t, gy, £ и П породит также равномерные сетки по х = (£-^0)/(Г gx) и у =П/(У Ту ) для плотности с .

Для повышения устойчивости вычисления преобразования Фурье предлагается использовать метод регуляризации Тихонова [11]. В этом случае формула (4) запишется в виде регуляризированно-го решения:

c

[-®0 Y gy *Y Tv

\

= F

S(t, gy )

-х

+ aM (t, gy)

X exp[- i [t + ng )] dt dg,

(5)

где a > 0 — параметр регуляризации, a M(t, gy) — функция вида

/- \2n f \2n

M (t, gy ) =

t

gy

g y max V ' J

(6)

P( x, y) = Y [ gy о y + AGy(x, y)] Ty T = t - 2T .

причем n = 1, 2, 3,... — порядок регуляризации

(обычно n = ^ tmax и gy max — максимальные

значения t и gy .

Существует ряд способов выбора параметра регуляризации a [5]: способ невязки (использует значение среднеквадратической погрешности измерения S ), способ подбора и др.

Решение ряда примеров показывает [5, с. 170172], [11], что использование метода

регуляризации Тихонова уменьшает погрешность вычисления ca согласно (5) по сравнению с вычислением с согласно (4) (а значит, увеличивает отношение сигнал/помеха) в 2-3 раза.

ЭКСПЕРИМЕНТ В СЛУЧАЕ НЕОДНОРОДНОСТИ ПОЛЕЙ И ЕГО МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ

В предыдущем пункте мы описали эксперимент, использующий высокооднородные магнитные поля. Теперь рассмотрим эксперимент, использующий вместо Во = const неоднородное поле В(x, y) = В0 + ДВ(x, y), а вместо градиентных полей Gx = gxx и Gy = gyy — градиентные поля вида Gx + ДGx (x, y) и Gy + ДGy (x, y). Здесь ДВ(x, y), ДGx (x, y) и ДGy (x, y) — неоднородности полей. Схема такого эксперимента аналогична схеме эксперимента с однородными полями (см. рисунок).

Математически задача описывается следующим двумерным интегральным уравнением [5, с. 54], [8, 9]:

АII с(x, y) exp{ i [(x, y) т + P(x, y) p]} dy =

= S (т, p), (7)

где

^(x, y) = Y [Во + ДВ(x, y) + gxx + ДGX (x, y)],

Здесь ДВ(x, y), ДGX (x, y), ДGy (x, y) — неоднородности полей. Имеются в виду неоднородности, обусловленные в первую очередь техническими особенностями. Под ДВ подразумевается отклонение от В = В0 = const, а под ДGX и ДGy —

отклонения от законов gxx и gyy . Например, если магнитное поле создано одной катушкой в виде соленоида с однородным распределением тока вдоль ее обмотки, то поле В будет падать от центра катушки к ее краям в «1.5 раза в типичных случаях [5, с. 59], а от оси катушки к ее обмотке возрастать. Ставится задача: в условиях высокой неоднородности полей выполнить реконструкцию МР-изображения, т. е. найти функцию

с(x,y) = cz (x, y) путем решения уравнения (7).

Отметим свойства функций ДВ( x, y),

ДGX (x, y), ДGy (x, y), которым они должны удовлетворять, чтобы задача реконструкции имела решение и притом единственное. Функции ДВ( x, y),

ДGX (x, y), ДGy (x, y) должны быть известными

гладкими детерминированными (неслучайными) функциями x и y . Существующими в настоящее время техническими средствами можно добиться выполнения данных свойств неоднородностей магнитных полей.

Рассмотрим ряд градиентных полей по у , равных p [gy 0y + Д<^у (X, y)], где p ^tPmn^ Pmax ], gy 0 ф 0, например, gy 0 = min | g, | ф 0, а p = = ±1, ± 2, ...

Остановимся на способе решения уравнения (7) относительно c(x, y). Введем новые переменные

[ =П(x, у) ,| n = P(x,y). I

(8)

Систему нелинейных уравнений (СНУ) (8) относительно х и у запишем подробнее:

Y В0 + Y ДВ(x, y) + Y gxX + Y ДGx (x,y) = [ ,|

Y gy 0 Уту +Y^Gy(X У) Ty =П. I

(9)

Запишем систему (9) в виде, удобном для решения методом итераций:

x=

y=

[-&0 ДВ( x, У) + ДGx (x, у)

Y g x

Y gy 0Ty

gx

ДGy (x, у)

gy0

c

a

+

П

д х дх

д! дп

д У ду

д! дп

где ю0 = у В0. При этом dх dу =| J(!,п)^! dn: где

J (!,п) =

— якобиан преобразования.

В результате решения СНУ (10) некоторым итерационным методом (Ньютона, хорд и т.д.) получим

X = х(! ,п), |

у = у(!,п) ,

где х(!,п) и у(!,п) — некоторые численные зависимости.

Уравнение (7) можно записать в виде двумерного НПФ:

А Цс(х(!,п), у(!,п))х

: О!,п)}}^(т, р) ехр[- і (!т + п p)] dт dp,

где

О(!,п) = ■

1

ЯЫ Цу-

. + аМ (т, р) х ехр[- і (!т + пp)]dтdp,

где (ср. (6))

М(т,p) =

ґ \2п т

с \2п

p

Параметр регуляризации а может быть выбран различными способами, например способом подбора контрастности изображения (чем меньше а , тем выше контраст томограммы, и наоборот) или способом невязки из уравнения [5, с. 55], [11]

1 1 |~(т, Р)

а М (т, p)

1 + а М (т, p)

dт dp = 8‘

при условии

1 1 | ~(т, p) dтdp >8

где

< 8:

х ехр[/ т + п р)]\т (^ ,п) ап =

= ^ (т р). (11)

Обращение (11) с помощью обратного преобразования Фурье по аналогии с решением уравнения (3) приведет к формуле

Ас(х(£,п),у(£,п))\ )(£,п)\ =

= 4ТТ7 ^ (т,Р)ехр[-[ (^т + ПР)Мт&р,

откуда искомое решение равно

с(х ,п), у(£,п))=

(12)

4п2 (£,п)\

Для повышения устойчивости решения (12) используем метод регуляризации Тихонова. Решение

(12) с регуляризацией будет иметь вид:

с а , п) , У(^,п)) =

~~ В (т, р)

(13)

т.е. 8 — среднеквадратическая погрешность измерения эхо-сигнала В(т, р), полагаемая известной (В — точный эхо-сигнал, а В — измеренный эхо-сигнал).

В работе [8] приведены результаты численного решения уравнения (7) (без регуляризации). В модельных примерах полагались следующие неоднородности полей:

ДВ/В0 - 3.6-10-5 , ^Gx/gx - 3.4мм,

ДСу/ёу - 41 мм или Д?хА?х - 28 •10"2,

Дgy /gy - 2.8-10-2.

Однако данная методика может быть использована и в случае, когда , ^х^х и Д?у^у

имеют заметно большие значения. Необходимо только знать функции ДВ(х,у), ДGx(х,у) и ДGy (х, у). Они могут быть рассчитаны численно

на основе заданной конфигурации катушек и распределения токов в них, например, по методике, изложенной в работах [6, 7]. Кроме того, возможно, что использование метода инвариантных аппроксимаций [12] позволит существенно уменьшить влияние неоднородностей магнитных полей при решении уравнения (7). Однако этот вопрос требует специального рассмотрения.

Применение метода регуляризации Тихонова для повышения устойчивости решения (12), как показывают примеры, рассмотренные в работе [11], понижает отношение \ Дса \/са в 2-3 раза

+

т

тах

тах

2

2

2

2

Ь

2

(где Аса — погрешность решения са ).

При численной реализации метода непрерывное преобразование Фурье (НПФ) (13) заменяется на дискретное преобразование Фурье (ДПФ) или на быстрое преобразование Фурье (БПФ) на равномерных сетках узлов по Т, р,£,П • В результате решение са будет получено на неравномерных сетках узлов по x, y . Поэтому в дальнейшем необходимо применение интерполяции.

СПИСОК ЛИТЕРАТУРЫ

1. Иванов В. А. Внутривидение (ЯМР-томо-графия). Л.: Знание, 1989. 32 с.

2. Эрнст Р., Боденхаузен Дж., Вокаун А. ЯМР в одном и двух измерениях. М.: Мир, 1990. 200 с.

3. Физика визуализации изображений в медицине, т. 2 / Ред. Уэбб С. М.: Мир, 1991. 408 с.

4. Cho Z.H., Jones J.P., Singh M. Foundations of medical imaging. NY: Wiley, 1993. 400 p.

5. Сизиков В. С. Математические методы обработки результатов измерений: Учебник для вузов. СПб.: Политехника, 2001. 240 с.

6. Галайдин П.А., Замятин А.И., Иванов В.А. Расчет и проектирование электромагнитных систем магниторезонансных томографов. СПб.: Изд-во ИТМО, 1998. 29 с.

7. Марусина М.Я. Коррекция неоднородности основного магнитного поля МР-томографа на

постоянных магнитах. Автореф. дис. ... канд. техн. наук. СПб., 1993. 18 с.

8. Kawanaka A., Takagi M. Estimation of static magnetic field and gradient fields from NMR image // J. Phys. E: Sci. Instrum. 1986. V. 19. P. 871-875.

9. Sizikov V.S. Integral equations in NMR-tomography: reconstruction of NMR images with a regularization // Proc. of the 5th Int. Conf. IMSE98 / B.S. Bertram (ed.). Houghton, USA, 1998. P.74-75.

10. Лукьянович И.К., Савицкий А.А. ЯМР-томо-графия в нестабильном и неоднородном поляризующем магнитном поле // Прикладная спектроскопия. 1999. Т. 66, № 2. С. 270-274.

11. Сизиков В. С. Использование регуляризации для устойчивого вычисления преобразования Фурье // Ж. вычисл. матем. и матем. физики. 1998. Т. 38, № 3. С. 376-386.

12. Математические методы построения и анализа алгоритмов / Отв. ред. А.О. Слисенко. Л.: Наука, 1990. 238 с.

Санкт-Петербургский государственный институт точной механики и оптики (технический университет)

Материал поступил в редакцию 10.04.2003.

RECONSTRUCTION OF NMR IMAGES WITH REGARD TO FIELD INHOMOGENEITIES

V. A. Ivanov, M. Ya. Marusina, N. G. Rushchenko, V. S. Sizikov

Saint-Petersburg State Institute of Fine Mechanics and Optics (Technical University)

The problem of reconstruction of NMR images subject to inhomogeneities of the ground magnetic field B and gradient fields Gx, Gy, Gz is considered. The inhomogeneities AB(x, y), AGx (x, y), AGy (x, y),

AGz (x, y) must be known deterministic monotone functions. The problem is solved by the modified Kummar-Welti-Ernst method and the Tikhonov regularization method.