Алгоритмы совмещения изображений - Реализация поступательного аддитивного алгоритма
ОГЛАВЛЕНИЕ
Реализация поступательного аддитивного алгоритма
Теперь будет подробно описано, как реализуется алгоритм Лукаса-Канаде. Как говорят Бейкер и Мэтьюс, поступательный аддитивный алгоритм состоит из следующих шагов:
Предварительно вычислить:
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);
}