Оператор Собеля-Фельдмана или Зачем нам так много фильтров14.01.2022 00:16

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

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

Пикселизация изображения

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

Итак, значение пиксела $I_{mn}$ — это НЕ значение функции I(x,y) в некоей дискретной точке (x_m, y_n) , а среднее значение интенсивности, взятое по площади пиксела:

$I_{mn} = \frac{1}{\Delta}\int\limits_{-\Delta/2}^{\Delta/2} dx \int\limits_{-\Delta/2}^{\Delta/2} dy\ I(x_m + x, y_n + y)\qquad(1)$

где $\Delta$ — это ширина (она же высота) пиксела.

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

Как видим, в формуле (1) есть две составляющих. Одна составляющая — это усреднение, когда вместо значения функции в точке (x,y) мы знаем её среднее значение в квадратной области вокруг этой точки. Вторая составляющая — это дискретизация, когда даже средние значения функции мы знаем не в любой точке, а только на дискретном наборе точек (x_m, y_n) . Задача вычисления градиента изображения состоит в том, чтобы по сеточной функции $I_{mn}$ вычислить производные оригинальной функции I_x(x,y) и I_y(x,y) .

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

Вычисление градиентов: влияние дискретизации

Чтобы посмотреть, какие эффекты вносит дискретизация, мы проигнорируем влияние усреднения. Другими словами мы хотим посмотреть, что получилось бы, если бы яркости пикселов $I_{mn}$ были бы равны значениям функции I(x,y) в точках сетки (x_m, y_n) . Эту, настоящую, интенсивность I(x,y) мы, как уже говорили, считаем хорошей функцией координат. Это значит, мы можем оценить точность численного вычисления производных с помощью разложения выражений в ряд Тейлора.

Рис. 1: Вычисление центральных разностей

Вычисление производной по оси в точке (Рис. 1) с помощью «правой разности» даёт оценку

$I_{x,A} \approx \frac{I_B - I_A}{\Delta} + O(\Delta) \qquad (2)$

где $O(\Delta)$ — это порядок ошибки численного приближения.

Более высокую точность даёт центральная разность

$I_{x,L} \approx \frac{I_B - I_A}{\Delta} + O(\Delta^2) \qquad (3)$

В этом случае ошибка приближения будет порядка $\Delta^2$ , однако есть проблема: у нас были значения интенсивности в точках сетки и , а значение производной по оси мы получили в точке , которая НЕ лежит на нашей сетке.

«Ну и что?» — скажете вы — «У нас просто появилась новая квадратная сетка со значениями производной по оси ».

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

$I_{y,K} \approx \frac{I_C - I_A}{\Delta} + O(\Delta^2) \qquad (4)$

которая не только не лежит на оригинальной сетке, но и не попадает на сетку значений производной по оси . То есть мы не можем вычислить градиент с точностью $O(\Delta^2)$ в одной и той же точке. В итоге мы имеем три разных сетки: значения интенсивности, значения производных по оси и значения производных по оси .

И тут пришёл Лоуренс Робертс (Lawrence Roberts), который решил проблему вычисления градиента в одной и той же точке. Решению дали весьма печальное название — «крест Робертса» (Roberts cross).

Рис. 2: Вычисление оператора Робертса

Лично я плАчу, когда читаю объяснение смысла креста Робертса в википедии. Там пишут про психофизические представления в теории изображений, хотя главное достижение Робертса в том, что он провёл новую систему координат x'y' , наклонённую к оригинальной под углом в $45^\circ$ (Рис. 2), что дало возможность вычислить с помощью центральных разностей обе производные в одной и той же точке

$I_{x',P} \approx \frac{I_D - I_A}{\sqrt{2}\Delta} + O(2\Delta^2) \ ,\qquad I_{y',P} \approx \frac{I_C - I_B}{\sqrt{2} \Delta} + O(2\Delta^2) \qquad (5)$

Длины отрезков между парами точек выросли: $\sqrt{2}\Delta$ вместо $\Delta$ в формулах (3) и (4), а значит и точность меньше, именно поэтому двойка в выражении $O(2\Delta^2)$ . Тем не менее по-прежнему второй порядок точности. Вычисление градиента по отношению к оригинальным осям координат можно выполнить по формулам пересчёта координат вектора при повороте системы координат на $-45^\circ$ (попутно забудем про общий множитель $1/\sqrt{2}$ , так как нам не очень интересны абсолютные
величины градиентов)

$I_{x,P} \propto I_{x',P} - I_{y',P} \approx \frac{I_B + I_D - (I_A + I_C)}{\sqrt{2}\Delta} + O(2\Delta^2) \qquad (6)$ $I_{y,P} \propto I_{x',P} + I_{y',P} \approx \frac{I_C + I_D - (I_A + I_B)}{\sqrt{2}\Delta} + O(2\Delta^2) \qquad (7)$

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

