Алгоритмы совмещения изображений - Реализация обратного композиционного алгоритма совмещения изображений
ОГЛАВЛЕНИЕ
Страница 4 из 5
Реализация обратного композиционного алгоритма совмещения изображений
В обратном композиционном алгоритме шаблон 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);
}