9.2. ДИСКРЕТИЗОВАННЫЙ ОПЕРАТОР СУПЕРПОЗИЦИИ

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

,            (9.2.1а)

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

.                    (9.2.1б)

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

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

,                (9.2.2)

где . Дельта-функцию можно ввести в подынтегральное выражение интеграла (9.2.1а); в результате получается соотношение

.             (9.2.3)

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

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

                    (9.2.4)

при  и . Тогда

.                       (9.2.5)

Следует отметить, что усечение импульсного отклика эквивалентно умножению его на весовую функцию («окно») , равную единице внутри квадрата ,  и нулю за его пределами. Согласно теореме о спектре свертки, спектр получаемого после усечения изображения равен свертке спектров изображения  и весовой функции (окна) , причем спектр последней есть двумерная sinc-функция. Такое искажение спектра  приводит к появлению паразитных компонент с высокими пространственными частотами, группирующихся вблизи частот, кратных  (явление Гиббса). Искажения, вносимые при усечении, можно уменьшить, используя более совершенные весовые функции, такие, как окна Бартлетта, Блэкмана, Хемминга или Ханна [3], которые позволяют сгладить нежелательные эффекты, связанные с наличием разрывов прямоугольного окна. Выбор подходящего окна имеет большое значение, поскольку некорректность интегрального преобразования (9.2.1а) может привести к резкому усилению искажений, вызванных усечением импульсного отклика в процессе восстановления изображения.

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

,                (9.2.6)

где , причем  и  обозначают минимальное и максимальное значения индекса .

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

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

,                       (9.2.7)

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

.            (9.2.8)

Тогда

.              (9.2.9)

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

,                     (9.2.10)

где символ  обозначает округление до ближайшего целого числа.

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

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

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

.                    (9.2.11)

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

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

Итак, построена дискретная конечная модель линейного интегрального преобразования, которая связывает отсчеты выходного изображения  со значениями элементов исходного изображения  посредством математической операции дискретного линейного преобразования. Такая операция является аппроксимацией непрерывного преобразования, поскольку импульсный отклик  подвергался усечению, а интеграл вычислялся по приближенным квадратурным формулам. Ясно, что, расширяя область определения импульсного отклика, ошибку усечения можно сделать сколь угодно малой, хотя это и связано с увеличением объема вычислений. Ошибки за счет приближенного вычисления интегралов можно уменьшить, применяя более точные квадратурные формулы, но также ценой усложнения вычислений. Следует, однако, отметить, что дискретный вариант линейного интегрального преобразования является точной аппроксимацией последнего, если все функции пространственных переменных, входящие в формулу (9.2.1), имеют спектры ограниченной ширины, а интервалы при дискретизации выходного изображения и численном представлении исходного изображения выбраны в соответствии с критерием Найквиста [4]. Вопрос о точности и устойчивости дискретного представления рассмотрен в гл. 14.

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

,                   (9.2.12а)

где , а

,               (9.2.12б)

где . Кроме того, определим импульсный отклик так, что

.             (9.2.12в)

На рис. 9.2.3 показаны геометрические соотношения между этими функциями.

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

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

,                      (9.2.13)

где , a

.

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

,                       (9.2.14)

где  - матрица размера , имеющая следующую структуру:

.                       (9.2.15)

Элементы матрицы  определяются равенством

                      (9.2.16)

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

,              (9.2.17)

то операция дискретной суперпозиции (9.2.13) превращается в дискретную свертку:

. (9.2.18)

Интересен частный случай, когда , где  - целое число. При этом

            (9.2.19)

Теперь элементы матрицы смазывания равны

,                        (9.2.20)

где  и . Далее, если , то элементы матрицы смазывания принимают вид

.                      (9.2.21)

Кроме того,

.                                                              (9.2.22)

Следовательно, все строки матрицы  получаются сдвигом первой строки. Тогда оператор  является дискретизованным оператором свертки, а представление суммы (9.2.19) упрощается:

.                (9.2.23)

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

Рис. 9.2.4. Примеры матриц дискретизованных операторов свертки.

а - общий случай, ; б - импульсный отклик гауссовой формы, .

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

,            (9.2.24)

где матрицы  и  размера  имеют следующую структуру:

.                      (9.2.25)

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

.              (9.2.26)

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

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