Цифровая обработка информации

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

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

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


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

(по другому, реперных или сопряженных) точек на изображениях, по которым можно осуществить координатную привязку снимков с одновременной геометрической коррекцией. (Точки на двух изображениях называются сопряженными, если они являются образами одной точки сцены [5.1, гл.13]). Например, аэрокосмический компьютерный мониторинг предполагает наличие дискретного по времени  наблюдения с небольшим временным интервалом, и поэтому, когда движущаяся камера фиксирует яркостный образ наблюдаемого объекта (оптическую поверхность) в виде последовательности изображений, то этот образ от снимка к снимку  деформируется вследствие перспективных искажений и изменения положения камеры. Геометрия соответствующих деформаций моделируется проективными преобразованиями, которые составляют более обширный класс, нежели известные  преобразования евклидовой геометрии (достаточно сказать, что длины и углы  в проективной геометрии не сохраняются, а параллельные линии могут пересекаться! [5.12]).

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

5.1. Геометрические преобразования на плоскости и в пространстве

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



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

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

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

                                       .                                                        (5.1)

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


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

5.1.1. Точки и прямые линии на плоскости - двойственность описаний

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

.

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

,

.

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

  .

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

,

.

Отсюда для координат точки пересечения  получаем соотношения, аналогичные вышеприведенным соотношениям для параметров линии :

  .

Здесь просматривается очень важная симметрия или двойственность

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

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



5.1.2. Однородные координаты

  Для преодоления отмеченных проблем описания геометрических объектов, а также для решения задач преобразования 3D-пространства и 2D-плоскости  в единообразном (матричном) виде вводится формализм так называемых однородных координат. Однородными координатами служат тройки чисел  (одновременно не равные нулю), связанные с обычными  координатами точек плоскости соотношением:

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

.

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

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

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

 или .

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

5.1.3. Евклидовы преобразования

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

               ,                                   (5.2)

с матрицей поворота на угол  вида и вектором трансляции (сдвига) .



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

                                                                (5.3)

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



Рис.5.1. Действие евклидова преобразования на пять точек плоскости

(сдвиг, поворот)

                   (5.4)

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

5.1.4. Аффинные преобразования



 Если матрицу вращения в (5.2) заменить  общей невырожденной матрицей , то получим преобразование

                                                        (5.5)

или, в матричном виде

                           ,                  



а)

              

б)

Рис.5.2.

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

б) исходное изображение (слева) и его аффинно-преобразованная копия  (справа);параметры аффинного преобразования:

и в однородных координатах

                                           (5.6)



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

Уравнения (5.5),(5.6) определяют общую форму записи хорошо известного аффинного преобразования (смотрите рис.5.2). Любое аффинное преобразование имеет обратное, которое также является аффинным. Произведение прямого и обратного преобразований дает единичное преобразование, оставляющее все на месте.  Аффинное преобразование является самым общим взаимно однозначным отображением плоскости на плоскость, при котором сохраняются прямые линии. Сохраняются также отношения длин отрезков, лежащих на одной прямой (или на параллельных прямых), и отношения площадей фигур. Параллельные прямые переходят в параллельные.

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

1.Растяжение (сжатие)

вдоль координатных осей, задаваемое в виде:



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

.

2.Поворот вокруг начальной точки на угол , описываемый формулой:



Матрица вращения (для однородных координат)

.

3.Перенос, задаваемый простейшими соотношениями:



Матрица переноса имеет вид

.

4.Отражение (относительно какой либо из осей, например оси абсцисс) задается при помощи формулы:



Матрица отражения, соответственно

.

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





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

евклидово преобразование  аффинное  проективное преобразование.

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



решив систему уравнений (5.5) относительно :



Если параметры таковы, что , то данное аффинное преобразование имеет неподвижную точку :



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



5.1.5. Проективные преобразования

Как выше уже было сказано, общими линейными преобразованиями  (в представлении однородными координатами)

                                                                                 (5.7)

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

.

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

коллинеацией.

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

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



=.                (5.8)



а)

 

б)

в)	г)

Рис.5.3.

а) действие проективного преобразования на пять точек плоскости;

б) исходное изображение; в) и г) проективно преобразованные образы; параметры проективного преобразования соответственно:



и



Однородные координаты векторов проективного пространства  и проективной плоскости соотносятся с неоднородными (евклидовыми) координатами векторов  и заданным выше образом:   и  .

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

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

что если расстояние от плоскости изображения до центра проекции равно , то координаты элементов изображения  соотносятся с пространственными координатами объекта следующим образом

                            .                                                 (5.9)

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


Заметим, что луч, проходящий через 2D точку , является направляющим вектором прямой, соединяющей точки   и . В то же время 3D точка  также лежит на этом луче (и представляет его). Припроизвольном (пробегающем всевозможные значения) получим координаты всех 3D точек на этом луче, соответствующих единственной точке на изображении в плоскости . По существу, каждая точка на изображении определяет луч, идущий от сцены в начало координат. Следовательно, однородные координаты точки (на снимке) в действительности представляют (и определяют) линию, проходящую через начало в евклидовом трехмерном пространстве . Такой набор всевозможных линий, проходящих через начало и формирует проективное 2D-пространство  в плоскости . (Можно, конечно, положить фокусное расстояние , поскольку различным значениямсоответствует разный масштаб изображения). В однородных координатах уравнения (5.9) естественно имеют вид (5.8)



Рис.5.4. Перспективная проекция

                                                  (5.9)

 Замечание. При , 3D точка  , в общем-то, определяет линию, параллельную плоскости  ,  не имеющую с ней точек пересечения и, следовательно, не имеет соответствия с какой-либо конечной точкой изображения.  Такие линии или однородные векторы могут, тем не менее, иметь смысл, если считать, что соответствующая им точка  удаляется на бесконечность в направлении, задаваемом этими координатами:

.

Мы можем добавить все такие точки к проективной плоскости. Эти точки называются «идеальными» или точками на бесконечности. На изображениях проективной плоскости добавленные точки на бесконечности формируют «линию горизонта» (см. рис.5.5). Существует разделение идеальных точек, обусловленное различными направлениями на плоскости; например, точки (1,0,0) и (0,1,0) связаны с горизонтальным и вертикальным направлениями (осями координат) соответственно. Можно также сказать, что все идеальные точки лежат на линии, называемой «идеальной линией» или линией на бесконечности, которая рассматривается, тем не менее, как и обычная линия.


Идеальная линия представляется в виде (0,0,1).

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



Рис.5.5. Проективная плоскость = аффинная плоскость + идеальные точки

(идеальная линия)

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

         

Таким образом, идеализация процесса формирования изображения камерой может быть представлена как перспективная проекция из в . Допустим, что 3D координаты точек объекта известны. Тогда, зная элементы матрицы , относящиеся к данному проективному преобразованию, точки пространственного объекта можно связать с соответствующими им координатами на снимке в виде (5.8)



         =.                       (5.10)

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

,

так что координаты изображения объекта имеют вид отношения

,

(5.11)

.

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


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

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

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

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

,

что влечет .   Матрица  аффинного  преобразования  имеет  вид

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



5.1.6. Полиномиальное преобразование

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


Поэтому привлекают полиномиальную аппроксимирующую функцию преобразования (рис.5.6)

                    ,

                    ,                  (5.12)   

где  - координаты точек эталонного снимка, - соответствующие им координаты на текущем (сопоставляемом) снимке.   

 

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

              





5.1.7. Оценивание параметров преобразования



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

,    .

Аналогично коэффициенты полиномов можно представить в векторной форме

,    .

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

,                        (5.13)

где                                       .                                   

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

                               .                                   (5.14)

Следовательно, искомые оценки имеют вид

                           ,.                             (5.15)



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

                                     (5.16)

                                     (5.17)

где -попарно ортогональные на заданном множестве точек ортогональные многочлены Чебышева, получаемые из последовательности методом ортогонализации Грама-Шмидта [5.3, гл.20]. Ортогональные многочлены до третьей степени имеют вид



Коэффициенты снова выбирают таким образом, чтобы минимизировать среднеквадратичную ошибку аппроксимации фактически «наблюдаемых» координат  их полиномиальной оценкой из (5.16),(5.17) для набора заданных узловых точек

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

            .                   (5.18)

5.2. Восстановление изображения в преобразованных

 координатах

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

