Алгоритмы совмещения изображений

Written on .

ОГЛАВЛЕНИЕ

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

•    Скачать двоичные файлы - 1.38 Мб
•    Скачать исходники - 81.68 Кб

Введение

Совмещение изображений – процесс сопоставления одного изображения под названием шаблон (обозначаемый как T) с другим изображением I (смотрите рисунок выше). У совмещения изображений много применений, типа отслеживания объектов на видео, анализа движения и многих других задач машинного зрения.

В 1981 Брюс Д. Лукас и Такео Канаде предложили новый метод, использовавший данные о градиенте интенсивности изображения для поиска лучшего совпадения между шаблоном T и другим изображением I. Предложенный алгоритм широко применяется в области машинного зрения в течение последних 20 лет и имеет много изменений и расширений. В их число входит алгоритм, предложенный Симоном Бейкером, Фрэнком Деллаертом и Иэном Мэтьюсом. Их алгоритм гораздо эффективней вычислительно, чем исходный алгоритм Лукаса-Канаде.

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

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

Зачем нужна эта статья? Она дает хорошее понимание алгоритмов совмещения изображений. Есть много статей с уймой информации о совмещении изображений. Большинство из них предназначено для опытных читателей, а не для начинающих. Среди этих статей есть несколько таких, которые описывают сложные вещи просто. Но в них не хватало примеров кода для испытания. Итак, остальные две задачи этой статьи – относительно простое объяснение принципов работы совмещения изображений и предоставление примеров исходного кода.

Кому будет интересно прочитать статью? Если вы помните математику из университетской программы, знаете C или C++, заинтересованы в том, чтобы узнать, как отслеживать объект на видео, и запаслись терпением, то можете прочитать статью.


Компиляция примера кода

Запустить пример очень легко – скачайте demo_binary.zip и запустите исполняемый файл.
Но для испытания кода вам придется скомпилировать его самостоятельно. Сначала установите Visual Studio .NET 2005. Экспресс-версия Visual C++ тоже подойдет.

Затем скачайте и установите библиотеку OpenCV отсюда. Скачайте OpenCV_1.0.exe и запустите его на вашем компьютере.

Наконец, настройте Visual Studio следующим образом:
•    В главном окне Visual Studio выберите пункт меню Инструменты -> Опции..., затем щелкните по элементу дерева Проекты и решения -> Каталоги VC++.
•    В комбинированном поле 'Показать каталоги для..' выберите 'Библиотечные файлы'. Добавляйте следующую строку к списку ниже:

C:\Program Files\OpenCV\lib

•    В комбинированном поле 'Показать каталоги для..' выберите 'Включаемые файлы'. Добавляйте следующие строки к списку ниже:

C:\Program Files\OpenCV\cv\include
C:\Program Files\OpenCV\cxcore\include
C:\Program Files\OpenCV\otherlibs\highgui

Вот и все! Теперь можно открыть align_demo.sln, скомпилировать и запустить его.

Несколько слов об устройстве примера проекта.

auxfunc.h и auxfunc.cpp содержат вспомогательные функции деформирования изображения. forwadditive.h и forwadditive.cppfiles содержат реализацию алгоритма Лукаса-Канаде. invcomp.h и invcomp.cpp содержат реализацию обратного композиционного алгоритма. И main.cpp содержит главную функцию.

Что делает программа?

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

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

Разрешенный набор деформаций определяется путем определения матрицы деформаций, . Матрица W зависит от вектора параметров, p=(wz, tx, ty). Здесь tx и ty являются компонентами переноса, а wz является компонентом поворота.

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

Логика приложения следующая.
1.    После компиляции и запуска демо выводится консольное окно, где будет отображаться диагностическая информация. Введите компоненты вектора p=(wz, tx, ty).
2.    Затем выведутся два окна, содержащие шаблон T и входное изображение I. Также отобразится белый прямоугольник на изображении I. Это начальное приближение положения бабочки. Просто предполагается, что бабочка расположена где-то вблизи этого участка.
3.    Нажмите любую клавишу, чтобы запустить алгоритм совмещения изображений Лукаса-Канаде. Этот алгоритм пытается найти местоположение бабочки, итерационно уменьшая разницу между шаблоном и изображением I. В консольном окне видно, как он работает. Видно текущее число итераций, увеличения параметров и среднее значение ошибки.
4.    После схождения процесса минимизации ошибки сводка выводится на консоль. Сводка включает в себя название алгоритма, время расчета, итоговую среднюю ошибку  и приближение параметров. Также виден белый прямоугольник, показывающий, как шаблон совмещен с изображением I.
5.    Теперь нажмите любую клавишу, чтобы запустить такой же процесс, но использующий обратный композиционный алгоритм совмещения изображений. Будут видны текущее число итераций и среднее значение ошибки.
6.    Отображается сводка для обратного композиционного алгоритма. Нажмите любую клавишу, чтобы выйти из программы.

