Предисловие
Доброго времени суток! Сегодня решил поделиться с Вами сокровенным — одним из своих любимых велосипедов.
Начну издалека — довольно долго я работал на одном радиозаводе в Челябинске, и был у нас (вообще и сейчас есть, просто я уже не там) один мега-проект: оптико-электронный модуль для охраны физических объектов. Это такая здоровая штука на поворотной установке, с тремя камерами на все случаи жизни (цветная — дневная, ЧБ светочувствительная — для сумерек, и тепловизор — для ночного наблюдения). Берётся такой модуль, ставится на вышку высотой метров 50 — и можно днём и ночью держать под наблюдением территорию в радиусе 4-5 километров. Подробности писать не стану, не о том пост. Кому интересно — сами найдут.
Разумеется, интересных задачек по обработке изображений было много. Об одной из таких я и хочу рассказать. А именно — как использовать массивно-парралельные вычисления для компенсации дрожания камеры в реальном времени, или почему SURF подходит не всегда. Добро пожаловать под кат.
Суть проблемы
Как я уже упоминал ранее — модули ставятся на довольно высокие вышки. Из-за сильного ветра на такой высоте их постоянно мелко трясёт, причём чем меньше угол обзора камеры — тем более это заметно. А за изображением с камер должен наблюдать оператор, причём довольно долго (смену). Это очень тяжело, может банально стошнить. Именно поэтому одним из основных требований к подобным системам является наличие механизма, позволяющего максимально снизить тряску картинки. Где-то для этого используются гиростабилизированные платформы, но далеко не везде. Чаще всего используется программная стабилизация, причём чаще всего — именно SURF. Но подходит он для этого плохо.
Почему SURF и аналогичные алгоритмы плохо для этого подходят
По порядку:
1 — SURF для этого не предназначен. Единственное что нам нужно для решения этой задачи — найти смещение нового кадра относительно предыдущего, в пикселях. Камера не поворачивается вдоль плоскости изображения, и не передвигается ближе или дальше к изображению. SURF-же больше подходит для отслеживания перемещающихся объектов в видеоряде, но не для нахождения смещения всего изображения целиком.
2 — Быстродействие. Не буду вдаваться в подробности, скажу лишь что ограничения по времени серьёзные — алгоритм должен находить смещение между двумя последними кадрами меньше, чем за 13 миллисекунд. Как можно меньше. От SURF нам такого даже близко добиться не удалось.
3 — Помехи. Анализируется смещение контрольных точек, а не всего кадра целиком. Из-за ошибок сопоставления возникают проблемы с отсеиванием ложных совпадений, которые довольно тяжело решить. На живой картинке это вызывало редкие «рывки» изображения в ту или иную сторону.
Мой велосипед
К тому моменту, когда я взялся за решение этой проблемы, я уже имел некоторый опыт с nVidia CUDA. Переписал пару фильтров обработки видео — яркость, контраст. Прирост скорости по сравнению с обработкой на CPU меня несказанно порадовал, именно с этого времени я и увлёкся массивно-парралельными вычислениями и GPGPU в частности.
Именно массивно-парралельные вычисления я и решил использовать для решения этой задачи. Поначалу искал готовые решения, пробовал. Искал реализации SURF на CUDA, но, насколько помню, так и не нашел ничего подходящего. И решил изобрести своё. Очень кстати мне попалась отличная статья товарища BigObfuscator, после которой я понял, в какую сторону мне копать.
Алгоритм я назвал «Ratel Deshake», и состоит он из 3 этапов:
— подготовка нового кадра
— перевод кадра в сегментно-интегральное представление
— поиск наилучшего совпадения
Все три этапа довольно хорошо подходят для массивно-парралельной обработки. Я реализовал на CUDA, но, насколько я понимаю, всё это можно довольно хорошо реализовать и на ПЛИС, например.
Подготовка нового кадра
На этом этапе нужно убрать из входных данных всё, что нас не интересует. Если точнее — цвет, и, по возможности, помехи. Я для этого использовал перевод изображения в чёрно-белое и выделение контуров с порогом. Если в виде кода, то что-то вроде этого:
struct pixel {
unsigned char R, G, B;
}
pixel Image[width, height];
unsigned int GSImage[width, height];
unsigned int processedImage[width, height];
for (int x = 0; x < width; x++) {
for (int y = 0; y < height; y++) {
GSImage[x, y] = Image[x, y].R + Image[x, y].G + Image[x, y].B;
}
}
for (int x = 1; x < width; x++) {
for (int y = 1; y < height; y++) {
preprocessedImage[x, y] = 0;
if (GSImage[x, y] - GSImage[x - 1, y] > threshold) preprocessedImage[x, y]++;
if (GSImage[x, y] - GSImage[x, y - 1] > threshold) preprocessedImage[x, y]++;
if (GSImage[x, y] - GSImage[x + 1, y] > threshold) preprocessedImage[x, y]++;
if (GSImage[x, y] - GSImage[x, y + 1] > threshold) preprocessedImage[x, y]++;
}
}
Перевод кадра в сегментно-интегральное представление
Пожалуй, самый интересный этап. Для начала необходимо разъяснить некоторые моменты.
Во-первых — что такое интегральное представление изображения. Это можно понять, прочитав вот эту статью, если вы не сделали этого по ссылке выше. Если коротко, то это двухмерный (в нашем случае) массив, в котором элемент с координатами X и Y равняется сумме всех элементов исходного массива с координатами [0..X, 0..Y].
Во-вторых — что я называю сегментно-интегральным представлением. Из-за особенностей моего алгоритма нам необходим двухмерный массив, в котором элемент с координатами X и Y равняется сумме всех элементов исходного массива с координатами [(X-N)..X, (Y-M)..Y], где N и M — произвольные числа, влияющие на соотношение быстродействие/точность алгоритма. Например я использовал значения 15 и 15, т.е. элемент в результирующем массиве представлял из себя сумму из 256 элементов (квадратика 16х16) исходного массива. Да, фактически это — некоего рода размытие изображения.
Но вернёмся к интегральному представлению изображения. Путём недолгих размышлений можно понять, что в один поток его можно оптимально посчитать следующим образом:
unsigned int Image[width, height];
unsigned int processedImage[width, height];
processedImage[0, 0] = Image[0, 0];
for (int x = 1; x < width; x++) {
processedImage[x, 0] = Image[x, 0] + processedImage[x-1, 0];
}
for (int y = 1; y < height; y++) {
processedImage[0, y] = Image[0, y] + processedImage[0, y-1];
}
for (int x = 1; x < width; x++) {
for (int y = 1; y < height; y++) {
processedImage[x, y] = Image[x,y]+processedImage[x-1,y]+processedImage[x,y-1]-processedImage[x-1,y-1];
}
}
Но есть одна проблема. Разбить эту часть алгоритма на большое количество потоков не получится. Или всё-таки получится? Разумеется. Достаточно разделить его на два отдельных этапа — по горизонтали и по вертикали. Тогда каждый из этих этапов отлично параллелится на количество потоков, равное количеству пикселей по горизонтали и вертикали соответственно. Примерно вот так:
unsigned int Image[width, height];
unsigned int horizontalImage[width, height];
unsigned int processedImage[width, height];
//по горизонтали
for (int y = 0; y < height; y++) {
horizontalImage[0, y] = Image[0, y];
for (int x = 1; x < width; x++) {
horizontalImage[x, y] = Image[x, y] + horizontalImage[x-1, y];
}
}
//по вертикали
for (int x = 0; x < width; x++) {
processedImage[x, 0] = horizontalImage[x, 0];
for (int y = 1; y < height; y++) {
processedImage[x, y] =horizontalImage[x, y] + processedImage[x, y-1];
}
}
А теперь то-же самое, но для сегментно-интегрального представления c N=M=15:
unsigned int Image[width, height];
unsigned int horizontalImage[width, height];
unsigned int processedImage[width, height];
//по горизонтали
for (int y = 0; y < height; y++) {
horizontalImage[0, y] = Image[0, y];
for (int x = 1; x < 16; x++) {
horizontalImage[x, y] = Image[x, y] + horizontalImage[x-1, y];
}
for (int x = 16; x < width; x++) {
horizontalImage[x, y] = Image[x, y] + horizontalImage[x-1, y] - Image[x-16, y];
}
}
//по вертикали
for (int x = 0; x < width; x++) {
processedImage[x, 0] = horizontalImage[x, 0]
for (int y = 1; y < 16; y++) {
processedImage[x, y] = horizontalImage[x, y] + processedImage[x, y-1];
}
for (int y = 16; y < height; y++) {
processedImage[x, y] = horizontalImage[x, y] + processedImage[x, y-1] - horizontalImage[x, y-16];
}
}
Такое вот небольшое колдунство. На самом деле распараллелить можно ещё больше, если делить изображение на несколько отдельных частей, но это уже отдельные пляски с бубном. В любом случае, после всего этого мы имеем сегментно-интегральное представление изображения. На этом этапе нужно сделать ещё одну важную вещь — кроме выбора N и M надо выбрать по какому количеству «блоков» будет идти сравнение.
Этот момент тоже стоит подробно пояснить. Дело в том что при поиске оптимального совпадения ищется именно совпадение предыдущего кадра с новым, а не наоборот. После прохода итерации алгоритма — для следующей итерации остаются данные о предыдущем кадре, но не все, а только значения определённого числа блоков (в моём случае, имеющих размер 16х16). Это ещё один фактор, влияющий на соотношение быстродействие/точность. Меньше блоков — больше быстродействие, больше блоков — больше точность.
Конкретно в моём случае размер изображения был 704x576 пикселей. Из сегментно-интегрального представления предыдущего кадра я оставлял значения только 1140 сегментов (38x30 блоков из центра изображения, взаимно не накладывающихся, 4560 байт). Понятней будет в виде картинки:
Здесь у нас гипотетическая картинка размером 256x256 пикселей. Для удобства она поделена на отдельные блоки по 16x16 пикселей (хотя в сегментно-интегральном представлении блоков значительно больше, и они накладываются друг на друга). В центре выделен массив блоков, 10 на 10 штук. Именно значение этих блоков необходимо запомнить для использования в следующей итерации алгоритма. В коде это будет выглядеть примерно так:
unsigned int IntegralImage[256, 256];
unsigned int NextIteration[10, 10];
for (int x = 0; x < 10; x++) {
for (int y = 0; y < 10; y++) {
NextIteration[x, y] = IntegralImage[(x+4)*16-1, (y+4)*16-1];
}
}
Зачем нам нужен этот массив — станет понятно на следующем этапе.
Поиск наилучшего совпадения
Почти всё, финишная прямая. Если через алгоритм был пропущен только один кадр — этот этап надо пропустить (нечего пока стабилизировать). Если-же это не первый кадр — от предыдущей итерации у Вас должен остаться массив NextIteration. Так как он должен остаться от предыдущей итерации — на этом этапе он будет называться PrevIteration, не перепутайте.
Итак, у нас есть сегментно-интегральное представление нового кадра и массив контрольных блоков от предыдущего. На этом этапе нужно всего-лишь найти положение, при котором наложение блоков предыдущего кадра на сегментно-интегральное представление нового кадра даст минимум различий.
Важное примечание — нет особого смысла проверять все возможные положения, лучше константой выставить максимальное ожидаемое отклонение в пикселях. Например — плюс-минус 30 пикселей по X и Y.
Первым делом нам нужен массив, отображающий разницу при том или ином смещении. Пример в коде для картинки с квадратиками:
unsigned int IntegralImage[256, 256];
unsigned int PrevIteration[10, 10];
unsigned int Differences[61, 61];
for (int X = 0; X < 61; X++) {
for (int Y = 0; Y < 61; Y++) {
Differences[X, Y] = 0;
}
}
//При максимальном ожидаемом отклонении в +-30 пикселей эта часть параллелится на 3721 поток
for (int X = 0; X < 61; X++) {
for (int Y = 0; Y < 61; Y++) {
for (int x = 0; x < 10; x++) {
for (int y = 0; y < 10; y++) {
Differences[X, Y] += abs(PrevIteration[x,y]-IntegralImage[(x+4)*16-31+X,(y+4)*16-31+Y];
}
}
}
}
Вот и всё. Осталось только найти минимальный элемент в массиве Differences. Если это будет [0, 30] — значит новый кадр сместился на 30 пикселей вправо. Если [60, 30] — то на 30 пикселей влево.
Итог
В этой статье я не рассказал ничего иновационного или очень уж неочевидного. Но это значительный кусок моего опыта по этой теме, и я искренне надеюсь что он кому-нибудь да пригодится. Статья вышла длинной и запутанной, так что наверняка где-нибудь есть явные ошибки, или я недостаточно внятно что-то изложил — буду рад отзывам, конструктивной критике и любой другой обратной связи. Огромное спасибо тем, кто осилил до конца.
Что до моей реализации — получилась более-менее рабочая тестовая софтина. На видеокарте уровня GeForce GTX560 позволяла выполнить возложенную задачу — обработать 75 кадров в секунду размером 704x576. И ещё оставалось прилично запаса по времени. Разумеется в статье я изложил только общие принципы — многое здесь можно оптимизировать. К сожалению было принято решение не использовать CUDA в проекте, поэтому я уже больше года его не трогал.
Ещё раз выражаю признательность товарищу BigObfuscator за отличную статью, без которой я бы врятли придумал хоть что-то путнее.
Сырцы
Как я уже писал выше — я уже больше года не занимался реализацией этого алгоритма. Исходник есть, но он настолько неприглядный, что мне стыдно его выкладывать в нынешнем виде. Если достаточному количеству людей будет интересно — могу стряхнуть пыль с файлов, привести всё в порядок и выложить.
Всем приятного кодинга и хороших идей!
Автор: EighthMayer