Будем считать, что заданы два снимка (и ) одной и той же местности, полученные с некоторыми отклонениями точек съемки и условий     



Рис.5.7. Вычисленные координаты  (выделены серым тоном), наложенные на исходную дискретную целочисленную решетку

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



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

                 ,                               (5.19)

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

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

                            ,                                     (5.20)

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

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


В одномерном случае это прямоугольные, треугольные, B-сплайн функции и т.п. [5.2, разд.4.3]. При выборе соответствующего ядра  исходят из соображений как необходимой точности интерполяции, так и вычислительной эффективности. Понятно, что здесь одномерные функции должны быть преобразованы в двумерные функции. Общий подход состоит во введении так называемых «сепарабельных» интерполяционных функций в виде произведения двух одномерных функций. Сепарабельность во многих отношениях достаточно привлекательна в приложениях хотя и влечет неизотропность (за исключением гауссовых функций). Однако данные на квадратной решетке дискретизованы также не изотропно.

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



Рис.5.8. Интерполятор по ближайшему соседу с прямоугольным ядром. Справа график модуля Фурье-образа ядра. Пунктирной линией показан идеальный низкочастотный фильтр с частотой среза  .



Рис.5.9. Линейная интерполяционная функция и модуль ее Фурье-образа

(на правом рисунке пунктиром отмечен идеальный низкочастотный фильтр)

Здесь интерполированный сигнал представляется в виде

            ,                 (5.21)

где



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



Наиболее подходящим для интерполяции изображений является кубический B-сплайн (рис.5.10), поскольку в результате его применения получается функция, непрерывная и гладкая в узлах интерполяции. Ядро кубической свертки составляется из кусков кубических полиномов, определенных на подинтервалах  (-2,-1), (-1,0), (0,1), (1,2) по каждой из координат. Вне интервала (-2,2) интерполяционное ядро равно нулю. Двумерный кубический B-сплайн может быть записан как произведение двух одномерных интерполяционных функций по каждой из координат

 где  , так что

                                                         (5.22)

и , то есть ядро симметричное.

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

                             .                         (5.23)



Рис.5.10. Кубический B-сплайн и модуль его Фурье-образа.

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

           .                       (5.24)



Когда константа отрицательная, ядро (5.24) положительное в интервале от 0 до 1 и отрицательное в интервале от 1 до 2. Когда  возрастает, глубина боковых лепестков в интервале от 1 до 2 также увеличивается. Таким образом, с отрицательным значением свободной константы  интерполяционное ядро имеет вид усеченной sinc-функции. Выяснилось[5.11], что эта функция имеет более предпочтительные высокочастотные свойства, нежели кубический B-сплайн, и было предложено называть данную функцию высокоразрешающим интерполяционным кубическим сплайном.   Варьируя  значением  параметра

в пределах от -1 (рис.5.11.а) до -1/2 (рис.5.11.б), в каждом конкретном случае можно добиться приемлемой точности при передискретизации.



а)



б)

Рис.5.11. Высокоразрешающий интерполяционный кубический сплайн  и его частотная характеристика :

a) - a=-1; б) - а= - 0.5

5.3. Привязка изображений

В практике обработки изображений  задача поиска соответствия получила большое распространение и известна как проблема «поиска по образцу». Формально ее  можно рассматривать как процесс отождествления  эталонного изображения (образа фрагмента) на первом снимке с одним из множества образов фрагментов, лежащих в некоторой (задаваемой) области (зоне поиска) второго снимка . Алгоритмы установления сходства в своих основополагающих вариантах в той или иной степени связаны с получением характеристик стохастической взаимосвязи сравниваемых фрагментов изображений [5.2, гл.19].  Все они основываются на идеях корреляционной и спектральной теории сигналов, и для соответствующих критериев  получены экспериментальные характеристики основных процедур поиска по образцу.

5.3.1. Корреляционный критерий сходства

Будем считать, что изображение эталонного фрагмента (выбранного на снимке A и представляемого матрицей размером ), сравнивается с изображениями фрагментов снимка B в «зоне поиска»  размером . Перекрытие между фрагментами определяется шагом дискретной решетки  (в плоскости ), на которой заданы наблюдаемые переменные на A или  на B.


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

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

                                (5.25)

где и - параметры контраста и средней освещенности; k, l - параметры относительного сдвига образца и его аналога на контролируемом снимке; - шум;



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

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



                                    .

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

                                          (5.26)

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

            .                       (5.27)

Первый член выражения (5.27) - «энергия» эталонного сигнала, является величиной постоянной, не зависящей от параметров сдвига . Поэтому точка экстремума не изменится, если мы нормируем среднеквадратичную ошибку к энергии эталона 

    ,

и вместо минимума нормированной среднеквадратичной ошибки будем искать максимум коэффициента корреляции текущего фрагмента с эталоном



                .                              (5.28)

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

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



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

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

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



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

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

=  ,   ,

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

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

5.3.2. Локальное уточнение сдвига

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


Чтобы найти соответствующий вектор сдвига  будем считать, что кросс-корреляционная функция  в окрестности точки экстремума  разлагается в ряд Тейлора:



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





В результате приходим к матричному уравнению для параметров :

  .                            (5.30)

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

  .                    (5.31)

5.3.3. Кросс-спектральная мера сходства

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





Поскольку при смещении «образца» в области поиска  меняется только его фазовый спектр  то будем минимизировать функционал

(5.32)

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

                               (5.33)

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

                                                                                  (5.34)



так как в этом случае обратное дискретное фурье-преобразование



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



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



          a)                                                                  б)

Рис.5.12. Меры эффективности алгоритмов идентификации сходства:

а - поведение кросс-спектральной меры сходства вблизи ее экстремума: 1 - отклик фазового фильтра, 2 - отклик корреляционного обнаружителя;

б - оценки вероятностей ложной привязки (ВЛП) в зависимости от отношения амплитуды шума  к амплитуде фонового сигнала (в %): 1 - ВЛП в фазовом алгоритме, 2 - в корреляционном алгоритме.

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

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


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

 

Рис.5.13.Радиолокационные изображения местности, снятые при двух пролетах самолета (разные ракурсы съемки)

5.3.4. Привязка по локальным неоднородностям

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

 

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

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


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

ВОПРОСЫ К ГЛАВЕ 5



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

5.2. Выстройте иерархию геометрических преобразований. Какие из них сохраняют параллельность прямых? 

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

5.4. Чем отличается аффинная плоскость от евклидовой плоскости и от проективной плоскости?

5.5. В каких ситуациях рекомендуется применять полиномиальную аппроксимацию (в том числе и полиномы Чебышева) для описания геометрических деформаций?

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

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

5.8.  Опишите процедуру привязки изображений и методы ее ускорения. 

Фотограмметрия и стереовидение

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

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

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

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

Определим в трехмерном пространстве ортогональную правую систему координат OXYZ, начало которой совпадает с оптическим центром, ось OZ - с оптической осью камеры. Такая система называется стандартной системой координат камеры. Пусть плоскость изображения находится на расстоянии  от оптического центра. В этой плоскости зададим систему координат oxy с началом в главной точке и осями ox и oy, параллельными осям OX и OY соответственно (рис. 6.1). Легко убедиться, что в стандартной системе координат проекцией точки трехмерного пространства  с координатами  является точка  в плоскости изображения с координатами , причем

, .                                            



Рис.6.1. Система координат проективной камеры

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


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

 ,  ,

где  - координаты главной точки относительно начала координат фотоприемника (в естественных координатах фотоприемника);

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

В новой системе координаты проекции точки примут вид

                                   ,   .                                                      (6.1)

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

,                                                         (6.2)

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

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

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

                                                      ,                                                       (6.3)

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





Рис.6.2. Переход от глобальной системы координат к стандартной системе координат камеры.

На рис. 6.2. схематически показано преобразование координат. Здесь - углы, образованные осью  с осями ,  и  соответственно. Элементы первой строки матрицы   [6.1, п.14.10] содержат  косинусы этих углов: ,  ,  . Аналогично, вторая и третья строки матрицы содержат косинусы углов, образованных соответственно осями  и  с осями глобальной системы координат.

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

                  ,  ,  ,  , ,  ,                         (6.4)

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

                           , , .                           (6.5)

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

 Смысл вектора  ясен непосредственно из рисунка.

6.3. Стереоскопическая система

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


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

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

,  , .

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

.                                                (6.6)