Содержимое main.cpp приведено ниже:

// План:
// 
// 1. Попросить пользователя ввести вектор параметров деформации изображения p=(wz, tx, ty)
// 2. Назначить прямоугольник шаблона ограничительным прямоугольником бабочки
// 3. Деформировать изображение с помощью матрицы деформации W(p)
// 4. Показать шаблон T и изображение I, ждать нажатия клавиши
// 5. Оценить параметры деформации с помощью метода Лукаса-Канаде, ждать нажатия клавиши
// 6. Оценить параметры деформации с помощью метода Бейкера-Деллаерта-Мэтьюса, ждать нажатия //клавиши
//

int main(int argc, char* argv[])
{
 // Попросить пользователя ввести вектор параметров деформации изображения.
 // p = (wz, tx, ty)
 
 float WZ=0, TX=0, TY=0;
 printf("Please enter WZ, TX and TY, separated by space.\n");
 printf("Example: -0.01 5 -3\n");
 printf(">");
 scanf("%f %f %f", &WZ, &TX, &TY);

В OpenCV изображения будут храниться в структуре IplImage, поэтому здесь определяются и вначале устанавливаются в нуль указатели на изображения.

 // Здесь будут храниться изображения.
 IplImage* pColorPhoto = 0; // Фото бабочки на цветке.
 IplImage* pGrayPhoto = 0; // Полутоновая копия фото.
 IplImage* pImgT = 0; // Шаблон T.
 IplImage* pImgI = 0; // Изображение I.

В OpenCV матрицы хранятся в структурах CvMat. Определяется и инициализируется нулем указатель на матрицу деформации W.

 // Здесь будет храниться матрица деформации.
 CvMat* W = 0; // Деформировать W(p)

Теперь создаются два окна для шаблона T и для изображения I. OpenCV поддерживает базовый GUI(графический пользовательский интерфейс). Поддержка GUI реализована в составе OpenCV под названием библиотека highgui. Функция cvLoadImage() позволяет загрузить изображение из файла JPG.
Теперь выделяется память для изображений с помощью функции cvCreateImage(). Первый параметр – размер изображения в пикселях, второй - глубина, а третий – число каналов. Изображения RGB имеют глубину IPL_DEPTH_8U и три канала. IPL_DEPTH_8U использует unsigned char как элемент изображения.IPL_DEPTH_16S использует short int как элемент изображения.

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

 // Создать два простых окна.
 cvNamedWindow("template"); // Здесь будет отображаться T(x).
 cvNamedWindow("image"); // Здесь будет отображаться I(x).

 // Загрузить фото.
 pColorPhoto = cvLoadImage("..\\data\\photo.jpg");

 // Создать другие изображения.
 CvSize photo_size = cvSize(pColorPhoto->width,pColorPhoto->height);
 pGrayPhoto = cvCreateImage(photo_size, IPL_DEPTH_8U, 1);
 pImgT = cvCreateImage(photo_size, IPL_DEPTH_8U, 1);
 pImgI = cvCreateImage(photo_size, IPL_DEPTH_8U, 1);

 // Преобразовать фото в полутоновое, потому что нужны значения интенсивности вместо RGB.
 cvCvtColor(pColorPhoto, pGrayPhoto, CV_RGB2GRAY);

Теперь выделяется память для матрицы деформации W с помощью функции cvCreateMat(). Первые два параметра определяют размер матрицы (3x3), а последний является типом элемента матрицы. CV_32F означает, что используется float в качестве типа элемента матрицы.

 // Создать матрицу деформации, используемую для создания изображения I.
 W = cvCreateMat(3, 3, CV_32F);

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

 // Создать шаблон T
 // задать нужный участок для ограничительного прямоугольника бабочки.
 cvCopy(pGrayPhoto, pImgT);
 CvRect omega = cvRect(110, 100, 200, 150);

Определяется, какое изображение будет использоваться, путем слабой деформации шаблона.

 // Создать I путем деформации фото. 
 init_warp(W, WZ, TX, TY);
 warp_image(pGrayPhoto, pImgI, W);
 
 // Нарисовать прямоугольник приближения начального положения шаблона.
 cvSetIdentity(W);
 draw_warped_rect(pImgI, omega, W);

 // Показать T и I и ждать нажатия клавиши.
 cvSetImageROI(pImgT, omega);
 cvShowImage("template", pImgT);
 cvShowImage("image", pImgI);
 printf("Press any key to start Lucas-Kanade algorithm.\n");
 cvWaitKey(0);
 cvResetImageROI(pImgT);

 // Напечатать параметры, введенные пользователем.
 printf("==========================================================\n");
 printf("Ground truth: WZ=%f TX=%f TY=%f\n", WZ, TX, TY);
 printf("==========================================================\n");

 init_warp(W, WZ, TX, TY);

Запускается алгоритм Лукаса-Канаде:

 // Восстановить изображение I
 warp_image(pGrayPhoto, pImgI, W);

 /* Лукас-Канаде */ 
 align_image_forwards_additive(pImgT, omega, pImgI);
Теперь ожидается нажатие клавиши, и запускается обратный композиционный алгоритм:

 // Восстановить изображение I
 warp_image(pGrayPhoto, pImgI, W);

 printf("Press any key to start Baker-Dellaert-Matthews algorithm.\n");
 cvWaitKey();

 /* Бейкер-Деллаерт-Мэтьюс*/
 align_image_inverse_compositional(pImgT, omega, pImgI);
Ожидается нажатие клавиши, освобождаются все используемые ресурсы, и завершается работа.

 printf("Press any key to exit the demo.\n");
 cvWaitKey();

 // освободить все используемые ресурсы и завершить работу
 cvDestroyWindow("template");
 cvDestroyWindow("image");
 cvReleaseImage(&pImgT);
 cvReleaseImage(&pImgI);
 cvReleaseMat(&W);
 
 return 0;
}

