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

О ВОССТАНОВЛЕНИИ СМАЗАННЫХ ТОМОГРАММ РАЗЛИЧНЫМИ МЕТОДАМИ

М.В. Римских

Научный руководитель - д.т.н., профессор В.С. Сизиков

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

Введение

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

• искажения, устранение которых не требует математической обработки — царапины, штрихи на томограмме, неудачная яркость и контрастность;

• искажения, требующие примитивной математической обработки - геометрические искажения, требующие изменения масштаба по вертикали и/или горизонтали или устранения нелинейности;

• искажения, требующие сложной математической обработки - смаз, дефокусировка, зашумленность томограммы [1, 3-8].

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

Математическое описание задачи реконструкции смазанной томограммы

Рассмотрим данную задачу на примере смазанной томограммы (см. рис. 1). Полагаем, что за время экспозиции носитель томографического изображения совершил прямолинейный и равномерный сдвиг (смещение) на величину А. Можно также считать, что сдвиг совершил объект (подвижный пациент). В результате на носителе (томограмме) зафиксируется смазанное изображение типа рис. 1.

Определим по штрихам на томограмме направление смаза и его величину А. Направим вдоль смаза ось x, а перпендикулярно ему - ось у. Математически задача смазывания изображения описывается соотношением [3-8]:

1 x+А

А JУ)ds = g(x у), (1)

x

где g ( x, у) - распределение интенсивности вдоль смазанной томограммы (измеренная функция), а w(s, у) - распределение интенсивности вдоль неискаженной томограммы, той томограммы, которая была бы получена в отсутствие сдвига, т.е. при А = 0 (искомая функция). При этом sOy - неподвижная система координат, а xOy - система координат, связанная с движущимся носителем изображения.

Отметим, что под функциями g (x, у) и w(s, у) подразумеваются только амплитуды излучения без учета фазы (как в голографии). Отметим также, что если на томограмме фиксируется серое изображение (gray image), то под g(x, у) и w(s, у) будем подразумевать интенсивности gg (x, у) и wg (s, у) в сером цвете. Если же фиксируется цветное изображение (RGB image), то можно преобразовать RGB-изображение в gray-

изображение (это особенно эффективно выполняется в Ма1ЬаЬ'е), или под записью (1) подразумеваем три соотношения для трех цветов - красного, зеленого и синего (Я, О, В), причем с помощью светофильтров нужно получить распределение интенсивностей §к(х,У), (х,У), §в(х,У), затем восстановить (5,у), (5,у), (5,у) и, наконец, вычислить суммарную интенсивность w(s, у) = wR (5, у) + wG (5, у) + wв (5, у).

Рис. 1. Смазанная томограмма. g (х, у) - распределение интенсивности по изображению, А - величина сдвига (смаза), ось х направлена вдоль смаза

Преобразование соотношения (1)

Соотношение (1) является неклассическим одномерным интегральным уравнением Вольтерра I рода относительно искомой функции w(s, у) при каждом фиксированном значении у, играющем роль параметра, другими словами, (1) есть множество одномерных интегральных уравнений.

Отметим, что в ряде работ [3, 4] рассмотрены более сложные задачи - неравномерный и/или непрямолинейный сдвиг изображения и т.д.

Запишем уравнение (1) в виде (опустив для простоты у):

1 х+А

А | w(s) ё5 = g(х), (2)

х

где ^5) = ^(5) = у), g(х)=gy(х) = g(х у).

Методы решения уравнения (2) недостаточно проработаны в виду того, что оно является неклассическим, нестандартным (оба предела интегрирования переменны, отсутствует в явном виде ядро) [9]. Однако уравнение (2) можно преобразовать к уравнению в стандартной форме, а именно, к одномерному интегральному уравнению Фред-гольма I рода типа свертки [3-8]:

да

| к ( х - 5) w(s)= g(х), - да < х < да, (3)

—да

где