Рис.6.3. Простейшая стереоскопическая система

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

, .                                         (6.7)

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



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

,                                              (6.8)

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



Рис.6.4. Система двух произвольно ориентированных камер

Используя (6.2) из (6.8) можно получить соотношение, связывающее координаты сопряженных точек (в координатах фотоприемника) :

,                                               (6.9)

где(см. (6.2))               ,  .                                (6.10)

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

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

,

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



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

,

разрешая которую относительно  и  получим

.            (6.11)

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

,    .                                (6.12)

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

6.4. Калибровка камеры

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

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

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

, .            

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



      ,     

,                                                   (6.13)

где ,  ,   ,

, ,  .

Очевидно, достаточно двух опорных точек, чтобы получить оценки внутренних параметров камеры:

,        ,

,   .

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

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

,       .

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

.                       (6.14)

Матрица  имеет размер  и называется калибровочной. Обозначим

,  ,

,  ,                      (6.15)

,                  .

Отметим два важных свойства этой матрицы, которые легко получить из (6.15), учитывая (6.4)  и (6.5):

 ,

.                                      (6.16)

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

,

,

,                                                     (6.17)

.         

Кроме того, непосредственно из (6.14) следует, что

,

,

,

,                                          (6.18)



,

.

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

Рассмотрим сначала линейный метод оценивания матрицы . Запишем матричное уравнение (6.14) как систему трех обычных уравнений

,

,

,

или, подставляя в два первых уравнения значение  из третьего,

,       

.

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

.   (6.19)

Представим эту систему в матрично-векторном виде:

,                                              (6.20)

где



,  .

Сначала рассмотрим некоторые общие особенности этой системы. Очевидно, что одним из решений этой системы является тривиальное , которое не имеет физического смысла. Известно [6.2, с.153], что если однородная линейная система имеет хотя бы  одно ненулевое решение, то она имеет бесконечное множество решений, причем, если  - решение, то и  , где - произвольное число, тоже является решением. Здесь необходимо различать два случая.

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


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

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



Рис.6.5. Калибровочный объект

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

,                                    

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

          В этой ситуации можно снова воспользоваться МНК, согласно которому в качестве оценки вектора  следует принять такой, который минимизирует значение функционала  при условии  (см.(6.16)). Поскольку функционал  представляет собой квадратичную форму с неотрицательно определенной симметричной матрицей , то  минимума по  при условии   он достигает, если  - собственный вектор матрицы , соответствующий ее минимальному собственному числу [6.4, с.64]. Заметим, что  .  (Собственным вектором матрицы  называется вектор , который удовлетворяет уравнению , где - скалярный множитель, называемый собственным числом матрицы  ).

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

,

которая представляет сумму квадратов этих расстояний, и минимизируем ее по  с учетом ограничений (6.16).


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

6.5. Взаимное ориентирование

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

Рассмотрим снова пару камер, внутренние параметры которых известны, но неизвестны внешние параметры (матрица  и вектор ). Умножив обе части выражения (6.8) слева сначала векторно на , а затем скалярно на , получим

.

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

                                     (6.21)

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

.                                  (6.22)

Эквивалентные соотношения (6.21) и (6.22) являются основой для  оценивания матрицы  и вектора . Предположим, что известны координаты  пар сопряженных точек и, соответственно,  пар векторов  и   .

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

,   где  .                     (6.23)



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

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

,

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

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

Значительно более сложной задачей является задача оценивания матрицы . Один из возможных алгоритмов состоит в следующем [6.6]. Известно [6.1, п.14.10], что матрица  может быть представлена в виде , где

 ,   ,

.

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


Затем в окрестности этих значений для уточнения положения минимума можно воспользоваться одним  из известных методов минимизации [6.5, гл.V] (например, наискорейшего спуска, Ньютона, Маркуардта).

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

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

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

6.6. Поиск сопряженных точек

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

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


Обратимся к соотношению (6.21). Известно [6.8, Приложение], что векторное произведение  можно представить как , где матрица  имеет вид:

.

Тогда (6.21) можно представить как

                                              (6.24)

или

                                                          ,

где

                                  .                                     (6.25)

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

                                                (6.26)

с вектором коэффициентов ,

а в плоскости изображения второй камеры - прямую

                                                (6.27)

с вектором коэффициентов .

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

Этот результат имеет наглядную геометрическую иллюстрацию (рис.6.6).



Рис. 6.6. Эпиполярные линии

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



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

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

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

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

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


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



Рис.6.7. Проективные и яркостные искажения

Одним из наиболее распространенных методов идентификации является корреляционный. Детальное описание корреляционного метода приведено в гл.6. При незначительных отличиях в ракурсах съемки и на достаточно гладких поверхностях от него можно ожидать хороших результатов. На рис.6.8 (а и б)  приведена пара вертолетных снимков района озера Карымское (Камчатка), на которые нанесена прореженная сетка сопряженных точек, найденных корреляционным методом, а на рис.6.8.в показана поверхность, «натянутая» на трехмерные точки, координаты которых вычислены методом, изложенным в п.6.2. Следует заметить, что построение поверхности, лежащей на заданном множестве трехмерных точек, представляет собой самостоятельную нетривиальную задачу, рассмотрение которой выходит за рамки темы, обсуждаемой в этой главе.

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

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

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


В приведенном примере для построения рельефа  найдено более 7000 точек.

а)	б)
в)
Рис.6.8. Восстановление поверхности трехмерной сцены по стереопаре (приведено с любезного согласия авторов [6.9])

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

Влияние яркостных искажений можно попытаться устранить, отказавшись от сопоставления фрагментов исходных изображений, и перейдя к сопоставлению фрагментов контурных изображений, где выделены края, т.е. переходы между областями с постоянной яркостью. На таких переходах, как известно,  градиент яркости имеет максимум, а лапласиан обращается в нуль. Это обстоятельство и используется для выделения краев [6.10, п.17.4]. Существенным моментом является то, что контуры на изображениях непрерывных поверхностей должны располагаться в одинаковом порядке, как это видно на рис. 6.9 (изображения точек ,  и ). Это ограничение позволяет в некоторой степени ослабить влияние масштабных искажений. Однако и здесь, как и в корреляционном методе, существует проблема ложного отождествления. Связана она, в частности, с тем, что видимая часть поверхности зависит от ракурса наблюдения (поэтому, например, на рис.6.9 изображение точки  в правой камере отсутствует), и на разных изображениях могут пропадать разные края.


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

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

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



Рис.6.9. Упорядочение сопряженных точек в случае непрерывных поверхностей



ВОПРОСЫ К ГЛАВЕ 6

6.1. Какой смысл имеет знак при  в выражении (6.1)?

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

6.3. Обратимся к рис. 6.3. Какой вид будут иметь матрицы , ,  и вектор  в выражении (6.9) для ситуации, приведенной на рисунке?

6.4. Пусть в ситуации, изображенной на рис. 6.3, векторы   и  внутренних координат проекций точки  в плоскостях изображений левой и правой камер известны. Получите оценки трехмерных координат точки  в системах координат правой и левой камер, пользуясь выражениями (6.11) и (6.12). Сравните полученный результат с (6.6) и (6.7). Объясните отличия.

6.5. Почему в предыдущем вопросе векторы   и  имеют одинаковые - компоненты?

6.6. Докажите справедливость соотношений (6.17).

Дискретизация и квантование непрерывных изображений

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

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



1.1. Дискретизация непрерывных изображений



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

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

,                                     (1.1)

где - соответственно вертикальный и горизонтальный шаги или интервалы дискретизации. Рис.1.1 иллюстрирует расположение отсчетов на плоскости  при прямоугольной дискретизации.

Рис. 1.1. Расположение отсчетов при прямоугольной дискретизации

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



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

,        (1.2)

которому отвечает двумерное обратное непрерывное преобразование Фурье:

.       (1.3)

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

.   (1.4)

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



Выполняя замену переменных по правилу  , добиваемся независимости области интегрирования от номеров  и :



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

     (1.5)

Теперь выражение (1.5) имеет форму обратного преобразования Фурье, следовательно стоящая под знаком интеграла функция

            (1.6)

является двумерным спектром дискретного изображения. В плоскости ненормированных частот выражение (1.6) имеет вид:

         (1.7)

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

а)	б)
Рис. 1.2. Частотные спектры непрерывного и дискретного изображений

<


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

  при ,  .     (1.8)

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

