Устойчивые методы восстановления зашумленных изображений

УДК 517.983

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

Т.И. Сережникова

STABLE METHODS FOR RECONSTRUCTION OF NOISY IMAGES

T.I. Serezhnikova

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

Ключевые слова: численные методы, алгоритмы, решение некорректных задач, обратные задачи, итеративная регуляризация, негладкая оптимизация.

We consider the numerical reconstruction of noisy images by using two regularization algorithms. The basis of algorithms is Tikhonov regularization with two special nondifferentiable stabilizers. To solve the problem of a nonsmooth minimizing the proximal method and the subgradient process are involved. The results of calculations on the supercomputer «Uranus» are presented.

Keywords: numerical methods, algorithms, ill-posed problems, inverse problems, iterative regularization, nonsmooth optimization.

Введение

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

Известно, что можно выделить группу методов, развитых советскими (российскими) учеными [5 - 10] - методы регуляризации Тихонова, итеративной, статистической, локальной, дискриптивной регуляризации, субоптимальной фильтрации, решения на компакте и др., а также - группу методов, развитых западными учеными [11 - 17] - методы оптимальной фильтрации Калмана - Бьюси и Винера, методы управляемой линейной фильтрации (Бей-куса - Гильберта) и др., которые являются в принципе более точными. Но методы первой группы (например, метод регуляризации Тихонова) требуют гораздо меньше дополнительной информации о решении и поэтому находят более широкое применение при решении обратных прикладных задач.

В настоящей работе рассматривается широко используемая постановка задачи восстановления изображения в виде решения двумерного интегрального уравнения типа свертки. В разделах 3 и 4 приведены результаты численного решения двух задач о восстановлении изображений на основе использования разрабатываемых численных алгоритмов для решения двумерного интегрального уравнения Фредгольма 1-го рода.

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

Известно, что изображение объекта, полученное с помощью оптического прибора, как правило, искажено вследствие дефокусировки, из-за влияния среды между объектом и прибором и в силу ряда других причин. В связи с этим, возникает задача реконструкции изображения. Постановка задачи восстановления изображения достаточно проста и наглядна. В то же время, решение этой задачи требует использования довольно разнообразных и сложных инструментов, в том числе математического моделирования, статистической теории оценивания, преобразования Фурье, методов оптимизации, численных методов линейной алгебры и др. В настоящее время в научной литературе можно найти решения нескольких задач реконструкции изображения. Например, в работе [8] приводятся решения трех задач обработки изображений: восстановление истинного изображения - характеристической функции перевернутых букв - по дефокусированному с помощью аппаратной функции и зашумленному изображению; решение тестовой задачи (хорошо известной в теории и практике обработки изображений) о восстановлении истинного изображения - фотографии искусственного спутника Земли; обработка изображений гравитационных линз.

Задача восстановления изображения сводится к решению двумерного интегрального уравнения типа свертки:

1 1

Аи = J ! К

о о

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

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

2. О регуляризирующих алгоритмах

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

Аи = /

на паре линейных нормированных пространств и, Р с разрывным, в общем случае, многозначным оператором А-1 и приближенными данными /<5, || / — /гЦ < 5 .

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

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

Для того, чтобы охватить случай разрывного решения и построить регуляризирующий алгоритм, в работе [1] предложен и теоретически исследован стабилизатор с нормой пространства Соболева И'^(П) с дробными производными (многомерный случай), и установлена сходимость приближенных решений в сильной топологии этого пространства и дискретная сходимость конечномерных аппроксимаций. В работе [2] был предложен другой параметрический класс недифференцируемых стабилизирующих функционалов на основе нормы пространства Липшица и обоснована равномерная сходимость тихоновских аппроксимаций к непрерывному необязательно дифференцируемому решению исходного уравнения. В данной работе для случая линейного непрерывного в С (П) оператора А приведены ниже теоремы 1 и 2, доказанные в [2] для более общего случая.

Пусть П — компакт в В,п . Определим множество функций, удовлетворяющих на П условию Липшица

3 й > 0 : |и(х\) — и(хг)| < й \х\ — х2\^ V 11,12 6 П,

и

где 0 < Ц < 1, \х\ — Х21 = \ХИ — х2'1\2)1^2  Обозначим через Н^ =

г=1

функций, удовлетворяющих условию (1.1) с фиксированным параметром ц,. По аналогии со случаем функции одной переменной определяется норма

и{х)\\нр = тах |гі(ж)| + вир