$I_x = \pmatrix{ -1 & 1 \cr -1 & 1 \cr } \qquad I_y = \pmatrix{ -1 & -1 \cr 1 & 1 \cr } \qquad (8)$

Наши матрицы подразумевают перемножение с матрицей изображения, а не свёртку!

Небольшое отклонение от темы. Раньше я думал, что вирусная идея использовать изображения из Плейбоя была запущена инженерами из Университетя Южной Калифорнии в 1973-м году фотографией Лены Шьёблом (боже… как же она, бедняжка, жила с такой фамилией?). Но нет, в 1961 году, то есть за 12 лет до Лены, Робертс в своей статье использовал картинки из Плейбоя (википедия утверждает, что и в своей диссертации тоже), правда тогда они были без обнажёнки. Многие из нас прошли период мастурбации на девочек из мужских журналов, но как-то никому из моих знакомых не приходило в голову вставлять эти картинки в научные статьи. Поэтому популярность изображения Лены среди образованной публики не перестаёт меня удивлять. И если вдруг кто-то не знает… Лену урезали, её полное изображение можно найти тут.

Но довольно невтемщины, подведём предварительные итоги. Выше мы показали, что прямолинейный подход (3) с центральными разностями даёт градиенты $I_{x}$ в одних точках, а градиенты $I_{y}$ — в других, в то время как подход Робертса хоть и даёт чуть худшую точность, но зато позволяет вычислять градиенты $I_{x}$ и $I_{y}$ в совпадающих точках сетки, а именно в углах пикселов.

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

Рис. 3: Центральные разности для фильтра размером 3 на 3

В соответствии с Рис. 3 мы можем написать

$I_{x,E} \approx \frac{I_F - I_D}{2\Delta} + O(4\Delta^2) \ , \qquad I_{y,E} \approx \frac{I_H - I_B}{2\Delta} + O(4\Delta^2) \qquad (9)$

Формулы (9) полностью решают поставленную задачу. Точность имеет порядок $O(4\Delta^2)$ , то есть как мы и хотели, квадратична по $\Delta$ , но множитель стал больше, ибо длина между вычитаемыми точками увеличилась.

Вычисление градиентов: усреднение

Пришло время вспомнить, что значение пиксела $I_{mn}$ — это не только дискретизация пространства, но и усреднение значений интенсивности в пределах пиксела (1).

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

Чтобы увидеть, какие проблемы оно несёт, вернёмся к подходу Робертса. Для проверки его формул мы создадим изображение, состоящее из двух областей с интенсивностями и , разделённых прямой линией (Рис. 4). Выше линии границы (по оси )
функция I(x,y) будет равна , а ниже линии границы — . Чтобы формулы давали наилучшую точность, мы проведём линию границы через точку (на рисунке она образует положительный угол $\alpha$ с осью ).

Теперь можно вычислить интенсивности пикселов $I_{mn}$ по тому, какая часть площади каждого пиксела будет освещена, и получить значение градиента Робертса по формулам (6) и (7). Если обозначить угол, который вектор градиента образует с осью как $\gamma$ , то понятно, что в идеале мы должны получить $\gamma \equiv \alpha$ . Давайте проверим.

При угле наклона световой границы $\alpha < 45^\circ$ :

Пиксел будет освещён полностью, $I_{C} = 1$
Пиксел будет полностью тёмным, $I_{B} = 0$ .
Пиксел будет освещена частично, $I_D = 1-\tan(\alpha)/2$
Пиксел будет освещен частично, $I_A = \tan(\alpha)/2$

Тогда по формулам (6) и (7) можно найти

$\tan(\gamma) = \frac{-I_{x,P}}{I_{y,P}} = \frac{I_A + I_C - (I_B + I_D)}{I_C + I_D -(I_A + I_B)}= \frac{\tan(\alpha)}{2-\tan(\alpha)}$

Откуда

$\gamma(\alpha) = \arctan \left( \frac{\tan(\alpha)}{2-\tan(\alpha)} \right) \qquad (10)$

Как видим, значение угла $\gamma$ (10), вычисленное по формулам Робертса, не даёт нам $\gamma\equiv\alpha$ . График функции (10) изображён на Рис. 5.

Рис. 5: График угла наклона вектора градиента в зависимости от угла наклона границы света для оператора Робертса.

Как видно из рисунка, формула Робертса даёт точные значения наклона вектора градиента только для углов $0^\circ$ , $45^\circ$ (а также $90^\circ$ , $135^\circ$ и так далее), для остальных углов ошибка вычисления доходит примерно до $8^\circ$ . Для описания данной проблемы используют термин «неизотропность оператора», подразумевая, что при использовании изотропного оператора вращение изображения должно приводить к такому же повороту вектора градиента. Именно эту проблему и пытался решить Ирвин Собель в своём докладе об операторе, получившем позже его имя (а не проблему усреднения шумов, как иногда пишут).

