Борискевич А.А. Цифровая обработка речи и изображений - файл тема10.doc

Реконструкция размытых изображений в MATLAB

Рассмотрим некоторые подходы к решению задачи восстановления изображений с использованием функции протяженности точки. При реализации алгоритмов в среде MATLAB будем использовать соответствующие функции (deblurring functions) пакета Image Processing Toolbox.

Определим некоторые понятия, которыми будем пользоваться в процессе изложения материала.

Деконволюция - процесс, обратный к свертке.

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

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

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

Рассмотрим коротко основы технологии обработки размытых изображений.

Причины размытости изображений.

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

Модель размытых изображений.

Размытые и искаженные изображения приближенно могут быть описаны выражением

где - размытое изображение;

Оператор искажения, определяемый функцией протяженности точки. Свертка этой функции с изображением, является причиной искажений;

Исходное изображение;

Аддитивный шум, который вносится во время захвата изображений и искажает их.

Функция протяженности точки

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

Проиллюстрируем это на примере. Для создания функции протяженности точки используем функцию fspecial, которая будет симулировать искажения, обусловленные движением, зададим также длину размытия (LEN=31), и угол степени размытия (THETA=11). После того, как создана функция протяженности точки, используем, например, функцию imfilter для свертки с исходным изображением I, и создания размытого изображения.

I = imread("peppers.png");
I = I(60+,222+,:);
figure; imshow(I); title("Original Image");
LEN = 31;
THETA = 11;
PSF = fspecial("motion",LEN,THETA); % создание функции протяженности точки
Blurred = imfilter(I,PSF,"circular","conv");
figure; imshow(Blurred); title("Blurred Image");

Исходное изображение Размытое изображение

Использование функций обработки размытых изображений

В пакете Image Processing Toolbox системы MATLAB существует четыре функции обработки размытых изображений:

1. deconvwnr - выполняет восстановление размытых изображений с использованием винероской фильтрации;
2. deconvreg - выполняет восстановление размытых изображений с использованием регуляризационных фильтров;
3. deconvlucy - выполняет восстановление размытых изображений с использованием алгоритма Лаки-Ричардсона (Lucy-Richardson);
4. deconvblind - выполняет восстановление размытых изображений с использованием алгоритма слепой деконволюции;

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

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

Функция deconvblind используется для реализации алгоритма слепой деконволюции, когда неизвестна функция протяженности точки, искажающая изображение. Результатом работы функции deconvblind является восстановленное изображение и функция протяженности точки. При выполнении алгоритма используются итерационные модели, аналогичные функции deconvlucy.

Восстановление размытых изображений методами винеровской фильтрации

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

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

1. Считываем изображение в рабочее пространство MATLAB.

I = imread("peppers.png");
I = I(10+,222+,:);
figure;imshow(I);title("Original Image");

Исходное изображение

2. Создание функции протяженности точки.

LEN = 31;
THETA = 11;
PSF = fspecial("motion",LEN,THETA);

3. Симуляция размытости на изображении.

Blurred = imfilter(I,PSF,"circular","conv");
figure; imshow(Blurred);title("Blurred Image");

Размытое изображение

4. Восстановление размытого изображения.

Изображение восстановленное методами винеровской фильтрации

Анализ результатов обработки

Результаты работы функции deconvolution можно использовать для определения оптимальных аргументов при реализации функции deconvwnr. Используя эти аргументы можно определить соотношение шум-сигнал и/или автокорреляционную функцию для усовершенствования результатов восстановления.

Мы полагаем, что на практике наиболее широко применимым метопом деконволюции слелуегг считать мультипликативный метод винеровской фильтрации (см. § 16). Рассматриваемый здесь пример иллюстрирует различные практические аспекты винеровской фильтрации с применением некоторых методов предварительной обработки, описанных в § 15.

