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

Written on 03 Ноября 2013.

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

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

Предварительно вычислить:
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);
}

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

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

ОГЛАВЛЕНИЕ

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