\и(хг) - и(х2)|

хЄП

,х,ЄП \Х!~Х2Ї

множество

При приближенно заданной правой части /г, 11 / — /г 11 < 5, для построения регуляризо-ванного семейства приближенных решений используется метод Тихонова в виде

Теорема 1. Пусть А-линейный непрерывный оператор в С(П), Кег(А) = {0}, ||/—/г|| < 5. Тогда задача (4) имеет решение иа, возможно неединственное, и при связи параметров (х(д) —> 0, ^ 0, 5 ^ 0 имеет место равномерная сходимость

Пт \\иа(-6) -й||с(п) = 0,

кроме того,

Замечание 1. Если вместо липшицевой нормы в (4) использовать стабилизатор

1,

П(и) = ЦгіЦям + — ||гі|, 1 < р < аз,

2

то решение иа задачи (4) единственно.

XI

Используя метод Рптца, можно построить аппроксимации регуляризованного приближенного решения иа элементами конечномерного пространства кусочно-линейных функций. Предположим, что конмакт П допускает разбиение симплексами n-го порядка с максимальным диаметром h = m&xdiam(Arn). Обозначим Uh С Н^ - конечномерное подпространство

ГП

кусочно-линейных (т.е. линейных на каждом симплексе) функций.

Теорема 2. Пусть выполнены предположения теоремы 1. Тогда конечномерная задача шт{\\Аи - /г||с[п] + «[1КН|Р[П] + IMI нА  uh G Uh},

имеет единственное решение uft, для которого выполнены соотношения

1^11 uh “ Л1с[п] = 0, lim ||г1^||ям[п] = \\иа\\н^[п],

где иа - решение, полученное методом Тихонова со стабилизатором (5).

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

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

ик = агдтт{Фа

где (3 > 0, || • ||я ~ некоторая гильбертова норма.

В отличие от выпуклого функционала Ф"(и), функционал = Ф"(и) + /3 \\и — г;|\2Н

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

Важно отметить, что итерации ик ргох-метода (7) могут вычисляться субградиентными методами с некоторой погрешностью с сохранением сходимости. А именно, как следует из результатов [5, 18], справедлива теорема 3.

Теорема 3. Пусть Т - ргох-отображение, т.е.

Т : v —> argrriin {Ф^

ОО

\\zk+1 ~ T(zk) || < 7к, z° G RN, < oo.

i=0

Тогда lim \\zk — = 0, где uN - решение дискретного аналога регуляризованной задачи

»оО

(6).

В работе [4] для рассматриваемого алгоритма, использующего норму Липшица и ргох-метод, были проведены расчеты по восстановлению трех модельных решений интегрального уравнения Фредгольма первого рода

2

/Ц

—^и(х)(1з = /(у), 0 < Н < 2, 0 < х,у < 2,

о

которое моделирует ситуацию, когда необходимо найти продолженное на глубину х = Н гравитационное поле и(х) по заданному на земной поверхности /(у) (плоский случай).

Для каждой модели рассчитывался вариант с точной правой частью и вариант с возмущенной правой частью.

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

Графики точных решений и восстановленных решений как для точных, так и для возмущенных данных для трех одномерных модельных решений из работы [4] приведены ниже на рис. 1. На каждом графике по горизонтальной оси расставлены номера узлов г сетки х-1 = (2/128) х г (г = 0, 2,..., 128) с шагом 20, начиная с 0.

Рис. 1. Графики точных и восстановленных решений для трех моделей

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

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

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

Пусть А — линейный ограниченный оператор, действующий из Рр{В) в Ьд(8), где 1 < р, д < (50, И С Мт, 5 С К4" — произвольные области с кусочно-гладкими границами.

Пусть уравнение (1) разрешимо в пространстве и = {и : и £ Ьр(0) : 7(и) < га} для точных данных А, /, которые заданы своими приближениями А^, /<5.

Рассматривается метод Тихонова в форме

где Л^и) - полная вариация, определенная формулой:

J

о

Теорема 4. Для любой пары А^ /г при условиях ||А — Аь\\ < Ь, ||/ — /$\\ < 5, для любых а > 0 и и0 £ и задача (9) имеет единственное решение иа, и при связи параметров а(5, К) —0; {5 + к)д/сх(5, /г.) —^ 0 имеет место сходимость

Нт \\иа(-6^ -й\\т = 0, Нт - и0) = Ли- и0),

где и — и0 — нормальное решение уравнения (2) относительно стабилизирующего функционала О,(и) = ||г4||^ + Ли)(т.е. решение, реализующее среди всех решений минимум величины функционала П(и — и0)).

Теорема сохраняет силу, если пространство и заменить на выпуклое замкнутое множество к а и.

Регуляризация - необходимый этап при построении устойчивого алгоритма решения некорректно поставленной задачи. Для практического нахождения регуляризованного решения требуется второй этап - подходящий итерационный метод решения задачи минимизации со строго выпуклым целевым функционалом О,. В данном случае исследование субдифференцируемости функционала полной вариации 7(и) оказалось крайне затруднительным. Поэтому основная идея заключается в том, чтобы строить минимизирующую последовательность не в пространстве Ьр(И) , а в более гладком пространстве ЩР^(И), в частности, Щ2,(0), а также на подмножестве С2(И) пространства Ь2(^). Тогда для функций

и Е Ж1 (О) функционал полной вариации 7(и) принимает вид, который обеспечивает его субдифференцируемость на Ж1 (О) и позволяет использовать технику гильбертова пространства при доказательстве сходимости целого класса методов субградиентного типа.

Принципиальная возможность построения в задаче (9) минимизирующей последовательности из более гладких функций и Е Ж^(О) вытекает из следующей леммы, доказанной в [3].

Лемма 1. Для любой функции конечной полной вариации и Е и = {и : V Е ЬР(И), 7(и) < га} существует последовательность функций ио Е С“(О), для которой

Пт ||ил — и|| г =0 (р > 1), Пт 7(и,-) = 7(и).

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

Пусть 7) — т-мерная область прямоугольных очертаний, для определенности, например, единичный куб. Вводится равномерная сетка с шагом к = 1/п по каждой переменной, рассматривается сеточный аналог М™ пространства М™

Ж™ = {ж е Мт : же (л/г,,.. л, = 0, ±1, ±2,...}

и сеточные функции ип : М™ —> М, здесь М™ — пт-мерное пространство векторов ип ;

индекс п означает, что функция задана на сетке с шагом к = 1/п. В соответствии с предложенной методикой для нахождения приближенного решения уравнения (1) на первом этапе проводится вариационная регуляризация в форме (9) при р = 2 . На втором этапе после дискретизации регуляризованная задача решается субградиентным методом.

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

3. Численные эксперименты для решения задачи восстановления изображения характеристической функции объекта

Рассмотрим решение задачи о восстановлении характеристической функции для множества точек прямоугольной области на плоскости. Для решения поставленной задачи в качестве стабилизатора используется норма Липшица, как это изложено в разделе 2. Опишем необходимые изменения для двумерного случая.

Для дискретной аппроксимации регуляризованной задачи в двумерном случае в единичном квадрате О строится равномерная сетка с шагом к = 1/п по каждой переменной. Упорядочим N = п2 узлов двумерной сетки естественным образом по столбцам, а внутри столбца - по строкам, и пронумеруем с помощью j. Получим: {х3 ,у3} - координаты узла с номером j, и^ = и(хЭ, у^); К1- = К(уЭ — г)г, ж-5, — £*), = 1, ...,п2. Тогда дискретный аналог

регуляризованной задачи в двумерном случае принимает вид

N

тах

1 <І<N

\^Ь2К)и3-П\ + 0=1

~\~СХ

N

> ЫиА2 + тах + тах

0=1

1

2

Как показано в [2], кусочно-линейные восполнения вектора , реализующего минимум в (11), равномерно сходятся к решению регуляризованной задачи (4) со стабилизатором (5). Дискретный аналог ргох-метода принимает форму

ик = агдтгп {Ф^{и) +/3 \\и — ик~г\\2яМ : и Є RN}.

Так как целевая функция является сильно выпуклой, то применим субградиентный

метод вида

ик'и

ик^+1 = ик,и _ Хк___^ г/ = 0,1, 2,и*;,

где ьк'и Є дФ^к(ик'и), Ф^д.—целевая функция в задаче (12), дФ— субдифференциал.

Напомним, что установлено, как это указано в разделе 2, сходимость по аргументу, т.е. Пт 11ик'и — гіЛ' ||=0.

1/—ЮО

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

1СГ3

Рис. 2. Восстановление структуры поверхности характеристической функции: а) точная поверхность; б) наблюдаемая; в) восстановленная

