Численное моделирование скалярной лучевой томографии с помощью преобразования Гильберта Текст научной статьи по специальности «Математика»

ВЕСТНИК ЮГОРСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2011 г. Выпуск 2 (21). С. 29-34

УДК 519.642

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СКАЛЯРНОЙ ЛУЧЕВОЙ ТОМОГРАФИИ С ПОМОЩЬЮ ПРЕОБРАЗОВАНИЯ ГИЛЬБЕРТА

В. А. Гафт

1. Теоретические сведения

1.1 Постановка задачи

Пусть D - ограниченная строго выпуклая область в евклидовой плоскости R2 с границей Г. Пусть о = {4 |,| 41= 1} - единичная окружность. Введем расслоение единичных окружностей над D:

Q = D хо = {(х, у, 4)1, (х, y) е D,4 = (cos р, sin р)}.

Его граница состоит из объединения «входящих» и «выходящих» векторов:

дО = дО U дО д±0 = {(х ^ ф) 1 ±(4, v(х у)) > 0},

где v - вектор внутренней нормали к Г.

Пусть f - гладкая функция на О . Рассмотрим интеграл вдоль луча, выходящего из точки (х,у) в направлении вектора 4 = (cos р, sin р)

т( х, у ,р)

и (х, у ,р) = J f (х +1 cos р, у + sin p,p)dt, (1)

0

где т(х, у, р) - длина отрезка луча, выходящего из точки (х,у) в направлении вектора 4 до пересечения с границей Г.

Рисунок 1. Функция и

Естественно положить, что

Т !я о= 013 _ о

Тогда и

и |я о= 013 _ о

Функция

1/ = и !а+о

называется лучевым преобразованием функции/. В случае, когда/ = Ах,у) будем писать IVи называть оператор 10 скалярным лучевым преобразованием.

Задача обращения скалярного лучевого преобразования состоит в нахождении функции Лх,у) по известному I'$.

Определим оператор переноса:

H = cosp------+ sinp—.

дх ду

Тогда задачу обращения If можно записать следующим образом: требуется найти правую часть уравнения (не зависящую от угловой переменной)

Hu = — f, если

u \д_0= 0 4+0= 10 f.

В такой форме задачу решать удобнее, чем находить f непосредственно из равенства (1) (при (х, у) е Г ).

Ключевую роль в решении задачи обращении скалярного лучевого преобразования играет преобразование Гильберта.

1.2 Преобразование Гильберта

Рассмотрим разложение гладкой 2п — периодической функции u (р) в ряд Фурье: u(р) = — + ^°°_1 (an cosnp+ bn sinnp),

Тогда функция Sn1(—bn cosпр + —n sinпр), ре [0,2п]

называется преобразованием Гильберта от функции u, которое будем обозначать Hu. Перечислим основные свойства этого преобразования:

1. H: C00 (о) ^ C00 (о),

2. H: L (о) ^ L (о),

3. H1 = 0,

4. H2u = -u + щ,

5. (Hu,v) = - (Hv,u).

1.3 Формула обращения

Для получения формулы обращения, сначала разложим функцию u (х, у, 4) на четную и нечетную часть относительно вектора 4;

u = u+ + u-,

u+(х у,4) = 2(u( х у,4) + u( х у,—4)),

u—(х у,4) = 2(u (х у,4) — u( х у,—4)).

Тогда Hu+ = 0, Hu = —f.

Имеет место следующая формула коммутации преобразования Гильберта и оператора переноса [1, 2]:

[H H]u = HHu - HHu = (H±u)0 + H±u0,

где ноль означает усреднение по окружности:

f2n

J0 v(х, ур^р

Тогда применяя преобразование Гильберта к равенству Hu— = — f , получим HHu- = —(H]U— )0 — H± (u— )0 — Hf = —(H±u)0.

Можно показать, также, что (H±u)0 = 0 [1, 2]. Таким образом, приходим к равенству HHu- = 0.

Это равенство означает, что функция Hu- постоянна вдоль луча х +1 cos р, у +1 sin р. Применяя еще раз преобразование Гильберта, получим 0 = HHHu- = HH2 u- + H± Hu-)0 = -Hu- + (H± Hu-)0 = f + (H± Hu-)0,

f = - (H ± Hu-)0.

Рассмотрим функцию u— (х + г(х, у, р) cos р, у + г(х, у, р) sin р, р).

Точка (х', у') = (х + т(х, у, p)cosp, у + т(х, у, p)sinp) (точка выхода луча) принадлежит границе Г. Но при (х, у) е Г

u—(х у,4) = <

1[u(х, у,4) — u(х, У, —4)] = 2 l0f (х У,4), при (х, у,4) ед+0 1[u(х, у,4) — u(х, ^ —4)] = — 2l0f (х, y,4), при (^ у,4) ед—0,

Т. е. u— \д0 = 2(I{)f) , где (I0 f) - нечетное продолжение функции If на все д0.

Таким образом, получаем формулу обращения:

1 (*2п д д

f (х, у) = —— l0 (cos^-—sin^)w(х, y,ф)dф, (2)

2П0 дх ду

w(х У, р) =1 Hlf (х + т(х, у, р) cos р, у + т(х, у, р) sin р, р). (3)

1.4 Особенности алгоритма

Отличие описанного алгоритма от существующих алгоритмов [3], решающих аналогичные задачи (обратное преобразование Радона, обратное верно-лучевое преобразование), состоит в использовании для обращения преобразования Гильберта, вместо преобразования Фурье по пространственным переменным, которое в них используется. Преобразование Гильберта здесь «работает» в касательном пространстве единичных окружностей, и поэтому, данный метод может применяться в случае интегрирования по геодезическим кривым двумерной римановой метрики общего вида. При этом формула обращения, практически, та же самая, но функция w определяется более сложным образом [1, 2].

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

2. Особенности численного моделирования

В решении задачи использовался метод конечных элементов. То есть в области (брался единичный круг) вводилась триангуляционная сетка. Для построения триангуляции (Делоне) использовались встроенные функции Matlab. Использование триангуляции для решения дает ряд преимуществ:

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

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

3. Простое получение матрицы масс для триангуляционной сетки, в связи, с чем не требуется отдельная реализация ее расчета. Матрица масс используется при вычислении дивергенции векторного поля.

Таким образом, использование метода конечных элементов упрощает задачу, чем, например, использование представления области по пикселям.

Как видно из формул (2) и (3), основная сложность реализации заключается в численном дифференцировании. Точнее, в расчете дивергенции векторного поля. При решении обратной задачи, основной процент ошибки будет появляться именно в этом месте. Поэтому, для выявления наилучшего метода, были реализованы несколько методов численного дифференцирования.

2.1 Вычисление с использованием формулы Грина (обобщенная производная)

Интегрируя тождество ) = (Ур, V) + (рсИуУ

при помощи формулы Гаусса-Остроградского, имеем:

\pidivVdD = р (V, у)СГ-1 (Ур ,У)СВ, (4)

В Г Б

где и СГ - элемент объема области В и ее границы Г. Разложим функцию divV по стан-

дартным базисным функциям метода конечных элементов:

N

=Е ср ,

к=1

где N - количество узлов конечноэлементной сетки.

Подставляя последнее равенство в формулу (4), получим:

]Тск | ррСВ = | (7, у)рСГ -1(у рз,V )СВ.

к=1 В Г В

Здесь, первый интеграл под знаком суммы - определяет матрицу масс. Все входящие сюда интегралы легко вычисляются подобно тому, как это описано в справке МаЛаЬ относительно метода конечных элементов.

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

2.2 Вычисление с интерполяцией на равномерную сетку

В данном случае дивергенция вычисляется следующим образом: полученное векторное поле (на триангуляционной сетке) интерполируется на равномерную сетку. Затем частные производные вычисляются как отношение приращения функции к приращению аргумента. И дивергенция - как сумма частных производных.

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

К достоинствам относится очевидность и простота реализации.

2.3 Вычисление с использованием дифференцирующих матриц

В данном случае снова используется метод конечных элементов. Разложим функцию / базисным функциям:

N

/(х)=Ё ^Рк(х).

к=1

Здесь х означает точку круга. Заметим, что в силу свойств базисных функций метода конечных элементов коэффициенты разложения совпадают со значениями функции в узлах сетки: /к = /хк). Разлагая и производную (штрих здесь означает дифференцирование по одной из переменных), будем иметь:

1'(х) = Е Р (х) = Е скРк (х X

к=1 к=0

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

N. N .

Е ск 1 РкР1Сх=Е 1 \р'к рА.

к=0 к=0

Обозначим Ыы = 1 ррСх и Лк1 = ^р\ р^х.

Таким образом, приходим к уравнению Мс = ЛГ

откуда с = М Л/(заметим, что матрица масс М положительна).

Матрицу В = М-1Л будем называть дифференцирующей (по той переменной, по которой вычисляется производная в матрице).

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

Замечание. Задача численного дифференцирования - слабо некорректная задача и требует регуляризации. В случае использования матриц В естественная регуляризация состоит в «обрезке» по большим сингулярным числам в ее Буё-разложении. В настоящей работе такая регуляризация не проводилась, поскольку задачи с зашумленными данными не рассматривались.

Примечание. Система Ма^аЬ позволяет обращать матрицы без использования дополнительных функций, однако, в зависимости от вычислительных возможностей компьютера, может возникнуть проблема нехватки памяти при обращении больших матриц. В этом случае предлагается вместо обращения матрицы использовать функции решения систем линейных алгебраических уравнений (Ьи- и ^К-разложения). Это позволяет сократить затраты памяти вдвое.

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

3. Результаты численных экспериментов

Ниже представлены примеры выполнения реализованной программы: исходные изображения / (слева) и результаты их реконструкции по /</ (справа).

Рисунок 2. Модель № 1 33

Рисунок 3. Модель № 2

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

Заключение

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

Небольшие изменения, внесенные в программу, позволяют перейти от скалярного лучевого преобразования к векторному.

Планируется дальнейшее проведение работ по улучшению качества и увеличению скорости выполнения алгоритма.

Как уже говорилось выше, данный метод решения может быть перенесен на случай криволинейных лучей - геодезических кривых, что и планируется сделать.

Автор благодарит научного руководителя Л. Н. Пестова за постановку задачи и внимание к работе.

ЛИТЕРАТУРА

1. Пестов, Л. Н. Вопросы корректности задач лучевой томографии [Текст] / Л. Н. Пестов // Сибирское науч. изд-во. - Новосибирск, 2003. - 110 с.

2. Pestov, L. On characterization of the Range and Inversion Formulas for the Geodesic X-Ray Transform [Text] / L. Pestov, G. Uhlmann. International Mathematical Research Notices. -2004. - p. 4331-4347.

3. Тихонов, А. Н. Математические задачи компьютерной томографии [Текст] / А. Н. Тихонов, В. Я. Арсенин, А. А. Тимонов. - М. : Наука, Гл. ред. физ.-мат. лит., 1987. - 160 с.