Реализация поступательного аддитивного алгоритма

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

Предварительно вычислить:
1.    Вычислить градиент  изображения I.

Повторять:
1.    Деформировать I с помощью W(x;p) для вычисления I(W(x,p)).
2.    Вычислить изображение ошибки T(x) – I(W(x;p)).
3.    Деформировать градиент  с помощью W(x;p).
4.    Вычислить якобиан   в (x;p).
5.    Вычислить изображения быстрейшего спуска  .
6.    Вычислить матрицу Гессе  
7.    Вычислить 
8.    Вычислить 
9.    Обновить параметры:  
до тех пор, пока  .

Далее реализуется данный алгоритм (смотрите forwadditive.h и forwadditive.cpp).

Сначала включаются заголовки для используемых функций.

#include <span class="code-keyword"><stdio.h> </span>
#include <span class="code-keyword"><time.h> </span>
#include <span class="code-keyword"><cv.h> // Включить заголовок для части машинного зрения OpenCV.</span>
#include <span class="code-keyword"><highgui.h> // Включить заголовок для части GUI OpenCV. </span>
#include <span class="code-string">"auxfunc.h" // Заголовок для функций деформации. </span>

Определяется функция, реализующая поступательный аддитивный алгоритм (Лукаса-Канаде). Изображение шаблона pImgT, ограничительный прямоугольник шаблона omega и другое изображение pImgI являются параметрами для алгоритма. В OpenCV изображения хранятся как структуры IplImage, и прямоугольник можно определить с помощью структуры CvRect.

// Метод Лукаса-Канаде
// @param[in] pImgT изображение шаблона T
// @param[in] omega прямоугольный участок шаблона
// @param[in] pImgI другое изображение I