,     ,                                 (1.9)

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

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

              (1.10)

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

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


Двумерная импульсная характеристика восстанавливающего фильтра, которую легко получить при помощи обратного преобразования Фурье от (1.10), имеет вид:

.

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

,

после выполнения свертки находим:

.   (1.11)

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

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

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


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

а)	б)
в)	г)
Рис. 1.3. Влияние интервала дискретизации на восстановление изображения «Отпечаток пальца»

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

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

а)	б)
в)	г)
Рис. 1.4. Влияние интервала дискретизации на восстановление изображения «Портрет»

<


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

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

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



Рис.1.5.Функция, описывающая квантование

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

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


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

,                                  (1.12)

в предположении, что яркость  - случайная величина с известной плотностью вероятности .

Cреднеквадратическая ошибка квантования  (1.12) равна

.        (1.13)

Дифференцируя (1.13) по переменным  ,   и приравнивая производные нулю, получаем   нелинейных уравнений

,            

(1.14)

.                               

Следует отметить, что крайние пороги  и  определяются динамическим диапазоном яркости. Уравнения (1.14) нетрудно привести к виду

,                                  (1.15)

.                                       

Из (1.15) следует, что пороги  должны располагаться по середине между двумя соседними уровнями  и . Решение этих уравнений можно найти итеративным способом. Оптимальный квантователь, удовлетворяющий критерию (1.12), называется квантователем Ллойда-Макса [1.3, 1.5], а среднеквадратическая ошибка для такого квантователя равна

                        (1.16)

При равномерном распределении яркости нелинейные уравнения (1.15) можно представить в виде [1.3]

,

,      

                          

а среднеквадратическая ошибка равна    .

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

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


На рис.1.7.а и 1.7. б приведены результаты равномерного квантования изображения «Портрет» соответственно на 256 и 14 уровней квантования.



Рис.1.6. К механизму возникновения ложных контуров

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

Рис.1.8. Результат неравномерного квантования	Рис.1.9. Гистограмма изображения “Портрет”

В темных частях изображения на рис. 1.7.б заметны ложные контуры. Использование квантователя Ллойда-Макса позволяет существенно снизить их уровень (см. рис. 1.8, где число уровней квантования также равно 14). На рис. 1.9 приведена гистограмма яркости изображения «Портрет» при 256 уровнях квантования и отмечены пороги  при . Из рисунка следует, что чаще квантуются те области динамического диапазона, в которых сгруппированы значения яркости отсчетов.

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



                   

Рис.1.10. Квантование с предварительным нелинейным преобразованием

Для разрушения ложных контуров Робертс предложил перед равномерным квантованием к отсчетам яркости добавлять шум с равномерной плотностью распределения вероятностей. Добавленный шум переводит одни отсчеты изображения на уровень выше, а другие на уровень ниже. Тем самым разрушаются ложные контуры. Дисперсия добавляемого шума должна быть небольшой, чтобы не привести к искажениям, воспринимаемым как «снег» на изображении, и в то же время достаточной для разрушения ложных контуров. Обычно используют равномерно распределенный шум на интервале . Результаты равномерного квантования на 14 и 8 уровней изображения «Портрет» с предварительным добавлением    шума    приведены    на рис.1.11.а    и    1.11.б.


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

а)	б)
Рис.1.11. Результаты равномерного квантования с предварительным добавлением шума

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

.

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



а)

б)	в)
Рис.1.12.Растрирование изображений

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

Улучшение визуального качества изображений путем поэлементного преобразования

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

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


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

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

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

,                                         (2.1)



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





,                                                    (2.2)



одинаковой для всех точек кадра.

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



2.1. Линейное контрастирование изображения

          Задача контрастирования связана с улучшением согласования динамического диапазона изображения и экрана, на котором выполняется визуализация. Если для цифрового представления каждого отсчета изображения отводится 1 байт (8 бит) запоминающего устройства, то входной или выходной сигналы могут принимать одно из 256 значений. Обычно в качестве рабочего используется диапазон 0...255; при этом значение 0 соответствует при визуализации уровню черного, а значение 255 - уровню белого. Предположим, что минимальная и максимальная яркости исходного изображения равны   и   соответственно. Если эти параметры или один из них существенно отличаются от граничных значений яркостного диапазона, то визуализированная картина выглядит как ненасыщенная, неудобная, утомляющая при наблюдении. Пример такого неудачного представления приведен на рис. 2.1.а, где диапазон яркостей имеет границы  , .

          При линейном контрастировании используется линейное поэлементное преобразование вида:

,                                                  (2.3)



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

  ,



относительно параметров преобразования   и  , нетрудно привести (2.3) к виду: 

.



Результат линейного контрастирования исходного изображения, представленного на рис. 2.1.а, приведен на рис.2.1.б при   и

а)	б)

Рис 2.1  Пример линейного контрастирования

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



2.2. Соляризация изображения



          При данном виде обработки преобразование  (2.2) имеет вид [2.1]:



,

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

Рис. 2.2. Функция, описывающая соляризацию

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

а)	б)
Рис. 2.3. Пример соляризации

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

2.3. Препарирование изображения

          Препарирование представляет собой целый класс поэлементных преобразований изображений. Характеристики применяемых на практике процедур препарирования приведены на рис.2.4. Остановимся на описании некоторых из них.

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





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

а)	б)	в)
г)	д)	е)
ж)	з)	и)
	к)
Рис.2.4.Примеры преобразований, используемых при препарировании

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

Рис.2.5. К выбору порога бинарного квантования

называемой бинарной сегментации. Здесь же ограничимся обсуждением частного, но практически важного случая. Иногда при обработке приходится иметь дело с изображениями, хранимыми как полутоновые, но по своему содержанию мало отличающимися от бинарных. К ним относятся текст, штриховые рисунки, чертежи, изображение отпечатка пальца, пример которого приведен на рис.2.6.а. Плотность вероятности , описывающая распределение яркости такого изображения, может содержать два хорошо разделяющихся пика. Интуитивно понятно, что порог бинарного квантования следует выбирать посредине провала между этими пиками, как это показано на рис.2.5. Замена исходного полутонового изображения бинарным препаратом решает две основные задачи. Во-первых, достигается бóльшая наглядность при визуальном восприятии, чем у исходного изображения. Во-вторых, ощутимо сокращается объем запоминающего устройства для хранения изображения, поскольку бинарный препарат для записи каждой точки бинарного изображения требует лишь 1 бит памяти, в то время как полутоновое изображение для решения той же задачи при наиболее часто применяемом формате представления - 8 бит.


Пример бинаризации изображения отпечатка пальца приведен на рис.2.6.б.

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

Рис.2.6. Пример бинаризации изображения

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

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

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

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


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

а)	б)
в)
Рис.2.7. Примеры препарирования изображения

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

         

2.4. Преобразование гистограмм, эквализация

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

Рис.2.8. Нелинейное преобразование случайной величины

имеющей однозначную обратную функцию . Предположим, что





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

,



где модули учитывают зависимость вероятностей от абсолютных длин интервалов (и независимость от знаков

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



.                                       (2.4)



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



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

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

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

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

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



                   (2.5)

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



.



Подставляя это выражение в условие вероятностной эквивалентности



=,



после простых преобразований получаем соотношение



,                                   (2.6)

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

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

           (2.7)



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

,                                         (2.8)



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

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

          На первом этапе измеряется гистограмма

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


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

.                                               



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

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

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

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

а) исходное изображение	б) результат обработки
Рис. 2.9. Пример эквализации изображения

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


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



2.5. Применение табличного метода

при поэлементных преобразованиях изображений

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

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

табл.2.1. 

			.	.	.	.	.
			.	.	.	.	.
Табл.2.1. Табличное задание функции

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



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

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





ВОПРОСЫ К ГЛАВЕ 2



2.1. В чем состоит сущность поэлементной обработки изображений ?

2.2. Как определяются параметры преобразования изображения при его линейном контрастировании ?

2.3. Поясните механизм действия поэлементных преобразований, применяемых при препарировании изображения (рис.2.4).

2.4. Каков механизм появления ложных контуров при применении                  пилообразного контрастирования? При каких еще процедурах могут возникать ложные контуры?

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

2.6. Докажите справедливость преобразования  (2.7) для получения гиперболического распределения  (2.6).

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

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

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

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

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


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

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

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

3.1. Оптимальная линейная фильтрация.