[1/ А, х е[—А,0],

к (х) = <! ' ' 1 ' ь (4)

[0, х е[-А,0].

Задачи решения уравнений (2) и (3) являются некорректными (в первую очередь, неустойчивыми), причем «степень неустойчивости» уравнения (3) выше, чем уравнения (2) [6, 7, 9]. Однако устойчивые методы решения уравнений типа (3) достаточно подробно разработаны [6, 7, 10, 11]. Основным (устойчивым) методом решения уравнения (3) считается метод преобразования Фурье (ПФ) с регуляризацией Тихонова [3, 4, 6, 7, 10, 11]. Назовем его первым методом.

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

Два метода решения интегрального уравнения

Регуляризованное (устойчивое) решение уравнения (3) методом ПФ с регуляризацией Тихонова (первым методом) имеет вид [3-7, 10, 11]:

1 м

^ (®) = — | (®) ^-'Ю^ , (5)

—<Х)

где регуляризованный Фурье-спектр решения равен

»а(«) = Г)С(шР, (6)

Дш) + аш у

причем

<Х)

К(ш) = Г к(х) Л = + с°5(шА) — 1 г, (7)

шА шА

—<Х)

<Х)

О(ш) =| g(х) егшх Сх (8)

— да

- Фурье-спектры ядра к(х) и правой части g(х) уравнения (3), а > 0- параметр регуляризации, р = 1,2,3, к - порядок регуляризации (обычно р = 1).

При практической реализации первого метода вычисления одномерных обратного и прямого преобразований Фурье (5), (7), (8) выполняются в виде дискретного преобразования Фурье (ДПФ) или быстрого преобразования Фурье (БПФ).

Выбор а можно осуществлять, например, способом невязки или обобщенным принципом невязки [4, 10, 11]. Однако для задачи реконструкции изображений, как показала практика, более эффективен способ подбора (подробности см. в [5-8]).

Теперь рассмотрим решение уравнения (3) методом квадратур с регуляризацией Тихонова (вторым методом). Будем полагать, что при некотором фиксированном у правая часть g(х) задана при х е [с,С], а функция ^(5) ищется при £ е[а, Ь] (обычно [а,Ь] с [с,С]), т.е. вместо уравнения типа свертки (3) рассматривается уравнение общего типа:

Ь

А*м = Г к (х, 5) ^(5) С5 = g(х), с < х < С , (9)

где

