Расширение на другое количество измерений

Оператор Собеля.

Введение

Паровая машина в цвете.

Применение оператора Собеля к изображению

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

2. Упрощённое описание

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

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

Строго говоря..

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

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

Расширение на другое количество измерений

Оператор Собеля состоит из двух отдельный операций ^[1]:

Сглаживание треугольным фильтром в перпендикулярном к производной направлении: h( 1) = 1,h(0) = 2,h(1) = 1
Нахождение простого центрального изменения в направлении производной: h'( 1) = 1,h'(0) = 0,h'(1) = 1

Фильтры Собеля для производных изображения в разных пространствах для :

1D: h_x'(x) = h'(x);

2D: h_x'(x,y) = h'(x)h(y)

3D: h_x'(x,y,z) = h'(x)h(y)h(z)

4D: h_x'(x,y,z,t) = h'(x)h(y)h(z)h(t)

Вот пример трёхмерного ядра Собеля для оси z:

Технические детали

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

и две производные, G_x и G_y, теперь можно вычислить как

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

Применение свертки K к группе пикселей P можно представить псевдокодом:

N(x, y) = Сумма { K(i, j).P(x-i, y-j)}, для i, j от 1 до 1.

N(x, y) представляет собой результат применения матрицы свёртки K к P.

Программная реализация оператора Собела может эффективно использовать SIMD-расширения системы команд современных процессоров (т. н. векторизация кода), при этом выигрыш в скорости вычисления оператора может составлять до 5 раз по сравнению с высокоуровневой реализацией ^[2]. Ручное кодирование на языке ассемблера позволяет обогнать по скорости такие компиляторы как Microsoft Visual C++ и Intel C++ Compiler. Вычисление оператора Собела элементарно распараллеливается на произвольное число потоков (в пределе каждую точку результирующего изображения можно вычислять независимо от соседних). Например, при наличии двух процессоров (ядер) верхний полукадр изображения может быть обработан одним из них, а нижний — другим.

Примеры

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

Полутоновое изображение кирпичной стены и стойки для велосипеда	Нормализованный Собелев градиент изображения кирпичной стены и стойки велосипеда
Нормализованный Собелев градиент изображения по x	Нормализованный Собелев градиент изображения по y

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

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

Изменим шейдер следующим образом:

float4 DiffuseColor = float4(0,0,1,1);

float kd = 1;

// поставим для kd значение побольше, чем было раньше, чтобы получше осветить
// модель

float4 PixelShaderFunction(VertexShaderOutput input) : COLOR0

{

// TODO: add your pixel shader code here.

float3 worldPosition = input.WorldPosition;

float3 worldNormal = normalize(input.Normal);

float4 Ambient = ka * AmbientColor;

float3 lightDirection = normalize(LightPosition – worldPosition);

float Diffuse = kd * max(0, dot(worldNormal, lightDirection));

float2 pos= float2(0,0);

pos.x = Diffuse;

float4 dColor = tex2D(LigthMaskSampler, pos) * DiffuseColor;

return Color + Ambient + dColor;

}

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

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

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

Упрощённое описание

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

Формализация

Строго говоря, оператор использует ядра 3×3, с которыми свёртывают исходное изображение для вычисления приближенных значений производных по горизонтали и по вертикали. Пусть A исходное изображение, а G_x и G_y — два изображения, где каждая точка содержит приближенные производные по x и по y. Они вычисляются следующим образом:

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

Координата x здесь возрастает «направо», а y — «вниз». В каждой точке изображения приближенное значение величины градиента можно вычислить, используя полученные приближенные значения производных:

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

где, к примеру, угол равен нулю для вертикальной границы, у которой тёмная сторона слева.

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

Так выглядит кусок шейдера, отвечающий за получение границ:

float color00 = dot(tex2D(Sampler, pos-deltaX-deltaY), luminance);

float color10 = dot(tex2D(Sampler, pos-deltaY), luminance);

float color20 = dot(tex2D(Sampler, pos+deltaX-deltaY), luminance);

float color01 = dot(tex2D(Sampler, pos-deltaX), luminance);

float color21 = dot(tex2D(Sampler, pos+deltaX), luminance);

float color02 = dot(tex2D(Sampler, pos-deltaX+deltaY), luminance);

float color12 = dot(tex2D(Sampler, pos+deltaY), luminance);

float color22 = dot(tex2D(Sampler, pos+deltaX+deltaY), luminance);

float gx = color00*1 + color10*2 + color20*1 +

color01*0 + 0 + color21*0 +

color02*-1 + color12*-2 + color22*-1;

float gy = color00*1 + color10*0 + color20*-1 +

color01*2 + 0 + color21*-2 +

color02*1 + color12*0 + color22*-1;

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

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

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

3. Расчёт локального оптического потока.

Вычисление оптического потока методом Лукаса-Канаде.

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

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

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

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

Метод Лукаса-Канаде

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

На математическом языке это допущение можно записать так: . Где I — это функция яркости пикселей от положения на кадре и времени. Другими словами x и y — это координаты пикселя в плоскости кадра, и — это смещение, а t — это номер кадра в последовательности. Условимся, что между двумя соседними кадрами проходит единичный отрезок времени.

Одномерный случай

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

Обратите внимание, что и при желании можно записать и в общем виде: ; где y и t зафиксированы и равны нулю.

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

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

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

Мы почти у цели. Смещение d — это наша искомая величина, поэтому надо что-то сделать с . Как мы условились ранее, , поэтому просто перепишем:

То есть:

Двумерный случай

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

Где — вектор смещения.
В соответствии со сделанным допущением . Обратите внимание, что это выражение эквивалентно . Это то, что нам нужно. Перепишем:

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

Перепишем ещё раз, раскрыв градиент:

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

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

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

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

Перепишем в более компактной форме и приравняем к нулю:

Перепишем эти два уравнения в матричной форме:

Где

Если матрица М обратима (имеет ранг 2), можем вычислить и , которые минимизируют ошибку E:

Вот собственно и все. Мы знаем приблизительное смещение пикселей между двумя соседними кадрами.

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

Недостатки метода

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

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

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

Ещё одна проблема связана с тем, что некоторые текстуры в изображении дают вырожденную матрицу М, для которой не может быть найдена обратная матрица. Соответственно, для таких текстур мы не сможем определить смещение. То есть движение вроде есть, но непонятно в какую сторону. В общем-то от этой проблемы страдает не только рассмотренный метод. Даже глаз человека воспринимает такое движение не однозначно (Barber pole).

Заключение

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

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

Источники и ссылки:

1. Optical Flow Estimation. David J. Fleet, Yair Weiss. — Обстоятельная статья, которая содержит детальное описание не только метода Лукаса-Канаде, но и других дифференциальных методов.

2. Image and Video Compression for Multimedia Engineering: Fundamentals, Algorithms, and Standards. Yun Q. Shi, Huifang Sun — Учебник по сжатию видео и изображений, но содержит также неплохой обзор истории вопроса.

3. Image Processing On Line — Большое количество актуальных алгоритмов обработки изображений. Что самое приятное, алгоритмы снабжены онлайн демо.

4. www.wikipedia – Свободная инциклопедия.

5. Ватутин Э.И., Мирошниченко С.Ю., Титов В.С. – Программная оптимизация Собеля с использованием SIMD – расширений процессов.