Уравнение Винера-Хопфа



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

.                 (3.1)

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

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


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

а)	б)	в)
Рис. 3.1  Примеры окрестностей различных видов

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

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

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

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


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

          При линейной фильтрации выходной эффект определяется линейной комбинацией входных данных:

.                            (3.2)

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

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

,             (3.3)

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

          Оптимизационную задачу (3.3) нетрудно свести к решению уравнения или системы уравнений.


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

.          (3.4)

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

,     .

С их учетом (3.4) примет более компактный вид:

                           (3.5)

Считая автокорреляционную   и взаимно корреляционную   функции известными, замечаем, что (3.5) представляет собой линейное относительно искомых коэффициентов  алгебраическое уравнение. Число неизвестных в этом уравнении равняется числу точек    в окрестности   и является конечным в случае КИХ-фильтра и бесконечным при БИХ-фильтрации. Ограничимся в данном параграфе рассмотрением КИХ-фильтрации. Линейное алгебраическое уравнение со многими неизвестными имеет бесконечное множество решений. Если повторить дифференцирование  (3.3)  по остальным   неизвестным, то получим  еще    уравнений, отличающихся друг от друга левыми частями  и коэффициентами   в правых частях, т.к. определяющие их корреляции вычисляются каждый раз в различных точках. В результате образуется система   линейных алгебраических уравнений с  неизвестными, называемая в теории фильтрации уравнением Винера-Хопфа:

                            (3.6)

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

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

,                                (3.7)

где  - средний квадрат ошибок фильтрации.

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


Вычислив математическое ожидание от обеих частей (3.2), находим:

,                                     (3.8)

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

,                                          (3.9)

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

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

3.2. Масочная фильтрация изображений

при наличии аддитивного белого шума



          Распространенным видом помехи является белый шум, аддитивно воздействующий на изображение. Наблюдаемое в этом случае изображение (3.1) имеет вид:

,                            (3.10)

а корреляционная функция шума     описывается выражением:

.

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

,

где  - дисперсия, а   - нормированная корреляционная функция полезного сигнала. Нетрудно видеть, что в этих условиях взаимная корреляционная функция   совпадает с корреляционной функцией полезного сигнала  . Поэтому уравнение Винера -Хопфа  (3.6) приводится к виду:



,             (3.11)

где   - отношение дисперсий сигнала и шума.

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

. Сравнивая это соотношение с (3.7), окончательно получаем:

,

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

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

.

          В практике цифровой обработки изображений широко используется масочная фильтрация. Ее линейная разновидность является одним из вариантов двумерной КИХ-фильтрации. В качестве маски используется множество весовых коэффициентов, заданных во всех точках окрестности  , обычно симметрично окружающих рабочую точку кадра. Распространенным видом окрестности, часто применяемым на практике, является квадрат 33 с рабочим элементом в центре, изображенный на рис. 3.1.б. Применяют различные разновидности масок, одним из эвристических вариантов является равномерная маска, все девять весовых коэффициентов которой равны 1/9. Такой выбор коэффициентов отвечает условию сохранения средней яркости  (3.9) и поэтому в процессе обработки центрировать изображение не требуется.

          Визуально эффективность фильтрации можно оценить с помощью рис.3.2. На рис. 3.2.а показан зашумленный портрет (изображение без шума приведено на рис. 1.3.а) при отношении сигнал/шум равном -5дБ.


Результат масочной фильтрации при оптимальном виде ИХ, найденной из (3.11), приведен на рис.3.2.б. Результат фильтрации, выполненной равномерным масочным оператором не приводится, поскольку с визуальной точки зрения он мало отличается от рис.3.2.б. При этом, однако, с количественной точки зрения различия достаточно заметны: если при оптимальной КИХ относительная ошибка , то при равномерной КИХ она возрастает почти на 30% и составляет  . Различие резко возрастает при более высоком уровне шума. Так, например, при отношении сигнал/шум равном -10дБ  имеем   и    , т.е. применение равномерной КИХ вместо оптимальной приводит в этом случае к увеличению ошибок более чем вдвое.

а)	б)
Рис. 3.2. П.ример масочной фильтрации при

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

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

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


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

3.3. Рекуррентная каузальная фильтрация изображений



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

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

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

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

.

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


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

.                              (3.12)

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

,                     (3.13)

где  - параметр, лежащий в пределах  , и получивший название коэффициента усиления фильтра Калмана. Первое слагаемое в алгоритме (3.12) определяет вклад в оценку сигнала на текущем  -м шаге фильтрации, вносимый оценкой предыдущего шага, и называется одношаговым прогнозом. Второе учитывает влияние текущего наблюдения   и называется новой информацией. Коэффициент усиления  определяет чувствительность фильтра к этой новой информации. При высоком уровне шума параметр    приближается к нулю. При этом, кроме общего снижения ИХ, увеличивается параметр  , приближаясь к единице. Это означает удлинение импульсной характеристики и, следовательно, сужение полосы пропускаемых фильтром частот. Очевидно, эти свойства ИХ способствуют эффективной фильтрации шума. При снижении шума, наоборот,  стремится к единице, - к нулю, что соответствует расширению полосы частот до бесконечности. Здесь уместно отметить, что фильтрация не только ослабляет действие шума, но, к сожалению, и вносит так называемые динамические искажения в полезный сигнал. Механизм этих искажений очень прост и заключается в неравной передаче на выход фильтра различных спектральных компонент сигнала. Фильтр Калмана “ведет себя” вполне разумно, когда при исчезновении шума на входе расширяет полосу пропускаемых частот до бесконечности, поскольку именно при этом условии исчезают и динамические ошибки фильтрации.

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

  .      (3.14)

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


Для определения параметра   (или  ), остающегося единственным неизвестным параметром двумерного фильтра, воспользуемся уравнением Винера-Хопфа в форме (3.4), переписав его в виде: 

.

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

.                         (3.15)

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

,                                             (3.16)

для чего достаточно вычислить левую часть этого выражения с учетом (3.2) и (3.15).

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

.                                       (3.17)

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

             (3.18)

относительно параметра  , численное решение которого не представляет проблемы. Анализируя (3.18), легко заметить, что при   это уравнение

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

          Подставив в (3.7) выражения ИХ (3.14) и корреляционной функции (3.17), можно получить следующую формулу для среднего квадрата ошибок фильтрации:



.

          Подставив далее выражение ИХ (3.14) в (3.2), можно привести выражение отклика фильтра к виду :

.       (3.19)

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

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

          Пример применения описанного двумерного фильтра показан на рис. 3.3, где представлен результат эксперимента с тем же портретом и при том же отношении сигнал/шум  -5 дБ, что и при испытании масочного фильтра.

а)	б)
Рис. 3.3. Пример двумерной рекуррентной фильтрации

Поэтому здесь не приводится показанное на рис.3.2.а входное изображение с шумом. Результат двумерной рекуррентной фильтрации представлен на рис.3.3.а, а на рис.3.3.б для сравнения повторен результат оптимальной масочной фильтрации (рис.3.2.б). Визуальная оценка говорит в пользу двумерного рекуррентного фильтра, поскольку уровень остаточного шума на рис.3.3.а ниже. Сравнение по среднему квадрату ошибок совпадает с субъективной оценкой: величина  при масочной фильтрации составляет, как говорилось ранее, 0.309, а при двумерной рекуррентной - 0.29.


Различие заметно усиливается при более высоком уровне шума. Так, при отношении сигнал/шум -10 дБ имеем соответственно   равным 0.57 и 0.43.

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

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

3.4. Применение фильтра Винера для некаузальной

двумерной фильтрации



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

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



3.4.1. Двумерное дискретное преобразование Фурье

          Обозначим через

                           (3.20)

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

.                            (3.21)

Если сигнал   существует только внутри прямоугольника  со сторонами элементов

а)	б)
Рис. 3.4. Реальное (а) и периодически продолженное (б) изображения

(рис. 3.4.а), то сигнал   определен на всей плоскости  и является на ней прямоугольно-периодическим (рис. 3.4.б).

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

.  (3.22)

Базисные функции этого двумерного представления - двумерные комплексные экспоненты (иногда называемые комплексными синусоидами)