[1/ А, х — 5 е[—А,0], к(х, 5) = к (х — 5) = Г ' А ' (10)

[0, х — 5 ё[—А,0],

A - оператор.

Применение метода регуляризации Тихонова к уравнению (9) порождает интегральное уравнение IIрода [6, 7, 10, 11]:

и

a wа (t) + J R(t, s) wa (s) ds = F(t), a < t < b,

a

с новым (симметричным, положительно определенным) ядром

d

R(t, s) = R(s, t) = J k (x, t) k (x, s) dx

(11)

(12)

и с новой правой частью

d

F (t) = J k ( x, t) g ( x) dx.

(13)

Для численного решения уравнения (11) используем метод квадратур. Полагаем, что при некотором y правая часть g(x) задана, а решение w(s) ищется на равномерных

сетках узлов (соответственно): xi = c + h • i, i = 0,1,...,m , sj = tj = a + h • j , j = 0,1,...,и ,

где h = Ax = As = At = const - шаг дискретизации по x, s и t (более общие случаи неравномерных сеток рассмотрены в [6, 7, 10, 11]). Интегралы в (11)—(13) аппроксимируем конечными суммами (обычно по формуле трапеций).

В результате вместо интегрального уравнения (11) получим систему линейных алгебраических уравнений (СЛАУ) [6, 7, 10, 11]:

(aC + ATA)wa= ATg, (14)

где C — некоторая диагональная матрица (в частности, единичная матрица), A — матрица M х N, A — транспонированная матрица N х M, где M = m +1, N = п +1 (подробности см. в [6, 7, 10, 11]).

T

У системы уравнений (14) матрица a C + A A является квадратной, симметричной, положительно определенной, поэтому для решения СЛАУ (14) обычно используется метод Холецкого.

Численная иллюстрация методов

Первые и второй методы проиллюстрируем численным примером. Был решен следующий модельный пример (ср. [6, стр. 198], [7, p. 162]). При некотором фиксированном значении y распределение интенсивности вдоль s на неискаженном томографическом изображении описывается следующей функцией (суммой 5 гауссиан):

w(s) = 6.5exp

+ 14exp

—

—

s + 0.66 0.085

s - 0.41 0.095

+ 9exp

+ 9exp

—

—

s + 0.41 0.075

s - 0.67 0.065

+12exp

—

s - 0.14 0.084

+

(15)

e [a,b], a = -0.85, b = 0.85, ядро k(x) аппроксимировалось также гауссианой:

(16)

k (x) = J2. e -

- q- 2

П

q = 47. Такое ядро примерно соответствует ядру (4) при А^ 0.3. Сначала решалась прямая задача: путем численного интегрирования вычислялась функция

2

2

2

2

2

ь

Е( х) = | к (х — я) ^(я) ds, с < х < d, (17)

а

на дискретной сетке узлов: х = с,с + И,...,d, причем к = 0.0125, с = —1, d = 1. Затем к значениям g(х) были добавлены погрешности 5 g, распределенные по нормальному закону с нулевым математическим ожиданием и среднеквадратическим отклонением а = 0.05, что соответствует относительной погрешности « 1 % .

14 12 10 8 6 4 2 0

~2-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 х,и 1

Рис. 2. Реконструкция изображения первым методом при некотором фиксированном значении у. 1 - интенсивность ^(я) вдоль неискаженного изображения, 2 - интенсивность g(х) вдоль смазанного изображения, 3 - ядро к(х) интегрального уравнения, 4 -реконструированное распределение интенсивности (5) после перестановки частей

решения на интервалах [—1,0] и [0,1]

На рис. 2 представлены: интенсивность ^(я) вдоль неискаженного изображения (кривая 1), интенсивность g(х) вдоль смазанного изображения (кривая 2) и ядро к(х) (кривая 3). Видно, что флуктуации, присутствующие в функции м>(я), в значительной степени заглажены в функции g(х).

Затем решалась обратная задача - восстановление ^(я) по измеренной g( х) = g (х) + 5g с помощью первого и второго методов.

Решение первым методом (методом ПФ с регуляризацией Тихонова) уравнения (3) с правой частью х) осуществлялось по формулам (5)-(8), причем непрерывные преобразования Фурье (НПФ) заменялись дискретными преобразованиями Фурье (ДПФ). При этом использовались одинаковые носители для функций, входящих в (3): зиррк(х) = зиррх) < [с,d], где с = —1, d = 1. Шаг дискретизации к = Ах = Ая = 0.0125, число шагов т = ^ — с) /к = 160, а число узлов М = т +1 = 161. На рис. 2 представлено регуляризованное решение (реконструированное распределение интенсивности при некотором фиксированном значении у) (я), рассчитанное первым методом при р = 1, —8

а = 10 (кривая 4).

,1

М

1

Г 1 / I У! \ч

/л\ \ ^ Л ]/ 1 ч I \Г

! / 5у; г 1"* к 2 1

■ / г/ 4 и\ V > / / * / 1 * ' \ и <

• V \/

Рис. 3. Реконструкция изображения вторым методом при некотором фиксированном значении у. 1 - интенсивность ^(5) вдоль неискаженного изображения, 2 - интенсивность g(х) вдоль смазанного изображения, 3 - реконструированное распределение

интенсивности (5)

Решение вторым методом (методом квадратур с регуляризацией Тихонова) уравнения (9) с правой частью х) осуществлялось по формулам (11)—(14) на сетках узлов: х = с,с + И,...,С , 5 = а,а + И,...,Ь, t = а,а + И,...,Ь , причем а = —0.85, Ь = 0.85, число узлов по х равно М = 161, а по 5 и t равно N = 137. На рис. 3 представлено регуляризо-

—3 5

ванное решение (5), полученное вторым методом при а = 10 ' (кривая 3).

Сравнение методов и выводы

Сравнение двух методов позволяет сделать следующие выводы.

1. Оба метода при удачно выбранных направлении смаза изображения, величине смаза А и параметре регуляризации а позволяют достаточно эффективно восстанавливать смазанные томографические изображения (см. рис. 2 и 3, а также [5], [6, стр. 65], [7, р. 193-197]).

2. Точность восстановления обоих методов приблизительно одинакова. Причем оба метода дают на краях, а также в местах широкой нулевой области решения ложные знакопеременные флуктуации - проявление эффекта Гиббса. Этот эффект можно исключить, полагая (5) = 0, если (5) < 0 . Но можно использовать более совершенный подход, а именно, метод регуляризации Тихонова с ограничением на решение (5) в

