Как найти спектр изображения

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

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

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

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

Пространственная область

Вспомним, что картинка представляется в виде функции от x и y. Если мы говорим о полноценном изображении, то каждое значение этой функции – трёхзначное число, которое представляет собой значение для каждого из цветовых каналов.

Представление изображения в пространственной области – это то, как мы привыкли понимать и видеть изображение. Есть x и y, и в каждой точке у нас есть какое-то значение интенсивности или значение цветового канала.

На этом слайде Лена и логотип библиотеки OpenCV разложены на три цветовых канала — красный, зелёный и синий.

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

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

Чтобы упростить задачу, представим, что наше изображение одномерное. Линия, которая идёт слева-направо, отображает все изменения яркости. Всплеск в начале соответствует паре пикселов белого цвета, потом идёт серая область, потом — снова белая. На черном цвете мы проваливаемся вниз. То есть по этой линии можно проследить, как изменяется яркость.

Посмотрим еще один пример.

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

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

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

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

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

Что мы можем сказать о коэффициентах базисных функций? Если у нас большой коэффициент перед базисной функцией с высокой частотой, это означает, что яркость изменяется достаточно часто. На картинке очень много перепадов яркости в маленьких локальных регионах. Если картинка описывается плавными синусоидами, с низкой частотой, то это означает, что на картинке много однородных областей, яркость изменяется плавно, или картинка, например, была «забурена».

Таким образом, можно использовать отображение в частотной области для описания изображений.

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

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

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

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

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

Здесь гармоники — это те самые синусы и косинусы, которые нарисованы на предыдущем слайде. Для каждой фиксированной частоты есть некая функция от x.

Надеюсь, что с одноименным случаем более или менее понятно. Теперь посмотрим на двумерный случай, потому что картинка у нас двухмерна.

Здесь мы тоже можем построить двухмерные гармоники, которые уже будут зависеть от четырёх параметров: x, y (от двух направлений) и от двух частот в направлении x и y соответственно.

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

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

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

Дальше, если каждую ячейку начнём складывать значение параметра F, коэффициента, который мы получили при разложении. Чем больше коэффициент, тем ярче он будет отображаться на спектре. То есть, мы хотим визуализировать спектр Фурье. Если коэффициент F=0, будем отображать его чёрным цветом. Чем он больше, тем светлее цвет.

На двумерном спектре всего один пиксел будет соответствовать u=0, v=0, следующий — U= -1, V=0. Значение пиксела будет равняться коэффициенту, полученному из преобразования. Важно то, что центральные коэффициенты соответствуют важности этой гармоники в представлении изображения, в центре стоят гармоники с нулевыми частотами. То есть, на картинке, если у нас будет очень большой отклик вот здесь, это значит, что картинка практически не меняет своей яркости. Если у картинки будет здесь тёмное пятно, а по периферии — светлое, то она пёстрая — там перепад яркости в каждой точке и в каждом направлении.

Картинка — это не спектр, это визуализация двумерной синусоиды.

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

Обработка в пространственной области

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

Инвертирование картинки

Если у нас яркость изменяется от 0 до 255, то для каждого пиксела пишем 255 минус его прежнее значение. Чёрное становится белым, а белое — чёрным.

Одним из наиболее простых и практичных видов представления картинки является гистограмма. Для любого изображения можно посчитать, сколько пикселов с нулевой яркостью, сколько с яркостью 50, 100 и так далее, и получить некое распределение частотности.

На этой картинке у Лены отрезаны все точки яркости, которые больше 125. Получаем гистограмму, смещённую влево, и тёмную картинку. На второй картинке все наоборот — есть только светлые пикселы, гистограмма смещена вправо.

Следующая картинка размазана, у неё нет чётких контуров. Для такой картинки характерен столбик яркости где-то посередине и отсутствие энергии по краям спектра.

Для второй Лены с этого слайда была сделана обработка, чтобы подчеркнуть каждый из контуров, то есть сделать сильнее перепады яркости. Здесь гистограмма занимает весь диапазон яркости.

Что можно делать с гистограммами? По форме гистограммы уже можно многое рассказать о её свойствах. Дальше можно что-то сделать с гистограммой, изменить её форму, чтобы картинка стала выглядеть лучше.

Для тёмной картинки гистограмму можно растянуть вправо, и тогда картинка станет светлее, а для светлой — наоборот. Это верно для любой формы. Если гистограмма картинки не покрывает весь частотный диапазон, то при равномерном растягивании картинка станет более контрастной и будет лучше видно детали.

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

Делается это так. В формуле xk — это некий уровень яркости, nk — количество пикселов такой яркости, n- общее количество пикселов. Мы выбираем какой-то пиксел и получаем вероятность, с которой он будет цвета xk. То есть, количество пикселов nk делить на общее количество пикселов. На самом деле — доля.

Попытаемся получить равномерную плотность вероятности, то есть, чтобы для каждого цвета была одинаковая вероятность его получить. Это достигается следующим преобразованием. Если мы по пикселю цвета xk вычислим новые значения яркости для него yk по формуле, то есть возьмём и пронумеруем все вероятности oт i до этого цвета, гистограмма получится более равномерная. То есть, исходная затемнённая картинка выглядела вот так, а если мы к яркостям этой картинки применим экранизацию, то в результате получим гистограмму вот такой формы. Как вы видите, они гораздо более равномерно распределены по всем возможным значениям, и картинка будет выглядеть вот так.

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

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

Бинаризация

Мы уже отказались от цветных изображений и говорим только о чёрно-белых, где есть различные градации серого. Бинаризация — это продолжение издевательства над изображением, когда мы отказываемся еще и от них. В итоге получаем картинку, где есть только чёрный и белый. Нам нужно понять, какие пикселы мы хотим сделать чёрными, а какие — белыми.

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

Поговорим о простейшем типе бинаризации — пороговой. Этот тип вообще-то не очень применим для фотографий, но, тем не менее, иногда используется.

Если мы применяем пороговую бинаризацию к гистограмме, то видим, что у нас в изображении есть два типа пикселов: потемнее и посветлее. Обычно предполагают, что большее количество пикселов соответствует фону. Из этого делаем вывод, что вот тут у нас кучка побольше, она тёмная, соответственно у нас фон — тёмный, и на нём есть один или несколько светлых объектов. Объект может быть составной. Здесь у нас два светлых объекта разного цвета.

Это очень красивые гистограммы, в жизни вы вряд ли такие увидите. Но по ним легко понять, где нужно провести порог, чтобы отделить фон от объекта. Здесь, если мы возьмём пороговое значение ровно между ними, и все пикселы, которые ярче порога, сделаем белыми, а те, которые темнее — чёрными, то мы превосходно отделим объект от фона. Выбрав нужный нам диапазон яркости, мы вырезаем из картинки тот или иной объект.

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

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

На этой картинке представлено преимущество локальной бинаризации.

Если бы бинаризация была глобальной, был бы выбран глобальный порог для всей картинки, то результат бы выглядел так, на картинке b. Пикселы здесь в одной области примерно одного уровня яркости, поэтому вся эта часть получается белой, засвеченной. В том случае, если мы применим простейший вариант локальной бинаризации — то есть поделим изображение на сегменты и будем выбирать порог для каждого сегмента отдельно, — результат бинаризации уже будет выглядеть получше .

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

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

В этих двух областях всё плохо, потому что тут есть совсем чёрные пикселы и есть — посветлее. Для каждого квадратика находим порог и применяем его внутри этого квадратика.

Выделение компонент связности

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

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

Простейший алгоритм — двухпроходный, он работает следующим образом. Мы начинаем из левого верхнего угла и первому пикселу присваиваем номер. Дальше двигаемся направо и смотрим, совпадает ли цвет пиксела с уже размеченным. Если совпадает, то присваиваем ему ту же самую метку. Если этот пиксел был помечен ноликом, то и этот получает тоже метку нолик. И так доходим до конца строки, потому что тут у нас все нолики. Дальше, цвет третьего пиксела на второй строчке не совпадает с цветом уже размеченных соседей. Мы увеличиваем счётчик областей и присваиваем этому пикселу номер следующей области — 1. Переходим сюда, у этого уже ест сосед такого же цвета, у которого есть метка, ему присваиваем такую же метку. Дальше мы увеличиваем увеличиваем счётчик и присваиваем номера 2, 3, 4, 5, 6, 7. Здесь представлен результат первого прохода. На втором проходе нам останется только объединить области, в которых соседи одного цвета носят разные метки. В результате получается изображение такого вида.

Здесь фон одного цвета, одна связная компонента другого, и вторая связная компонента третьего цвета. Красная и белая компоненты на чёрном фоне.

Уже в первом проходе накапливается информация о том, что одной компоненте соответствует несколько меток. Понятно, что все они одного цвета и 1 и 2 относятся к одной компоненте. Эта информация записывается и во второй проход перенумеруется.

Где используется? Это один из способов сегментации картинки. Простейший пример — когда картинка бинарная, можно выделять компоненты связности на цветном изображении. Только там критерием присоединения к одной компоненте будет не полное соответствие цвета, а то, что они отличаются на какое-то пороговое значение.

Алгоритмы фильтрации

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

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

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

То есть, если у нас есть вот такая картинка и такой фильтр (3 на 3), и если мы ничего не делаем с картинкой, то получается отфильтрованной только вот эта часть. Теряем по одному пикселу с каждой стороны. Мы можем нарастить исходную картинку еще на один пиксел, вставив по краям чёрные либо белые пикселы, либо зеркально отразив край. Если у нас фильтр большего размера, то мы соответственно увеличиваем наращиваемую по краю область.

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

Дистрибутивность по сложению — можем сложить два фильтра, просто сложить их веса и получить новый фильтр, или применить сначала один фильтр, потом второй и сложить результаты. Скаляр можем выносить за скобки. То есть если здесь есть коэффициент, можем все веса поделить на какое-то число, применить фильтр, а потом результат умножить на эту константу. Что произойдет, если мы картинку свернём вот с таким фильтром? Сдвиг на один пиксел.

Я обещала рассказать, почему теорема Фурье так хорошо работает в частотной области. Существует теорема свёртки. Что это такое? Теорема свёртки говорит о том, что если мы производим свёртку в пространственной области, это то же самое, что перемножить результаты преобразования Фурье в частотной области. Если есть картинка F и фильтр H, и мы их попытаемся свернуть, это будет то же, как сделать перевод картинки и фильтра в частотную область. То есть найти для картинки и фильтра коэффициенты преобразования Фурье и их между собой перемножить. И в обратную сторону. Если у нас есть перемножение в пространственной области, это то же самое, что и операция свёртки над соответствующими коэффициентами преобразования Фурье.

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

Какие могут быть фильтры/ядра и что они делают с изображением

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

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

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

Этот фильтр уже будет нелинейным. Медианные фильтры более устойчивы к выбросам, сильным отклонениям от средних значений, поэтому медианные фильтры значительно лучше справляются с подавлением шума типа “соль и перец”.

Еще один очень распространённый фильтр — фильтр Гаусса. Это свёртка с функцией Гаусса, которая выглядит таким образом. Здесь есть зависимость от радиуса. Вот представление фильтра в пространственной области и представление его же в частотной области. Можно ли по виду характеристик фильтра определить, что он сделает с картинкой?

Вспомним теорему о свёртке. Мы собираемся перемножать образ нашего изображения вот с таким вот. Что получится? Мы существенно обрезаем высокие частоты, изображение получится более гладким. Вот результат сглаживания.

Исходное изображение, картинка сглаженная фильтром Гаусса с =1,4 размера 5 и =2,8 с размером 10. Чем больше радиус, тем более размытая получается картинка, хотя явно нет зависимости между размером фильтра(количеством отсчётов) и на практике не имеет смысла выбирать большую с малым размером.

Здесь у нас размер фильтра 3х3 (маленький), если большая, это означает, что мы не сможем отразить все значения. Негласное правило — размер должен составлять порядка 5-6. Сглаживание следует использовать, когда есть шум и мы хотим от него избавиться. То есть, пытаемся усреднить значение в каждой точке, посмотрев на значения соседей, и если это какой-то выброс, то мы это значение уберём.

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

Почему тут практически не осталось вот этих контуров? На картинке их очень хорошо видно. На свёртке точка стала хорошо видна, а контуров практически не видно. Потому что точка изолированная и значение в точке над которой расположена маска гораздо выше, чем всё вокруг.

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

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

Как найти перепад яркости на изображении. Что есть перепад яркости? Это когда у нас были пикселы одного цвета, а потом стали пикселы другого цвета. Если мы посмотрим на профиль картинки, сначала было всё чёрное, потом — раз, стало серое. Или оно плавно изменялось. То есть, у нас есть скачок или перепад яркости. Для того чтобы его выделить, можно просто взять производную. В плоских значениях (там, где яркость не меняется) производная равна 0. Там где яркость меняется, там чем выше перепад, тем больше будет значение производной. Первая производная позволит выделить участок плавного изменения яркости. Есть такое понятие как zero crossing line. Если мы посчитаем значения производных и их соединим, то можем предположить, что контур расположен где-то вот в этой точке, там, где мнимая прямая пересекает вот этот уровень.

Как можно посчитать производную на картинке и что из себя представляет градиент изображения? Мы смотрим, насколько изменяется вдоль x и y наша яркость.

Подсчитав градиент, мы можем получить направление градиента, которое задаётся углом и величину, которая задаётся изменением по каждому из направлений.

Вычисление градиента. В дискретном случае дискретная производная по x будет выглядеть как-то вот так. Считаем разницу в значении пиксель сдвинутого в одну сторону на 1 по сравнению с текущим пикселям и можно это тоже задать маской. Существуют разные маски для вычисления градиента.

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

С помощью определённых масок можно посчитать вторую производную, по x и по y. Помимо того что мы можем посчитать контуры, иногда нужно повысить резкость. Когда мы сглаживаем картинку, мы теряем детали. Если мы возьмём оригинал и вычтем из него результат сглаживания, то мы получим некие контуры.

Для того чтобы получить более резкую и контрастную картинку, нужно добавить эти контуры. Один из способов повышения контраста — лапласиан Гаусса. Проделаем ровно то, что здесь нарисовано, берём исходное изображение, и прибавляем к нему разность исходного изображения и сглаженного, предположим, что это у нас гуанина. Путём некоторых математических преобразований получаем вот такое выражение, то есть нам нужно свернуть исходное изображение вот с таким вот фильтром, в котором — это константа, e — единичный фильтр, который выглядит вот таким образом, g — это лапласиан.

То есть, из единичного фильтра вычитаем Гаусса и получаем вот такую бяку, называемую лапласиан Гаусса, или иначе — Mexican Hat.

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

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

# -*- coding: utf-8 -*-
import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

# Читать изображение
img = cv.imread('test.png', 0)

 # Алгоритм быстрого преобразования Фурье для получения частотного распределения
f = np.fft.fft2(img)

 # По умолчанию центральная точка результата находится в верхнем левом углу,
 # Вызов функции fftshift () для перехода в среднее положение
fshift = np.fft.fftshift(f)       

 #fft Результат - комплексное число, а его абсолютное значение - амплитуда.
fimg = np.log(np.abs(fshift))

 # Показать результаты
plt.subplot(121), plt.imshow(img, 'gray'), plt.title('Original Fourier')
plt.axis('off')
plt.subplot(122), plt.imshow(fimg, 'gray'), plt.title('Fourier Fourier')
plt.axis('off')
plt.show()

https://www.cnblogs.com/xh6300/p/5956503.html

Исследование характеристик двумерной спектрограммы с преобразованием Фурье изображения

1. Сначала сделайте несколько важных выводов:

1. Амплитуда преобразования Фурье называется спектром Фурье или частотным спектром.

2. Нулевое положение F (u) обратно пропорционально ширине W «коробчатой» функции.

3. Теорема о свертке: преобразование Фурье свертки двух функций в пространственной области равно произведению преобразования Фурье двух функций в частотной области. f (t) * h (t) <=> H (u) F (u)

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

5. Сильное замешательство может даже вызвать полное непонимание.

6. Наиболее медленно изменяющаяся частотная составляющая (u = v = 0) пропорциональна среднему уровню серого изображения. Термин DC определяет средний уровень серого изображения.

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

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

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

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

11. Инвертируйте изображение в оттенках серого, и “узор” его спектра не изменится. (Личное понимание: инверсия только переворачивает черный и белый, но не меняет контраст при изменении уровня серого)

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

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

Используемый инструмент анализа преобразования Фурье – Halcon, код выглядит следующим образом:

read_image (Image, 'C:/Users/xiahui/Desktop/1.jpg')
fft_image (Image, ImageFFT)