,                     (3.23)

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

 Коэффициенты Фурье ряда  (3.22) образуют двумерный частотный спектр сигнала   и определяются формулой прямого преобразования Фурье:

                       (3.24)

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

          Заметим, что для точного представления дискретного сигнала с двумерным периодом   элементов согласно формулам БПФ достаточно конечного  числа базисных функций (3.23) - ряд (3.22) является конечным.


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

          Остановимся на наиболее существенных свойствах двумерного дискретного спектра Фурье. Вычислим спектральные коэффициенты (3.24) в частотных точках  :

                                 (3.25)

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

,

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

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

,

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


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

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

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

3.4.2. Циклическая свертка

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


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

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

.                           (3.26)

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

,

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

.                             (3.27)

также является периодическим. Покажем, что именно его спектр  определяется выражением:

.                                     (3.28)

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

.

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



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

,                                            (3.29)

                                            (3.30)

в обычном и циклическом вариантах соответственно. Рис. 3.5 поясняет процесс

а)	б)	в)
Рис. 3.5. Сравнение обычной и циклической свертки

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



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

3.4.3. Решение уравнения Винера-Хопфа в циклическом приближении

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

                   (3.31)

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

                     (3.32)

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

.                           (3.33)

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

.                                       (3.34)

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


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

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

                                (3.35)

                                   (3.36)

Рис. 3.6 иллюстрирует формирование сумм, входящих в правые части этих равенств при некотором произвольном значении сдвига   и достаточно большом значении интервала наблюдения  . Здесь показаны некаузальная ИХ  и корреляционная функция входного сигнала, а штриховкой условно отмечена область суммирования. Из сравнения рис. 3.6.а, соответствующего (3.35), и рис. 3.6.б, соответствующего (3.36), видно, что, хотя во втором случае область суммирования и является двухсвязной, это не приводит к различию результатов

а)	б)
Рис. 3.6. Сравнение обычного и циклического уравнений Винера-Хопфа

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

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


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

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

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

а)	б)
Рис. 3.7. Пример винеровской фильтрации шума при

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

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



          Анализ эффективности метода будет неполным, если не сделать оценки вычислительной эффективности реализующей его процедуры. Для вычисления ДПФ разработаны эффективные вычислительные методы, воплощенные в процедурах быстрого преобразования Фурье (БПФ). Количество комплексных умножений, составляющих основную трудоемкость двумерного БПФ, оценивают числом   [3.2]. Поскольку полный цикл обработки предполагает выполнение прямого и обратного БПФ, то это число следует удвоить. По отношению к одному элементу кадра число умножений, таким образом, составляет  . При  число умножений в каждой точке кадра равно 32. Для сравнения напомним, что, например, рекуррентный двумерный фильтр, описанный выше, реализуется всего тремя вещественными умножениями в каждой точке кадра (при различных значениях одношагового коэффициента корреляции изображения по строкам и по столбцам - четырьмя умножениями).

3.5. Байесовская фильтрация изображений



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

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


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

3.5.1. Сущность байесовской фильтрации

          Полагаем, что на входе фильтра действует сигнал

,            (3.37)

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

,                          (3.38)

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

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



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

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

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

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

                                 (3.39)

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


На пояснении последнего остановимся подробнее. Последовательность является марковской, если ее совместное распределение вероятностей может быть представлено в виде:

.                               (3.40)

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

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

.                         (3.41)

В соотношение (3.41) входит распределение последнего элемента и цепочка одношаговых распределений перехода в обратном времени , не совпадающих с .

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

           (3.42)

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


Это является результатом того, что апостериорное распределение вероятностей для произвольного элемента последовательности может быть представлено в виде [3.7]:

.                         (3.43)

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

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

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

.            (3.44)

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


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

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

3.5.3. Двухэтапная марковская фильтрация изображений

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

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


Рис. 3.8 иллюстрирует геометрию задачи.

Рис. 3.8. Геометрия использования данных при двухэтапной фильтрации

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

Будем, кроме того, рассматривать такие случайные поля , которые обладают свойством условной независимости. Это означает, что совместное распределение всех его элементов , расположенных на “кресте”  (рис. 3.8), можно представить в виде:

,                (3.45)

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

.          (3.46)

Используя эту математическую модель изображения в случае независимой помехи , можно одноточечное апостериорное распределение представить в виде [3.6]:

.                                    (3.47)

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

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


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

          Соотношение (3.47) дает возможность выполнить двумерную обработку изображения в виде некоторой совокупности одномерных процедур. Весь цикл вычислений можно представить следующим образом. Выполняется обработка всех строк изображения в прямом направлении (слева направо), в результате чего в каждой точке образуется распределение . При этом используются одномерные рекуррентные процедуры, описанные выше. Далее происходит повторное сканирование строк, но в “обратном времени” - справа налево, в процессе которого вычисляются распределения . Затем изображение аналогично дважды обрабатывается по столбцам - сверху вниз и снизу вверх, в результате чего определяются частные АРВ  и . Вычислением одноточечного АРВ  завершается первый этап обработки. На втором этапе происходит объединение всех частных АРВ в каждой точке кадра в окончательное АРВ, а также на его основе вычисляются точечные оценки изображения .

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



Структура распределений очень сильно влияет на требуемые объем вычислений и ресурс памяти. Имеются очень “удобные” в этом смысле виды распределений. Например, если для описания изображения применима модель случайного поля с гауссовским распределением, то для представления каждого из частных и финального АРВ в (3.47) требуется наличие всего двух параметров - математического ожидания и дисперсии. Именно это и определяет конкретный характер и количество вычислений в процессе фильтрации, а также объем необходимой памяти.

Другим примером такого рода может служить математическая модель бинарного случайного поля, которое в различных точках принимает значения  или . Такое описание также является очень экономичным, поскольку АРВ содержит всего две вероятности  и , непосредственное вычисление которых и выполняется при помощи (3.47).

Существует отдельный вопрос, связанный с применимостью марковских двумерных моделей (3.45), (3.46), позволяющих  построить эффективные двухэтапные процедуры. Его изучение является достаточно непростой теоретической задачей. В частности, в работах [3.6.,3.8] установлено, что и для гауссовских, и для бинарных случайных полей необходимым и достаточным условием применимости (3.45) является возможность представления двумерных корреляционных функций этих полей в разделимом виде, т.е. в виде произведения двух множителей, один из которых описывает корреляцию изображения по строке, а второй - по столбцу. Дополнительные требования, вытекающие из (3.46), сводятся к существованию марковских свойств у одномерных последовательностей в горизонтальном и вертикальном сечениях изображения. В двух указанных примерах наличие таких свойств связано с экспоненциальным видом корреляционных функций этих одномерных сечений изображения.

На рис. 3.9 приведены результаты экспериментальной проверки двумерных двухэтапных алгоритмов фильтрации изображения. На рис. 3.9.а показано тестовое бинарное изображение “острова”, на рис. 3.9.б - изображение, искаженное белым гауссовским шумом (отношение сигнал/шум дБ).


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

а)	б)	в)
г)	д)	е)
Рис. 3.9. Двухэтапная марковская фильтрация изображения

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

3.6. Медианная фильтрация



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


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

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

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

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

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



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

а)	б)
Рис. 3.10. Примеры окон при медианной фильтрации

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

.                               (3.48)

Рассмотрим пример. Предположим, что выборка имеет вид: , а элемент 250, расположенный в ее центре, соответствует текущей точке фильтрации  (рис. 3.10). Большое значение яркости в этой точке кадра может быть результатом воздействия импульсной (точечной) помехи. Упорядоченная по возрастанию выборка имеет при этом вид {45,55,75,99,104,110,136,158,250}, следовательно, в соответствии с процедурой (3.48), получаем .


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

          Из (3.48) следует, что действие МФ состоит в “игнорировании” экстремальных значений входной выборки - как положительных, так и отрицательных выбросов. Такой принцип подавления помехи может быть применен и для ослабления шума на изображении. Однако исследование подавления шума при помощи медианной фильтрации показывает, что ее эффективность при решении этой задачи ниже, чем у линейной фильтрации [3.9].

Результаты экспериментов, иллюстрирующие работу МФ, приведены на рис. 3.11. В экспериментах применялся МФ, имеющий квадратную апертуру со