На рис. 3, а показан фрагмент изображения а из примера 1. Здесь мы принимаем изображение а за истинное изображение (см. § 4, первый абзац). Предположим теперь, что этот фрагмент фотографируется несфокусированной камерой, так что получается искаженное изображение с которой свертывается изображение а для получения изображения имеет форму однородного кружка (см. табл. 1.1) диаметром 15 элементов, что соответствует изображению а на рис. 2 в примере 2. Допустим, что, если не считать несфокусированности объектива, получаемые нами фотографии имеют столь высокое качество, что на искаженном изображении в имеется только шум квантования, связанный с восьмиразрядной точностью дискретизации. Далее мы покажем, что такое практически незашумленное исходное идеальное искаженное изображение [см. определения (3.4) и Изображение - это полный искаженный вариант изображения а. Последнее содержит 242 х 242 элементов, тогда как изображение состоит из 256 х 256 элементов. Поскольку диаметр ФРТ равен 15 элементам, максимальное горизонтальное (или вертикальное) размытие каждой точки на изображении а составляет 7 элементов слева и справа (или нверх и вниз), чем и объясняется, что изображение имеет на 14 элементов больше, чем изображение а, как в горизонтальном, так и в вертикальном направлениях. Поскольку изображение является «незашумленным», можно без труда практически идеально восстановить изображение а методом винеровской фильтрации. Отмстим, что «полнота» изображения имеет даже более существенное значение, чем отсутствие зашумленности. Различные виды искаженных изображений, которые встречается в практических приложениях, всегда записаны в кадрах конечных размеров, так что информация, содержащаяся в «заданном» искаженном изображении, всегда как-то «обрезана» вблизи границ. При попытках восстановить истинное изображение, исходя из усеченного искаженного изображения, мы обычно встречаемся со значительными трудностями.

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

(кликните для просмотра скана)

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

Изображения и демонстрируют, что происходит, когда винерозекан фильтрации осуществляется без отмеченных мер предосторожности (т. е. без предварительной обработки. рассматриваемо» в § 15). Величина введенная в формуле (16.5), предполагается постоянной, т. е. равной константе фильтра которая введена после формулы (16.8). При получении изображений гид были взяты жачення константы равные 0.0009 и 0,01. Артефакты в виде «ряби» и низкая четкость «восстановленных» деталей типичны для искаженных изображений, подвергнутых винеровской фильтрации без предварительной обработки. Эти два изображения позволяют сделать еще один важный практический вывол: качество восстановленного изображения в какой-то мере зависит от взятого значения константы

Отмстим олин вычислительный момент практического характера, связанный к изображениями Поскольку винсровская фильтрация оенпнанз на преобразовании Фурье, для экономии машинного времени приходится прибегать к алгоритму БПФ (см. § 12 и пример 2). Для этого нужно, чтобы изображения, которые трубуется фильтровать. состояли из элементов, где целое положительное число. Мы могли бы «дополнить нулями» изображение в [см. § 12, перед абэаием, содержащим формулы (12.14) - (12.17)]. Но мы предпочли «участить» его методом билинейной интерполяции (см. § 45 и текст, относящийся к формуле (47.1)]. гак чтобы и нем было элементов. После вннеровскон фильтрации «частота» каждого восстановленного изображения (т. е. была снопа уменьшена до элементов.

Прежде чем демонстрировать выгоды предварительной обработки, рассмотрим кадра которые определены в § 15 и представлены на рис. 3. о. Кадр это кадр записи, в пределах которого существует фактически записанное изображение (см. § 4). В нашем примере кадр из 228 х 228 элементов, содержащий изображение в. Иногда оказывается удобным считать, что имеется некое «исходное изображение». лишь частью которого является рассматриваемое истинное изображение (т. е. в нашем примере изображение а). Но если фактически записанное изображение есть все, что нам известно, то мы можем получить какую-либо информацию только о той части исходного изображения, которая непосредственно влияет на фактически записанное изображение. Действие состоит в том, что каждая точка исходного изображения размывается к некую область плоскости изображения с центром в Исходной точке, по площади равную кадру ФРТ. Последний представляет собой минимальный прямоугольный кадр, содержащий все точки х, в которых величина заметно отлична от нуля (в данном примере это квадрат размером элементов). Вообще говоря, ширина ФРТ неодинакова слева (скажем, в пределах элементов), справа (в пределах элементов), вверху (в пределах элементов) и внизу (в пределах элементов); поэтому кадр состоит из элементов, так как нулевое искажение характеризуется кадром ФРТ, содержащим один элемент (а не нуль элементов!) в рассматриваемом примере Таким образом, ФРТ размывает часть исходною изображения, лежащую в пределах области в кадр на рис. 3, а), размеры которого равны размерам области плюс элементов в горизонтальном направлении и элементов в вертикальном направлении (в нашем примере кадр состоит их элементов). Обратно, ФРТ размывает на область по

крайней мере частично, каждую точку исходного изображения, которая лежит в пределах калра И. Поэтому часть исходного изображения, лежащая в пределах кадра и есть та часть, которая оказывает непосредственное влияние на фактически записанное изображение. Назовем эту часть исходною изображения истинным изображением (в нашем примере это изображение а). Кроме того, размывает истинное изображение на больший кадр на рис. 3. а), размеры которою равны размерам кадра увеличенным на элементов изображения в горизонтальном направлении и элементов в вертикальном направлении; в нашем примере кадр (состоит из 256 х 256 элементов (это кадр, ограничивающий изображение 6).

Желательно иметь какой-то объективный критерий качества восстановления. Мы здесь выбрали критерий средней абсолютной ошибки

Элемент восстановленного изображения, истинного изображения (обозначение показывает, что все элементы изображения лежат в пределах калра . На рис. 3. и представлена зависимость величины (5 от Кривая А подтверждает вывод, сделанный выше о том, что изображение а можно более или менее хорошо восстановить по изображению если величина достаточно мала, то величина (2 становится ничтожно малой. Кривая И показывает, что наилучшее восстановление изображения в (без предварительной обработки) соответствует весьма значительному уровню зашумленности (эквивалентный среднеквадратичный «уровень шума» составляет что существенно выше, чем в большей части важных практических приложений, и этим подтверждается сделанный ранее вывод о том. что эффект усечения изображения в может доминировать над влиянием любого реального уровня зашумленности).

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

Теперь мы можем проиллюстрировать улучшение качества восстановленного изображения, которое достигается методом предварительной обработки - расширением границ, описанным в абзаце, содержащем формулу (15.34), и в шести последующих абзацах. Изображение это результат применения метода расширения границ с целью получить из изображения в оценку идеального искаженного изображения (которое существует на кадре Подчеркнем, что не может быть и речи о том, чтобы достоверно восстановить часть изображения лежащую между границами кадров . От применения метола расширения границ можно ожидать только компенсации эффектов усечения - о восстановлении истинного изображения вне области нечего и думать. Изображение есть наилучший результат восстановления по изображению (соответствует минимуму кривой С на рис. 3, н). По сравнению с изображением здесь заметное улучшение, но многие детали изображения остаются еще несколько искаженными и имеется остаточная «рябь», которая обусловлена тем, что простое расширение границ не позволяет получить искаженное изображение, достаточно близкое к свертке истинною изображения с в этом проявляется принципиальная несогласованность конечных сверток (см. § 14).

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

детали, имеющиеся на изображении а, не четко восстанавливаются на изображении и. шссь значительно меньше «ряби», которая сильно портит изображения и в меньшей степени Это проявление принципиальной несогласованности периодических сверток, отмеченное в следующем абзаце после формулы (14.8). Здесь видно, насколько важно с чисто практической точки зрения следить за тем, чтобы расчеты как можно более согласовались с формулами математической физики, на которых оснеован тот или иной процесс обработки, если это, конечно, можно сделать без больших усложнения. Изображение к - результат применения винеровской фильтрации к изображению при Изображение к лишь немного хуже изображения и, чем подтверждается наш вывод о том, что качество восстановленного изображения почти не зависит от величины если она достаточно близка к оптимальному значению.

Теперь покажем, к чему приводит в случае винеровского фильтра применение неверной ФРТ. Изображение восстановлено по изображению з с применением гауссовской ФРТ (см. табл. 1.1.) вместо ФРТ, отвечающей расфокусировке. Эффективный диаметр гауссовского «колокола» был таким же, как и диаметр кружка ФРТ в случае расфокусировки (этот кружок представлен на рис. 2 а, а случай гауссовского размытия - на рис. 1, е). Изображение получено в результате восстановления с ФРТ. отвечающей расфокусировке, но при диаметре кружка, равном 19. а не 15 элементам изображения. Константа в случае изображений была взята такой же, как и в случае изображения и. Изображения заметно хуже изображения и. чем подтверждается наш вывод о том, что правильная оценка формы ФРТ важнее, чем точное задание оптимального значения константы фильтра

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

Спектральная плотность сигнала определяется соотношением:

где – автокорреляционная функция.

Взаимная спектральная плотность сигнала определяется соотношением:

, (14)

где – функция взаимной корреляции.

При построении фильтра Винера ставится задача минимизации среднеквадратического отклонения обработанного изображения от предмета:

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

.

Дальнейший анализ показывает, что восстановление изображения, формирование которого описывается выражением должно осуществляться с использованием следующего ОПФ восстанавливающего преобразователя:

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

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

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

Рис. 1. Примеры фильтров

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

Фильтр Винера использует оптимальный критерий в виде минимизации среднеквадратичного отклонения между фильтрованным изображением и истинным объектным изображениемf (i ,j ). В частотном домене фильтр Винера имеет вид :

(5.18)

где |N (u ,v )| 2 и |F (u ,v )| 2 – спектр мощности щума n (i ,j ) и объекта f (i ,j ) (см. уравнение 5.13).

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

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

Рис. 5.16. Клинические изображения костного скелета, получающиеся после разных видов фильтрации: верх слева – без фильтрации; верх справа – фильтр Баттеруорта четвертого порядка с пороговой частотой 0,4; низ слева – фильтр Метца; низ справа – фильтр Винера

Контрольные вопросы

Опишите структуру цифрового изображения гамма-камеры.

Какие факторы влияют на размер пикселя изображения?

В чем отличия фреймового способа запоминания данных от листингового и байт-моды от слово-моды?

Что такое формат DICOM и для чего он применяется?

Какие задачи выполняет PACS ?

Назовите физические факторы, влияющие на качество изображения и на пространственное разрешение гамма-камеры.

На какие параметры изображения влияет комптоновское рассеяние фотонов?

Почему возникает шум в изображении?

Как определяется информационная плотность (ID ) изображения?

Что такое контраст изображения и какая его величина требуется для визуального обнаружения патологических очагов в организме пациента?

С какой целью и каким образом производится преобразование изображения в частотное пространство?

Как создается выборочная версия непрерывной функции?

Какой критерий должен выполняться, чтобы непрерывная функция однозначно определялась из N выборочных значений?

Опишите математическую модель процесса визуализации.

С какой целью проводится фильтрация изображения.

На какие группы подразделяются фильтры?

Для чего применяются низкочастотные фильтры?

С какой целью применяется восстановительная фильтрация?

Как зависит пороговая частота восстановительного фильтра от уровня шума?

Охарактеризуйте особенности фильтров Винера и Метца.

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

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

Преобразуем передаточную функцию фильтра Винера (4.28) следующим образом:

(4.29)

Анализируя соотношения (4.28) и (4.29), можно отметить следующее:

1. При отсутствии шума фильтр Винера переходит в инверсный фильтр. Следовательно, в области низких частот, где, как правило, отношение сигнал/шум велико, передаточные функции инверсного и винеровского фильтров практически совпадают.

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

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

На рис. 4.18 приведены одномерные сечения типичных передаточных функций винеровских фильтров (сплошная линия). Здесь же для сравнения приведены сечения передаточных функций инверсных фильтров (4.15) и (4.16), которые обозначены штриховой линией.

Рассмотрим результаты моделирования винеровского алгоритма восстановления. На рис. 4.19.а и 4.21.а приведены результаты искажения изображений «Сатурн» и «Часы» сверткой с гауссовской ФРТ () с последующим «обрезанием» краев и добавлением аддитивного дельта-коррелированного шума (). На рис. 4.20.а и 4.22.б приведены изображения, полученные в результате смаза () изображений «Сатурн» (рис. 4.6) и «Часы» (рис. 4.22.а) () также с последующим «обрезанием» краев и добавлением аддитивного дельта-коррелированного шума ().

Размеры всех наблюдаемых и восстановленных изображений равны элементов. Результаты восстановления винеровским фильтром изображения «Сатурн» (рис. 4.19.б и рис.4.20.б) свидетельствуют о том, что фильтр Винера значительно лучше подавляет шумы. Осциллирующая помеха на результатах восстановления изображения «Часы» (рис. 4.21.б и рис.4.22.в) вызвана краевыми эффектами. Очевидно, что ее уровень существенно меньше, чем при инверсной фильтрации (см. рис.4.13.в). Однако винеровский фильтр лишь частично компенсирует краевые эффекты, которые делают качество восстановления неудовлетворительным.

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