2. Анализ спектрограмм различных изображений.

Слева – исходное изображение, а справа – спектр после преобразования Фурье.

1. Все черное изображение-спектрограмма тоже вся черная (Разрешение изображения 240 * 240

2. Серое изображение – в центре спектрограммы есть белый маленький квадрат размером с один пиксель, координаты (120,120), значение (30480,0).

3. Все белое изображение – в центре спектрограммы есть небольшой белый квадрат размером в один пиксель. Координаты (120, 120) и значение (61200, 0).

4. Нарисуйте круг на картинке – спектрограмма имеет форму концентрических кругов, значение центра (координаты 120, 120) равно (3852,64, 0), значения в других местах положительные или отрицательные, а тенденция ближе к центральному значению. Чем больше абсолютное значение.

5. Нарисуйте квадрат на картинке – центральное значение (5143,03, 0), а значения в других местах – положительные или отрицательные. Тенденция такова, что чем ближе к центральному значению, тем больше абсолютное значение.

6. Поверните квадрат на приведенном выше рисунке на 45 градусов – самое центральное значение (5140,22, 0), которое можно рассматривать как 5140,22≈5143,03; значения в других местах положительные или отрицательные, и тенденция такова, что чем ближе к центральному значению, тем больше абсолютное значение.

7. Нарисуйте два круга – центральное значение (7704,13,0).

8. Нарисуйте белый и серый кружки – центральное значение (5772,82,0).

9. Нарисуйте четыре круга – самое центральное значение (15402,8,0).

10. Нарисуйте 2 квадрата – центральное значение (10061,8,0).

11. Нарисуйте четыре квадрата, повернутых на 45 градусов – самое центральное значение (20114.3, 0)

12. Проведите прямую линию – центральное значение (766,0).

13. Проведите линию под углом 15 ° – самое центральное значение будет (876,571,0).

14. Нарисуйте пару пересекающихся линий – самое центральное значение (1194,55,0).

15. Нарисуйте плавный круг – хотя он также концентрический, он не так очевиден, как раньше; самое центральное значение – (1849,6, 0)

16. Заполните все изображение градиентом – центральное значение (30470, 0), которое можно считать равным значению серого изображения (30480).

17. Нарисуйте полосу градиента с уклоном 15 ° – самое центральное значение (12051,9,0).

18. Сделайте снимок и проведите различную обработку, а затем посмотрите на их спектрограммы соответственно.

Исходное изображение, обратное изображение, слабое размытие по Гауссу, сильное размытие по Гауссу, среднее серое изображение, обратное среднее серое изображение.

Обработанный код выглядит следующим образом:

read_image (Image1, 'C:/Users/xiahui/Desktop/1.jpg')
read_image (Image2, 'C:/Users/xiahui/Desktop/2.jpg')
read_image (Image3, 'C:/Users/xiahui/Desktop/3.jpg')
read_image (Image4, 'C:/Users/xiahui/Desktop/4.jpg')
read_image (Image5, 'C:/Users/xiahui/Desktop/5.jpg')
read_image (Image6, 'C:/Users/xiahui/Desktop/6.jpg')
fft_image (Image1, ImageFFT1)
fft_image (Image2, ImageFFT2)
fft_image (Image3, ImageFFT3)
fft_image (Image4, ImageFFT4)
fft_image (Image5, ImageFFT5)
fft_image (Image6, ImageFFT6)

Из вышесказанного видно, что после инверсии «узор» спектрограммы изображения не изменяется. Но на самом деле значение меняется. Значения центральных точек (120, 120) на 6 изображениях следующие:

(47235,0)、(13965,0)、(47169.5,0)、(47130.4,0)、(47280,0)、(13920,0)

По приведенному выше примеру можно сделать два вывода:

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

2. Значение центральной точки (120, 120) изображения должно быть связано со средним уровнем серого изображения. В примере 3 (полностью белое изображение) уровень серого не изменяется, поэтому нет изменения частоты. Вся энергия сосредоточена в спектрограмме. В блоке одного белого пикселя в центре его значение составляет 61200 (чем выше разрешение изображения, тем больше предельное значение). В примере 18 значение центральной точки среднего спектра серого изображения составляет 47280, а среднее значение серого инвертировано. Спектр фигуры 13920, а 47280 + 13920 = 61200. Уровень серого для среднего изображения с серой шкалой составляет 198, а уровень серого для инвертированного среднего изображения с серой шкалой составляет 57.

Информация в данной статье относится к релизам программы MATLAB ранее 2016 года, и поэтому может содержать устаревшую информацию в связи с изменением функционала инструментов. С более актуальной информацией вы можете ознакомиться в разделе документация MATLAB на русском языке.

Краткий курс теории обработки изображений. Автор – И.М.Журавель.

Содержание:

  • Свойства зрительной системы человека
  • Возможности цифровой обработки изображений в Matlab
  • Типы изображений
    – Бинарные изображения: геометрические характеристики
    – Бинарные изображения: топологические характеристики
  • Формирование и обработка цифровых изображений
  • Локально-адаптивная обработка изображений
    – Адаптивное повышение контрастности изображений
    – Энтропия изображения. Использование среднеквадратического отклонения значений яркостей элементов окрестности в методах контрастирования изображений. Нелинейное растяжение локальных контрастов.
    – Анализ некоторых характеристик локальных окрестностей
    – Статистическое определение локального контраста
    – Локально-адаптивное улучшение качества изображений
  • Фильтрация изображений
    – Алгоритмы сглаживания изображений
    – Двумерное сглаживание изображений
    – Обобщенная линейная фильтрация
    – Градиентный метод выделение контуров объектов на цветных изображениях
    – Пространственная фильтрация
  • Деконволюция
    – Предварительная обработка изображений
  • Расширение границ изображений. Сверхразрешение
  • – Реконструкция размытых изображений в MATLABСтруктурное распознавание на основе меры схожести символьных строк
  • Границы изображений
    – Края и их обнаружение
  • Оптимизация палитры изображений
  • Кодирование и сжатие изображений
    – задачи кодирования
    – основные методы кодирования
  • Некоторые области практического применения методов обработки изображений и распознавания образов (геофизические наблюдения, применение в биологии, применение на транспорте)
    – Распознавание рукописных знаков
    – Коррекция неравномерной засветки изображения.
    – Сегментация цветных изображений на основе кластеризации по методу k-средних.
    – Сегментация цветных изображений на основе цветового пространства L*a*b*
    – Уменьшение количества градаций цветных изображений
    – Обнаружение вращений и масштабных искажений на изображении
    – Регистрация изображений с помощью нормированной кросс-корреляции
    – Наложение двух изображений 
    – Технология повышения контрастности изображений.
    – Улучшение мультиспектральных цветных изображений
    – Регистрация аэрофотографий на ортофотоснимках
    – Пространственные преобразования изображений
    – Исследование конформных преобразований
    – Создание обивочных материалов с использованием изображений
    – Извлечение данных из трехмерных магниторезонансных изображений 
    – Поиск длины маятника в движении
    – Гранулометрия
    – Идентификация округлых предметов
    – Измерение углов пересечения
    – Измерение радиуса части мотка ленты
    – Обнаружение объектов с помощью сегментации изображений
    – Реконструкция изображений по их проекционным данным
    – Сегментация методом управляемого водораздела
    – Восстановление изображений методом слепой деконволюции
    – Реконструкция изображений с использованием регуляризационного фильтра
    – Восстановление изображений с использованием метода Лаки-Ричардсона
    – Некоторые подходы к улучшению визуального качества изображений с затемненными участками
    – Реализация некоторых методов видоизменения гистограмм в системе Matlab
    – Подавление шумов на изображениях
    – Текстурная сегментация с использованием текстурных фильтров
    – Поиск растительности на мультиспектральных изображениях
    – Распознавание объектов на основе вычисления их признаков
    – Распознавание объектов на основе вычисления коэффициента корреляции
    – Анализ признаков объектов
    – Некоторые аспекты задачи распознавания номерных знаков автомобилей
    – Анализ серии изображений с распределенной обработкой данных
    – Сглаживание цветных изображений
    – Построение гистограмм
    – Видоизменение гистограмм
    Визуализация объектов
  • Применение методов улучшения изображений при разработке системы видеонаблюдения
  • Улучшение изображений с яркостными искажениями
  • Некоторые алгоритмы повышения визуального качества изображений
  • Пороговая обработка цветных изображений
  • Формирование ночного изображения на основе дневного и наоборот
  • Обнаружение лиц на основе цвета
  • Метод управления яркостью изображения

Наверх

Свойства зрительной системы человека

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

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

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

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

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

Зрительная система по-разному реагирует на излучения, которые равны по мощности, но излучаемые из различных диапазонов спектра. Такая чувствительность называется спектральной.

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

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

Цветоощущение характеризуется тремя основными характеристиками – светлота, цветовой тон и насыщенность. Для классификации цветов используются цветовые пространства.

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

Наверх

Возможности цифровой обработки изображений в Matlab

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

Геометрические преобразования изображений

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

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

L=imread('cameraman.tif');
imshow(L);
imcrop;

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

В пакете Image Processing Toolbox существует функция imrotate, которая осуществляет поворот изображения на заданный угол.

L1=imrotate(L,35,'bicubic');
figure,imshow(L1)

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

Анализ изображений

Для работы с отдельными элементами изображений используются такие функции как imhist, impixel, mean2, corr2 и другие.

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

L=imread('cameraman.tif');
figure, imshow(L);
figure, imhist(L);

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

>> impixel

ans =

   173   173   173
   169   169   169
   163   163   163
    39    39    39

Следует отметить, что функция impixel по своим возможностям в некоторой степени повторяет опцию Data Cursor, пример использования которой при веден на изображении внизу.

Еще одной широко применяемой функцией является функция mean2 – она вычисляет среднее значение элементов матрицы.

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

Улучшение изображений

Среди встроенных функций, которые реализуют наиболее известные методы улучшения изображений, выделим следующие – histeq, imadjust та imfilter(fspecial).

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

L=imread('cameraman.tif');
figure, imshow(L);
L1=histeq(L);
Figure, imsow(L1);


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


Изображение после эквализации гистограммы

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

L=imread('cameraman.tif');
figure, imshow(L);
L1=imadjust(L);
figure, imshow(L1);
figure, imhist(L);
figure, imhist(L1);

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

L=imread('cameraman.tif');
figure, imshow(L);
H = fspecial('unsharp');
L1 = imfilter(L,H,'replicate');
figure, imshow(L1);

Фильтрация изображений

Пакет Image Processing Toolbox обладает очень мощным инструментарием по фильтрации изображений. Среди множества встроенных функций, которые решают задачи фильтрации изображений, особо следует выделить fspecial, ordfilt2, medfilt2.

Функция fspecial является функцией задания маски предопределенного фильтра. Эта функция позволяет формировать маски:

  1. высокочастотного фильтра Лапласа;
  2. фильтра, аналогичного последовательному применению фильтров Гаусса и Лапласа, так называемого лапласиана-гауссиана;
  3. усредняющего низкочастотного фильтра;
  4. фильтра, повышающего резкость изображения.

Рассмотрим примеры применения названных выше фильтров:

L=imread('cameraman.tif');
figure, imshow(L);
h=fspecial('laplasian',.5);
L1 = imfilter(L,h,'replicate');
figure, imshow(L1);

h=fspecial('log', 3, .5);
L1 = imfilter(L,h,'replicate');
figure, imshow(L1);

h=fspecial('average', 3);
L1 = imfilter(L,h,'replicate');
figure, imshow(L1);

h=fspecial('unsharp', .5);
L1 = imfilter(L,h,'replicate');
figure, imshow(L1);

Сегментация изображений

Среди встроенных функций пакета Image Processing Toolbox, которые применяются при решении задач сегментации изображений, следует выделить qtdecomp, edge и roicolor.

Функция qtdecomp выполняет сегментацию изображения методом разделения и анализа однородности не перекрывающихся блоков изображения.

I = imread('cameraman.tif');
S = qtdecomp(I,.27);
blocks = repmat(uint8(0),size(S));
 
for dim = [512 256 128 64 32 16 8 4 2 1];    
  numblocks = length(find(S==dim));    
  if (numblocks > 0)        
    values = repmat(uint8(1),[dim dim numblocks]);
    values(2:dim,2:dim,:) = 0;
       blocks = qtsetblk(blocks,S,dim,values);
  end
end
 
blocks(end,1:end) = 1;
blocks(1:end,end) = 1;
 
imshow(I), figure, imshow(blocks,[])

Одной из наиболее часто применяемых является функция выделения границ edge, которая реализует такие встроенные методы – Собела, Превит, Робертса, лапласиан-гауссиана, Канни и др.

Приведем примеры реализации функции edge с использованием различных фильтров.

clear;
I = imread('cameraman.tif');
BW1=edge(I,'sobel');
figure,imshow(BW1);title('sobel');
BW2=edge(I,'prewitt');
figure,imshow(BW2);title('prewitt');
BW3=edge(I,'roberts');
figure,imshow(BW3);title('roberts');
BW4=edge(I,'log');
figure,imshow(BW4);title('log');
BW5=edge(I,'zerocross');
figure,imshow(BW5);title('zerocross');
BW6=edge(I,'canny');
figure,imshow(BW6);title('canny');

Еще одной функцией, которая применяется для бинаризации по заданным цветам, является функция roicolor. Приведем пример ее использования.

I = imread('cameraman.tif');
figure, imshow(I);
BW = roicolor(I,128,255);
imshow(I);
figure, imshow(BW);

Морфологические операции над бинарными изображениями

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

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

clear;
L=imread('circles.bmp');
L=L(:,:,1);
imshow(L);
BW1=L<150;
figure,imshow(BW1);
BW2=bwmorph(BW1,'erode',12);
figure,imshow(BW2);
BW2=bwmorph(BW2,'thicken',Inf);
figure,imshow(BW2);
BW1=BW1&BW2;
figure,imshow(BW1);

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

Кроме перечисленных возможностей, пакет прикладных программ Image Processing Toolbox владеет широкими возможностями при решении задач анализа изображений, в частности, при поиске объектов и вычислении их признаков.

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

Наверх

Типы изображений

Интегрированные среды для моделирования и исполнения программ цифровой обработки изображений и сигналов содержат мощные средства для инженерно–научных расчетов и визуализации данных. Большинство современных пакетов поддерживает визуальное программирование на основе блок–схем. Это позволяет создавать программы специалистам, не владеющим техникой программирования. К таким пакетам относится Image Processing Toolbox системы MATLAB, разработанный фирмой MathWorks. Этот пакет владеет мощными средствами для обработки изображений. Они имеют открытую архитектуру и позволяют организовывать взаимодействие с аппаратурой цифровой обработки сигналов, а также подключать стандартные драйвера.

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

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

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

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

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

Полутоновое

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

В палитровых

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

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

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

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

Тип изображения double uint8
Бинарное 0 и 1 0 и 1
Полутоновое [0, 1] [0, 255]
Палитровое [1, размер палитры],

где 1 – первая строка палитры

[0, 255],

где 0 – первая строка палитры.*

Полноцветное [0, 1] [0, 255]

*

Примечание. При программной реализации лучше избегать использования такой индексации строк. MATLAB корректно воспринимает индексацию с первой, а не нулевой строки.

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

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

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

Наверх

Бинарные изображения: геометрические характеристики

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

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

Бинарные изображения

Начнем со случая, когда в поле зрения находится объект, а все остальное считается “фоном”. Если объект оказывается заметно темнее (или светлее), чем фон, то легко определить характеристическую функцию , которая равна нулю для всех точек изображения, соответствующих фону, и единице для точек на объекте (рис.1) или наоборот.

Рис. 1. Бинарное изображение, определяемое характеристической функцией , которая принимает значение “нуль” и “единица”.

Часто бинарное изображение получают пороговым разделением обычного изображения. К нему также можно прийти путем порогового разделения расстояния на “изображении”, полученном на основе измерений расстояний.

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

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

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

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

Пример:

Этой операции соответствует функция BWEULER – вычисление чисел Эйлера в пакете Image Processing Toolbox:
L=imread(‘test.bmp’);
L=double(L);
imshow(L);

e=bweuler(L(:,:,1), 4)
e =
1; % На объекте действительно одно отверстие.

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

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

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

Простые геометрические характеристики

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

,

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

Пример:

В системе Matlab этой операции соответствует функция BWAREA – вычисление площади объектов.
L=imread(‘test.bmp’);
L=double(L);
imshow(L);

S=bwarea(L(:,:,1))
e =
24926; %Площадь объекта в пикселах (размер изображения 236х236).

Площадь и положение

Как определить положение объекта на изображении? Поскольку объект, как правило, состоит не из одной единственной точки, мы должны четко определить смысл термина “положение”. Обычно в качестве характерной точки объекта выбирают его геометрический центр (рис. 2).

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

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

,

а относительно оси по формуле

,

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

Ориентация

Мы также хотим определить, как расположен объект в поле зрения, т. е. его ориентацию. Сделать это несколько сложнее. Допустим, что объект немного вытянут вдоль некоторой оси; тогда ее ориентацию можно принять за ориентацию объекта. Как точно определить ось, вдоль которой вытянут объект? Обычно выбирают ось минимального второго момента. Она представляет собой двумерный аналог оси наименьшей инерции. Нам необходимо найти прямую, для которой интеграл от квадратов расстояний до точек объекта минимален; этот интеграл имеет вид

,

где — расстояние вдоль перпендикуляра от точки с координатами до искомой прямой.

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

(В пакете Matlab операции определения центра масс, ориентации а также другие морфометрические признаки вычисляются с помощью функции IMFEATURE.)

Проекции

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

Дискретные бинарные изображения

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

,

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

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

Кодирование с переменной длиной

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

0

1

1

1

1

0

0

0

1

1

0

0

0

0

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

Другой подход к обработке изображений описан в книге [2]. Книга [3] содержит сведения о некоторых работах в области обработки бинарных изображений. Много интересных бинарных изображений, полученных художниками, можно найти в книге [4].

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

В системе Matlab также рассматривается один из видов кодирования, который содержится в описании функции BWPACK.

Литература

  • Хорн Б.К.П. Зрение роботов: Пер. с англ. – М.: Мир, 1989. – 487 с., ил. ISBN 5–03–000570–6.
  • Rosenfeld A., Kak A.C., Digital Picture Processing, Vols. 1, 2, Second Edition, Academic Press, New York, 1982.
  • Stoffel J.C. (ed.), Graphical and Binary Image Processing and Applications, Artech House, Inc., Massachusetts, 1982.
  • Grafton C.B. (ed.), Silhouettes – A Pictorial Archive of Varied Illustrations, Dover Publications, New York, 1979.
  • Mitchell J.L., Goertzel G., Two–Dimensional Facsimile Coding Scheme, IBM Reserch Report RC 7499, Jan., 1979.

Наверх

Бинарные изображения: топологические характеристики

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

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

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

Сложные объекты

Иногда в поле зрения попадает более одного объекта (рис. 1). В этом случае вычисление площади, геометрического центра и ориентации приведет к значениям, “усредненным” по всем компонентам бинарного изображения. Как правило, это не то, что требуется. Желательно как-то пометить отдельные компоненты изображения и вычислить значения площади, первых и вторых моментов для каждой компоненты в отдельности.

Разметка компонент

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

Рис. 1. Изображение, состоящее из нескольких областей, для каждой из которых необходимо проводить расчет положения и ориентации.

Элементы изображения необходимо пометить таким образом, чтобы принадлежащие одной области были отличимы от остальных Для этого нам необходимо решить, какие точки принадлежат одной и той же области. На рисунке 1 точка А считается связанной с точкой В, поскольку мы можем найти непрерывную кривую, целиком лежащую в затененной области и соединяющую указанные точки. Ясно, что точка А не связана с точкой С, так как в этом случае подобной кривой найти нельзя.

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

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

Связность

Теперь нужно аккуратно рассмотреть смысл термина сосед. Если мы имеем дело с квадратным растром, то, по-видимому, соседями следует считать четыре элемента изображения, касающиеся сторон данного элемента. Но как быть с теми, которые касаются его в углах? Существуют две возможности:

— четырехсвязность – соседями считаются только элементы, примыкающие к сторонам;

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

Оказывается, ни одна из этих схем не является полностью удовлетворительной. В этом можно убедиться, если вспомнить, что фон также можно разбить на несколько связных компонент. Здесь нам хотелось бы применить наши интуитивные представления о связности областей на непрерывном бинарном изображении. Так, например, простая замкнутая кривая должна разделять изображение на две связные области (рис. 2). Это так называемая теорема Жордана о кривой.

Рис. 2. Простая замкнутая кривая, разбивающая плоскость на две связные области.

Теперь рассмотрим простое изображение, содержащее четыре элемента со значением “единица”, которые примыкают к центральному элементу со значением “нуль”:

Это — крест с выброшенным центром. Если принять соглашение о четырехсвязности, то на изображении окажутся четыре различные компоненты (,, и ):

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

Итак, мы получили замкнутую кривую и только одну связную компоненту фона.

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

или

Для обеспечения симметричности отношения связности два угловых элемента должны находиться на одной и той же диагонали: если элемент А — сосед элемента В, то элемент В должен быть соседом элемента А. В дальнейшем мы будем пользоваться первым из двух возможных вариантов, приведенных выше, считая соседями элементы в направлениях N, E, SE, S, W и NW. С помощью шестисвязности как объект, так и фон можно трактовать единообразно без каких-либо дальнейших неувязок. Для изображений на квадратном растре мы примем именно такое соглашение.

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

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

Локальные вычисления и итеративная модификация

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

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

Рис. 5. Методы локальных вычислений. Комбинируются выходы отдельных вычислительных ячеек, каждая из которых соединена с несколькими элементами изображения, лежащими вблизи нее.

Локальные вычисления.

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

Какие другие характеристики представимы в виде суммы результатов локальных операций? Например, периметр: достаточно просто подсчитать количество участков на изображении, где рядом с нулями стоят единицы. Имеются два типа локальных операторов (рис. 6): операторы одного типа просматривают два соседних элемента, расположенных в одной строке, а операторы другого типа — два соседних элемента, расположенных в одном столбце. В обоих случаях результат есть ИСКЛЮЧАЮЩЕЕ ИЛИ (аÄb) двух значений на входе. Сумма всех получаемых выходов представляет собой оценку периметра.

Рис. 6. Возможность использования операции ИСКЛЮЧАЮЩЕЕ ИЛИ к двум соседним элементам изображения для выделения участков, находящихся на границе областей.

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

Вычисленный периметр представляет собой лишь приблизительную оценку, поскольку, как правило, дискретное бинарное изображение строится на основе непрерывного, и при этом границы объектов становятся более изрезанными. Например, оценка длины диагональной прямой в раз больше “истинной”:

Усреднение по всем углам наклона дает среднее значение коэффициента, показывающего, во сколько раз увеличено полученное значение. Оно составляет 4/=1,273…. Разделив на это число, можно улучшить оценку периметра.

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

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

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

Итеративная модификация

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

В пакете Image Processing Toolbox системы Matlab существует много функций, осуществляющих обработку бинарных изображений, в частности, морфологические операции. Среди них – BWMORPH, DILATE, ERODE, BWPERIM, MAKELUT, BWFILL, BWSELECT, IMFEATURE и другие.

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

  • Arcelli C., Pattern Thinning by Contour Tracing, Computer Graphics and Image Processing, 17, № 3, 130 – 144 (1981).
  • Dyer C.R., Rosenfeld A., Thinning Algorithms for Gray–Scale Picture, IEEE Trans. on Pattern Analysis and Machine Intelligence, 1, №1, 88 – 89 (1979).
  • Stefanelli R., Some Parallel Thinning Algorithms for Digital Pictures, Jornal of the ACM, 18, № 2, 255–264 (1971).
  • Хорн Б.К.П. Зрение роботов: Пер. с англ. – М.: Мир, 1989. – 487 с., ил. ISBN 5–03–000570–6.

Наверх

Формирование и обработка цифровых изображений

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

Методы получения цифровых изображений

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

Системы получения рентгенографических изображений

Рентгеновское излучение активно используется для получения изображений с момента его открытия в 1895 г. Изображение формируется в результате взаимодействия квантов рентгеновского излучения с приемником и представляет собой распределение квантов, которые прошли через объект диагностики и были зарегистрированы детектором (рис. 1). Последние делятся на первичные, т.е. те, которые прошли через объект

s6.gif (3816 bytes)

Рис. 1. Компоненты системы для получения рентгеновских изображений. B и E – кванты, которые прошли через исследуемый объект без взаимодействия; C и D – рассеянные кванты. Квант D отсеивается сеткой, которая препятствует рассеянному излучению, а квант A – поглощается в объекте.

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

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

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

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

Блок-схема типичной цифровой рентгенографической системы представлена на рис. 2.

s7.gif (4391 bytes)

Рис. 2. Элементы цифровой системы получения рентгеновских изображений.

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

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

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

К преимуществам цифровых рентгенографических систем относятся следующие факторы: цифровое отображение информации; низкая доза облучения; цифровая обработка изображений и улучшения качества. Рассмотрим эти преимущества более подробно [1-3].

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

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

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

Получение изображений с помощью радиоизотопов

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

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

Ультразвуковая диагностика

Ультразвуковые методы визуализации широко применяются при разных диапазонах частот – от подводной локации и биоэхолокации (частоты до 300 КГц) до акустической микроскопии (от 12 МГц до 1ГГц и выше). Промежуточное расположение по частотам занимают ультразвуковая диагностика и терапия, а также неразрушающий контроль в промышленности. Информация о структуре исследуемого объекта закодирована в лучах, которые прошли через него и в рассеянном излучении. Задача системы визуализации состоит в расшифровке этой информации. В отличие от рентгеновских лучей, ультразвуковые волны преломляются и отбиваются на границах раздела сред с разными акустическими показателями преломления. Эти эффекты могут быть довольно заметными, что разрешает создать фокусирующие системы.

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

Использование эффекта ядерного магнитного резонанса (ЯМР) для получения изображений

Несмотря на то, что во многих больших исследовательских центрах ЯМР-визуализация является одним из важных диагностических средств, сам метод еще находится на относительно ранней стадии своего развития. Само явление ядерного магнитного резонанса было открыто в 1946 году независимо Блохом и Парселлом с Паундом. Этот метод с помощью небольших изменений резонансной частоты (через наличие околомолекулярной электронной тучи) позволяет идентифицировать ядра в разном химическом окружении. Сначала ЯМР-методы с высокой разрешающей способностью разрабатывались как универсальное средство изучения химического состава и структуры твердого тела и жидкостей, а далее нашли свое применение и в других областях, в частности, медицине. Рядом с развитием ЯМР-спектроскопии развивались и методы визуализации – это и точечные методы, методы “быстрой” визуализации и прочие. Роль центрального процессора в современных ЯМР-системах выполняет мощный миникомпьютер, который обеспечивает канал связи с оператором и контроль функций узлов системы. Компьютер также обеспечивает запоминание и архивирование информации, отображение результатов исследований и во многих случаях соединяется с устройством быстрой обработки типа матричного процессора.

Пример обработки рентгеновских биомедицинских изображений с использованием системы MATLAB

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

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

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

%Пример программы обработки биомедицинских изображений в среде MATLAB
L=imread(‘cardial.bmp’);
figure, imshow(L);
L1=imadjust(L,[min(min(L)) max(max(L))]/255,[],1);
figure, imshow(L1);
L=L1(:,:,1);
L=double(L);
Filter=[1 1 1,1 1 1,1 1 1];
Lser=filter2(Filter,L(:,:,1),’same’)./9;
C=abs(L(:,:,1)-Lser)./(L(:,:,1)+Lser);
C=C.^.55;
[N M]=size(L);
for i=1:N;
disp(i);
for j=1:M;
if L(i,j,1)>Lser(i,j);
Lvyh(i,j)=Lser(i,j)*(1+C(i,j))/(1-C(i,j));
else
Lvyh(i,j)=Lser(i,j)*(1-C(i,j))/(1+C(i,j));
end;
end;
end;
figure, imshow(Lvyh/255);

Рис. 3.

Литература

  • Физика визуализации изображений в медицине: в 2–х томах. Т. 2: Пер. С англ. / Под ред. С. Уэбба. – М.: Мир, 1991. – 408 с., ил.
  • Беликова Т.П. Моделирование линейных фильтров для обработки рентгеновских изображений в задачах медицинской диагностики // Цифровая оптика. Обработка изображений и полей в экспериментальных исследованиях / Под ред. В.И.Сифорова и Л.П.Ярославского. – М.: Наука, 1990. – 176 с.
  • Н.Н. Блинов, Е.М. Жуков, Э.Б. Козловский, А.И. Мазуров. Телевизионные методы обработки рентгеновских и гамма–изображений. М.: Энергоатоиздат, 1982. – 200 с.

Наверх

Адаптивное повышение контрастности изображений

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

– исходное изображение и его элемент с координатами соответственно;

– контраст элемента изображения с координатами ;

– преобразованное значение контраста ;

– характеристики локальных окрестностей ( – энтропия, – среднеквадратичное отклонение, – функция протяженности гистограммы);

– элемент обработанного изображения с координатами .

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

Основные шаги реализации методов адаптивного преобразования локальных контрастов такие:

Шаг 1. Для каждого элемента изображения вычисляют значение локального контраста в текущей окрестности с центром в элементе с координатами .

Шаг 2. Вычисляют локальную статистику для текущей скользящей окрестности .

Шаг 3. Преобразуют (усиливают) локальный контраст , употребляя для этого нелинейные функции и учитывая локальную статистику текущей скользящей окрестности .

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

Шаги 1 и 2 могут выполняться в различной последовательности или параллельно.

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

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

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

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

Основные шаги реализации этого метода такие.

Шаг 1. Вычисляем локальный контраст элемента.

Шаг 2. Определяем характеристику локальной скользящей окрестности, используя функцию протяженности гистограммы

, (1)

где , – соответственно максимальное и минимальное значения яркостей элементов скользящей окрестности с центром в элементе с координатами ; – максимальное значение гистограммы уровней яркости элементов окрестности с центром в элементе с координатами .

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

, (2)

где

,

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

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

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

Первый тип – это однородный участок изображения, который характеризуется примерно одинаковыми уровнями яркостей элементов; гистограмма такой окрестности показана на рис. 2.

С рис. 2 видно, что , а следовательно , согласно выражению (1), функция протяженности гистограммы локальной окрестности будет равна нулю.

Рис. 2. Гистограмма распределения яркостей элементов однородной окрестности.

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

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

Рис. 3. Гистограмма распределения яркостей элементов бинарной окрестности.

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

, (3)

где и – размеры скользящей окрестности , выражение (1) будет иметь вид

, (4)

Если , , а размеры локальной окрестности такие, что допускают присутствие элементов со всеми возможными уровнями яркостей [0,255], например элементов, тогда функция протяженности гистограммы в соответствии с выражением (4) примет значение .

Третьим характерным типом возможной локальной окрестности является такая окрестность, где в примерно одинаковой мере присутствуют элементы со всеми возможными яркостями с диапазона [0,255]. Такие окрестности характеризуются гистограммой равномерного распределения яркостей, которая показана на рис. 4. Тогда согласно изложенных предположений относительно размера локальной окрестности и характера его гистограммы получим, что , . В этом случае функция протяженности гистограммы примет значение . Для такой окрестности будем считать, что она высококонтрастна и не нуждается в усилении контраста.

Рис. 4. Гистограмма скользящей окрестности с равномерно распределенными яркостями элементов.

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

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

Рис. 5. Зависимость показателя степени преобразования локального контраста от функции протяженности гистограммы : 1 – в известном подходе [1], 2 – в предложенном методе.

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

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

Выражение для определения (рис. 5 , кривая 2) такое:

, (5)

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

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

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

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

%=======Считывание данных======
clear;
L=imread('test.bmp');%Исходное изображение полутоновое, поэтому  L(:,:,1)=L(:,:,2)=L(:,:,3);
L=L(:,:,1);
L=im2double(L);
m=15;n=m;n1=fix(n/2);m1=fix(m/2); %Определение размеров локальных окрестностей

%=======Преобразование матрицы яркостей изображения для устранения краевого эффекта=======
%=======В новых версиях системы Matlab существуют функции, которые реализуют эту процедуру=======
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;
  for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
  end;
end;
   L2=L(1,1:M);   L02=L2;
    for i=1:n1-1;
      L2=[L2;L02];
    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;
          L7=[L7;L07];
        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;
            L4=[L4;L04];
          end;
       L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;
      L5=[L5;L05];
    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';  L2=[L2;L];  L2=[L2;L7];  L2=L2'; 
 L3=[L3;L5]; L3=[L3;L8];  L3=L3';  L1=[L1;L2];  L1=[L1;L3];
  Lr=L1';
clear L2;clear L3;clear L4;clear L5;clear L6;
clear L7;clear L8;clear L02;clear L04;clear L05;
clear L07;clear L1;clear L;

%=======Определение параметров локальной окрестности (функции протяженности гистограммы)=======
HP=zeros(N+2*n1,M+2*m1);
 for i=1+n1:N+n1;
 disp(i)
 for j=1+m1:M+m1;
                 if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=Lr(i+a,j+b);
                        end;
                        end;
                 end;
           if j>1+m1;
            for a=-n1:n1;
              D(n1+1+a,m+1)=Lr(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
          end;            
               LMIN=min(min(D));               LMAX=max(max(D));
               H_lokal=hist(D(:)+1,LMAX-LMIN+1);
               H_lokal_max=max(H_lokal);
               clear H_lokal;
               HP(i,j)=(LMAX-LMIN)/H_lokal_max;
               clear LMIN;               clear LMAX;               clear H_lokal_max;
 end;
 end;
n_filter=3;m_filter=n_filter;
F=ones(n_filter,m_filter);
Lser=filter2(F,Lroshyrena,'same')/(n_filter*m_filter);
clear n_filter;clear m_filter;
amax=.7;amin=.5;

%=======Определение и преобразование локального контраста с учетом локальных характеристик=======
C=(Lr-Lser)./(Lr+Lser+eps);
C=abs(C);
for i=1+n1:N+n1;
 disp(i)
for j=1+m1:M+m1; 
                 if j==1+m1;
                        TM=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           TM(n1+1+a,m1+1+b)=HP(i+a,j+b);
                        end;
                        end;
                 end;
         if j>1+m1;
            for a=-n1:n1;
              TM(n1+1+a,m+1)=HP(i+a,j+m1);
            end;
             TM=TM(1:n,2:m+1);
         end;    
    HP_MIN=min(min(TM));
    HP_MAX=max(max(TM));
            C(i,j)=C(i,j)^(amin+(amax-amin)*(HP(i,j)-HP_MIN)/(HP_MAX-HP_MIN));
 if Lroshyrena(i,j)>Lser(i,j);    
      Lvyh(i,j)=Lser(i,j)*(1+C(i,j))/(1-C(i,j));
 else
      Lvyh(i,j)=Lser(i,j)*(1-C(i,j))/(1+C(i,j));
 end; 
       if Lvyh(i,j)>=255;
          Lvyh(i,j)=255;
       end;
        if Lvyh(i,j)<=0;
          Lvyh(i,j)=0;
        end;
end;
end;
Lvyh=round(Lvyh);
Lvyh=Lvyh(1+n1:N+n1,1+m1:M+m1);
L=Lr(1+n1:N+n1,1+m1:M+m1);

%=======Визуализация=======
colormap(gray(255));
subplot(221);image(L');axis('image');
subplot(222);image(Lvyh');axis('image');   

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

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

Литература.

  1. Dhawan A.P., Buelloni G., Gordon R. Enhancement of mammographic features by optimal adaptive neighbourhvod image processing // IEEE Trans. Med. Imaging. – 1986. – v.5. – P.8-15.
  2. 2. Gordon R., Rangayyan R.M. Feature enhancement of film mammograms using fixed and adaptive neighbourhood // Applied optics. – 1984. – v.23. – P. 560-564.

Наверх

Энтропия изображения. Использование среднеквадратического отклонения значений яркостей элементов окрестности в методах контрастирования изображений. Нелинейное растяжение локальных контрастов.

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

Шаг 1. Вычисляем локальный контраст элемента изображения .

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

, (1)

где

. (2)

Шаг 3. Вычисляем степенное преобразование локального контраста, которое в связи с использованием локальной энтропии приобретает адаптивный характер:

, (3)

где , – максимальное и минимальное значения энтропии скользящей окрестности размером элементов.

Шаг 4. Восстанавливаем изображение за выражением, которое определяется из выражения определения локального контраста.

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

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

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

, (4)

где – значение гистограммы для элемента с значением яркости .

Кроме этого нами предложено выражение преобразования локального контраста для модификации шага 3 известного подхода

, (5)

где – параметр нелинейного усиления контраста.

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

Использование среднеквадратического отклонения значений яркостей элементов окрестности в методах контрастирования изображений

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

Определим величину показателя степени так:

, (6)

где – нормирующий коэффициент, , – среднеарифметическое значение яркости исходного изображения

, (7)

, – размеры изображения ( , ),

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

. (8)

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

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

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

Нелинейное растяжение локальных контрастов

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

Шаг 1. Определение количественной меры локального контраста.

Шаг 2. Увеличение по определенному закону некоторой количественной меры локального контраста.

Шаг 3. Восстановление элемента преобразованного изображения с усиленным локальным контрастом.

Такие преобразования выполняются для каждого элемента изображения.

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

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

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

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

Рис. 1. Типичная гистограмма распределения значений локальных контрастов исходного изображения.

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

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

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

(9)

где – значение локального контраста элемента исходного изображения с координатами ,

– усиленное значение локального контраста элемента изображения с координатами ;

– максимально возможное значение локального контраста ;

, – максимальное и минимальное значения локального контраста исходного изображения,

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

, – коэффициенты постоянного смещения;

– показатель степени .

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

Рис. 3. Графики функций преобразования локального контраста (1 – степенная функция, 2 – за выражением (9)).

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

%Программа, реализующая метод повышения контрастности изображений
%с использованием энтропии
clear;
%=====Считывание исходных данных=====
L=imread('image1.tif');
[N M]=size(L);
m=15;n=15;n1=fix(n/2);m1=fix(m/2);
%=====Преобразование матрицы элементов изображения для избежания краевого эффекта=====
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;  for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
  end;end;
   L2=L(1,1:M);   L02=L2;
    for i=1:n1-1;
      L2=[L2;L02];
    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;
          L7=[L7;L07];
        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;
            L4=[L4;L04];
          end;
       L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;
      L5=[L5;L05];
    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';  L2=[L2;L];  L2=[L2;L7];  L2=L2';
  L3=[L3;L5];  L3=[L3;L8];  L3=L3';  L1=[L1;L2];  L1=[L1;L3];
  Lr=L1';
clear L2;clear L3;clear L4;clear L5;clear L6;clear L7;clear L8;clear L02;
clear L04;clear L05;clear L07;clear L1;clear L;
%=====Определение параметров локальной окрестности=====
Entropia=ones(N+2*n1,M+2*m1);
      for i=1+n1:N+n1;
           disp(i)
      for j=1+n1:M+m1;
                 if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=Lr(i+a,j+b);
                        end;
                        end;
                 end;
         if j>1+m1;
            for a=-n1:n1;
              D(n1+1+a,m+1)=Lr(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
         end; 
              Pk_vektor=hist(D(:),max(max(D))-min(min(D))+1)/(n*m);
               for ind=1:length(Pk_vektor);
                   if Pk_vektor(ind)==0;
                      Pk_vektor(ind)=1;
                   end;
               end;
              Entropia(i,j)=-sum(Pk_vektor.*log(Pk_vektor));        
     end;
     end;
n_filter=3;m_filter=n_filter;F=ones(n_filter,m_filter);
Lser=filter2(F,Lr,'same')/(n_filter*m_filter);
clear n_filter;clear m_filter;clear F;
amax=.7;amin=.5;
%=====Определение локальных контрастов и функции их преобразования=====
C=abs(Lr-Lser)./(Lr+Lser+eps);
alfa=amin+.2.*(1-exp(-((Entropia./max(max(Entropia))-.5).^2) ./(2*.14^2))).^3;
C=C.^alfa;
for i=1+n1:N+n1;
for j=1+n1:M+m1;
                if   Lr(i,j)>=Lser(i,j);
                     Lvyh(i,j)=Lser(i,j)*(1+C(i,j))/(1-C(i,j)+eps);
                else
                     Lvyh(i,j)=Lser(i,j)*(1-C(i,j))/(1+C(i,j));
                end;             
%=====Проверка корректности диапазона=====
   if Lvyh(i,j)>=255;      Lvyh(i,j)=255;   end;
   if Lvyh(i,j)<=0;      Lvyh(i,j)=0;   end;
end;
end;
Lvyh=Lvyh(n1+1:N+n1,m1+1:M+m1);
Lvyh=round(Lvyh);
L=Lr(n1+1:N+n1,m1+1:M+m1);
%=====Визуализация результатов обработки=====
colormap(gray(255));
subplot(221);image(L');axis('image');  
subplot(222);image(Lvyh');axis('image');
Результаты работы программы представлены на рис.

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

%Программа, реализующая метод контрастирования изображений с использованием
%среднеквадратических отклонений значений яркостей элементов локальной окрестности.
clear;
%=====Считывание исходных данных=====
L=imread('image2.tif');
[N M]=size(L);
m=15;n=15;n1=fix(n/2);m1=fix(m/2);
%=====Преобразование матрицы элементов изображения для избежания краевого эффекта=====
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;  for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
  end;end;
   L2=L(1,1:M);   L02=L2;
    for i=1:n1-1;
      L2=[L2;L02];
    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;
          L7=[L7;L07];
        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;
            L4=[L4;L04];
          end;
       L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;
      L5=[L5;L05];
    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';
  L2=[L2;L];  L2=[L2;L7];  L2=L2';
  L3=[L3;L5];  L3=[L3;L8];  L3=L3';  L1=[L1;L2];
  L1=[L1;L3]; Lr=L1';
clear L2;clear L3;clear L4;clear L5;clear L6;clear L7;clear L8;
clear L02;clear L04;clear L05;clear L07;clear L1;clear L;
%====Определение среднеквадратических отклонений значений яркостей элементов локальной окрестности====
SV=zeros(N+2*n1,M+2*m1);
for i=1+n1:N+n1;
for j=1+m1:M+m1;
                 if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=Lr(i+a,j+b);
                        end;
                        end;
                 end;
           if j>1+m1;
            for a=-n1:n1;
              D(n1+1+a,m+1)=Lr(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
          end;            
               D=D(:);
               SV(i,j)=std(D);
 end;
 end;
%=====Определение массива усредненных значений элементов изображения=====
n_filter=3;m_filter=n_filter;
F=ones(n_filter,m_filter);
Lser=filter2(F,Lr,'same')/(n_filter*m_filter);
amin=.35;amax=.95;
C=(Lr-Lser)./(Lr+Lser+eps);
C=abs(C);
for i=1+n1:N+n1;
 disp(i)
for j=1+m1:M+m1;
      if SV(i,j)<=1;
        SV(i,j)=1;
      end;
                 if j==1+m1;
                        TM=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           TM(n1+1+a,m1+1+b)=SV(i+a,j+b);
                        end;
                        end;
                 end;
         if j>1+m1;
            for a=-n1:n1;
              TM(n1+1+a,m+1)=SV(i+a,j+m1);
            end;
             TM=TM(1:n,2:m+1);
         end;    
    SV_MIN=min(min(TM));
    SV_MAX=max(max(TM));
            C(i,j)=C(i,j)^(amin+(amax-amin)*(SV(i,j)-SV_MIN)/(SV_MAX-SV_MIN));
    if Lroshyrena(i,j)>=Lser(i,j);
      Lvyh(i,j)=Lser(i,j)*(1+C(i,j))/(1-C(i,j));
    else
      Lvyh(i,j)=Lser(i,j)*(1-C(i,j))/(1+C(i,j));
    end;
%=====Проверка корректности диапазона=====
   if Lvyh(i,j)>=255;      Lvyh(i,j)=255;   end;
   if Lvyh(i,j)<=0;      Lvyh(i,j)=0;   end;
end;
end;
Lvyh=Lvyh(n1+1:N+n1,m1+1:M+m1);
Lvyh=round(Lvyh);
L=Lroshyrena(n1+1:N+n1,m1+1:M+m1);
%=====Визуализация результатов=====
colormap(gray(255));
subplot(221);image(L);axis('image');
subplot(222);image(Lvyh);axis('image');
Результаты работы программы представлены на рис.

Рис. Локально-адаптивний метод улучшения визуального качества изображений с использованием среднеквадратических отклонений значений яркостей элементов локальной окрестности: а) исходное изображение; б) изображение а, обработанное известным методом; в) изображение а, обработанное предложенным методом.

%Программа, реализующая метод нелинейного (степенного) преобразования (растяжения)
%локальных контрастов (глобальная реализация)
clear;
%=====Считывание исходных данных=====
L=imread('im.tif');
[N M]=size(L);
%=====Задание параметров локальных окрестностей=====
m=5;n=m;n1=fix(n/2);m1=fix(m/2);
R=255;
%=====Преобразование матрицы интенсивностей с целью избежания краевого эффекта=====
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;
  for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
  end;
end;
   L2=L(1,1:M);   L02=L2;
    for i=1:n1-1;
      L2=[L2;L02];
    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;
          L7=[L7;L07];
        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;
            L4=[L4;L04];
          end;
       L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;
      L5=[L5;L05];
    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';
  L2=[L2;L];  L2=[L2;L7];  L2=L2';
  L3=[L3;L5];  L3=[L3;L8];  L3=L3';
  L1=[L1;L2];  L1=[L1;L3];
  Lr=L1';
clear L2;clear L3;clear L4;clear L5;clear L6;
clear L7;clear L8;clear L02;clear L04;
clear L05;clear L07;clear L1;clear L;
%=====Определение усредненного локального среднего изображения=====
F=ones(n,m);
Lser=filter2(F,Lr,'same')/(n*m);
С=abs(Lr-Lser)/(Lr+Lser);
CMIN=min(min(C));
CMAX=max(max(C));
alfa=.9; 
%=====Параметр нелинейного преобразования=====
%=====Нелинейное (при a>1 или a<1) растяжение локальных контрастов=====
C=((C-CMIN)./(CMAX-CMIN)).^alfa;
%=====Восстановление яркостей результирующего изображения=====
for i=1+n1:N+n1;
 disp(i)
for j=1+m1:M+m1;
    if Lr(i,j)>=Lser(i,j);
      Lvyh(i,j)=Lser(i,j)*(1+C(i,j))/(1-C(i,j));
    else
      Lvyh(i,j)=Lser(i,j)*(1-C(i,j))/(1+C(i,j));
    end;
   if Lvyh(i,j)>=R;      Lvyh(i,j)=R;   end;
   if Lvyh(i,j)<=0;      Lvyh(i,j)=0;   end;
end;
end;
Lvyh=Lvyh(n1+1:N+n1,m1+1:M+m1);
Lvyh=round(Lvyh);
L=Lroshyrena(n1+1:N+n1,m1+1:M+m1);
%=====Визуализация результатов=====
subplot(221);imshow(L');axis('image');
subplot(222);imshow(Lvyh');axis('image');

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

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

%Программа, реализующая метод нелинейного (степенного) преобразования (растяжения)
%локальных контрастов (скользящая реализация)
clear;
%=====Считывание исходных данных=====
L=imread('im1.bmp');
L=L(:,:,1);
L=double(L)./255;
[N M]=size(L);
%=====Задание параметров локальных окрестностей=====
m=17;n=m;n1=fix(n/2);m1=fix(m/2);
R=255/255;
%=====Преобразование матрицы интенсивностей с целью избежания краевого эффекта=====
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;
  for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
  end;
end;
   L2=L(1,1:M);   L02=L2;
    for i=1:n1-1;
      L2=[L2;L02];
    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;
          L7=[L7;L07];
        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;
            L4=[L4;L04];
          end;
       L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;
      L5=[L5;L05];
    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';
  L2=[L2;L];  L2=[L2;L7];  L2=L2';
  L3=[L3;L5];  L3=[L3;L8];  L3=L3';
  L1=[L1;L2];  L1=[L1;L3];
  Lr=L1';
clear L2;clear L3;clear L4;clear L5;clear L6;
clear L7;clear L8;clear L02;clear L04;
clear L05;clear L07;clear L1;clear L;
alfa=5; 
%=====Параметр нелинейного преобразования=====
%=====Восстановление яркостей результирующего изображения=====
C=zeros(N+n+1,M+m+1);
C1=zeros(N+n+1,M+m+1);
for i=1+n1:N+n1;
 disp(i)
for j=1+m1:M+m1;
                 if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=Lr(i+a,j+b);
                        end;
                        end;
                 end;
           if j>(1+m1);
            for a=-n1:n1;
              D(n1+1+a,m+1)=Lr(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
           end;            
                               Lser(i,j)=sum(sum(D))/(n*m);
                               C(i,j)=abs((Lr(i,j)-Lser(i,j))/(Lr(i,j)+Lser(i,j)));
 end;
 end;
for i=1+n1:N+n1;
 disp(i)
for j=1+m1:M+m1;
                 if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=C(i+a,j+b);
                        end;
                        end;
                 end;
           if j>(1+m1);
            for a=-n1:n1;
              D(n1+1+a,m+1)=C(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
           end;            
                               CMIN(i,j)=min(min(D));
                               CMAX(i,j)=max(max(D));
                                C1(i,j)=((C(i,j)-CMIN(i,j))/(CMAX(i,j)-CMIN(i,j)))^alfa;
   if C1(i,j)>1;      C1(i,j)=1;   end;
   if C1(i,j)<0;      C1(i,j)=0;   end;                                
    if Lr(i,j)>=Lser(i,j);
      Lvyh(i,j)=Lser(i,j)*(1+C1(i,j))/(1-C1(i,j)+eps);
    else
      Lvyh(i,j)=Lser(i,j)*(1-C1(i,j))/(1+C1(i,j));
    end;   
   if Lvyh(i,j)>R;      Lvyh(i,j)=R;  end;
   if Lvyh(i,j)<0;      Lvyh(i,j)=0;  end;
end;
end;
Lvyh=Lvyh(n1+1:N+n1,m1+1:M+m1);
L=Lr(n1+1:N+n1,m1+1:M+m1);
%=====Визуализация результатов=====
subplot(221);imshow(L');axis('image');
subplot(222);imshow(Lvyh');axis('image');

Литература.

  1. Dash L., Chatterji B.N. Adaptive contrast enhancement and de-enhancement // Pattern Recognition, 1992. – V. 24, № 4. – P.289 – 302.
  2. Воробель Р.А. Цифровая обработка изображений на основе теории контрастности: Дис… докт. техн. наук: 05.13.06. – Львов, 1999. – 369 с.

Наверх

Анализ некоторых характеристик локальных окрестностей

Аналитическим выражением, описывающим количественное определение реакции зрительной системы на световое возбуждение, является его контраст. Вид этого выражения определяется свойствами конкретной зрительной системы восприятия [1, 2]. То есть изменение выражения определения контраста отвечает изменению типа зрительной системы или ее параметров. Это создает возможности адаптации зрительной системы путем изменения выражения для определения локального контраста. Разумеется, что при этом аналитическое выражение должно обеспечивать сохранение основных предельных свойств. Эти свойства заключаются в том, что локальный контраст приобретает максимальное значение только тогда, когда его компоненты имеют значения, которые лежат на противоположных краях диапазона, и равный нулю – в случае равенства этих компонентов по величине. Критерием оценки выражений контраста является эффективность их применения при цифровой обработке изображений. Следовательно, удачный выбор того или иного выражения определения контраста существенно влияет на дальнейшее применение метода.

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

, , .

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

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

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

(1)

где – вычисляется за выражением

,

– значение гистограммы локальной окрестности (количество элементов с яркостью в окрестности ) для величины яркости элемента с координатами .

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

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

, (2)

где , – минимальное и максимальное значения яркости в скользящей окрестности с центром в элементе с координатами ;

– максимальное значение гистограммы скользящей локальной окрестности с центром в элементе с координатами .

Эта характеристика локальной окрестности приобретает минимальные значения на однородных участках, и максимальные – на бинарных участках.

Следующей характеристикой гладкости локальных окрестностей является среднеквадратическое отклонение значений яркостей элементов скользящей окрестности , которое вычисляется за выражением:

, (3)

где – среднеарифметическое значение яркостей элементов локальной окрестности с центром в элементе с координатами .

Выражение (3) равно нулю для однородных окрестностей и возрастает с увеличением неоднородности. Более наглядно характер изменения значений локальных характеристик в зависимости от типа окрестности демонстрирует рис. 1.

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

Рис. 1. Схематическое отображение изменения значений локальных характеристик в зависимости от типа скользящей окрестности (стрелками указано направление возрастания значения локальной статистики для различных типов скользящих окрестностей).

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

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

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

Следует отметить, что должно описываться относительно простыми выражениями, поскольку это влияет на вычислительную сложность метода.

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

Таблица 1.


п/п

Степень преобразования локальных контрастов , где – некоторая локальная статистика.

Типичные виды зависимостей

1

где .

2

,

где .

3

где .

4

,

где – среднее значение яркостей элементов изображения;

; – константа.

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

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

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

Литература

  1. Dash L., Chatterji B.N. Adaptive contrast enhancement and de-enhancement // Pattern Recognition, 1992. – V. 24, № 4. – P.289 – 302.
  2. Воробель Р.А. Цифровая обработка изображений на основе теории контрастности: Дис… докт. техн. наук: 05.13.06. – Львов, 1999. – 369 с.

Наверх

Статистическое определение локального контраста

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

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

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

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

, (1)

где – среднее значение яркостей элементов локальной окрестности .

Из выражения (1) следует, что , а . Второй момент, который называется дисперсией и обозначается как , служит для описания текстуры. Он является также мерой контраста интенсивности и применяется для описания однородности поверхностей. В некоторых работах в качестве меры контраста текстуры предложено использовать выражение

, (2)

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

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

На втором этапе осуществляем нелинейное преобразование локального контраста .

На третьем этапе восстанавливаем изображение путем определения нового значения яркости элемента с координатами . Для этого используем выражение, которое определяется из формулы (2):

. (3)

Описанную процедуру повторяем для каждого элемента изображения.

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

%Программа, реализующая метод контрастирования изображений
%с использованием статистического определения локальных контрастов 
clear;
%Считывание исходного изображения
L=imread('lena.tif');
L=double(L);
L=L(:,:,1)./255; %Поскольку изображение полутоновое, то L(:,:,1)=L(:,:,2)=L(:,:,3) 
[N M]=size(L);
m=5;n=m;n1=fix(n/2);m1=fix(m/2); %Размеры локальной области
Lmax=max(L(:));Lmin=min(L(:));
H=imhist(L,256);
figure, imshow(L);
%Преобразование матрицы значений яркостей исходного изображения 
%с целью избежания краевого эффекта
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
end;end;
   L2=L(1,1:M); L02=L2;
    for i=1:n1-1;
      L2=[L2;L02];
    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;
          L7=[L7;L07];
        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;
            L4=[L4;L04];
          end;
L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;
      L5=[L5;L05];
    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';  L2=[L2;L];
  L2=[L2;L7];  L2=L2';  L3=[L3;L5];  L3=[L3;L8];  L3=L3';
  L1=[L1;L2];  L1=[L1;L3];  L1=L1';
clear L2;clear L3;clear L4;clear L5;clear L6;
clear L7;clear L8;clear L02; clear L04;clear L05;clear L07;clear L;
F=ones(n,m);
Lser=filter2(F,L1,'same')/(n*m);
alfa=1;
P=3; %Параметр нормирования 
for i=n1+1:n1+N;
disp(i);
for j=m1+1:m1+M;
           if j==1+m1;
             D=0;
                 for a=-n1:n1;
                 for b=-m1:m1;
                     D(n1+1+a,m1+1+b)=(Lser(i+a,j+b)-L1(i+a,j+b))^2*H(round(255*L1(i,j))+1);
                 end;
                 end;
           end;
        if j>1+m1;
            for a=-n1:n1;
             D(n1+1+a,m+1)=(Lser(i+a,j+m1)-L1(i+a,j+m1))^2*H(round(255*L1(i,j))+1);
            end;
              D=D(1:n,2:m+1);
        end;
    Dyspers(i,j)=(1/(n*m))*sum(sum(D));
     C(i,j)=1-1/(Dyspers(i,j)/P+1);
   C(i,j)=C(i,j)^.67;
suma(i,j)=sum(sum(D))-D(n1+1,m1+1);
DRD(i,j)=sqrt(n*m*C(i,j)/(1-C(i,j)+eps)-suma(i,j));
if L1(i,j)>=Lser(i,j);
Lvyh(i,j)=Lser(i,j)+alfa*DRD(i,j);
else
Lvyh(i,j)=Lser(i,j)-alfa*DRD(i,j);
end;
if Lvyh(i,j)>1;
   Lvyh(i,j)=1;
end;
if Lvyh(i,j)<0;
   Lvyh(i,j)=0;
end;
end;
end;
Lvyh=Lvyh(1+n1:n1+N,1+m1:m1+M);
figure, imshow(abs(Lvyh)); %Визуализация результата 

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

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

в) изображение а), обработанное методом контрастирования с использованием предложенного выражения для статистического определения локального контраста

Рис. 1. Обработка изображений методом степенного преобразования локальных контрастов с использованием известного ( б) ) и статистического ( в) ) выражений их определения ( – количественная оценка качества изображений).

Наверх

Локально-адаптивное улучшение качества изображений

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

Типичное локальное преобразование, основанное на этих параметрах, переводит интенсивность исходного изображения Lin в интенсивность нового изображения Lout путем осуществления следующей операции над расположением (i,j) каждого пикселя:

где

– среднее значение интенсивностей элементов всего изображения Lin;

– среднеквадратическое отклонение интенсивностей элементов локальной окрестности изображения в точке с координатами (i,j);

– среднее значение интенсивности для окрестности с центром в точке (i,j);

k – некоторая константа, 0 < k < 1.

Следует отметить, что значения параметров и    зависят от заданной окрестности точки   , что делает

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

интенсивности пикселя исходного изображения Lin(i,j) и локальным средним    на

.

Среднеквадратическое отклонение будет принимать меньшие значения в малоконтрастных окрестностях и более

высокие значения в окрестностях с более высокой контрастностью. Учитывая это, а также то, что    находится

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

во избежание больших отклонений интенсивностей на отдельных участках.

%	Метод повышения визуального качества изображений.
clear;
L=imread('krepost.bmp');          % считывание исходных данных изображения
L=L(:,:,:);
L=double(L)./(256);
figure, imshow(L);  title('Input image'); 
[N M s]=size(L); %визуализация исходных данных
m=15; % размер локального окна
n=m;n1=fix(n/2); m1=fix(m/2);
Ltemp=L;
for r=1:s;
L=Ltemp(:,:,r);
% Предварительна обработка матрицы исходного изображения
% для устранения краевого эффекта
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;  for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
  end;end;
   L2=L(1,1:M);   L02=L2;
    for i=1:n1-1;      L2=[L2;L02];    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;          L7=[L7;L07];        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;            L4=[L4;L04];          end;
       L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;      L5=[L5;L05];    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';  L2=[L2;L];  L2=[L2;L7];  L2=L2';
  L3=[L3;L5];  L3=[L3;L8];  L3=L3';  L1=[L1;L2];  L1=[L1;L3];  Lr=L1';
clear L2;clear L3;clear L4;clear L5;clear L6;clear L7;clear L8;
clear L02;clear L04;clear L05;clear L07;clear L1;clear L;
Lser=mean(mean(Lr));
k=.4;
for i=1+n1:N+n1; 
    disp(i)
    for j=1+m1:M+m1;
                if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=Lr(i+a,j+b);
                        end;
                        end;
                 end;
           if j>1+m1;
            for a=-n1:n1;
              D(n1+1+a,m+1)=Lr(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
           end;      
LS=mean(mean(D));      
sigma=std2(D)+eps;
         Lvyh(i,j)=(k*Lser/sigma)*(Lr(i,j)-LS)+LS; 
end;
end;
Lvyh(Lvyh>1)=1;
Lvyh(Lvyh<0)=0;
Lvyh=Lvyh(n1+1:N+n1,m1+1:M+m1);
Ltemp(:,:,s)=Lvyh;
end;
figure, imshow(Ltemp);

Результаты компьютерного моделирования рассмотренного метода представлены на рис. 1-4.


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


Рис. 2. Обработанное изображение при m=35 и k=0.7.


Рис. 3. Обработанное изображение при m=35 и k=0.3.


Рис. 4. Обработанное изображение при m=15 и k=0.7.

Проанализируем результаты моделирования. В методе есть два основных параметра, которые существенно влияют на результат обработки – размер локальной окрестности m и коэффициент k. Рассмотрим два изображения, которые представлены соответственно на рис. 2 и рис. 4, которые представляют результат обработки при одинаковом коэффициенте k, но разных размерах локальных окон m. Из этих изображений видно, что уменьшение размера локального окна приводит к увеличению детальности обработки. Для анализа влияния коэффициента k при одинаковых размерах локальной окрестности m рассмотрим два других изображения, которые представлены на рис. 2 и рис. 3. Уменьшение коэффициента k приводит к устранению резких перепадов на изображении и понижению его контрастности. Таким образом, используя различные значения параметров m и k, можно управлять уровнем контрастности и детальности обработки изображений.

Литература:

  1. Прэтт У. Цифровая обработка изображений – М.: Мир, 1982. – 790 с.

Наверх

Фильтрация изображений: алгоритмы сглаживания изображений

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

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

Понятие сглаживания всегда подразумевает некоторое представление об “идеально гладком” сигнале. Такой сигнал – цель сглаживания.

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

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

Сглаживание для аддитивной модели

Аддитивная модель шума предполагает, что наблюдаемый сигнал представляет собой сумму полезного сигнала и шума. Ранговые алгоритмы сглаживания аддитивного шума проще всего обосновывать с позиций кусочно-постоянной модели изображения. При таком подходе, сглаживание можно определить как оценку параметра кластера, к которому принадлежит данный элемент. Для того чтобы найти эту оценку, необходимо определить границы кластера. Можно предложить два способа определения границ кластера: адаптивное квантование мод [1, 2] и “выращивание” кластера.

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

image67.gif (15797 bytes)

Рис. 1. Адаптивное квантирование мод: а – исходная гистограмма распределения значений видеосигнала; б – гистограмма после адаптивного квантирования.

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

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

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

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

Сглаживание для модели импульсных помех

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

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

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

Если

,

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

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

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

,

являющийся, как известно, устойчивой к распределению оценкой разброса значений в выборке.

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

Таким образом, алгоритм сглаживания импульсного шума может быть представлен в виде:

image215.gif (1832 bytes)

где SMTH(M) означает сглаживание по некоторой окрестности М, из которой исключены точки, подлежащие исправлению.

Характерными при сглаживании импульсных помех являются ошибки ложного обнаружения, которые приводят к нежелательному сглаживанию деталей изображения, и ошибки пропуска, из–за которых на изображении могут остаться несглаженные выбросы помех. Доля этих ошибок зависит от порогов: с увеличением порогов доля ложных обнаружений падает, а доля пропусков возрастает. Следует учитывать, что число ложных обнаружений и пропусков возрастает также из–за возможного наличия в S–окрестности не одного, а нескольких больших выбросов помех. Поэтому для повышения качества сглаживания импульсных помех его целесообразно проводить итеративно, начиная с больших значений порогов, и, по мере удаления больших выбросов помех, понижая пороги на каждой итерации.

Увеличение детальности изображений

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

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

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

image216.gif (1112 bytes),

где сумма элементов локальной окрестности, взятых с некоторыми весами; g – коэффициент усиления.

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

image219.gif (1156 bytes),

где .

Преимущества нерезкого маскирования с ранговым сглаживанием вместо линейного – адаптивность и меньшая пространственная инерционность – вытекают из преимуществ рангового сглаживания.

Обнаружение деталей и их границ

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

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

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

Само по себе обнаружение состоит в сравнении измеренной степени соответствия с порогом. При препарировании изображений имеет смысл также предъявлять для визуализации саму величину соответствия, а не только бинарный результат сравнения с порогом. При этом обнаружение осуществляется оператором визуально [2].

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

Применения ранговых алгоритмов

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

Автоматическая диагностика параметров помех и искажений видеосигнала.

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

Стандартизация изображений.

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

Определение статистических характеристик видеосигнала и измерение текстурных признаков.

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

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

Кодирование изображений.

Возможность применения ранговых алгоритмов для кодирования изображений связана с использованием алгоритмов адаптивного квантования мод в режиме пофрагментной обработки. В этом случае анализируется гистограмма распределения значений элементов изображения в пределах фрагмента (или, как принято говорить в кодировании, блока), находятся границы кластеров, которые выбираются в качестве границ интервалов квантования, и производится квантование всех отсчетов фрагмента в соответствии с найденными границами. Как правило, если размеры фрагмента не слишком велики, количество уровней квантования Qs отсчетов фрагмента намного меньше количества Q уровней квантования, выбираемого из условия качественного воспроизведения всего изображения. Нетрудно подсчитать, что количество бит, требуемых для передачи значений NB отсчетов фрагмента, будет равно сумме Qs log2Q бит на передачу таблицы квантования и Nslog2Qs бит на передачу номера уровня квантования, т.е. на один отсчет изображения требуется в среднем Iog2Qs+(Qslog2Q)/Ns бит вместо log2Q безадаптивного квантования по фрагментам. Отсюда вытекает, что площадь фрагментов целесообразно увеличивать до тех пор, пока количество уровней квантования Qs не превысит нескольких единиц. Опыты, проведенные по пофрагментному квантованию мод, показывают, что это возможно при размерах фрагмента до 30х30 элементов. Следовательно, оценкой потенциальных возможностей кодирования изображений этим методом является величина порядка 1–2 бит на элемент.

Литература

  • Rosenfeld A., Troy E.B. Visual Texture Analysis // Conference Record of the Symposium on Feature Extraction and Selection in Pattern Recognition. – IEEE Publ. – 1970. – 70C51 – C.
  • Ярославский Л.П. Цифровая обработка сигналов в оптике и голографии: Введение в цифровую оптику – М.: Радио и связь. – 1987. – 296 с.: ил.
  • Беликова Т.П., Ярославский Л.П. Использование адаптивных амплитудных преобразований для препарирования изображений // Вопросы радиоэлектроники. Сер. Общетехн. – 1974. – Вып. 14.

Наверх

Двумерное сглаживание изображений

Довольно часто развитие вычислительной техники и информационных технологий приводит к необходимости модификации известных алгоритмов обработки изображений. Рассмотрим один из таких методов – двумерное сглаживание изображений.

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

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

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

(1)

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

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

Как указал Д. Габор [3], спектральное описание статистических свойств изображения оказывается недостаточным, так как оно не отражает его локально-анизотропную (но изотропную в целом) структуру.

Для описания такой структуры в статье предлагается “составная” модель фрагмента изображения. Эта модель используется для синтеза алгоритма сглаживания, дающего оценку элементов изображения с минимальным СКО.

“Составная” модель фрагмента изображения. Одноцветное неподвижное изображение может быть описано как распределение яркости , где пара чисел – координаты точек плоскости изображения. Будем рассматривать далее только дискретизированные изображения с целочисленными координатами.

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

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

Одни классы образованы фрагментами с изотропной структурой (без преобладания связей в каком-либо направлении на плоскости изображения), другие – фрагментами с той или иной анизотропией. Пусть – -мерная плотность вероятности фрагмента при условии, что фрагмент принадлежит классу (, а – распределение вероятностей классов . Тогда

. (2)

Выражение (2) есть разложение плотности по системе плотностей 1,…,M .Такое представление особенно полезно, когда хорошо аппроксимируется с помощью небольшого набора гауссовских распределений:

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

Результаты статистических измерений фрагментов реальных изображений показывают, что хорошее описание их может быть получено с помощью модели, имеющей всего пять классов. Четыре из них соответствуют преобладающим корреляционным связям в одном из четырех направлений, составляющих углы 0°, 45°, 90° и 135° с горизонталью (выбор направлений обусловлен наличием квадратной решетки, на которой задано изображение). Пятый класс описывает фрагменты с “изотропной” структурой. Для нескольких изображений были измерены матрицы и распределение вероятностей ,

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

Известно, что оптимальной оценкой элемента является апостериорное условное среднее значение

, (4)

где

. (5)

В этих выражениях обозначает ; и есть условная плотность наблюдаемого фрагмента при заданном фрагменте . Используя (2), перепишем (5) в следующем виде:

, (6)

где

(7)

и

. (8)

Подставив (6) в (4), получим:

, (9)

где

. (10)

Величина есть условная оценка при заданном классе фрагмента. Оценка есть взвешенная сумма условных оценок. Вес каждой условной оценки есть апостериорная вероятность (8) класса при данном фрагменте .

Пусть шум имеет гауссовское распределение с известной ковариационной матрицей N. Тогда

. (11)

Подставив (11) и (3) в (7), получим затем из (10) известную формулу Винера [2]

, (12)

где – элемент матрицы .

В этом случае апостериорная вероятность состояния (8) переходит в

. (13)

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

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

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

Реализация алгоритма. На основе теоретического материала было проведено компьютерное моделирование предложенного в [1] алгоритма.

Для моделирования зашумленного изображения на оригинал накладывался белый гауссовский шум [4]. Наблюдаемый фрагмент содержал 5×5 элементов, причем оцениваемый элемент находился в центре фрагмента. Предполагалось, что имеется пять классов фрагментов, и использовались пять матриц , одна из которых соответствовала “изотропным” корреляционным связям и четыре были анизотропными. Эти матрицы соответствовали связям только вдоль одной из четырех прямых, проходящих через центр (вертикальной, горизонтальной и составляющей угол + 45° с горизонталью). Таким образом, имелось пять режимов сглаживания: по всему фрагменту из 25 элементов и по отрезкам прямых, содержащих по пять элементов. Пять оценок, получаемых в результате пяти процедур сглаживания, суммировались с весами, равными вычисленным апостериорным вероятностям. На рис. 1 показан оригинал с наложенным на него шумом, среднеквадратическая величина которого составляла 5% от диапазона яркостей изображения. На рис. 2 приведено изображение, полученное в результате сглаживания с зашумленного изображения в соответствии с алгоритмом (9). На рис. 3 показан результат винеровского сглаживания. Сравнение рис. 2 и 3 показывает, что описанный выше алгоритм (9) приводит к меньшей нерезкости изображения, чем алгоритм Винера, и, следовательно, использованная модель более адекватна структуре изображения, чем гауссовская.

Литература

  1. Лебедев Д.С., Миркин Л.И. Двумерное сглаживание с использованием “составной” модели фрагмента. Сб. “Иконика. Цифровая голография. Обработка изображений”. М., “Наука”, 1975.
  2. W. Wiener. The Interpolation, Extrapolation and Smoothing of Stationary Time Series. N.Y., 1949.
  3. D.Gabor. The Smoothing and Filtering of Two-Dimensional Images. – “Progress in Biocybernetics”, v. 2. Amsterdam, 1965.
  4. D.Gabor. The Smoothing and Filtering of Two-Dimensional Images. – “Progress in Biocybernetics”, v. 2. Amsterdam, 1965.
  5. R.E.Graham. Snow Removal – A noise -Stripping Process for picture Signals. – IRE Transaction on Information Theiry, 1962. V.IT – IT-8, № 2.
  6. Г. Г. Вайнштейн. Пространственная фильтрация изображений средствами аналоговой вычислительной техники.- В сб. “Иконика. Пространственная фильтрация изображений. Фотографические системы”. М., “Наука”, 1970.

Наверх

Обобщенная линейная фильтрация

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

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

Задача линейной фильтрации как это констатируется, связана с применением линейной системы для извлечения сигнала из суммы сигнала и шума. С точки зрения векторного пространства задачей линейной фильтрации можно считать определение такого линейного преобразования в векторном пространстве, которое сводит длину или норму вектора ошибки к минимуму. Норма для данного векторного пространства определяет используемый критерий ошибки. Во многих случаях, когда сигнал суммируется с шумом, линейная система не является лучшей системой. Рассмотрим, например, квантованный сигнал с уровнями квантования 1, 2, 3,…, и допустим, что к нему добавились шумы с пиковыми значениями ±0,25. Ясно, что сигнал может быть точно восстановлен с помощью квантизатора, хотя его нельзя формально обосновать как оптимальный нелинейный фильтр. В менее очевидных случаях могут существовать одновременно формальные обоснования как для “лучшего” линейного фильтра, так и для «лучшего» нелинейного фильтра из некоторого класса, но при этом не всегда может быть проведено полное и точное сравнение этих фильтров, хотя бы из-за того, что они часто используют различную информацию о входных сигналах.

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

Если x1 и х2 обозначают два сигнала, которые объединяются с помощью операции О, то каноническая форма для класса гомоморфных фильтров, которые можно было бы использовать для восстановления x1 или х2, имеет вид, приведенный на рис. 1.

Рис. 1. Каноническая форма класса гомоморфных фильтров, используемых для разделения сигналов, объединенных с помощью операции О.

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

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

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

Применение гомоморфной фильтрации

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

Гомоморфная обработка изображений

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

,

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

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

image204.gif (963 bytes)       (1)

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

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

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

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

,

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

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

 

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

Целесообразно предположить, что  содержит главным образом низкие пространственные частоты. Быстрые изменения в  способствуют появлению высоких пространственных частот в , хотя отражательная способность вносит некоторый вклад и в низкие пространственные частоты. Таким образом, возможна только частично независимая обработка. Тем не менее на практике оказалось полезным связать низкие пространственные частоты с , а высокие пространственные частоты — с . При таком предположении линейный фильтр (рис. 3) выбирался так, чтобы он производил умножение низких пространственных частот на  и высоких пространственных частот на .

s3.gif (2633 bytes)

Рис. 3. Частотная характеристика линейного фильтра (рис. 2) для одновременного сжатия динамического диапазона и усиления контрастности.

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

Рис. 4. Изображения а) и в) после обработки с целью одновременного изменения (растяжения или сжатия) динамического диапазона и усиления контрастности. Результаты обработки представлены на рис. 4 б) и г) соответственно.

Литература

  1. Blackman R.B., Tukey Y.W., “The Measurement of Power Spectra from the Point of View of Communication Engineering”, Dover, 1959.
  2. Голд Б., Рэйдер Ч. Цифровая обработка сигналов. Пер. с англ., под ред. А.М.Трахтмана. М., “Сов. радио”, 1973, 368 с.

Наверх

Градиентный метод выделение контуров объектов на цветных изображениях

Если изображение представить двумерной функцией I = ƒ(x,y), то модуль градиента вычисляется за выражением

Одним из свойств градиента является то, что он всегда направлен в сторону возрастания функции в точке с координатами (x,y).

При практической реализации этого подхода производные аппроксимируются разностями значений пикселей. Представим некоторую локальную окрестность в виде

z1 z2 z3
z4 z5 z6
z7 z8 z9

Тогда производная по представляется в виде

Тогда производная по представляется в виде

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

Пример. Представим реализацию рассмотренного выше подхода для выделения выделение контуров объектов. Считаем и визуализируем некоторое изображение.

I=imread(‘zhaba.bmp’);

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

I=I(:,:,1);

I=double(I)./255;

figure,imshow(I);

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

Сформируем маску фильтра для выделения контуров объектов изображения.

f=fspecial(‘sobel’)

f =

     1     2     1

     0     0     0

    -1    -2    -1

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

I_filtered_x=imfilter(I,f,’replicate’);

figure,imshow(I_filtered_x);

Визуализация градиента по направлению без коррекции динамического диапазона

Полученное изображение I_filtered_x выходит за пределы динамического диапазона, поэтому его можно скорректировать.

I_filtered_x=I_filtered_x+min(min(I_filtered_x));

MN=min(min(I_filtered_x));

MX=max(max(I_filtered_x));

I_filtered_x=(I_filtered_x-MN)/(MX-MN);

figure,imshow(I_filtered_x);

Визуализация градиента по направлению y c коррекцией динамического диапазона

Полученное изображение I_filtered_y выходит за пределы динамического диапазона, поэтому его можно скорректировать.

MN=min(min(I_filtered_y));

MX=max(max(I_filtered_y));

I_filtered_y=(I_filtered_y-MN)/(MX-MN);

figure,imshow(I_filtered_y);

Визуализация градиента по направлению c коррекцией динамического диапазона

Далее сформируем градиент изображения на основе вычисленных градиентов по направлению x и y.

I_filtered=sqrt(I_filtered_x.^2+I_filtered_y.^2);

figure,imshow(I_filtered);

Градиент изображения без коррекции динамического диапазона

Полученное изображение I_filtered также выходит за пределы динамического диапазона, поэтому его можно скорректировать.

I_filtered=I_filtered+min(min(I_filtered));

MN=min(min(I_filtered));

MX=max(max(I_filtered));

I_filtered=(I_filtered-MN)/(MX-MN);

figure,imshow(I_filtered);

Градиент изображения с коррекцией динамического диапазона

Рассмотренный нами метод выделения контуров объектов может применяться для обработки полутоновых изображений.

Рассмотрим вопрос вычисления градиента изображений в цветовом пространстве RGB. Считаем и визуализируем исходное изображение.

L=imread(‘peppers.bmp’);

L=double(L)./255;

figure,imshow(L);

Выполним фильтрацию по отдельным цветовым компонентам и сформируем из них результат в виде RGB изображения.

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

L_R=L(:,:,1); figure, imshow(L_R);

L_G=L(:,:,2);

figure, imshow(L_G);

L_B=L(:,:,3);

figure, imshow(L_B);

Градиент изображения без коррекции динамического диапазона

Полученное изображение I_filtered также выходит за пределы динамического диапазона, поэтому его можно скорректировать.

I_filtered=I_filtered+min(min(I_filtered));

MN=min(min(I_filtered));

MX=max(max(I_filtered));

I_filtered=(I_filtered-MN)/(MX-MN);

figure,imshow(I_filtered);

Градиент изображения с коррекцией динамического диапазона

Рассмотренный нами метод выделения контуров объектов может применяться для обработки полутоновых изображений.

Рассмотрим вопрос вычисления градиента изображений в цветовом пространстве RGB. Считаем и визуализируем исходное изображение

L=imread(‘peppers.bmp’);

L=double(L)./255;

figure,imshow(L);

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

L_R=L(:,:,1);

figure, imshow(L_R);

L_G=L(:,:,2);

figure, imshow(L_G);

L_B=L(:,:,3);

figure, imshow(L_B);

Фильтрация по цветовой компоненте B изображения

L_B_filtered_x=imfilter(L_B,f,’replicate’);%фильтрация по направлению х

L_B_filtered_y=imfilter(L_B,f’,’replicate’);%фильтрация по направлению у

L_B_filtered=sqrt(L_B_filtered_x.^2+L_B_filtered_y.^2);

figure,imshow(L_B_filtered);%визуализация отфильтрированной цветовой

%компоненты изображения

Результат объединения отдельных отфильтрированных цветовых компонент в одно RGB-изображение представлен на изображении внизу.

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

[VG, A, PPG]=colorgrad(L, T);

где L – исходное RGB изображение; Т – порог; VG – модуль RGB градиента; A – матрица углов; PPG – градиент, который получен сложением двух одномерных градиентов по отдельным цветовым составляющим.

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

Сегментация RGB изображений

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

Рассмотрим решение этой задачи на конкретном примере.

Для этого считаем и визуализируем некоторое исходное изображение.

L=imread(‘peppers.bmp’);

figure,imshow(L);

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

mask=roipoly(L);

red=immultiply(mask,L(:,:,1));

green=immultiply(mask,L(:,:,2));

blue=immultiply(mask,L(:,:,3));

g=cat(3,red, green,blue);

figure,imshow(g);

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

[M, N, K]=size(g);

I=reshape(g, M*N, 3);

idx=find(mask);

I=double(I(idx,1:3));

[C,m]=covmatrix(I);

d=diag(C);

sd=mean(sqrt(d));

После этого можно приступить к выполнению сегментации с помощью встроенной функции colorseg при разных порогах – sd,3*sd и при 5*sd.

E25=colorseg(‘euclidean’,L,sd,m);

figure,imshow(E25);

E25=colorseg(‘euclidean’,L,3*sd,m);

figure,imshow(E25);

E25=colorseg(‘euclidean’,L,5*sd,m); figure,imshow(E25);

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

E25=colorseg(‘mahalanobis’,L,sd,m);

figure,imshow(E25);

E25=colorseg(‘mahalanobis’,L,3*sd,m);

figure,imshow(E25);

E25=colorseg(‘mahalanobis’,L,5*sd,m);

figure,imshow(E25);

Выбор различных мер сходства цветов (евклидово расстояние или расстояние Махаланобиса) определяет методики отслеживания цветовых данных. А выбор порога (sd,3*sd и 5*sd) влияет на уровень выделения областей.

Литература.

Гонсалес Р., Вудс Р., Эддинс С. Цифровая обработка изображений в среде Matlab. Москва: Техносфера, 206. – 616 с.

Наверх

Пространственная фильтрация

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

В этом материале рассмотрим некоторые методы пространственной фильтрации. Большинство из них относится к локальным методам. Реализация этих методов состоит из четырех основных этапов:

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

Далее приведем примеры применения некоторых видов пространственных фильтров.

Сначала считаем исходное изображение

L=imread('leo1.bmp');

Для лучшего понимания сути приведем реализацию метода для одной цветовой составляющей.

L=L(:,:,1);
figure, imshow(L);


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

Испортим исходное изображение шумом типа ”соль и перец”.

L=imnoise(L,'salt & pepper',0.01);
figure, imshow(L);


Рис. 2.

Далее применим к этому зашумленному изображению различные виды пространственных фильтров и проанализируем результат обработки.

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

– ‘P’ – границы изображения расширяются значением Р. По умолчанию значение Р=0;
– ‘replicate’ – размер изображения увеличивается повторением величин на его боковых границах;
– ‘symmetric’ – размер изображения увеличивается путем зеркального отражения через границы;
– ‘circular’ – размер изображения увеличивается периодическим повторением двумерной функции;

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

%n, m – размеры локальной окрестности
n1=fix(n/2);m1=fix(m/2);
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;
  for j=1:m1;
    L1(i,j)=a;
    L3(i,j)=b;
    L6(i,j)=c;
    L8(i,j)=d;
  end;
end;
   L2=L(1,1:M);
   L02=L2;
    for i=1:n1-1;
      L2=[L2;L02];
    end;
    L7=L(N,1:M);
    L07=L7;
        for i=1:n1-1;
          L7=[L7;L07];
        end;
     L4=L(1:N,1);
     L4=L4';
     L04=L4;
          for i=1:m1-1;
            L4=[L4;L04];
          end;
  L4=L4';
  L5=L(1:N,M);
  L5=L5';
  L05=L5;
    for i=1:m1-1;
      L5=[L5;L05];
    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';  L2=[L2;L];  L2=[L2;L7];
  L2=L2';  L3=[L3;L5];  L3=[L3;L8];  L3=L3';  L1=[L1;L2];  L1=[L1;L3];
  Lr=L1';

Приведем программную реализацию непосредственно самих методов

for i=1+n1:N+n1;
    disp(i);
    for j=1+m1:M+m1;

% Формирование локальной окрестности D
                 if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=Lr(i+a,j+b);
                        end;
                        end;
                 end;
           if j>1+m1;
            for a=-n1:n1;
              D(n1+1+a,m+1)=Lr(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
          end;    

% Формирование результата обработки

% Контргармоническое среднее
               Lout1(i,j)=sum(sum(D.^(q+1)))/sum(sum(D.^q+eps));
% Арифметическое среднее
               Lout2(i,j)=sum(sum(D))./(m*n); 
% Геометрическое среднее               
               Lout3(i,j)=prod(prod(D)).^(1/(m*n));    
 % Гармоническое среднее           
               Lout4(i,j)=(m*n)/sum(sum(1./(D+eps))); 
  % Медиана           
               Lout5(i,j)=median(median(D)); 
 % Максимум по окрестности 
               Lout6(i,j)=max(D(:)); 
  % Минимум по окрестности 
               Lout7(i,j)=min(D(:));  
  % Срединная точка
               Lout8(i,j)=0.5*(max(D(:))+min(D(:)));         
          end;
  end;
  % а–усеченное среднее
d=6;
               Dus=sort(D(:));
               Dus=Dus(d/2+1:end-d/2);
               Lout9(i,j)=mean(Dus);

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

Lout1=Lout1(n1+1:N+n1,m1+1:M+m1);
Lout2=Lout2(n1+1:N+n1,m1+1:M+m1);
Lout3=Lout3(n1+1:N+n1,m1+1:M+m1);
Lout4=Lout4(n1+1:N+n1,m1+1:M+m1);
Lout5=Lout5(n1+1:N+n1,m1+1:M+m1);
Lout6=Lout6(n1+1:N+n1,m1+1:M+m1);
Lout7=Lout7(n1+1:N+n1,m1+1:M+m1);
Lout8=Lout8(n1+1:N+n1,m1+1:M+m1);
Lout9=Lout9(n1+1:N+n1,m1+1:M+m1);

Визуализируем полученные результаты.

figure, imshow(Lout1);
title('Контргармоническое среднее');

Контргармоническое среднее

Рис. 3. Контргармоническое среднее при q=1,5.

Контргармоническое среднее

Рис. 4. Контргармоническое среднее при q=-1,5.

Уравнение контргармонического фильтра можно представить следующим выражением:

          (1)

В зависимости от значения параметра q в выражении (1) метод приводит к удалению белых или темных импульсных выбросов (см. рис. 3).

figure, imshow(Lout2);
title('Арифметическое среднее');

Арифметическое среднее

Рис. 5. Арифметическое среднее при размере локального окна 33.

Арифметическое среднее

Рис. 6. Арифметическое среднее при размере локального окна 99.

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

figure, imshow(Lout3);
title('Геометрическое среднее');

Геометрическое среднее

Рис. 7. Геометрическое среднее при размере локального окна 33.

Геометрическое среднее

Рис. 8. Геометрическое среднее при размере локального окна 99.

Фильтр типа геометрическое среднее можно представить следующим выражением:

          (2)

На основе выражения (2) можно предположить, что результат обработки будет существенно зависеть от размеров локальной апертуры. Это подтверждают обработанные этим методом изображения, которые представлены на рис. 7 и 8.

figure, imshow(Lout4);
title('Гармоническое среднее');

Гармоническое среднее

Рис. 9. Гармоническое среднее при размере локального окна 33.

Гармоническое среднее

Рис. 10. Гармоническое среднее при размере локального окна 99.

Исходя из выражения фильтра гармонического среднего

          (3)

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

figure, imshow(Lout5);
title('Медиана');

Медиана

Рис. 11. Медиана при размере локального окна 33.

Медиана

Рис. 12. Медиана при размере локального окна 99.

Для устранения данного вида шума (‘соль и перец’) медианная фильтрация является одним из наиболее простых и эффективных методов. Размер локальной апертуры влияет на результат обработки, но не в такой степени как в других рассмотренных нами методах. Одно из основных преимуществ медианной фильтрации состоит в том, что она приводит к устранению импульсных выбросов, не размывая границ объектов (см. рис. 11, 12).

figure, imshow(Lout6);
title('Максимум по окрестности');

Максимум по окрестности

Рис. 13. Максимум по окрестности при размере локального окна 33.

Максимум по окрестности

Рис. 14. Максимум по окрестности при размере локального окна 99.

figure, imshow(Lout7);
title('Минимум по окрестности');

Минимум по окрестности

Рис. 15. Минимум по окрестности при размере локального окна 33.

Минимум по окрестности

Рис. 16. Минимум по окрестности при размере локального окна 99.

На рисунках 13–16 приведены результаты обработки изображений с помощью фильтра типа максимум (рис. 13, 14) и минимум (рис. 15, 16) по окрестности. Конечно, результаты такой обработки существенно зависят от размеров локальной апертуры.

figure, imshow(Lout8);
title('Срединная точка');

Срединная точка

Рис. 17. Срединная точка при размере локального окна 33.

Срединная точка

Рис. 18. Срединная точка при размере локального окна 99.

Фильтр типа срединная точка может быть представлен выражением

          (4)

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

figure, imshow(Lout9);
title('a-усеченное среднее');

a-усеченное среднее

Рис. 19. a –усеченное среднее при размере локального окна 33.

a-усеченное среднее

Рис. 20. a –усеченное среднее при размере локального окна 99.

Фильтр типа a –усеченное среднее можно представить выражением

(5)

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

Из приведенных на рис. 19 и 20 изображений видно, что фильтр типа α–усеченное среднее является эффективным инструментом для удаления импульсных шумов. Конечно, нужно подобрать соответствующий размер локальной апертуры.

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

Литература.

  1. Прэтт У. Цифровая обработка изображений – М.: Мир, 1982. – 790 с.
  2. Гонсалес Р., Вудс Р., Эддтис С. Цифровая обработка изображений в среде MATLAB. Москва: Техносфера, 2006. – 616 с.

Программная реализация приведенных в материале методов (m-файл).

clear;
L=imread('leo1.bmp');
L=L(:,:,1);
figure, imshow(L);
[N M]=size(L);
L=imnoise(L,'salt & pepper',0.01);
figure, imshow(L);
L=double(L)./255;
q=1.5;n=3; m=n;
 
n1=fix(n/2);m1=fix(m/2);
a=L(1,1);b=L(1,M);c=L(N,1);d=L(N,M);
for i=1:n1;  for j=1:m1;
    L1(i,j)=a;    L3(i,j)=b;    L6(i,j)=c;    L8(i,j)=d;
end;end;
   L2=L(1,1:M);   L02=L2;
    for i=1:n1-1;      L2=[L2;L02];    end;
    L7=L(N,1:M);    L07=L7;
        for i=1:n1-1;          L7=[L7;L07];        end;
     L4=L(1:N,1);     L4=L4';     L04=L4;
          for i=1:m1-1;            L4=[L4;L04];          end;
       L4=L4';  L5=L(1:N,M);  L5=L5';   L05=L5;
    for i=1:m1-1;      L5=[L5;L05];    end;
     L5=L5';  L1=[L1;L4];  L1=[L1;L6];  L1=L1';  L2=[L2;L];  L2=[L2;L7];
  L2=L2';  L3=[L3;L5];  L3=[L3;L8];  L3=L3';  L1=[L1;L2];  L1=[L1;L3];
  Lr=L1';
  clear L1;  clear L2;  clear L3;  clear L4;  
    clear L5;  clear L6;  clear L7;  clear L8;
  
for i=1+n1:N+n1;
    disp(i);
    for j=1+m1:M+m1;
                 if j==1+m1;
                        D=0;
                        for a=-n1:n1;
                        for b=-m1:m1;
                           D(n1+1+a,m1+1+b)=Lr(i+a,j+b);
                        end;
                        end;
                 end;
           if j>1+m1;
            for a=-n1:n1;
              D(n1+1+a,m+1)=Lr(i+a,j+m1);
            end;
             D=D(1:n,2:m+1);
          end;            

Lout1(i,j)=sum(sum(D.^(q+1)))/sum(sum(D.^q+eps));	%Контргармоническое среднее
Lout2(i,j)=sum(sum(D))./(m*n);				%Арифметическое среднее
Lout3(i,j)=prod(prod(D)).^(1/(m*n));			%Геометрическое среднее
Lout4(i,j)=(m*n)/sum(sum(1./(D+eps)));			%Гармоническое среднее
Lout5(i,j)=median(median(D));				%медиана
Lout6(i,j)=max(D(:));					%max
Lout7(i,j)=min(D(:));					%min
Lout8(i,j)=0.5*(max(D(:))+min(D(:)));			%срединная точка 
   d=6;
   Dus=sort(D(:));
   Dus=Dus(d/2+1:end-d/2);
Lout9(i,j)=mean(Dus);					%а–усеченное среднее
          end;
  end;
          
Lout1=Lout1(n1+1:N+n1,m1+1:M+m1);
Lout2=Lout2(n1+1:N+n1,m1+1:M+m1);
Lout3=Lout3(n1+1:N+n1,m1+1:M+m1);
Lout4=Lout4(n1+1:N+n1,m1+1:M+m1);
Lout5=Lout5(n1+1:N+n1,m1+1:M+m1);
Lout6=Lout6(n1+1:N+n1,m1+1:M+m1);
Lout7=Lout7(n1+1:N+n1,m1+1:M+m1);
Lout8=Lout8(n1+1:N+n1,m1+1:M+m1);
Lout9=Lout9(n1+1:N+n1,m1+1:M+m1);

figure, imshow(Lout1);title('Контргармоническое среднее');
figure, imshow(Lout2);title('Арифметическое среднее');
figure, imshow(Lout3);title('Геометрическое среднее');
figure, imshow(Lout4);title('Гармоническое среднее');
figure, imshow(Lout5);title('Медиана');
figure, imshow(Lout6);title('Максимум по окрестности');
figure, imshow(Lout7);title('Минимум по окрестности');
figure, imshow(Lout8);title('Срединная точка');
figure, imshow(Lout9);title('a-усеченное среднее');

Наверх

Деконволюция

В пакете Image Processing Toolbox системы MATLAB существует довольно много функций для решения тех или иных задач обработки изображений, которые оперируют такими понятиями как функция распространения точки, обращение двумерной свертки (например, deconvblind, deconvlucy, deconvreg, deconvwnr [2] и др.). Рассмотрим это более детально с точки зрения теории.

Деконволюция

Решение задачи деконволюции заключается в обращении двумерной свертки. Термин “деконволюция” охватывает наиболее важные и широко используемые методы обработки изображений. Необходимость в такой операции возникает во всех областях науки, связанных с измерениями. По методах деконволюции существует большое число работ [1].

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

Какой бы метод не использовался, почти всегда необходимо провести предварительную обработку заданного искаженного изображения для преобразования его в форму, удобную для выполнения процедуры деконволюции. Предварительную обработку целесообразно разделить на пять категорий: сглаживание, разбиение на фрагменты, аподизацию (взвешивание обрабатываемого отрезка сигнала весовой функцией), расширение границ и сверхразрешение. Под сглаживанием изображения здесь понимается уменьшение зашумленности. Разбиение на фрагменты включает разделение изображения с пространственно-зависимой ФРТ на фрагменты, в каждом из которых ФРТ может приближенно рассматриваться как пространственно-инвариантная. Аподизация – это метод, позволяющий уменьшить влияние кадрового окна (записывающего устройства), которое производит усечение изображения. Однако этот метод может быть менее эффективным, нежели метод расширения границ, который мы удачно применяли в ряде случаев. Известны две модификации метода расширения границ – простое расширение и расширение с перекрыванием. Второй метод, как правило, более предпочтителен, поскольку в нем используются преимущества условия согласованности периодических сверток. Это еще один пример того, как повышается эффективность численного метода, когда более полно учитываются особенности исследуемой задачи в плане математической физики. Сверхразрешение рассматривается как процедура предварительной обработки, поскольку в конечном счете она позволяет уменьшить зашумленность.

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

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

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

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

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

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

Интеграл свертки представляется выражением

(1)

где h(x) – функция, задающая искажение; f(x) – функция, которую необходимо восстановить.

Согласно теореме о свертке фурье-образ величины (1) равен

(2)

где F(u) – функция, связанная с функцией f(x) двумерным преобразованием Фурье; H(u) – фурье-образ оптической передаточной функции.

Идеализированная задача конечной деконволюции такова: заданы функции b(x) и h(x) , требуется восстановить функцию f(x) при условии, что все три величины имеют конечную протяженность.

Из соотношения (2) следует, что эту задачу можно решить следующим образом

(3)

Операция деления внутри фигурных скобок в выражении (3) называется простой инверсной фильтрацией. Термин “фильтрация” здесь употребляется по аналогии с классической теорией цепей и современной теорией обработки сигналов. Классический фильтр представляет собой устройство, которое изменяет спектр временных частот сигнала. Спектр B(u) есть функция пространственной частоты.

Оптическая передаточная функция H(u) изменяет спектр пространственных частот B(u) в результате применения указанной выше операции деления.

Поскольку обработанные изображения обычно хранятся в памяти ЭВМ в виде квантованных значений, в технике обработки изображений, как правило, используются цифровые, а не классические аналоговые фильтры. Цифровой фильтр определяется дискретным массивом, вообще говоря, комплексных чисел, который изменяет в процессе некоторой операции обработки спектр пространственных частот. Следовательно, обе функции, h(x) в формуле (1) и H(u) в формуле (2), могут рассматриваться как фильтры (и в большей части приложений они реализуются в цифровом виде). Общепринятая классификация цифровых фильтров возникла в теории обработки сигналов как функций времени, и этой классификацией можно пользоваться в теории обработки одномерных изображений, т. е. сигналов как функций (одной) пространственной переменной. Мы перенесем соответствующую терминологию на двумерный случай. Понятие “отсчета” в теории обработки сигналов переходит в понятие “элемента изображения” в теории обработки изображений. Как отсчеты, так и элементы изображений должны квантоваться по амплитуде до их цифровой обработки. Изображение, к которому должна быть применена операция фильтрации, называется заданным изображением, и о нем говорят как о состоящем из заданных элементов изображения. Элементы профильтрованного изображения называются выходными элементами изображения. В случае нерекурсивного цифрового фильтра каждый выходной элемент изображения представляет собой взвешенную сумму заданных элементов изображения. В случае же рекурсивного цифрового фильтра каждый выходной элемент изображения есть взвешенная сумма заданных элементов изображения и рассчитанных ранее выходных элементов изображения. Все практически реализуемые цифровые фильтры, конечно, описываются массивами конечных размеров (в одномерном случае конечный фильтр часто называют коротким). Цифровой фильтр называется прямым, если он применяется в плоскости изображения, и спектральным – если он применяется в частотной плоскости. Каузальный фильтр является односторонним в том смысле, что его отклик всегда отстает от входного воздействия (это несколько искусственно в двумерном случае, но, конечно, имеет очень важное значение для операций одномерной фильтрации, которые лежат в основе обработки сигналов как функций времени). Каузальные фильтры почти всегда реализуются как прямые. Мультипликативный цифровой фильтр представляет собой спектральный фильтр, в котором каждый выходной отсчет получается как произведение заданного элемента входного сигнала на один элемент массива фильтра.

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

Прежде чем ставить практическую задачу деконволюции, исследуем некоторые свойства согласованности сверток.

В одномерном случае соотношение (2) представляется в виде

(4)

где вещественная переменная u заменена комплексной переменной w . Если функции f(x) и h(x) имеют конечную протяженность, так что протяженность функции b(x) тоже конечна, то их спектры характеризуются множествами нулей в комплексной w -плоскости.

Если заданное множество Zg представить в виде множества вещественных нулей Zgr и нулей, которые могут быть комплексными Zgc , то можно записать

(5)

Это означает, что одномерная задача деконволюции является согласованной только в том случае, если все нули функции H(w) будут также и нулями функции B(w) . Следовательно, величины b(x) и h(x) нельзя задавать независимо; заранее должно быть известно, что они удовлетворяют соотношению (1). То же относится и к двумерным сверткам.

Теперь вернемся к периодически продолженному (с перекрыванием) идеальному искаженному изображению imb(x) и к его спектру IMb(x) . Последний можно записать в виде

(6)

где (·) – дельта-функция; Fi,m – коэффициенты Фурье истинного изображения f(x) , являющиеся также отсчетами функции F(u) , которые рассматриваются в теореме отсчетов. Эти отсчеты берутся в точках растра (l/L1 , m/L2 ) в частотной плоскости. Величины Hl,m входящие в выражение (6) – это отсчеты оптической передаточной функции H(u) в тех же самых точках растра:

(7)

где l и m – произвольные целые числа.

Теперь мы можем поставить идеализированную задачу периодической деконволюции: заданы функции imb(x) и h(x), требуется найти функцию f(x) [зная, что f(x) и h(x) – функции конечной протяженности, а imb(x) – периодическая функция].

По заданной функции b(x) можно рассчитать функцию B(u) и сразу же найти, что

(8)

Аналогичным образом вычисляются отсчеты оптической передаточной функции Hl,m . Из выражения (6) видно, что каждое значение Fl,m дается операцией деления , которая всегда может быть выполнена, если значения Hl,m отличны от нуля. Такой простой подход адекватен в случае функций b(x) и h(x) , выбранных достаточно независимо, поскольку функция Imb(u) в соответствии с выражением (6) фактически существует только в вышеупомянутых точках растра. Но подобный подход неприемлем в идеализированной задаче в случае конечной свертки, так как тогда B(u) – непрерывная функция переменной u.

Поэтому удивительно, что единственным условием согласованности для периодических сверток оказывается требование, чтобы величины Hl,m могли быть нулевыми только при тех значениях l и m , при которыхBp,l,m=0 . Это условие называется условием согласованности периодических сверток. Подчеркнем, что ни одна величина .Hl,m не может быть точно равна нулю при реальном измерении функции h(x) , или, что эквивалентно, функции H(u) , так что периодические свертки всегда на практике являются согласованными (они, конечно, очень сильно зашумлены, когда большое число величин Hl,m “малы” при значениях l и m , отвечающих существенно отличным от нуля значениям величин Bp,l,m ).

Практическая задача деконволюции ставится следующим образом: заданы функции b(x) и h(x) , требуется найти функцию f(x) , зная, что – усеченный вариант функции записываемого изображения r(x).

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

Любой вид предварительной обработки может, конечно, вносить свой шум в добавление к искажению изображения f(x) , уже имеющемуся в записываемом изображении r(x) . Но если разрывы не устранены, то соответствующие артефакты, как правило, преобладают над любым дополнительным шумом, вносимым предварительной обработкой. “Выровненную” форму изображения обозначим здесь через a(x) и будем называть предварительно обработанным записанным изображением. Хотя в результате проведения предварительной обработки должны изменяться все три величины, редко имеется какой-либо способ оценить, насколько именно, а потому обычно не имеет смысла говорить о различии между изображениями a(x) и r(x). Далее мы будем рассматривать эти два изображения как идентичные, по крайней мере на том кадре (т. е. в той области плоскости изображения), где умещается предварительно обработанный вариант изображения смысла . Поэтому будем считать, что

(9)

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

Теперь мы введем понятие “восстановимого истинного изображения” . Это оценка изображения f(x) , которую можно получить, исходя из изображения .

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

(10)

Коэффициенты Фурье функции удобно обозначить через , а для обозначения спектров функций a(x) , c(x) , и использовать соответствующие заглавные буквы со “шляпкой” или без нее.

Если есть опасение, что различия между и f(x) сильно увеличатся из-за отсутствия согласованности между функциями a(x) и h(x) , взятыми явно конечными, то можно обратиться к формуле

для периодического изображения imb(x) с заменой b на a . Тогда спектр IMb(u) свертки дается выражением (6), но с заменой величин и величинами Fl,m и Hl,m соответственно. Напомним, что на периодические свертки не оказывает влияния несогласованность, которая, как уже говорилось, может искажать свертки величин, имеющих конечные протяженности.

Литература:

1.Бейтс Р., Мак-Доннелл М. Восстановление и реконструкция изображений: Пер. с англ. – М.: Мир, 1989. – 336 с.
2. Дьяконов В. MATLAB. Обработка сигналов и изображений. Специальный справочник. – СПб.: Питер, 2002. – 608 с.

Наверх

Предварительная обработка изображений

Успешность восстановления изображений сильно зависит от качества предварительной обработки, в результате которой из записанного изображения получают изображение a(x) . Мы разделяем предварительную обработку на пять категорий: сглаживание, разбиение на фрагменты, аподизацию, расширение границ и сверхразрешение [1].

Обычно для более полного уменьшения эффектов зашумления проводят сглаживание изображения . Хотя эта процедура часто носит главным образом косметический характер, она может иметь и более важное практическое значение. Напомним, что величина c(x) (см. “Деконволюция”, формула (9)) учитывает эффекты, связанные с нелинейностями записи, шумом записи изображения, ошибками в передаче битов, отсутствием некоторой информации (т. е. отсутствием отдельных элементов изображений или целых групп их), насыщением, а также с загрязнением и царапинами, которые искажают фотографии. Сглаживание можно рассматривать как двумерный аналог простейшей обработки сигналов, имеющей целью исключить весь шум, спектральные составляющие которого лежат вне полосы временных частот, соответствующей сигналу, передаваемому рассматриваемым каналом связи. Большинство видов помех, перечисленных выше, можно считать помехами с независимыми отсчетами, тогда как характерные детали изображений обычно коррелированны в пределах нескольких соседних элементов изображения. Иначе говоря, спектр пространственных частот шума существенно шире, чем спектр изображения, и в этом случае весьма эффективна пространственная фильтрация изображения , оставляющая только те спектральные составляющие шума c(x) , которые разрешены в той же степени, что и деталь в истинном изображении.

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

Методы деконволюции прямо применимы только в случае пространственно-инвариантной функции распространения точки (ФРТ).

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

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

Обозначим область плоскости изображения, занимаемую заданным (но сглаженным, как описано выше) записанным изображением, через . Это согласуется с определениями, которые связывают записываемое и фактически записанное изображения r(x) и . Однако записываемое изображение r(x) в большей части практических приложений фактически усекается до вида , чем и объясняется, почему эти два изображения, вообще говоря, различны. Напомним, что – зашумленный вариант изображения b(x) . Последнее является фактически интересующим нас изображением, поскольку к нему мы хотели бы применить операцию деконволюции. Поскольку изображение b(x) выходит за пределы вышеуказанного кадра, логично предположить, что то же самое имеет место и для изображения . Заметим, что изображение может все же отличаться от изображения r(x) . Поэтому для этого кадра мы введем другой символ Г . Имеет смысл пытаться восстанавливать только те части изображения f(x) , которые оказывают влияние на вид изображения b(x) в пределах кадра Г . Это части изображения f(x) , которые в исходном состоянии находятся в кадре Г , а также части этого изображения, которые вносятся в результате действия ФРГ в пределы кадра Г извне его. Обозначим через кадр, содержащий сумму этих частей. Кадр может быть построен путем размещения центра кадра ФРТ на внешней границе кадра Г и перемещения ее по этой границе (см. “Деконволюция”, формула (1)). Тогда кадр и будет представлять объединение всех точек в кадре Г и всех точек, охватываемых кадром при его прохождении по области Г .

Поскольку части истинного изображения, лежащие вне кадра , полностью теряются, можно предположить, что изображение f(x) лежит в пределах кадра . Поэтому далее будем считать, что

(1)

откуда, естественно, следует равенство

(2)

Так как изображение f(x) существует на кадре , из-за “размывающего” действия ФРТ h(x) оно должно сказываться в пределах большего кадра. Этот больший кадр содержит все части изображения b(x) , которые теряются при усечении изображения r(x) , а потому мы обозначим его через . Поскольку же есть зашумленный вариант изображения b(x) , область идентична области .

Хотя мы знаем о существовании изображения в пределах кадра , оно задается только на кадре Г . Обычно целесообразно провести дальнейшую обработку (кроме сглаживания) для более полной компенсации эффектов усечения, а также несогласованности операции свертки. Итак, нужно, чтобы предварительно обработанное записанное изображение (см. “Деконволюция”, формула (9)) удовлетворяло условию

(3)

где через pre{} обозначены операции, описываемые ниже

Усечение изображения имеет столь важное значение с практической точки зрения, что нужно остановиться на его последствиях. Сначала конкретизируем форму кадра Г . Мы будем рассматривать только прямоугольные и круговые кадры, поскольку они чаще встречаются в приложениях. Таким образом, если L1 и L2 есть x и y -протяженности прямоугольного кадра г или если R – радиус кругового кадра Г , следует, что

(4)

или

(5)

где за начало координат взят центр кадра Г . Эти два варианта изображения удобно исследовать раздельно. Из определений

(6)

а также выражения (4) и теоремы о свертке следует, что

(7)

где sinc(u) – фурье-образ rect(x) [2].

Взяв теперь фурье-образ функции (5) и вспомнив первое из двух определений (6), увидим, что

(8)

где – радиальная координата в частотной плоскости и

(9)

причем J1 – функция Бесселя первого рода первого порядка.

Отметим, что sinc(t) – осциллирующая функция, имеющая центральный пик (часто называемый основным лепестком ) приблизительно единичной ширины и бесконечную последовательность меньших пиков (иногда называемых боковыми лепестками), каждый из которых имеет эффективную ширину, равную 1/2, и амплитуду, которая уменьшается сравнительно медленно (по закону ). Эти боковые лепестки могут привести к неприемлемым артефактам, если изображение подвергается операции фильтрации без соответствующей предварительной обработки. Хотя это относится в первую очередь к изображению , определенному выражением (4), то же самое справедливо и для изображения, заданного выражением (5). Функция jinc , введенная в формуле (8), аналогична функции sinc . Она фактически эквивалентна двум функциям sinc , входящим в формулу (7). Отметим, что типичная фильтрация может быть описана соотношением

(10)

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

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

Если мы знаем, что более интересные для нас части изображения f(x) лежат ближе к центру кадра Г , то в тех случаях, когда размер последнего существенно больше размера кадра , предварительная обработка может состоять в аподизации. Она заключается в умножении функции на функцию окна m=m(x) , которая плавно уменьшается до нуля на внешней границе кадра Г и равна нулю везде вне кадра Г . Вследствие этого область оказывается равной кадру Г . Обращаясь теперь к формуле (3), можно получить, что предварительно обработанное записанное изображение принимает вид

(11)

где – сглаженное изображение, полученное из фактически записанного изображения.

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

Поскольку здесь рассматриваются только изображения , описываемые выражениями (4) и (5), по-видимому, не имеет особого смысла изучать функции окна, которые не обладают свойством круговой симметрии или не разделяются на сомножители, зависящие от переменных x и y по отдельности. Поэтому достаточно исследовать одномерные функции окна, например m(x) (через x обозначены переменные x , y или r ). В качестве величины L , удобно взять размер (усеченного) фактически записанного изображения в x-направлении. Приняв обозначение M(u)=M=F(m) и предположив, что m – функция, аналитическая на интервале <L/2 (т. е. “непрерывно гладкая”, иначе говоря, бесконечно дифференцируемая), мы видим, что интеграл Фурье, определяющий величину M , можно взять по частям и получить следующее выражение:

где введено обозначение

(13)

Если взять интеграл (12), то появляется возможность анализировать боковые лепестки функции M(u) . Чтобы боковые лепестки быстрее уменьшались с возрастанием величины , необходимо обеспечить выполнение условия

(14)

где – положительное целое число.

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

Литература.

1. Бейтс Р., Мак-Доннелл М. Восстановление и реконструкция изображений: Пер. с англ.- М.: Мир, 1989. – 336.
2. Bracewell R.N. The Fourier Transform and its Applications. – N.Y.: McGraw-Hill, 1978.

Наверх

Расширение границ изображений. Сверхразрешение

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

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

а) свободно от скачкообразных изменений вблизи своей границы;

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

в) содержит всю записанную информацию.

Такой вид предварительной обработки будем называть расширением границ. Эта процедура описывается соотношением

при .

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

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

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

при ,

при .

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

Отметим, что изображением , полученным в результате расширения границ с перекрыванием в пределах каждой основной ячейки, аппроксимируется периодически продолженное идеальное искаженное изображение. Практически успешность процедуры расширения границ сильно зависит от выполнения следующего требования “гладкости”. Изображение в пределах границы произвольной ячейки должно быть по крайней мере столь же гладким, как и изображение в пределах любой внутренней ячейки. Оценку выполнения этого требования вполне допустимо осуществлять визуально. Данным требованием предотвращается появление в изображении ложных составляющих с высокими пространственными частотами, благодаря чему повышается общая устойчивость процесса восстановления изображений. Еще одно преимущество процедуры расширения границ с перекрыванием перед процедурой простого расширения состоит в том, что покрывается менее половины области и допускается меньше произвола. Эта процедура требует минимума экстраполяции для заполнения границы каждой ячейки и поэтому оказывается, вообще говоря, более точной. Если вспомнить основную задачу деконволюции, то методы простого расширения границ и расширения с перекрыванием можно рассматривать как средства сведения практической задачи деконволюции к задачам, соответственно, идеализированной конечной и периодической деконволюции, но с минимизацией вредного влияния искажений, которые неизбежно зашумляют записанные изображения. Очевидно, что процедуру расширения границ с перекрыванием можно реализовать столь же непосредственно, как и простое расширение границ, но практически она часто дает гораздо лучшие результат ы. Тем не менее процедура простого расширения границ тоже нередко применяется, особенно в тех случаях, когда метод с перекрыванием по какой-либо технической причине нельзя использовать.

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

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

Для этого необходимо рассмотреть одномерное изображение конечной протяженности и его фурье-образ .

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

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

. (1)

Неотрицательная вещественная константа обычно полагается равной нулю или единице. Но, по-видимому, нет особых оснований для такого ограничения, если исходить из термодинамических аналогий. Единственный критерий приемлемости метода обработки изображений – это качество получаемых результатов, а потому можно испробовать и другие значения константы . Рассматриваемый метод заключается в продолжении функции , обеспечивающем максимум энтропии . Вся суть этого метода, возможно, в том, что логарифмы в определении (1) исключают получение отрицательных значений для изображения с вещественными значениями. Данный метод частот приводит к впечатляющим результатам, но они очень сильно за висят от некоторых деталей истинного изображения, например спектр может резко измениться, если к изображению добавить небольшой фон. В настоящее время метод максимальной энтропии изучен недостаточно, чтобы о нем можно было судить вполне объективно. Однако это один из немногих методов, при которых в сверхразрешенном изображении автоматически учитывается требование неотрицательности.

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

Пусть и – верхняя и нижняя границы величины , и – соответствующие границы величины

; (2)
. (3)

Предположим, что существует такая функция , для которой неравенства (2) и (3) действительно выполняются. Для этого, конечно, нужно лишь, чтобы ограничения были не слишком жесткими. Практика показывает, что при итерационном применении условий (2) и (3) начальная оценка для функции сходится к некоторой функции, скажем , удовлетворяющей условиям (2) и (3). Кроме того, величина уменьшается при сужении указанных ограничений, если это производится постепенно в ходе выполнения итераций. На начальных итерациях важно предусмотреть “запас” при выборе нижних границ для изображения и его спектра. Слишком жесткие границы, при которых функция не удовлетворяет условиям (2) и (3), приводят к ложным результатам.

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