Массивно-параллельная стабилизация изображения

    image
    Предисловие

    Доброго времени суток! Сегодня решил поделиться с Вами сокровенным — одним из своих любимых велосипедов.

    Начну издалека — довольно долго я работал на одном радиозаводе в Челябинске, и был у нас (вообще и сейчас есть, просто я уже не там) один мега-проект: оптико-электронный модуль для охраны физических объектов. Это такая здоровая штука на поворотной установке, с тремя камерами на все случаи жизни (цветная — дневная, ЧБ светочувствительная — для сумерек, и тепловизор — для ночного наблюдения). Берётся такой модуль, ставится на вышку высотой метров 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 байт). Понятней будет в виде картинки:
    image
    Здесь у нас гипотетическая картинка размером 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 за отличную статью, без которой я бы врятли придумал хоть что-то путнее.

    Сырцы

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

    Всем приятного кодинга и хороших идей!
    AdBlock похитил этот баннер, но баннеры не зубы — отрастут

    Подробнее
    Реклама

    Комментарии 19

    • НЛО прилетело и опубликовало эту надпись здесь
        0
        Ну, для VirtualDub и подобного софта этот алгоритм мало пригодится — его целью была риалтайм работа. То есть быстродействие в ущерб точности, но чтобы точность при этом сохранялась на приемлемом уровне. Уже записанное видео можно обработать гораздо лучше более серьёзными алгоритмами.

        Выложу, но за сроки не ручаюсь (возможно на следующей неделе, сейчас у меня есть другая интересная идея, для следующей статьи). Но мне всё равно совесть не позволит в очень уж сыром виде это выложить, надо подштукатурить. =)
        +5
        а про wavelet image stabilization вы что-нибудь читали?
          0
          Возможно, я просто уже не помню. Почитаю, спасибо. =)
          0
          Было бы неплохо привести примеры видео без стабилизации, со стабилизацией с использованием SURF (пусть и не в режиме реального времени) и с Вашей стабилизацией.
            +1
            Ну, сравнительного примера делать не буду, т.к. я не профи в применении SURF. Ну и, как я уже говорил, SURF немного не для этого. А вот видео-демонстрацию — пожалуйста. Выложу вместе с исходниками.
            0
            Да, для стационарных камер проблем меньше, там можно сравнивать только маленький кусочек изображения с высококонтрасными границами.
            Но у меня в обычной бытовой видеокамере 10-летней давности на ленте (Digit8) цифровая стабилизация работала очень хорошо. Очень плавно видео снималось с рук даже когда я поворачивал камеру вслед за движением. Т.е. как-то там более дешевым(менее числодробительным) способом делали. Может просто 2-4 области всего использовалось)
              +6
              Свое время писал стабилизатор правда для стационарных камер. Нахождения смещения основывалось на простой корреляции методом последовательного приближения на разных масштабах. После нахождения смещения, его значение уточнялось с субпиксельной точностью. Оптимизировал нахождение корреляции и выполнение смещения под SSE2. Общее время работы алгоритма было менее 1 мс для мегапиксельного изображения.
                +2
                Было-бы очень интересно почитать подробности по Вашему алгоритму. Очень уж заинтриговали производительностью. =)
                  0
                  Тоже очень интересно, жду статью)
                  0
                  А подробности? Что это за «простая корреляция»? И каков самый большой масштаб в пикселях?
                    0
                    И мне бы тоже почитать…
                    +1
                    Есть ещё метод стабилизации, основанный на фазовой корреляции двумерных преобразований Фурье кадров (или блоков в кадре). Работает достаточно быстро, есть есть быстрая реализация двумерного БПФ (сорри за каламбур...). Не пробовали?
                    Подробнее по теории вот тут.
                      0
                      Пробовали этот метод в схожих обстоятельствах: камера на вышке сотового оператора в лесу, высота — 60 метров, башню качает ветром. Фазовая корреляция дает довольно плохие результаты. Причин две:
                      1) В условиях лесного массива в ветреную погоду довольно много «локальных», внутрикадровых движений — ветром качаются не только вышки, но и деревья. Из-за этого локализация корреляционного пятна нарушается — оно скачет туда сюда. Результата два: либо невозможность вывода о локации корреляционного пятна, либо неправильный вывод и, как следствие, дестабилизация вместо стабилизации.
                      2) В сумерках (даже ранних) шумы матрицы становятся настолько большими, что фильтрация либо не может с ней справиться, либо вместе шумами «вымывает» из картинки значительную часть полезной информации (лес на не слишком больших зумах — довольно мелкогранулярная структура, ее всякими сглаживателями можно отутюжить до плоского состояния и останется вместо леса одна бурая масса, на которой не за что зацепиться). Результат — пятно корреляции настолько большое, что можно даже не пытаться понять куда уехала картинка.

                      Хотя в охранке, возможно, таких проблем нет, картинка может быть довольно гладкой, с крупными деталями.
                        0
                        Не обязательно считать именно фазовые корреляции, можно и амплитудные
                      0
                      Как показывает опыт, в таких системах проще установить ВОГи и стабилизировать картинку по ним. Единственно, что при узком поле зрения уход гироскопа довольно заметен. Но таким образом мы отвязываемся от картинки.
                        0
                        Здесь может быть нюанс инсталляционного характера. Я с железками-стабилизаторами не работал, но классически при установке на вышки есть три проблемы:
                        1) Масса того, что устанавливается. Камера (по-крайней мере которые мы ставим) имеет массу 5 кг. Лишний килограмм стоит денег при аренде установочной площадки. К тому же непонятно как вкручивать стабилизатор между кронштейном камеры и самой камерой.
                        2) Питание. Протащить лишний провод на высоту 60м — проблема.
                        3) Защита от окружающей среды. Молнии, снег, дождь, мороз: железяка должна работать. Сколько будет стоить стабилизатор, который защищен таким образом — тайна.
                        Так что когда кто-то говорит «проще» всегда задумываешься какой именно аспект имеется ввиду. «Проще» — это «дешевле»? Вряд ли. Даже если сам гироскоп стоит копейки, кожух + система подогрева + спец крепления + бог знает еще что может стоить больших денег, если вообще найдется на рынке. «Проще» — это «лучше результат»? Не знаю, возможно, но тут эксперимент — критерий истины. «Проще» — «меньше усилий на внедрение»? Если готовое решение есть на рынке — да, но интуиция говорит, что наружнее видеонаблюдение — довольно экзотическая сфера применения электромеханических стабилизаторов такого рода, так что готовое решение вряд ли есть.
                          0
                          1) масса волоконно оптического гироскопа 30грамм, для стабилизации картинки нужно 2 шт, устанавливается на саму камеру.
                          2) провод тащить не надо, достаточно одной платы с ПЛИС которая будет обрабатывать данные с гироскопа и электронным контуром стабилизировать картинку.
                          3) Если честно думал что объем, в котором у вас установлены камеры уже защищён от воздействий внешней среды.

                          Да, ВОГи удовольствие не дешёвое, но есть более дешёвые мемс гироскопы, которые сейчас практически не уступают по точностям ещё и по размерам/весу значительно меньше. Про электромеханические стабилизаторы я ничего не говорил. Про лучше результат — тут есть один аспект. В начале статьи говорится про поворотную установку, что подразумевает управление ею с пульта. Если мы говорим о чисто программной стабилизации — она не будет работать во время поворота камер с пульта или допустим изменении трансфокации, как следствие будет значительно сложнее навестись на интересующий объект на узком поле. Да и вообще корреляторы сдвига изображений — очень нестабильные алгоритмы, сложные для реализации и требующие адоптации под конкретные условия. Допустим было бы интересно посмотреть как ваша система будет работать в дождь или снегопад, или допустим когда через весь кадр движется дымка с большой скоростью. Единственное преимущество корреляторов — возможность компенсации атмосферных явлений.
                            +1
                            Лучше всего комплексное решение (из своего опыта) — сначала механика, потом программная обработка

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

                      Самое читаемое