void align_image_forwards_additive(IplImage* pImgT, CvRect omega, IplImage* pImgI)
{
 // Несколько констант для процесса итерационного уменьшения.
 const float EPS = 1E-5f; // Пороговое значение для условия завершения.
 const int MAX_ITER = 100; // Максимальное число итераций.

В OpenCV изображения будут храниться в структуре IplImage, поэтому здесь определяются и инициализируются нулями указатели на изображения.

 // Здесь будут храниться градиентные изображения. 
 IplImage* pGradIx = 0; // Градиент I в направлении X.
 IplImage* pGradIy = 0; // Градиент I в направлении Y.

Используются матрицы. В OpenCV матрицы хранятся в структурах CvMat.

 // Здесь будут храниться матрицы.
 CvMat* W = 0; // Текущее значение деформации W(x,p)
 CvMat* H = 0; // Гессиан
 CvMat* iH = 0; // Обратный гессиан
 CvMat* b = 0; // Вектор в правой части системы линейных уравнений.
 CvMat* delta_p = 0; // Значение обновления параметра.
 CvMat* X = 0; // Точка в системе координат T.
 CvMat* Z = 0; // Точка в системе координат I.

Теперь выделяется память для матриц с помощью функции cvCreateMat(). Первые два параметра определяют размер матрицы (3x3 или 3x1), а последний является типом элемента матрицы. CV_32F значит, что в качестве типа элемента матрицы используется float.

 // Создать матрицы.
 W = cvCreateMat(3, 3, CV_32F);
 H = cvCreateMat(3, 3, CV_32F);
 iH = cvCreateMat(3, 3, CV_32F);
 b = cvCreateMat(3, 1, CV_32F);
 delta_p = cvCreateMat(3, 1, CV_32F);
 X = cvCreateMat(3, 1, CV_32F);
 Z = cvCreateMat(3, 1, CV_32F);

Выделяется память для изображений с помощью функции cvCreateImage(). Первый параметр – размер изображения в пикселях, второй - глубина, а третий – число каналов. Полутоновые изображения имеют глубину IPL_DEPTH_8U и один канал. IPL_DEPTH_16S использует short int в качестве элемента изображения.

 // Создать полутоновые изображения.
 CvSize image_size = cvSize(pImgI->width, pImgI->height);
 pGradIx = cvCreateImage(image_size, IPL_DEPTH_16S, 1);
 pGradIy = cvCreateImage(image_size, IPL_DEPTH_16S, 1);

Сохранить текущее время:

 // Получить текущее время. Оно используется позже для получения полного времени расчета.
 clock_t start_time = clock();

Предварительно вычисляется градиент изображения  . Заметьте, что функция cvSobel() выдает расширенный градиент изображения (потому что используется гауссово ядро 3x1 или 1x3), поэтому приходится вызывать cvConvertScale() для нормализации изображения градиента.

 /*
 * Стадия предварительного вычисления.
 */
 
 // Вычислить градиент I.
 cvSobel(pImgI, pGradIx, 1, 0); // Градиент в направлении X
 cvConvertScale(pGradIx, pGradIx, 0.125); // Нормализовать
 cvSobel(pImgI, pGradIy, 0, 1); // Градиент в направлении Y
 cvConvertScale(pGradIy, pGradIy, 0.125); // Нормализовать

Войти в цикл итераций. Повторять до схождения или до достижения максимального числа итераций:

 /*
 * Стадия итераций.
 */

 // Здесь будет храниться приближение параметра.
 float wz_a=0, tx_a=0, ty_a=0;
 
 // Здесь будет храниться значение средней ошибки.
 float mean_error=0;

 // Повторять
 int iter=0; // номер текущей итерации
 while(iter < MAX_ITER)
 {
   iter++; // Увеличить счетчик итераций на единицу

   mean_error = 0; // установить значение средней ошибки в нуль.

   int pixel_count=0; // Число обработанных пикселей

Инициализировать деформацию W(x,p) и установить матрицу Гессе в нуль:

   init_warp(W, wz_a, tx_a, ty_a); // Инициализировать деформацию W(x, p)
   cvSet(H, cvScalar(0)); // Установить гессиан в нуль
   cvSet(b, cvScalar(0)); // Установить матрицу b в нуль

   // (u,v) – координаты пикселя в системе координат T.
   int u, v;
 
   // Обойти все пиксели в шаблоне T.
   int i, j;
   for(i=0; i<omega.width; i++)
   {
     u = i + omega.x;

     for(j=0; j<omega.height; j++)
     {
       v = j + omega.y;

       // Установить вектор X в координаты пикселя (u,v,1)
       SET_VECTOR(X, u, v);

       // Деформировать Z=W*X
       cvGEMM(W, X, 1, 0, 0, Z);

       // (u2,v2) - координаты пикселя в системе координат I.
       float u2, v2;
    
       // Получить координаты деформированного пикселя в системе координат I.
       GET_VECTOR(Z, u2, v2);

       // Получить ближайшие целые координаты пикселя (u2i;v2i).
       int u2i = cvFloor(u2);
       int v2i = cvFloor(v2);

       if(u2i>=0 && u2i<pImgI->width && // проверить, находится ли пиксель внутри I.
          v2i>=0 && v2i<pImgI->height)
       {
         pixel_count++;

         // Вычислить градиент I в W(x,p) с подпиксельной точностью
         // с помощью двухлинейной интерполяции.
         float Ix = interpolate<short>(pGradIx, u2, v2);
         float Iy = interpolate<short>(pGradIy, u2, v2);
 
         // Вычислить элемент изображения быстрейшего спуска.
         float stdesc[3]; // элемент изображения быстрейшего спуска
         stdesc[0] = (float)(-v*Ix+u*Iy);
         stdesc[1] = (float)Ix;
         stdesc[2] = (float)Iy;

     // Вычислить интенсивность преобразованного пикселя с подпиксельной точностью
         // с помощью двухлинейной интерполяции.
     float I2 = interpolate<uchar>(pImgI, u2, v2);
    
     // Вычислить разность изображений D = T(x)-I(W(x,p)).
     float D = CV_IMAGE_ELEM(pImgT, uchar, v, u) - I2;
    
     // Обновить значение средней ошибки.
     mean_error += fabs(D);

         // Добавить член в матрицу b.
         float* pb = &CV_MAT_ELEM(*b, float, 0, 0);
         pb[0] += stdesc[0] * D;
         pb[1] += stdesc[1] * D;
         pb[2] += stdesc[2] * D;
 
        // Добавить член в гессиан.
        int l,m;
        for(l=0;l<3;l++)
        {
          for(m=0;m<3;m++)
          {
            CV_MAT_ELEM(*H, float, l, m) += stdesc[l]*stdesc[m];
          } 
        }
      } 
    }
  }

  // Наконец, вычислить итоговую среднюю ошибку.
  if(pixel_count!=0)
     mean_error /= pixel_count;

Инвертировать матрицу Гессе:

  // Инвертировать гессиан.
  double inv_res = cvInvert(H, iH);
  if(inv_res==0)
  {
    printf("Error: Hessian is singular.\n");
    return;
  }

Обновить параметры:  :

  // Найти увеличение параметров. 
  cvGEMM(iH, b, 1, 0, 0, delta_p);
  float delta_wz = CV_MAT_ELEM(*delta_p, float, 0, 0);
  float delta_tx = CV_MAT_ELEM(*delta_p, float, 1, 0);
  float delta_ty = CV_MAT_ELEM(*delta_p, float, 2, 0);
 
  // Обновить приближение вектора параметров.
  wz_a += delta_wz;
  tx_a += delta_tx;
  ty_a += delta_ty;
 
  // Вывести диагностическую информацию на экран.
  printf("iter=%d dwz=%f dtx=%f dty=%f mean_error=%f\n", 
    iter, delta_wz, delta_tx, delta_ty, mean_error); 
 
  // Проверить условие завершения.
  if(fabs(delta_wz)<EPS && fabs(delta_tx)<EPS && fabs(delta_ty)<EPS) break;
 }

Показать результаты совмещения:

 // Получить текущее время и получить полное время расчета.
 clock_t finish_time = clock();
 double total_time = (double)(finish_time-start_time)/CLOCKS_PER_SEC;

 // Распечатать сводку.
 printf("===============================================\n");
 printf("Algorithm: forward additive.\n");
 printf("Caclulation time: %g sec.\n", total_time);
 printf("Iteration count: %d\n", iter);
 printf("Approximation: wz_a=%f tx_a=%f ty_a=%f\n", wz_a, tx_a, ty_a);
 printf("Epsilon: %f\n", EPS);
 printf("Resulting mean error: %f\n", mean_error);
 printf("===============================================\n");
 
 // Показать результат совмещения изображений.
 init_warp(W, wz_a, tx_a, ty_a);
 draw_warped_rect(pImgI, omega, W);
 cvSetImageROI(pImgT, omega);
 cvShowImage("template",pImgT);
 cvShowImage("image",pImgI);
 cvResetImageROI(pImgT);

Освободить все используемые ресурсы:

 // Освободить используемые ресурсы и завершить работу.
 cvReleaseMat(&W);
 cvReleaseMat(&H);
 cvReleaseMat(&iH);
 cvReleaseMat(&b);
 cvReleaseMat(&delta_p);
 cvReleaseMat(&X);
 cvReleaseMat(&Z);
}

