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

УДК 621.397.3:519.642.3

В. С. Сизиков, К. А. Кирьянов, Р. А. Экземпляров

ДВА БЫСТРЫХ АЛГОРИТМА ВОССТАНОВЛЕНИЯ СМАЗАННЫХ ИЗОБРАЖЕНИЙ

Представлены два быстрых алгоритма восстановления смазанного изображения, позволяющих численно решать некорректное интегральное уравнение Фредгольма I рода методом регуляризации Тихонова с использованием заранее рассчитанной („заготовленной") матрицы. Алгоритмы сводятся к умножению этой матрицы на каждую вектор-строку смазанного изображения. Выполнено сравнение с другими быстрыми алгоритмами.

Ключевые слова: смазывание изображения, восстановление изображения, интегральное уравнение, быстрые алгоритмы, „ заготовленная " матрица.

Введение. Задача восстановления искаженных, в частности смазанных, изображений является актуальной [1—5]. Смазывание изображения объекта может произойти из-за сдвига цифрового фотоаппарата за время экспозиции, перемещения самого объекта (самолета, автомобиля), несовпадения угловой скорости вращения телескопа и небесной сферы и т.д. Восстановление смазанного изображения сводится к решению множества одномерных интегральных уравнений (ИУ) Фредгольма I рода. Задача решения таких ИУ является некорректной (существенно неустойчивой) [6], обычно используются метод регуляризации Тихонова [1, 2, 4—6], метод фильтрации Винера [3—6] и другие устойчивые методы [4].

Метод регуляризации Тихонова обычно реализуется с использованием преобразования Фурье (ПФ) в случае ИУ типа свертки, когда функция рассеяния точки (ФРТ) является пространственно-инвариантной [1—6], или с использованием квадратур в случае ИУ общего типа (когда ФРТ является пространственно-неинвариантной) [3—6]. При этом иногда не требуется высокая скорость восстановления (восстановление старой смазанной фотографии, изображения астрономического объекта). Если требуется восстановление в режиме реального времени (восстановление опознавательных знаков на самолете-нарушителе границы, повышение качества теле- или кинокадров, томограмм динамической работы органов и др.) [4, 5], то для решения ИУ необходимо использовать быстрые алгоритмы.

В настоящей работе представлены два модифицированных быстрых алгоритма и выполняется их сравнение с существующими алгоритмами.

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

| И(х-£) ^ (£) = gy (х) + 5g, -ю< х <ю, (1)

где

(1А, -А< х < 0,

К х) = Г ' ' (2)

[0, иначе,

— ФРТ, обычно пространственно-инвариантная (разностная), н и g — распределение интенсивности по неискаженному (искомому) и искаженному (измеренному) изображению соот-

ветственно, А — величина смаза, 8g — помеха. В (1) ось х направлена вдоль смаза, а у — перпендикулярно смазу и играет роль параметра. Фактически (1) есть множество одномерных интегральных уравнений. Пространственная инвариантность функции к имеет место в случае, когда величина А не зависит от х (но может зависеть от у, в этом, довольно редком, случае А = Ау, к = ку ).

Если значение А зависит от х (пример: смазы различны на снимке, сделанном неподвижным фотоаппаратом, где представлены легкоатлеты, бегущие с разной скоростью), то уравнение задачи о смазывании изображения будет иметь вид:

ь

\к(х, = gy(х) + 5^, с < х < й, (3)

а

где ФРТ равна (ср. с (2))

fV A r, -А г < х— £< 0 или х <£< x + А x, h( х, Q = Г x' х ^ ^ - (4)

[0, иначе.

Как в случае уравнения (1), так и (3), необходимо решить столько самостоятельных (не образующих систему) одномерных уравнений, сколько строк содержит изображение.

Решение ИУ (1) и (3). Чтобы решение wy (£) было устойчивым, воспользуемся методом

регуляризации Тихонова. Применительно к ИУ (1) решение методом преобразования Фурье с регуляризацией Тихонова имеет вид [1, 2, 4—6]:

1 ^ H*(ш)Gy(ш) d (5)

Way (О = — J --,2 y 2 p е ^ d(D , (5)

2п—ГО H(ш)|2 +аш2p

где

ГО ГО

Gy (ш) =J gy(x)ei&xdx, H(ш) =J h(x)eiaxdx (6)

—ГО —ГО

— преобразования Фурье от gy(x) и h(x) соответственно, a > 0 — параметр регуляризации, а p > 0 — порядок регуляризации (обычно 1 или 2). При компьютерной реализации формул (5) и (6) все непрерывные преобразования Фурье заменяются на дискретные, а также быстрые ПФ (БПФ) [1, 4, 6]. Число отсчетов по x и по ш полагаем одинаковым и равным п.

Для решения ИУ (3) и (1) применим метод конечных сумм (квадратур) с регуляризацией Тихонова, согласно которому интеграл при каждом x заменяется конечной суммой с целым шагом дискретизации (в пикселах) и получается система линейных алгебраических уравнений (СЛАУ) при некотором фиксированном у:

Awy = gy , (7)

где A — (пхп)-матрица, связанная с ФРТ (4): Aik = h(xi, ^) , причем xi =^i = 1,2,..., п — дискретные целочисленные отсчеты по пикселам вдоль y-й строки (п — число столбцов в изображении). Полагаем, что векторы Wy и gy имеют размер (длину) п (в [5, с. 89—90] рассмотрены другие варианты размеров A, Wy и gy ).

Согласно методу регуляризации Тихонова вместо неустойчивой СЛАУ (7) решается:

(aE + A7 A) Way = ATgy, (8)

где E — единичная (пхп)-матрица, таким образом

Way = (aE + ATA)—1 ATgy. (9)

Способы выбора параметра регуляризации а приведены, например, в [1, 3, 6, 7], в настоящей работе они подробно не рассматриваются. Ниже рассмотрение ведется при одном значении а (выбранном некоторым способом).

Быстрые алгоритмы вычисления решения (5). Запишем уравнение (5) в виде

1 да

WayG) = — J ба (ш) Gy(ш) ed® , (10)

-го

где

*

еа («) = . н,(ш) (ii)

|Я (ш)|

2 2 v + аш2 p

Как правило, уравнения (5), (6), (10), (11) в дискретном виде решаются с использованием алгоритма БПФ [1], считающегося самым быстрым. Оценим, сколько требуется операций умножения/деления для вычисления решения (6), (10), (11) в дискретном виде с использованием БПФ (операции сложения и сопряжения не учитываем). Полагаем, что n — целая степень числа 2.

1. Один раз для всех строк вычисляем: H(ш) (потребуется n log2 n операций), | H(ш) | (n операций), аш p (2n операций при p = 1) и

H *(ш)/| H (ш)|2 +аш p (n операций) — всего

n log2 n + 4n операций для вычисления ба (ш) .

2. На вычисление Gy (ш) при некоторомy потребуется n log2 n операций.

3. На умножение ба(ш) на Gy (ш) при некоторомy потребуется n операций.

4. На вычисление w^ (£) при некоторомy потребуется nlog2 n операций.

5. В результате на вычисление w(X (£) для всех m строк потребуется n(log2 n + 4) + +n(2 log2 n + 1)m « 2mn log2 n операций.

Если m = n , то нужно « 2n2 log2 n операций для вычисления w(X (£) . Таким образом, для H(ш) , ба (ш), Gy (ш), ба(ш) Gy (ш) и way (£) требуется по n ячеек (всего 5n) памяти, при некотором y, а для всех m строк потребуется 2n + 3mn « 3mn ячеек памяти. Если m = n, то нужно 3n2 ячеек.

Предлагаем другой быстрый алгоритм, основанный также на методе ПФ с регуляризацией (см. (5), (6)). Изменив порядок интегрирования (ср. [6, с. 261]); решение (5) можно представить в виде

Way Ф = J

JL J_H (ш)_e-*»«-*) dш

2п j |h(ш)|2 +аш2p

—да

gv(x) dx (12)

или

где

way

£) = J Яа (£-X) gy (X) dx , (13)

1 ^

Ка (О = — I & (®)е"®^® . (14)

-го

Суть предлагаемого алгоритма состоит в том, что заранее может быть вычислен „заготовленный" вектор Ка (строка теплицевой матрицы при дискретизации) и восстановление

изображения сведется к вычислению way (£) для каждой у-строки путем умножения Ra на gy. Дискретизируем задачу (12)—(14). Для этого введем вдоль каждой у-й строки дискретные целочисленные отсчеты по пикселам = xi = 1,2,..., п . Заменив в (13) интеграл конечной суммой, получим:

п

way ) = Е ^ - ч ) gy (xk ). (15)

k=1

2

Вычисление (15) при некотором у потребует п операций умножения, а для вычисления всего изображения потребуется mп . При этом не учитываем затраты на вычисление Ra, так

как эти вычисления могут быть сделаны заранее. Если m = п , то потребуется п3 умножений, а предыдущий алгоритм (10) (с использованием БПФ) требует 2п п , т.е. алгоритм (15) требует больше операций умножения, чем алгоритм (10).

Что касается требуемой памяти, то алгоритм (15) требует 2п ячеек для Ra, а также по

2

mп ячеек для g и для wa, т.е. всего « 2mп ячеек. Если m = п , то требуется 2п ячеек, т.е.

чуть меньше, чем для алгоритма (10).

Хотя алгоритм (15) медленнее алгоритма (10), тем не менее, например, в системе Ма1:-ЬаЬ чрезвычайно быстро выполняется операция умножения матрицы Ra на вектор gy. Кроме

того, в алгоритме БПФ обычно не учитывают время на вычисление экспонент („поворачивающих коэффициентов"). Если это учитывать, а алгоритм (15) реализовывать на „матричном" языке или системе (типа МаЛаЬ), то эти два алгоритма могут сблизиться по скорости. Кроме того, алгоритм (15) требует меньше памяти и он проще алгоритма (10).

Рассмотрим также алгоритмы, основанные на соотношениях (3), (4), (7)—(9). Если ре-

T 3

шать СЛАУ (8), то на однократное вычисление матрицы ae + a a потребуется « п опера-

T 2

ций умножения, на вычисление a gy — п , а на решение СЛАУ (8) методом Холецкого (он наиболее подходит для решения СЛАУ с положительно определенной и симметричной матрицей

aE + aa [8,

с. 22]) потребуется « п . В результате на получение всего изображения

3 2 2 2

wa потребуется п + m(п + п ) = (2m + п)п операций.

3 T 2

Если m = п , то потребуется 3п . Для матрицы aE + a a требуется п ячеек памяти,

для вектора a gy — п, а для решения waу — п, всего п(2m + п) ячеек. Если m = п , то потребуется 3п2 ячеек памяти.

Если искомое решение получать согласно (9), то на вычисление матрицы aE + a a по-

3 T —1

требуется п операций умножения, на вычисление обратной матрицы (aE + a a) методом

3 / T —1 T 3

Холецкого [8, с. 20—22] — п /2, на умножение

{ож + AJ A)

на a — п , т.е. на вычисление матрицы

ba= (ож + ATA)—1at (16)

потребуется 2п3 + п3/2 = 2,5 п3. Далее, на умножение (пхп)-матрицы Ba на gy потребуется

22 п операций, а на все строки — mп . В результате на вычисление wa согласно (9) потребу-

3 2 2 3

ется 2, 5п + mп = (ш + 2, 5п)п операций, если m = п , то — 3,5п .

T — 1

Требуемая память: если матрицу (аЕ + А Л)'1 сохранять в тех же ячейках памяти

Т 2

аЕ + Л Л, то потребуется п(2т + п) ячеек; если т = п, то — 3п .

Представим второй алгоритм, основанный на соотношении (9), записанном в виде (для каждой у-й строки изображения, ср. [5, с. 137, 6, с. 248, 251]):

ж

ау

= Ва Е

а о у >

(17)

Характерная особенность решения (17) состоит в том, что матрица Ва может быть рассчитана заранее (несущественно, каким методом), и восстановление изображения сведется к умножению „заготовленной" матрицы Ва на каждую у-ю строку изображения е. Для этого с

учетом всех т строк изображения потребуется тп2 операций. Требуемая память: для матри-

22 цы Ва потребуется п ячеек, для е и для жа — по тп, т.е. всего п + 2тп = п(2т + п). Если

т = п , то потребуется 3п2 ячеек.

Видно, что алгоритм (17) требует такого же объема памяти, как и (8) и (9), но он требует в 3 и 3,5 раза соответственно меньше операций умножения, т.е. в 3—3,5 раза быстрее, чем (8)—(9). Кроме того, он проще в программировании.

Техническая реализация алгоритмов. Если, например, т = п = 400, то алгоритм (17)

потребует п = 64 млн операций умножения. Если скорость компьютера ~ 1 млрд оп/с, то восстановление изображения займет ~ 0,1 с. Можно время реализации алгоритма еще сократить, если использовать параллельные вычисления по строкам.

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

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

и высоту полета самолета и по этим данным служба обработки изображения сможет оценить А и выбрать из памяти несколько матриц Ва с некоторыми типичными значениями А и а.

В качестве примера на рис. 1 представлено смазанное изображение самолета (т = 510, п = 640, А^ 20). На несколько экранов могут быть выданы несколько восстановленных изображений жа, и опытные

диспетчеры выберут значения А и а (рис. 2)

Видно, что при А = 20, а = 10-4 получается вполне удовлетворительное восстановление изображения самолета: из-за смаза на самолете не были видны опознавательные знаки (рис. 1), а после восстановления алгоритмом „заготовленной" матрицы, согласно (17), на самолете стали видны опознавательные знаки. Восстановление одного кадра на рис. 2 потребовало < 1 с машинного времени.

Рис. 1

Рис. 2

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

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

список литературы

1. Арефьева М. В., Сысоев А. Ф. Быстрые регуляризирующие алгоритмы цифрового восстановления изображений // Вычислит. методы и программирование. 1983. Вып. 39. С. 40—55.

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

3. Грузман И. С., Киричук В. С., Косых В. П., Перетягин Г. И., Спектор А. А. Цифровая обработка изображений в информационных системах. Новосибирск: Изд-во НГТУ, 2002. 352 с.

4. Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера, 2006. 1072 с.

5. Сизиков В. С. Обратные прикладные задачи и Ма1ЬаЬ. СПб: Лань, 2011. 256 с.

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

7. Воскобойников Ю. Е., Мухина И. Н. Локальный регуляризирующий алгоритм восстановления контрастных сигналов и изображений // Автометрия. 2000. № 3. С. 45—53.

8. Уилкинсон, Райнш. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра. М.: Машиностроение, 1976. 389 с.

Валерий Сергеевич Сизиков

Константин Александрович Кирьянов

Роман Алексеевич Экземпляров

Сведения об авторах

— д-р техн. наук, профессор; Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, кафедра измерительных технологий и компьютерной томографии; E-mail: sizikov2000@mail.ru

— аспирант; Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, кафедра измерительных технологий и компьютерной томографии; E-mail: kiryancon@front.ru

— аспирант; Санкт-Петербургский государственный политехнический университет; E-mail: rexe@yandex.ru

Рекомендована кафедрой измерительных технологий и компьютерной томографии

Поступила в редакцию 13.04.13 г.

УДК 53.082.5

Г. Д. Фефилов

ГЕОМЕТРИЧЕСКОЕ ПРЕДСТАВЛЕНИЕ СИГНАЛА, ОСНОВАННОЕ НА ПОНЯТИИ О ФАЗОВОМ ПРОСТРАНСТВЕ, В ЛАЗЕРНОЙ ДИФРАКТОМЕТРИИ МИКРООБЪЕКТОВ

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

Ключевые слова: лазерная дифрактометрия микрообъектов, информативный параметр сигнала.

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

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