Отметим, что в работе другого автора [8] приведены похожие результаты восстановления истинного изображения - характеристической функции двух букв (положенных на плоскость) с помощью регуляризующих алгоритмов, основанных на двумерном преобразовании Фурье и на использовании в одном из алгоритмов гладких функций из пространства Соболева, а в другом - класса функций двух переменных с ограниченными вариациями.

4. Задача восстановления изображения искусственного спутника земли

Рассматривается известная тестовая задача о восстановлении изображения искусственного спутника земли. Численные решения этой задачи можно найти в работе [8], в которой изображения реконструировались по методу регуляризации в пространствах Соболева с использованием преобразования Фурье и на трех классах функций с ограниченными вариациями. В работе [8] отмечается, что полученное качество решения указывает на необходимость конструирования и других регуляризирующих алгоритмов.

В работе [3] проведена численная апробация алгоритма на основе пространства функций с ограниченной полной вариацией, представленного в разделе 2. Результат восстановления изображения спутника на сетке 128x128 узлов воспроизведен на рис. 3.

Рис. 3. Реконструкция изображения на классе функций с ограниченной полной вариацией: а) наблюдаемое; б) восстановленное

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

2, основанного на использовании пространства непрерывных функций и введении нормы Липшица. В разделе 2 приведены примеры успешной апробации алгоритма для решения одномерных интегральных уравнений Фредгольма 1-го рода.