Объяснения логики вывода, которые опубликовал Собель через 46 лет после своего доклада, можно кратко описать так: он хотел сделать изотропный оператор, поэтому усреднил 4 возможных направления, а для проверки точности они с Гарри Фельдманом зрительно (!) проверили, что всё выглядит хорошо. Несмотря на такой фривольный подход, оператор получился весьма приличный, что мы покажем ниже. Но сначала объясним логику вывода.

Вспомним, что для вычисления градиента по формулам Робертса мы использовали 4 точки, через которые можно провести линии, под углами $0^\circ$ , $45^\circ$ , $90^\circ$ , и так далее. То есть как раз под углами, при которых формулы Робертса дали нам правильный угол наклона градиента. Поэтому можно предположить, что если точек, используемых для вычисления будет больше, чтобы появилось больше возможных направлений, то и точность вычисления угла будет больше. То есть для начала нужно взять сетку большего размера.

На Рис. 6 изображена сетка размером 3 на 3.

Рис. 6: Вычисление градиента по формулам Собеля

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

$I_{x,E} \approx \frac{I_F - I_D}{2} + O(4\Delta^2) \ , \qquad I_{y,E} \approx \frac{I_H - I_B}{2} + O(4\Delta^2) \qquad (11)$

Однако можно точно также сказать, что давайте проведём оси координат (x'y') через точки и . Тогда можно вычислить центральные разности в новой координатной системе как

$I_{x',E} \approx \frac{I_I - I_A}{2\sqrt{2}} + O(8\Delta^2) \ , \qquad I_{y',E} \approx \frac{I_G - I_C}{2\sqrt{2}} + O(8\Delta^2) \qquad (12)$

Формулы (11) и (12) — это два различных приближённых вычисления одной и той же величины. Чтобы сделать вычисления более изотропными, мы просто усредним результаты (11) и (12). Единственное, что усреднение нужно сделать правильно, так как нам нужно найти среднее от двух векторов, заданных в разных системах координат. Поэтому сначала нужно пересчитать вектор~(12) в системе (xy) с помощью формул
поворота координатной системы на $-45^\circ$ , а потом уже сложить их. Подробный вывод приводить не буду, чтобы не раздувать текст, а результирующие матрицы вы наверняка видели в википедии (в них снова избавились от общих множителей)

$I_x = \pmatrix{ -1 & 0 & 1 \cr -2 & 0 & 2 \cr -1 & 0 & 1 \cr } \qquad I_y = \pmatrix{ -1 & -2 & -1 \cr 0 & 0 & 0 \cr 1 & 2 & 1 \cr } \qquad (13)$

Снова заметим, что форма матриц у нас отличается от википедии из-за того, что предполагается умножение на матрицу изображения, а не свёртка.

Интересно, что формулы (13) можно получить, используя логику вывода отличную от логики Собеля. Для этого рассмотрим правую часть Рис. 6. Мы можем вычислить градиент в угловой точке по формулам Робертса (8), используя интенсивности
в точках , , и . Аналогично по формулам Робертса можно вычислить градиенты в угловых точках , , . После этого градиент в точке можно вычислить как среднее от градиентов в 4-х угловых точках , , , и , окружающих её. Результат будет идентичен формуле (13).

Чтобы проверить изотропность оператора Собеля-Фельдмана, мы повторим схему, описанную для оператора Робертса (смотри Рис. 4 и описание к нему). А именно, мы будем проводить линии границы света через точку , вычислять освещённости пикселов, далее по значениям освещённости пикселов вычислять угол наклона градиента, который даст оператор Собеля-Фельдмана, и сравнивать его с углом наклона линии границы света.

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

Рис. 7: Точность вычисления угла наклона градиента для различных фильтров.

Для оператора Собеля-Фельдмана можно было бы написать аналитические формулы, аналогичные формуле (10), но я немного облегчил себе работу и вычислял освещённости пикселов, генерируя случайные наборы точек, поэтому графики получились немного шумными.

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

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

$I_x = \pmatrix{ -3 & 0 & 3 \cr -10 & 0 & 10 \cr -3 & 0 & 3 \cr } \qquad I_y = \pmatrix{ -3 & -10 & -3 \cr 0 & 0 & 0 \cr 3 & 10 & 3 \cr } \qquad (14)$

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

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

$I_x = \pmatrix{ -1 & 0 & 1 \cr -1 & 0 & 1 \cr -1 & 0 & 1 \cr } \qquad I_y = \pmatrix{ -1 & -1 & -1 \cr 0 & 0 & 0 \cr 1 & 1 & 1 \cr } \qquad (15)$

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

Наша модель вычисления ошибки (Рис. 7) показывает, что этот подход почти так же плох, как оператор Робертса.

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