виде его положительности [11, стр. 118]. Однако такой подход наиболее эффективен при использовании только квадратур, но не ПФ.

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

4. Первый метод вследствие замены НПФ на ДПФ отягощен такими особенностями ДПФ, как периодичность, эффект наложения, равенство длин носителей функций и др. [6, стр. 169], [11, стр. 39, 123]. В результате этого решение (5) первым методом сначала получилось таким, как на рис. 4, а окончательное решение на рис. 2 (кривая 4) получено путем смещения решения (5), 5 < 0 вправо на период С — с = 2, другими словами, пу-

тем перестановки частей решения на интервалах [-1,0] и [0,1]. Второй метод такими особенностями не отягощен. 14

12 10 8

-2

1 V

А • * / 1 1 1

/ V | / 1 1 ( 1

| | / | / 1л 1 /' 1 ^ 1 1 1 N 1

\ [ 1 \ Я \ ' 1 4 П V ** 1/ 1 к 1 / 1 1/ V п

у ^ г А/ 1 \ К1 V; \ I V \ 1 1 1 /

1 :' / (I / / / / I \ ' ч / / \\ V ч л,

• « ' V/ V ч / * V

1 -0,8 -0,6 -0.4 -0.2 0 0.2 0.4 0,6 0Лх,5 1

Рис. 4. Реконструкция изображения первым методом при некотором фиксированном значении у. 1, 2, 3 - как на рис. 2, 4 - реконструированное распределение интенсивности м>а (без перестановки частей решения

5. Если в первом методе использовать БПФ для ускорения обработки, то это потребует, чтобы при каждом значении у число отсчетов М вдоль х (см. рис. 1) было целой степенью числа 2, иначе измерения функции g нужно искусственно добавлять нулями или же использовать обычное ДПФ, что может понизить скорость обработки. Во втором методе числа узлов М и N могут быть произвольными (вообще говоря, зависящими от у). Если же требуется высокая скорость обработки, то второй метод может быть модифицирован следующим образом.

Запишем решение СЛАУ (14) в виде:

Уа = Ва g , (18)

где

Ва =

(а С + ЛГЛ)_1 Л1

(19)

- матрица N х М, которая может быть заранее рассчитана (может быть одна и та же для всех значений у, а может быть для каждого у своя). В результате реконструкция изображения сведется к умножению матрицы В на вектор g при каждом значении у. Это

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

Работа выполнена при поддержке РФФИ (грант № 05-08-01304-а).

Литература

1. Бейтс Р., Мак-Доннелл М. Восстановление и реконструкция изображений. М.: Мир, 1989. 336 с.

2. Физический энциклопедический словарь / Гл. ред. Прохоров А.М. М.: Сов. энциклопедия, 1984. 944 с.

3. Тихонов А.Н., Гончарский А.В, Степанов В.В., Ягола А.Г. Обратные задачи обработки фотоизображений. Некоторые задачи естествознания. / Под ред. Тихонова АН., Гончарского А.В. М.: Изд. МГУ, 1987. С. 185-195.

4. Бакушинский А.Б., Гончарский А.В. Некорректные задачи. Численные методы и приложения. М.: Изд. МГУ, 1989. 199 с.

5. Сизиков В.С., Белов И.А. Реконструкция смазанных и дефокусированных изображений методом регуляризации. // Оптический журнал. 2000. Т. 67. № 4. С. 60-63.

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

7. Petrov Yu.P., Sizikov V.S. Well-posed, ill-posed, and intermediate problems with applications. Leiden-Boston: VSP, 2005. 234 p.

8. Сизиков В.С., Российская М.В., Козаченко А.В. Обработка смазанного изображения методами дифференцирования, преобразования Хартли и регуляризации Тихонова. // Изв. вузов. Приборостроение. 1999. Т. 42, № 7. С. 11-15.

9. Апарцин А.С. Неклассические уравнения Вольтерра I рода: теория и численные методы. Новосибирск: Наука, 1999. 193 с.

10. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. Киев: Наукова думка, 1986. 544 с.

11. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. М.: Наука, 1990. 232 с.