Дискретный аналог регуляризованной задачи для решения двумерного интегрального уравнения (1) уже представлен уравнениями (11), (12), (13) в разделе 3.

Результаты проведенных расчетов с использованием алгоритма для двумерного случая изображены на рис. 4 а, б, в.

Величина параметра регуляризации а из формулы

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

5. Заключение

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

Рис. 4. Реконструкция изображения в пространстве функций с нормой Липшица: а) точное изображение; б) наблюдаемое; в) восстановленное

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

Работа поддержана Российским фондом фундаментальных исследований, грант №09-01-00053.

Литература

1. Vasin, V.V. Regularization and iterative approximation for linear ill-posed problems in the space of functions of bounded variation / V.V. Vasin // Proc. Steclov Inst. Math. Supl. -

2. Васин, В.В. Аппроксимация негладких решений линейных некорректных задач /

3. Васин, В.В. Двухэтапный метод аппроксимации негладких решений и восстановление зашумленного изображения / В.В. Васин, Т.И. Сережникова // Автоматика и телеме-

4. Васин, В.В. Регулярный алгоритм аппроксимации негладких решений для интегральных уравнений Фредгольма первого рода / В.В. Васин, Т.И. Сережникова // Вычисли-

5. Васин, В.В. Проксимальный алгоритм с проектированием в задачах выпуклого промеханики УНЦ АН СССР).

6. Тихонов, А.Н. Методы решения некорректных задач / А.Н. Тихонов, В.Я. Арсенин. -М.: Наука, 1976.

7. Сизиков, В.С. Математические методы обработки результатов измерений / В.С. Сизи-

8. Леонов, А.С. Решение некорректно поставленных обратных задач: очерк теории, практические алгоритмы и демонстрации в МАТЛАБ / А.С. Леонов. - М.: Книж. дом «Либ-роком>, 2010. - 336 с.

9. Восстановление изображения гравитационной QSO линзы 2237+0305 «Крест Эйнштейна:» / В.А. Белокуров, Е.В. Шимановская, М.В. Сажин и др. // Астроном, журнал. -

10. Бакушинский, А.Б. Некоторые нестандартные регуляризующие алгоритмы и их численная реализация / А.Б. Бакушинский, В.С. Сизиков // Журн. вычисл. матем. и матем.

11. Vogel, C.R. Computational methods for inverse problems / C.R. Vogel. - Philadelphia: SIAM, 2002.

12. Физика визуализации изображений в медицине: в 2-х т.: пер. с англ. / под ред. С. Уэбба.

13. Чейссон, Э.Дж. Первые результаты с космического телескопа «Хаббл> / Э.Дж. Чейссон

14. Эндрюс, Г. Применение вычислительных машин для обработки изображений / Г. Энд-

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

16. Kawanaka, A. Estimation of static magnetic field and gradient fields from NMR image /

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

18. Rockafellar, R.T. Monotone operators and the proximal point algorithm / R.T. Rockafellar

Татьяна Ивановна Серeжникова, кандидат физико-математических наук, доцент, отдел «Некорректные задачи анализа и приложений:», Институт математики и механики УрО РАН, sti@imm.uran.ru.

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