Реализация обратного композиционного алгоритма совмещения изображений

В обратном композиционном алгоритме шаблон T и изображение I меняются ролями.

Предварительно вычислить:

1.    Вычислить градиент  изображения T(x).
2.    Вычислить якобиан   в (x;0).
3.    Вычислить изображения быстрейшего спуска,  .
4.    Вычислить матрицу Гессе,

Повторять:

1.    Деформировать I с помощью W(x;p) для вычисления I(W(x,p)).
2.    Вычислить изображение ошибки I(W(x;p))-T(x).
3.    Деформировать градиент  с помощью W(x;p).
4.    Вычислить  
5.    Вычислить  
6.    Обновить параметры

до тех пор, пока  .

Код этого метода находится в invcomp.h и invcomp.cpp. Ниже показано, что содержит invcomp.cpp.

Включить все требуемые заголовки:

#include <span class="code-keyword"><stdio.h></span>
#include <span class="code-keyword"><time.h></span>
#include <span class="code-keyword"><cv.h>        // Включить заголовок для части машинного зрения OpenCV.</span>
#include <span class="code-keyword"><highgui.h>    // Включить заголовок для части GUI OpenCV.</span>
#include <span class="code-string">"auxfunc.h"    // Заголовок для функций деформации.</span>

Реализация алгоритма Бейкера-Деллаерта-Мэтьюса (обратного композиционного) приведена ниже:

// Обратный композиционный метод Бейкера-Деллаерта-Мэтьюса
// @param[in] pImgT   Изображение шаблона T
// @param[in] omega   Прямоугольный участок шаблона
// @param[in] pImgI   Другое изображение I

void align_image_inverse_compositional(IplImage* pImgT, CvRect omega, IplImage* pImgI)
{
    // Несколько констант для процесса итеративного уменьшения.
    const float EPS = 1E-5f; // Пороговое значение для условия завершения.
    const int MAX_ITER = 100;  // Максимальное число итераций.

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

    // Здесь будут храниться используемые внутри изображения.
    IplImage* pGradTx = 0;       // Градиент I в направлении X.
    IplImage* pGradTy = 0;       // Градиент I в направлении Y.
    IplImage* pStDesc = 0;        // Изображения быстрейшего спуска.

    // Здесь будут храниться матрицы.
    CvMat* W = 0;  // Текущее значение деформации W(x,p)
    CvMat* dW = 0;  // Обновление деформации.
    CvMat* idW = 0; // Обратное обновление деформации.
    CvMat* X = 0;  // Точка в системе координат T.
    CvMat* Z = 0;  // Точка в системе координат I.

    CvMat* H = 0;  // Гессиан.
    CvMat* iH = 0; // Обратный гессиан.
    CvMat* b = 0;  // Вектор в правой части системы линейных уравнений.
    CvMat* delta_p = 0; // Значение обновления параметра.
    
    // Создать матрицы.
    W = cvCreateMat(3, 3, CV_32F);
    dW = cvCreateMat(3, 3, CV_32F);
    idW = cvCreateMat(3, 3, CV_32F);
    X = cvCreateMat(3, 1, CV_32F);
    Z = cvCreateMat(3, 1, CV_32F);

    H = cvCreateMat(3, 3, CV_32F);
    iH = cvCreateMat(3, 3, CV_32F);
    b = cvCreateMat(3, 1, CV_32F);
    delta_p = cvCreateMat(3, 1, CV_32F);

IPL_DEPTH_16S использует короткий тип в качестве элемента изображения градиента. Здесь используется короткий тип, потому что cvSobel выдает не нормализованное изображение, и может произойти переполнение uchar. Изображение градиента нормализуется посредством cvConvertScale().

    // Создать изображения.
    CvSize image_size = cvSize(pImgI->width, pImgI->height);
    pGradTx = cvCreateImage(image_size, IPL_DEPTH_16S, 1);
    pGradTy = cvCreateImage(image_size, IPL_DEPTH_16S, 1);
    pStDesc = cvCreateImage(image_size, IPL_DEPTH_32F, 3);

    // Получить текущее время. Оно используется позже для получения полного времени расчета.
    clock_t start_time = clock();

    /*
     *  Стадия предварительных вычислений.
     */

    // Вычислить градиент T.
    cvSobel(pImgT, pGradTx, 1, 0); // Градиент в направлении X
    cvConvertScale(pGradTx, pGradTx, 0.125); // Нормализовать
    cvSobel(pImgT, pGradTy, 0, 1); // Градиент в направлении Y
    cvConvertScale(pGradTy, pGradTy, 0.125); // Нормализовать
    

    // Вычислить изображения быстрейшего спуска и гессиан

    cvSet(H, cvScalar(0)); // Установить гессиан в ноль

Обойти пиксели изображения шаблона pImgT и собрать матрицу Гессе H и матрицу b:

    int u, v;    // (u,v) – координаты пикселя в системе координат T.
    float u2, v2; // (u2,v2) - координаты пикселя в системе координат I.
    
    // Обойти пиксели в шаблоне T.
    int i, j;
    for(i=0; i<omega.width; i++)
    {
        u = i + omega.x;

        for(j=0; j<omega.height; j++)
        {
            v = j + omega.y;

            // вычислить градиент T.
            short Tx = CV_IMAGE_ELEM(pGradTx, short, v, u);    
            short Ty = CV_IMAGE_ELEM(pGradTy, short, v, u);    
            
            // Вычислить элемент изображения быстрейшего спуска.
            // элемент изображения быстрейшего спуска
            float* stdesc = &CV_IMAGE_ELEM(pStDesc, float, v,u*3);
            stdesc[0] = (float)(-v*Tx+u*Ty);
            stdesc[1] = (float)Tx;
            stdesc[2] = (float)Ty;

            // Добавить член в гессиан.
            int l,m;
            for(l=0;l<3;l++)
            {
                for(m=0;m<3;m++)
                {
                    CV_MAT_ELEM(*H, float, l, m) += stdesc[l]*stdesc[m];
                }
            }
        }    
    }

    // Инвертировать гессиан.
    double inv_res = cvInvert(H, iH);
    if(inv_res==0)
    {
        printf("Error: Hessian is singular.\n");
        return;
    }

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

    /*
     *   Стадия итераций.
     */

    // Установить деформацию в единицу.
    cvSetIdentity(W);

    // Здесь будет храниться текущее значение средней ошибки.
    float mean_error=0;

    // Повторять
    int iter=0; // номер текущей итерации
    while(iter < MAX_ITER)
    {
        iter++; // Увеличить счетчик итераций на единицу

        mean_error = 0; // Установить значение средней ошибки в нуль

        int pixel_count=0; // Число обработанных пикселей
        
        cvSet(b, cvScalar(0)); // Установить матрицу b в ноль
            
        // Обойти пиксели в шаблоне T.
        int i, j;
        for(i=0; i<omega.width; i++)
        {
            int u = i + omega.x;

            for(j=0; j<omega.height; j++)
            {
                int v = j + omega.y;

                // Установить вектор X в координаты пикселя (u,v,1)
                SET_VECTOR(X, u, v);

                // Деформировать Z=W*X
                cvGEMM(W, X, 1, 0, 0, Z);

                // Получить координаты деформированного пикселя в системе координат I.
                GET_VECTOR(Z, u2, v2);

                // Получить ближайшие целые координаты пикселя (u2i;v2i).
                int u2i = cvFloor(u2);
                int v2i = cvFloor(v2);

                if(u2i>=0 && u2i<pImgI->width && // проверить, находится ли пиксель внутри I.
                    v2i>=0 && v2i<pImgI->height)
                {
                    pixel_count++;

                    // Вычислить интенсивность преобразованного пикселя с подпиксельной точностью
                    // с использованием двухлинейной интерполяции.
                    float I2 = interpolate<uchar>(pImgI, u2, v2);

                    // Вычислить разность изображений D = I(W(x,p))-T(x).
                    float D = I2 - CV_IMAGE_ELEM(pImgT, uchar, v, u);

                    // Обновить значение средней ошибки.
                    mean_error += fabs(D);

                    // Добавить член в матрицу b.
                    float* stdesc = &CV_IMAGE_ELEM(pStDesc, float, v, u*3);
                    float* pb = &CV_MAT_ELEM(*b, float, 0, 0);
                    pb[0] += stdesc[0] * D;
                    pb[1] += stdesc[1] * D;
                    pb[2] += stdesc[2] * D;                    
                }    
            }
        }

        // Наконец, вычислить итоговую среднюю ошибку.
        if(pixel_count!=0)
            mean_error /= pixel_count;

        // Найти увеличение параметров. 
        cvGEMM(iH, b, 1, 0, 0, delta_p);
        float delta_wz = CV_MAT_ELEM(*delta_p, float, 0, 0);
        float delta_tx = CV_MAT_ELEM(*delta_p, float, 1, 0);
        float delta_ty = CV_MAT_ELEM(*delta_p, float, 2, 0);

        init_warp(dW, delta_wz, delta_tx, delta_ty);
        // Инвертировать деформацию.
        inv_res = cvInvert(dW, idW);
        if(inv_res==0)
        {
            printf("Error: Warp matrix is singular.\n");
            return;
        }

        cvGEMM(idW, W, 1, 0, 0, dW);
        cvCopy(dW, W);

        // Вывести диагностическую информацию на экран.
        printf("iter=%d mean_error=%f\n", iter, mean_error);

        // Проверить условие завершения.
        if(fabs(delta_wz)<=EPS && fabs(delta_tx)<=EPS && fabs(delta_ty)<=EPS) break;
    }

    // Получить текущее время и получить полное время расчета.
    clock_t finish_time = clock();
    double total_time = (double)(finish_time-start_time)/CLOCKS_PER_SEC;

Распечатать информационную сводку в stdout:

    // Распечатать сводку.
    printf("===============================================\n");
    printf("Algorithm: inverse compositional.\n");
    printf("Caclulation time: %g sec.\n", total_time);
    printf("Iteration count: %d\n", iter);
    printf("Epsilon: %f\n", EPS);
    printf("Resulting mean error: %f\n", mean_error);
    printf("===============================================\n");

    // Показать результат совмещения изображений.
    draw_warped_rect(pImgI, omega, W);
    cvSetImageROI(pImgT, omega);
    cvShowImage("template",pImgT);
    cvShowImage("image",pImgI);
    cvResetImageROI(pImgT);

    // Освободить используемые ресурсы и завершить работу.
    cvReleaseMat(&W);
    cvReleaseMat(&dW);
    cvReleaseMat(&idW);
    cvReleaseMat(&H);
    cvReleaseMat(&iH);
    cvReleaseMat(&b);
    cvReleaseMat(&delta_p);
    cvReleaseMat(&X);
    cvReleaseMat(&Z);
    
}

