CC BY
66
В. А. Фурсов
ВОССТАНОВЛЕНИЕ ИЗОБРАЖЕНИЙ КИХ-ФИЛЬТРАМИ, ПОСТРОЕННЫМИ ПУТЕМ НЕПОСРЕДСТВЕННОЙ ИДЕНТИФИКАЦИИ ИНВЕРСНОГО ТРАКТА
Аннотация
Решается задача восстановления изображений, подвергшихся линейным искажениям типа дефокусировки или смаза по результатам непосредственной идентификации параметров инверсного восстанавливающего КИХ-фильтра. Показана принципиальная возможность достижения различного регуляризующего эффекта путем простого линейного преобразования данных. Предложена итерационная процедура, использующая набор моделей невысокого порядка.
Введение
Проблема восстановления изображений обычно сводится к задаче построения обратного оператора искажающей системы. При этом предполагается, что искажающий оператор Щ(.), действующий на входную последовательность х(п1, п2) каким-то
образом задан. Часто оператор Щ(.) определяют путем предварительного решения задачи идентификации характеристик искажающей системы по специально подобранным для этой цели тестовым (например мирам) изображениям.
Если задачи идентификации и восстановления решаются изолированно, то усилия, направленные на определение, по-возможности, более точной модели искажающей системы на этапе идентификации могут оказаться в некоторой степени напрасными. Дело в том, что на этапе восстановления изображений, обычно, все равно приходится в той или иной степени (в зависимости от свойств оператора) осуществлять регуляризацию решений, по существу означающую намеренное искажение обратного оператора с целью уменьшения чувствительности к ошибкам в исходных данных. Поэтому представляется целесообразным задачи идентификации моделей искажающих систем и восстановления изображений с самого начала рассматривать с учетом их взаимосвязи.
1. Постановка задачи
Ограничимся рассмотрением линейных искажений типа дефокусировки. В качестве модели искажающей системы будем использовать, обладающий конечной импульсной характеристикой, КИХ-фильтр с заданной опорной областью Б:
у(щ, «2) = X к(к1, к2) Х(П1 - к1, п2 - Ю +£(«1, n2), (1)
{к1,кг }еП
где к(к1, к2), х(п1 - к1, п2 - к2) - значения импульсной характеристики и входных отсчетов в точках опорной области Б, а £(п1, п2)- дискретная последовательность, включающая погрешности измерений и ошибки аппроксимации.
Из М различных отсчетов у(п1, п2) изображения, для которых существуют опорные области Б(.) составим вектор у размерности Мх1. Если в каждой из этих опорных областей Б(п1, п2) число отсчетов
одинаково и равно N, а соответствующий N^1 вектор к параметров импульсной характеристики инва-
риантен к сдвигу, можно записать матричное соотношение
У = Хк + £, (2)
где X - матрица MхN, каждая строка которой составлена из N отсчетов изображения в соответствующей области Б, а Мх1 - вектор £ составлен из соответствующих компонентам у(п1, п2) в (1) ошибок £(п1, п2).
По аналогии с (1) модель инверсного тракта, если она существует, можно представить в виде х(«1;п2) = X кг (к1,к2)у(п1 -к1,п2 -к2) + ц(п1,п2), (3)
{к1,к2}еК
где ц(п1, п2) - ошибки, связанные с наличием ошибок £(п1, п2) в исходных данных и ограничением порядка КИХ-фильтра. Для Ь опорных областей, каждая из которых включает Б отсчетов выходного (искаженного) изображения, по аналогии с (2), можно записать матричное соотношение
х = Укг + ц , (4)
где х, ц- векторы Ьх1, У -матрица ЬхБ, а кг - вектор Бх1.
Восстанавливающий фильтр также строится в классе КИХ- фильтров:
Х(п1,п2) = X кг(кх,к2)у(пх -кх,п2 -к2) (5)
{к1 ,к2 }еЯ
или в матричном виде
Х = Укг, (6)
где У - фигурирующая в (4) матрица, кг - вектор
размерности Бх1 оценок коэффициентов разностного уравнения, описывающего инверсный фильтр. Компоненты вектора Х размерности Ьх1 - суть оценки значений поля яркости восстановленного изображения. Они содержат ошибки, связанные, как с наличием в (3) ошибок ц/(пх, п2), так и с погрешностями оценивания вектора параметров кг.
Для построения вектора кг есть, по крайней мере, две возможности. Во-первых, можно по вектору у и матрице X, фигурирующим в (2), найти оценку к вектора к разностной модели искажающей системы, а затем с ее использованием одним из известных способов построить вектор кг - восстанавливающего фильтра. Во-вторых, можно (поменяв изображения местами) по х и У сразу найти
оценку кг - параметров инверсного фильтра, решив соответствующую задачу идентификации.
В настоящей работе рассматривается второй подход, основанный на непосредственной идентификации характеристик инверсного тракта. Этот подход исследуется с целью, с одной стороны, показать, что можно обойти трудности, связанные с реализацией оптимального "винеровского" фильтра, с другой стороны, изучить возможность достижения различного качества восстановления изображений (в т.ч. и по субъективным критериям) путем простой модификации процедуры оценивания сразу на этапе идентификации модели инверсного тракта. При этом, естественно, возникают другие проблемы. В работе выявляются эти специфические трудности и обсуждаются пути их преодоления. Развивается качественная теория идентификации, опираясь на выводы которой строится итерационная схема восстановления изображений. Для этого используется набор моделей полученных на этапе изучения путем последовательной идентификации моделей "фиктивных" инверсных трактов, в которых в качестве искаженных используются промежуточные восстановленные тестовые изображения.
2. Обсуждение проблемы
В соответствии с (2), МНК-оценка к вектора параметров к модели искажающей системы имеет вид [1]
к = [ Г ХТУ.
(7)
Алгоритм для получения этой оценки с автоматическим отбором наиболее информативных фрагментов изображений и устойчивый к аномальным ошибкам в наблюдениях рассмотрен в работе [2]. Если оцениваются непосредственно параметры инверсного фильтра, МНК-оценка кг вектора параметров кг инверсного тракта определяется как
кг = ]-1 уТх.
(8)
Важной особенностью задачи (8) является то, что матрица У в отличие от матрицы X составлена из искаженных и зашумленных отсчетов изображения. Поэтому соответствующая матрица [гТ7] может оказаться невырожденной даже в ситуациях, когда исходное неискаженное и незашумленное изображение не содержит никаких информативных элементов (отсутствие объектов на фоне постоянной яркости). Результаты оценивания параметров инверсного тракта на таких фрагментах изображения могут ввести в заблуждение.
Вторая важная особенность заключается в следующем. Для задачи (7) естественным является требование достижения наивысшей точности определения к в смысле близости к к. Для этого наряду с
оценками МНК, часто, строят также, например, непараметрические оценки, обладающие устойчивостью к ошибкам в исходных данных. В данном случае (8) является наилучшей в смысле критерия :
кг : д(кг) = шшб(к,) = шп([ - ¥кг ]Т [х - ¥кг ])
оценкой вектора параметров разностной модели инверсного тракта, а поскольку кг ищется в условиях присутствия в элементах матрицы У ошибок, то эта оценка фактически оказывается регуляризованной в среднеквадратическом смысле. Построение в данном случае другой, может быть более "хорошей" (в смысле близости к кг) оценки, может привести к большим ошибкам
оценивания вектора х при реализации процедуры восстановления (5),(6). Поэтому представляет интерес установление общих закономерностей формирования оценок, обладающих различной степенью близости к их истинным значениям.
3. Теоретические результаты
Рассмотрим вначале случай решения задачи (7) с использованием модели (2). Для оценки степени близости оценок к к их истинным значениям к будем использовать скалярную величину [2]
|Дк| Е = |к - к!' =#ТХ[ХТХ]-2 ХТ#, (9)
означает
величиной
где Дк = [ХТХ]-1 ХТ% , а [ХТХ]-[ ХТХ ]-1 [ ХТХ ]-1. Аналогичной ]
\дк\Е = \к - С = ^ XX [XТ XX]-2 XT| (10)
будем характеризовать точность при использовании набора данных, полученных из исходных путем линейного преобразования:
у = Оу, Х = ОХ, # = в£, (11)
где О вещественная невырожденная ограниченная матрица. Ясно, что при таком преобразовании равенство в (2) сохраняется, а Дк = [ХТХ ] .
Теперь поставим следующий вопрос: существует ли матрица О, для которой
Дк\ <|Дк|Е.
(12)
Указанная постановка является предметом качественной теории идентификации, сформулированной в работе [3]. Из соотношений (9),(10) ясно, что неравенство (12) выполняется при любых если имеет место неотрицательная определенность матрицы
в = х[хТ х]-2 хТ - оТ X [ХТ XX]-2 ХХТо. (13)
Введем в рассмотрение ортонормированную ихи матрицу Т =[ТЛ': Т0 ]. Матрица Тя размерно-
сти nхm составлена из нормированных собственных векторов, соответствующих ненулевым собственным значениям матрицы ХХТ , а Т0 - пх(п-ш) матрица ортонормированных собственных векторов, соответствующих нулевым собственным значениям той же матрицы. Такая матрица всегда существует, -1
причем Т = Х¥Л 2 [4], где F - mхm матрица ортогонального преобразования:
¥Т¥ = ¥¥Т = Ет , ¥ТХТХ¥ = Л ,
-1 -1 -1 -1 Л2 = diag (42Л2,...А2),
Л = diag (Л1 ,Я7,...,Лп )-диагональные матрицы составленные из собственных значений ^, i = 1, п матрицы ХТX .
Применяя к матрице В, заданной соотношением (13), ортогональное преобразование вида ТТВТ , с учетом того что ТтлХ [ХТX]-2 ХтТя = Л-1, а
Т0ТХ [ ХТХ ]-2 ХТТ0 = 0, получаем
к = T'BT =
Л-1 0 0 0
1-2
TT AT,
TTATo
TT AT, TTATo
(14)
что
(15)
где А = вТХ [ХТХ]" ХТв . Покажем,
Т1АТЛ =Л-1. Подставив вместо матриц A и Тл соответствующие им выражения, имеем
Т^АТЛ =Л 2¥ТХОТХ[ХТХ]" ХТОХ¥Л 2.
последнюю матрицу
Так как ¥T¥ = ¥¥т = E,
можно переписать в виде
-1 -1 Л2¥ТХОТОХ¥Я2¥ТХТвТОХ¥Л 2, (16)
где Я = ¥Т [ХТвТОХ]-1 ¥.
Далее, поскольку ¥Т = ¥-1 и, следовательно, ¥ = (¥Т )-1, можно показать, что
¥Т [ХТвТОХ]-1 ¥ = [¥ТХТвТ ОХ¥]-1.
С учетом последнего равенства из (16) непосредственно следует (15). Следовательно
K =
-T0TCTFЛ 2
-Л 2 ¥T CT0
TTCTCTÜ
(17)
где
C = хХ [хХTхХ]-1 хХTG . Можно
показать, что все
блоки матрицы K становятся нулевыми при G=E.
Матрица K симметрическая, поэтому все ее собственные значения вещественны. Для нас важным является тот факт, что она не может быть ни положительно, ни отрицательно определенной, ввиду наличия нулевого блока mхm в левом верхнем углу. Это означает, что среди собственных значений этой матрицы кроме нулевых могут быть, как положительные, так и отрицательные, следовательно не
существует матрицы преобразования G, которая приводила бы к повышению точности оценок регре-сии в смысле критерия (9) при любом векторе ошибок £. Вместе с тем, в n - мерном пространстве ошибок £ для любой невырожденной ограниченной матрицы G всегда существует множество
□ = £:\Ah\E <|М|E } . (18)
Другими словами, при конкретной реализации вектора £ всегда можно подобрать матрицу линейного преобразования (11) так, что при этом повышается (в худшем случае сохраняется достигнутая) точность оценивания.
Можно показать, что аналогичные выводы справедливы и для задачи (8). Для этого матрицу Y необходимо представить в виде суммы Y = Y + SY , где Y - матрица, составленная из отсчетов искаженного, но не зашумленного изображения, а SY из соответствующей помеховой последовательности. При этом в соответствии с (4) х = Y*hr + п , где П = ц + 8Yhr. Нетрудно заметить, что все приведенные выше соотношения и выводы остаются справедливыми, если сделать замены y на x, X на Y и £ на п.
4. Построение процедур идентификации и восстановления
Опираясь на выявленные выше закономерности, можно сформулировать общий подход к построению восстанавливающих фильтров, основанный на линейных преобразованиях (11). Важной особенностью такого подхода является возможность достижения различного регуляризующего эффекта простым изменением способа формирования матрицы преобразования G. В частности, в настоящей работе строится итерационная процедура восстановления изображений с использованием набора инверсных КИХ-фильтров невысокого порядка, построенных заранее на этапе идентификации искажающей системы по тестовым изображениям.
Общая схема формирования совокупности моделей следующая. Вначале по исходным тестовым изображениям в результате решения задачи идентификации строится модель инверсного тракта и с помощью полученного разностного соотношения осуществляется обработка искаженного тестового изображения. Затем исходное искаженное тестовое изображение заменяется обработанным и вновь решается задача идентификации. В результате определяется новый КИХ-фильтр невысокого порядка и т.д. Процесс формирования множества моделей останавливается, если очередная модель не дает существенного улучшения качества восстановления тестового изображения.
Указанный набор моделей строится так, чтобы степень близости оценок к их истинным значениям, в смысле критерия (9) изменялась от итерации к итерации. В частности, при определении параметров
моделей, предназначенных для реализации на завершающих этапах восстановления изображения необходимо, чтобы критерий был менее регуляри-зующим, т.е. обеспечивал получение более точной инверсной модели, возможно, игнорируя присутствующие в наблюдениях помехи.
Для придания КИХ-фильтру указанных свойств используется диагональная матрица преобразования О, процедура построения которой описана в работе [2]. Эта процедура может быть выведена из семейства критериев вида
Q(K) = Шп Q(к^■)=шп X "=1 Iе' г, (19)
для 0<г<2, где е - '-я компонента вектора е = х - Укг. Соответствующие этим критериям оценки кг вектора параметров кг называются Ьг-оценками [1]. В частном случае при г = 2 оценки совпадают с приведенными выше оценками МНК. При г = 1 минимизируется сумма абсолютных отклонений (метод наименьших модулей - МНМ). В настоящей работе исследовалась процедура, которая предусматривает построение набора из трех моделей инверсного тракта с помощью последовательности трех Ьу - оценок с изменяющимися от шага к шагу значениями параметра г:г = 2, г = 1, г ^ 0.
Задача восстановления решается путем последовательной обработки реальных искаженных изображений полученным набором инверсных фильтров в той же последовательности как они были оценены. Можно показать, что если опорная область квадрат со строной N отсчетов, то каждые два соседних шага обработки формально соответствуют одному шагу обработки фильтром, имеющим опорную область со стороной 2№1.
Следовательно, для реализации указанной итерационной схемы восстановления изображений, даже при сравнительно больших размерах размазывающих масок, может применяться набор КИХ-
фильтров фиксированного, и притом невысокого, порядка.
5. Результаты экспериментов
Для формирования набора моделей инверсных фильтров использовалась тестовые изображения 128x128 с диапазоном яркости 0-256, показанные на рис. 1: а - исходное, б - искаженное, в - искаженное и зашумленное. Линейные искажения формировались 3-х кратным "проходом" КИХ-фильтром, обладающим радиальной симметрией с опорной областью 5x5 (без угловых отсчетов), показанной на рис 2. Компоненты вектора, характеризующего импульсный отклик системы, задавались следующими: к0 = 0.359118, к1 = 0.395030, к2 = 0.190333, к3 = 0.049917, к4 = 0.005602
(корни соответствующего разностного уравнения: -0.4; -0.3; -0.2±]0.3). Аддитивный шум имитировался псевдослучайной последовательностью с дисперсией 36. По приведенным тестовым изображениям строился набор из трех инверсных восстанавливающих КИХ-фильтров с такой же опорной областью. Заметим, что в данном случае хорошая обусловленность задачи идентификации обеспечивалась специально подобранным тестовым изображением, поэтому процедура отбраковки неинформативных фрагментов [2] не использовалась. В общем случае, когда задача идентификации инверсных фильтров должна решаться по реальным изображениям, целесообразно на начальном этапе многошагового процесса формировать матрицу [ХтХ ] (по
неискаженному изображению) и с ее использованием строить вектор признаков для отбраковки неинформативных фрагментов изображения на последующих этапах идентификации моделей инверсных фильтров (уже по матрице У искаженного изображения), как это описано в работе [2].
в)
Рис. 1. Тестовые изображения: а) неискаженное б) искаженное, в) искаженное и зашумленное.
Рис. 2. Опорная область
Подлежавшие восстановлению изображения "часы" (размером 256x256) приведены на рис 3: б -искаженное и в - искаженное и зашумленное. Для искажения и зашумления использоваласть таже модель, что и для тестового изображения (различались лишь генерирующие числа псевдослучайной последовательности шумов). Рис. 3 а дает представление как выглядели "часы" до искажения и зашумления. Дисперсия разности исходного и расфокусирован-
ного изображений (сг2 =0) составляла а2^ = 465.06, а при добавлении шумов (ст? = 36) - стД. = 500.62.
На рис. 4 а и б соответственно приведены изображения, полученные путем последовательного применения к искаженному изображению (рис. 3 б) набора из трех КИХ-фильтров, построенных по методу наименьших квадратов и с использованием Ьу - оценок (г=2,1,г^0).
На рис 5. приведены аналогичные результаты по восстановлению искаженного (расфокусированного) и зашумленного изображения. Здесь отчетливо наблю-
дается эффект подчеркивания шумов инверсным фильтром, хотя контуры и элементы восстановленного изображения воспринимаются более четкими.
в)
Рис. 3. Изображение "часы": а) исходное, б) расфокусированное, в) расфокусированное и зашумленное.
а) б)
Рис. 4. Искаженное изображение "часы" после восстановления набором из 3-х КИХ-фильтров полученных: а) с использованием МНК, б) с использованием Ьу - оценок (у=2,1, у-—0)
Рис 5. Искаженное и зашумленное изображение "часы" после восстановленния набором из 3-х КИХ-фильтров полученных: а) с использованием МНК, б) с использованием Ьу - оценок (У=2,1, у-—0)
В таблице 1 для исследовавшихся методов приведены оценки параметров инверсных фильтров и дисперсии разности исходного искаженного и восстановленных изображений (в последней колонке).
Результаты, приведенные в таблице показывают, что применение Ьу - оценок, обеспечивающих большую близость оценок параметров к истинным значениям, дает выигрыш при отсутствии или низкой интенсивности помех в исход-
ных данных. При наличии заметных помех в исходных данных дисперсии разности исходного и восстановленного изображений для МНК и Ьу оценок приблизительно одинаковы. Интересно, что эта дисперсия более чем в два раза выше дисперсии разности исходного и искаженного и зашумленного изображений, хотя субъективно отдельные элементы (цифры на циферблате) восстановленного изображения воспринимаются как более четкие.
Таблица 1.
Оценки параметров восстанавливающих фильтров и дисперсии разностей восстановленных и исходных изображений.
«Л Тип № Параметры инверсных КИХ-фильтров „2 О ix
оц. шага hr 0 hn hr 2 hr3 hr 4
1 8.687516 -4.199416 -4.739308 -2.618809 3.870017 205.50
МНК 2 1.693985 -0.634953 -0.369007 -0.621824 0.931799 149.04
0 3 1.033900 -0.067391 0.013760 -0.020369 0.040100 145.36
1 8.687516 -4.199416 -4.739308 -2.618809 3.870017 205.50
L-v 2 1.830662 -0.630797 -0.749170 -0.669917 1.219222 131.75
оц. 3 1.098201 -0.080034 -0.153213 0.026276 0.108770 114.85
1 4.415751 -0.873153 -2.678178 0.747165 -0.611585 1018.90
МНК 2 0.740701 0.972306 -0.543510 0.110665 -0.280162 1047.19
36 3 0.624883 0.956389 -0.635098 -0.088430 0.142256 1100.91
1 4.415751 -0.873153 -2.678178 0.747165 -0.611585 1018.90
L-v 2 0.706808 1.019733 -0.572980 0.177245 -0.330806 1058.30
оц. 3 0.535312 1.017578 -0.616960 0.050619 0.013451 1091.90
6. Заключение
Восстанавливающие КИХ-фильтры всегда могут быть построены путем непосредственной идентификации инверсного тракта с использованием, как реальных, так и полученных моделированием, тестовых изображений. Это возможно даже в случаях, когда соответствующий ряд является расходящимся (качество восстановления при этом будет, конечно, ниже). При непосредственной идентификации инверсных КИХ-фильтров удается избежать трудностей, связанных с регуляризацией и устойчивостью.
Из приведенных результатов следует, что необходимо осторожно относиться к использованию дисперсии разности исходного и восстановленного изображений в качестве критерия качества восстановления. Возможно, более предпочтительным является путь непосредственного подбора подходящих процедур формирования матриц линейного преобразования вида (11) с использованием субъективных оценок качества восстановления.
Полученные теоретические результаты указывают на принципиальную возможность достижения различного регуляризующего эффекта за счет непосредственного изменения степени близости оценок параметров к их истинным значениям. Ясно, что качество восстановления, в конечном итоге, будет определяться достоверностью предположений о
соответствии моделей искажений тестовых и реальных изображений.
7. Благодарности
Автор выражает признательность д.т.н., проф. Сергееву В.В. за обсуждение работы и ценные замечания, способствовавшие ее улучшению, а также научному сотруднику лаборатории Анализа оптических сигналов Института систем обработки изображений РАН Устинову А.В. и студенту Баранову В.Г. за помощь в подборе и формировании тестовых изображений.
Литература
1. Демиденко Е.З. Линейная и нелинейная регрессии.- М.: Финансы и статистика, 1981, 303с.
2. Fursov V.A. Identification of optical distorting systems with selecting image informative fragments. Workshop on Digital Image Processing and Computer Graphics. Proceedings SPIE. - 1994. 2363.
3. Шамриков Б.М. Сравнительный анализ точности параметрической идентификации динамических объектов в разомкнутых и замкнутых автоматических системах. - Изв. АН СССР, Техн. кибернетика, 1986, №3, с. 143-150.
4. Фурсов В. А. Анализ точности и построение алгоритмов идентификации по малому числу на-блюдений.-Изв. АН СССР, Техн. кибернетика, N 6, 1991 г.