стороной равной 3. В левом ряду представлены изображения, искаженные помехой, в правом - результаты их медианной фильтрации. На рис. 3.11.а и рис. 3.11.в показано исходное изображение, искаженное импульсной помехой. При ее наложении использовался датчик случайных чисел с равномерным на интервале [0, 1] законом распределения, вырабатывающий во всех точках кадра независимые случайные числа. Интенсивность помехи задавалась вероятностью  ее возникновения в каждой точке. Если для случайного числа , сформированного в точке , выполнялось условие , то яркость изображения в этой точке  замещалась числом 255, соответствующим максимальной яркости (уровню белого). На рис. 3.11.а действием импульсной помехи искажено 5 % (=0.05), а на рис. 3.11.в - 10 % элементов изображения. Результаты обработки говорят о практически полном подавлении помехи в первом случае и о ее значительном ослаблении во втором.

а)	б)
в)	г)
д)	е)
Рис. 3.11. Примеры медианной фильтрации

<


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

	масочный фильтр с оптимальн. КИХ	масочный фильтр с равномерн. КИХ	двумерный рекуррентн. фильтр	двумерный фильтр Винера	медианный фильтр
	0.309	0.395	0.29	0.186	0.539
	10.2	8.0	10.9	17.0	5.86
Табл.3.1. Сравнение эффективности подавления шума при фильтрации изображений,  дБ

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

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

Восстановление изображений

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



и т.п. Устранение или ослабление искажений с целью повышения резкости относится к задаче восстановления изображений.

Наиболее общая схема формирования изображения представлена на рис. 4.1,



Рис.4.1. Схема формирования изображения



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

.                                          

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

4.1. Модели изображений и их линейных искажений

4.1.1. Формирование изображений

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

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

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



Рис.4.2. Линейная модель формирования изображения

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

,

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



         (4.1)

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

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

                     (4.2)

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



Рис.4.3. Относительные размеры изображения и ФРТ

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


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

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

,        (4.3)

где

.      (4.4)

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

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

Спектр линейно-искаженного изображения  равен произведению спектра  исходного изображения  и передаточной функции  искажающей системы:

,                              (4.5)



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

Рассмотрим импульсные и частотные характеристики формирующих систем при смазе и расфокусировке.

4.1.2. Размытие вследствие движения (смаз)

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

                            

,

где длина смаза  - равна произведению скорости движения камеры на время экспозиции.  Соответственно в дискретном случае ФРТ  смаза равна

                (4.6)

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

Рис.4.4. Взаимное расположение  изображения и ФРТ при смазе	Рис.4.5. Изображние модуля частотной характеристики искажающей системы

Дополняя ФРТ (4.6) нулями до размеров кадра исходного изображения и применяя двумерное ДПФ, получим частотную характеристику искажающей системы:

       (4.7)

Изображение модуля  приведено на рис.4.5 при   и размерах исходного изображения элементов.

На рис.4.7 приведен искаженный вариант исходного изображения «Сатурн» (рис.4.6). Горизонтальный смаз составляет 15 элементов. Исходное изображение содержит  элементов, а искаженное -  элементов.

Рис.4.6.Исходное изображение “Сатурн”	Рис.4.7. Смазанное изображение “Сатурн”

4.1.3. Расфокусировка

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

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

                            (4.8)

Из (4.8) следует, что размеры кадра . Взяв двумерное преобразование Фурье от  (4.8), получим передаточную функцию оптической системы

,                (4.9)

где  - функция Бесселя первого порядка.

В дискретном случае ФРТ (4.8) имеет вид:

                     (4.10)

На рис.4.8 и 4.9 показаны ФРТ для тонкой линзы (4.10) и модуль ее передаточной функции при радиусе кружка размытия  и размерах  кадра изображения  элементов.

Рис.4.8. ФРТ  тонкой линзы	Рис.4.9. Изображение модуля частотной характеристики тонкой линзы

<


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

,                                (4.11)

который в дискретном случае имеет вид

,                                (4.12)

где -  нормирующий коэффициент,  - коэффициент пространственной нерезкости. Передаточная функция, соответствующая ФРТ (4.11), определяется выражением

.                  (4.13)

Очевидно, что точки, для которых выполняется условие (4.2), образуют круг радиусом

 .                                               (4.14)

Следовательно, чем больше , тем меньше расфокусировка наблюдаемого изображения. ФРТ для земной атмосферы и соответствующая ей передаточная функция при  приведены на рис.4.10 и 4.11. Радиус кружка размытия примерно равен 15. Размеры пятна ФРТ на рис. 4.10 визуально кажутся меньше чем размеры пятна для тонкой линзы (рис. 4.8), т.к. гауссовский импульс является быстро убывающей функцией.

Рис.4.10. ФРТ атмосферы Земли	Рис.4.11. Изображение модуля частотной характеристики атмосферы Земли

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

Рис.4.12. Дефокусированное изображение “Сатурн”

искаженный вариант изображения «Сатурн» (рис.4.6). Свертка исходного изображения производилась с гауссовским импульсом при . Искаженное изображение содержит   элементов.

Таким образом,

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

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


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

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

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

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

4.2. Алгебраические методы восстановления изображений



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



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

,                                        (4.15)

где искаженное изображение

.                                            (4.16)

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

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

 .                                                  (4.17)

Если бы  была квадратной матрицей и существовала бы обратная матрица , то, очевидно, что решение системы имело бы вид

.                                                 (4.18)

Однако матричное уравнение (4.16) представляет собой недоопределенную

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


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

,                       (4.19)

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

,

где  - обобщенная обратная матрица. В общем случае норма ошибки не равна нулю.

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

Для  искаженных изображений, наблюдаемых в присутствии шумов, к элементам вектора-столбца  добавляются отсчеты вектора-столбца . Это делает систему уравнений, как правило, неразрешимой. Неразрешимость системы означает, что не существует оценки исходного изображения, при которой она перейдет в тождество. Можно найти лишь приближенное решение  неразрешимой системы, которое определяется из условия минимума нормы ошибки  [4.4, 4.5]

 .                  (4.21)

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

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


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

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

4.3. Методы восстановления изображений на основе

 пространственной фильтрации



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

При циклической свертке модель (4.4) формирования искаженного изображения определяется соотношением



      (4.22)

где  - кадр изображения, одинаковый для всех изображений и ФРТ, входящих в (4.22). Размеры кадра равны периоду повторения изображений и ФРТ.

Применяя к (4.22) ДПФ, получим

        (4.23)

 

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

                              (4.24)

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

         (4.25)

4.3.1. Инверсный фильтр

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

.                                     (4.26)

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

                (4.27)

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

На рис. 4.13. и 4.14 приведены результаты восстановления изображений «Часы» и «Сатурн» инверсным фильтром: а) исходные изображения размером  элементов; б) дефокусированные изображения, полученные в результате свертки с гауссовским импульсом при  с последующим «обрезанием» краев до размеров  элементов; в) изображения, восстановленные инверсным фильтром.


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

На рис. 4.15 и 4.16 приведены сечения типичных частотных характеристик ФРТ  и соответствующих им инверсных фильтров, из которых следует,  что модуль передаточной функции формирующей системы, как правило, стремится к нулю   на   высоких    частотах.   Кроме    того,   нули   в передаточной функции имеются в рабочей полосе частот при расфокусировке камеры (4.10) и смазе (4.6). В этом случае инверсный фильтр является сингулярным, т.к. модуль его передаточной функции становится бесконечно большим на пространственных частотах, соответствующих нулевым значениям модуля передаточной функции  искажающей системы. Причем наличие даже относительно слабого шума приводит к появлению интенсивных шумовых составляющих на выходе инверсного фильтра, полностью разрушающих изображение. Этот факт иллюстрируется рис.4.17.  К дефокусированному изображению «Сатурн» (рис. 4.14.б) был добавлен аддитивный дельта-коррелированный шум. Из восстановленного изображения видно, что даже при пренебрежимо малом уровне шума (отношение сигнал/шум ) метод инверсной фильтрации дает очень плохой результат.

а)	б)
в)
Рис.4.13. Результаты восстанвления изображения “Часы”

а)	б)
в)
Рис.4.14. Результаты восстанвления изображения “Сатурн”

	Рис.4.15. Частотные характеристики искажающей системы с цилиндрической ФРТ  и инверсного фильтра	Рис.4.16. Частотные характеристики искажающей системы с гауссовской ФРТ  и инверсного фильтра
Рис.4.17. Результат восстановления изображения “Сатурн” при

<


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

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

4.3.2. Фильтр Винера

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

 , (4.28)

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

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

                   (4.29)



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

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

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

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

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