Двухлинейная интерполяция

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

Используется шаблонная функция interpolate(), определенная в auxfunc.h:

// Интерполирует интенсивность пикселя с подпиксельной точностью.
// О двухлинейной интерполяции в Википедии:
// http://en.wikipedia.org/wiki/Bilinear_interpolation

template <class T>
float interpolate(IplImage* pImage, float x, float y)
{
    // Получить ближайшие целые координаты пикселя (xi;yi).
    int xi = cvFloor(x);
    int yi = cvFloor(y);

    float k1 = x-xi; // Коэффициенты для формулы интерполяции.
    float k2 = y-yi;

    int f1 = xi<pImage->width-1;  // Проверить, что пиксели справа
    int f2 = yi<pImage->height-1; // и ниже существуют.

    T* row1 = &CV_IMAGE_ELEM(pImage, T, yi, xi);
    T* row2 = &CV_IMAGE_ELEM(pImage, T, yi+1, xi);
                
    // Интерполировать интенсивность пикселя.
    float interpolated_value = (1.0f-k1)*(1.0f-k2)*(float)row1[0] +
                (f1 ? ( k1*(1.0f-k2)*(float)row1[1] ):0) +
                (f2 ? ( (1.0f-k1)*k2*(float)row2[0] ):0) +                        
                ((f1 && f2) ? ( k1*k2*(float)row2[1] ):0) ;

    return interpolated_value;

}

Заключение

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

Теперь сравнивается, какой алгоритм работает быстрее на Pentium IV 3 ГГц (конфигурация выпуска). Вводятся следующие компоненты вектора p:

-0.01 5 -3

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

Сводка для поступательного аддитивного алгоритма:

===============================================
Algorithm: forward additive.
Caclulation time: 0.157 sec.
Iteration count: 13
Approximation: wz_a=-0.007911 tx_a=4.893255 ty_a=-3.944573
Epsilon: 0.000010
Resulting mean error: 2.898076
===============================================

Сводка для обратного композиционного алгоритма:

===============================================
Algorithm: inverse compositional.
Caclulation time: 0.078 sec.
Iteration count: 11
Epsilon: 0.000010
Resulting mean error: 2.896847
===============================================

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

форум разработчиков