Рис.4.18. Частотный характеристики фильтра Винера при цилиндрической и гауссовской ФРТ

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

а)	б)
Рис.4.19. Восстановление дефокусированного изображения “Сатурн” при

а)	б)
Рис.4.20. Восстановление смазанного изображения “Сатурн” при

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


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

а)	б)
Рис.4.21. Восстановление дефокусированного изображения “Часы” при
а)	б)
в)
Рис.4.22. Восстановление смазанного изображения “Часы” при

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

4.3.3. Компенсация краевых эффектов при

восстановлении линейно-искаженных изображений

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

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


Некоторые из них будут рассмотрены в данном подразделе.

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

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

          Хорошие результаты дает функция окна [4.8]

,       (4.30)

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

           На рис 4.23 и 4.24 приведены результаты восстановления изображения «Часы» при горизонтальном смазе, где а - результаты умножения строк искаженного изображения, приведенного на рис.4.22.б, на окно Кайзера и окно (4.30); б - результаты восстановления фильтром Винера. Параметры окон подбирались, исходя из визуального качества восстанавливаемых изображений.

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

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


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

а)	б)
Рис. 4.23. Восстановление с использованием окна Кайзера
а)	б)
Рис. 4.24. Восстановление с использованием окна (4.30)

                 (4.31)

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

 ,                  (4.32)

где   - спектральная плотность окна (4.31). Следует подчеркнуть, что импульсная характеристика фильтра (4.32) не сводится к произведению импульсной характеристики фильтра Винера и  регуляризирующего двумерного треугольного окна (4.31).

На рис. 4.25 приведен результат восстановления изображения «Часы» фильтром (4.32), откуда следует, что фильтр (4.32) практически полностью компенсирует краевые эффекты. Это позволяет отказаться от предварительной обработки. Качество восстановления изображения в центре и на краях почти одинаковое. Параметры фильтра (4.32) полностью определяются исходными данными и не требуют дополнительной подстройки. При использовании быстрого преобразования Фурье для обработки изображений объем вычислений при реализации фильтра (4.32) такой же, как и для фильтра Винера (4.28).

Рис.4.25. Результат восстановления с компенсацией краевых эффектов

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



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

.

Коэффициенты  определяются исходя из требований, предъявляемых к свойствам функции яркости. Например, на границе кадра  функция яркости должна равняться нулю, она должна быть неотрицательной, максимальное значение экстраполирующей функции не должно превышать максимального значения наблюдаемого изображения и т.п. Метод экстраполяции иллюстрируется рис.4.25, где а - экстраполированное изображение; б - результат восстановления. Размер наблюдаемого изображения «Часы» (см. рис. 4.22.б)  равен  элементов, экстраполированного - элементов. В качестве экстраполирующей функции использовался полином первой степени.

а)	б)
Рис. 4.26. Восстановление с применением экстраполяции

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

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


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



Рис.4.27. Восстановление с использованием экстраполяции и

компенсации краевых эффектов

4.4. Итерационные методы восстановления изображений

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

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

.                           (4.33)

Представим передаточную функцию инверсного фильтра   в виде геометрической прогрессии:

.                     (4.34)

Подставляя (4.34) в (4.33), получим

          (4.35)

Соотношение (4.35) позволяет представить процедуру нахождения оценки  в виде последовательных приближений:

                                                                                (4.36)









где каждое последующее приближение вычисляется по предыдущему. Взяв преобразование Фурье от соотношений (4.36), получим итерационную процедуру Ван Циттера [4.11]:



(4.37)



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


В итерационном алгоритме (4.37) нахождение обратного оператора заменяется на многократное вычисление свертки.

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

 .                                           (4.38)

Условие (4.38) выполняется для гауссовской ФРТ.  При цилиндрической ФРТ и равномерном смазе соотношение (4.33) заменяют на эквивалентное соотношение

.

Тогда итерационный алгоритм  (4.37) имеет вид [4.6]



(4.39)



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

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



Рис. 4.28. Частотные характеристики итерационного фильтра на разных шагах



На рис 4.29  приведены результаты восстановления изображения «Часы», где а и б - повторно приведенные исходное (рис.4.22.а) и искаженное в результате смаза (рис 4.22.б) изображения; в - восстановленное изображение итерационным алгоритмом (4.37) ( число итераций ); г - результат восстановления по экстраполированному наблюдаемому изображению 4.26.а ( число итераций ). В качестве критерия остановки итеративного процесса использовался критерий минимума нормированной среднеквадратической ошибки оценивания:

,(4.40)

а)	б )
в)	г)
Рис. 4.29. Восстановление изображения “Часы” итерационным алгоритмом

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

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

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

Итерационный алгоритм, например (4.39), с ограничением имеет вид



(4.41)

,

где  оператор ограничения.

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


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

 

                      (4.42)

Для большинства цифровых изображений диапазон изменения яркости равен . Использование нелинейного алгоритма с ограничением диапазона для восстановления изображения «Часы» (рис. 4.29.б) обеспечивает уменьшение среднеквадратической ошибки до 4%. Особенно эффективен этот алгоритм при восстановлении изображений  с распределением яркости, близким к бинарному. На рис.4.30 приведены результаты восстановления изображения «Текст», где а - исходное изображение размером  эл.;  б - часть исходного изображения, попадающая в кадра  размером   эл.; в - наблюдаемое изображение размером  эл., полученное в результате   смаза исходного изображения (, ) ; г - экстраполированное  изображение  размером      эл.; д   и   е   - изображения,  восстановленные по экстраполированному изображению итерационным линейным алгоритмом () и итерационным алгоритмом с ограничением диапазона яркости (). Ошибка вычислялась по кадру   размером  эл. На рис.4.30 размеры изображений увеличены в полтора раза.

а)	б)
в)	г)
д)	е)
Рис. 4.30.Восстановление изображения “Текст” нелинейным итерационным алгоритмом

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

Цифровая обработка

ВВЕДЕНИЕ
          Многие отрасли техники, имеющие отношение к получению, обработке, хранению и передаче информации, в значительной степени ориентируются в настоящее время на развитие систем, в которых информация имеет характер изображений. Изображение, которое можно рассматривать как двумерный сигнал, является значительно более емким носителем информации, чем обычный одномерный (временной) сигнал. Вместе с тем, решение научных и инженерных задач при работе с визуальными данными требует особых усилий, опирающихся на знание специфических методов, поскольку традиционная идеология одномерных сигналов и систем мало пригодна в этих случаях. В особой мере это проявляется при создании новых типов информационных систем, решающих такие проблемы, которые до сих пор в науке и технике не решались, и которые решаются сейчас благодаря использованию информации визуального характера.
          В связи с этим, в вузовских программах появляются дисциплины, направленные на изучение принципов обработки изображений, причем, приоритетное внимание уделяется цифровым методам, привлекательным своей гибкостью. Отсутствие учебной литературы является сильным препятствием данному изучению, что и побудило авторов к написанию пособия. Следует отметить, что ограниченный объем не позволил охватить многие важные аспекты проблемы цифровой обработки изображений. Авторы пособия, читающие курс цифровой обработки изображений  в НГТУ и НГУ, исходили из своих представлений  о важности тех или иных разделов, а также опирались на многолетний научно-исследовательский и педагогический опыт.
          Работа над пособием распределилась следующим образом. Глава 1 написана совместно И.С. Грузманом и А.А. Спектором, главы 2 и 3 - А.А. Спектором, глава 4 - И.С. Грузманом, глава 5 - В.С. Киричуком и Г.И. Перетягиным, глава 6 - В.С. Киричуком и В.П. Косых. При подготовке пособия авторы использовали известные издания в области цифровых методов обработки изображений, а также имеющиеся у них представления об аналогичных дисциплинах, присутствующих в учебных планах ряда высших учебных заведений страны.



    Работа с информацией: Cистемы - Технологии - Рынок

• Анализ информационных систем
  
• Методы информационных систем
  
• Интернет как информационная система
  
• Искусственный интеллект в информационных системах
  
• Обработка информации информационными системами
  
  
• Информационные системы в офисе
  
• Управление информационными системами
  
• Технологии информационных систем
  
• Теория информационных систем
  
• Почта - информационная система
  
  
• Outlook и информационные системы
  
  
• Информационный рынок
  
• Информационный рынок - IT
  
• Технологии информационного рынка
  
• Безопасность на информационном рынке
  

---

---