前面我们讲到梯度下降法时，就有提到：人们经常要求解一个问题的最优解，通常做法是对该问题进行数学建模，转换成一个目标函数，然后通过一定的算法寻求该函数的最小值，最终寻求到最小值时的输入参数就是问题的最优解。
当我们有两张图像A和B，图像A与图像B形状相似，但具有位置偏移：平移、旋转、缩放、局部扭曲形变等。如果以图像A为基准图像，以图像B为浮动图像，对图像B进行配准，使其与图像A位置对齐，那么我们可以使用FFD自由形变作为坐标变换模型，对图像B进行形变，并计算图像A与形变之后的图像B的相似度，通过求解FFD形变的最优控制参数，使得两者相似度达到最大，从而使用最优控制参数对图像B进行形变，实现其配准。
我们把图像A与FFD形变之后图像B的相似度作为目标函数，然后使用优化算法来求解这个目标函数的最优解，本文我们使用的优化算法是梯度下降法。
目标函数的示意图所下图所示：
梯度下降法求最优控制参数的示意图如下图所示：
下面我们分点讲解所有的步骤，然后再上C++代码。
1. 梯度下降法原理。
其原理我们在之前的文章已经详细讲解，此处不再重复，读者可参考博主的以下博文：
2. FFD自由变换原理。
FFD的原理我们在之前的博文中也讲解过：
3. 图像的相似度衡量与目标函数。
常见的图像相似度衡量指标有峰值信噪比（PSNR）、结构相似度（SSIM）、归一化互相关（NCC）、归一化互信息（NMI）、均方误差（MSE）等。由于归一化互相关系数计算相对简单，且具有的良好凹凸特性利于求解最优参数，因此经过综合考虑，本文选择归一化互相关系数作为相似度的衡量指标。假设图像A与图像B都是m行n列的图像，那么归一化互相关系数可以按照下式计算：
上式中，归一化互相关NCC越大，说明图像A与图像B越相似，反之则两者差异越大。由于梯度下降法为求解目标函数的最小值，所以我们需要把以上函数取个倒数，使得NCC越小，图像的相似度越高：
需要注意，计算NCC的图像B是经过FFD形变的，因此目标函数的数学表达式如下，其中X为FFD形变模型的所有控制参数。
很多时候，图像之间存在很多局部的形变差异，所以通过分块来求解的NCC'更能表现两图的相似度。因此我们在以上基础上，在对图像A与图像B进行相同的分块，计算两图中对应位置块的NCC'，最后再取所有块的NCC'的平均值作为整张图像的NCC'。假设把图像的高平均分为r块，宽平均分为c块，那么最终的目标函数F的表达式如下：
4. C++代码实现。
(1) 分块计算归一化互相关代码
-
double cal_cc_block(Mat S1, Mat Si, int row, int col)
-
{
-
int ksize_row = S1.rows/row;
-
int ksize_col = S1.cols/col;
-
-
-
Mat tmp1, tmpi;
-
double sum =
0.0;
-
for(
int i =
0; i < row; i++)
-
{
-
int i_begin = i*ksize_row;
-
for(
int j =
0; j < col; j++)
-
{
-
double sum1 =
0.0;
-
double sum2 =
0.0;
-
double sum3 =
0.0;
-
int j_begin = j*ksize_col;
-
for (
int t1 = i_begin; t1 < i_begin+ksize_row; t1++)
-
{
-
uchar *S1_data = S1.ptr<uchar>(t1);
-
uchar *Si_data = Si.ptr<uchar>(t1);
-
for (
int t2 = j_begin; t2 < j_begin+ksize_col; t2++)
-
{
-
sum1 += S1_data[t2]*Si_data[t2];
-
sum2 += S1_data[t2]*S1_data[t2];
-
sum3 += Si_data[t2]*Si_data[t2];
-
}
-
}
-
sum +=
sqrt(sum2*sum3)/(sum1+
0.0000001);
-
-
-
}
-
}
-
-
-
sum /= (row*col);
-
-
-
return sum;
-
}
(2) 初始化控制参数代码
-
#define randf(a, b) (((rand()%10000+rand()%10000*10000)/100000000.0)*((b)-(a))+(a))
-
-
-
void init_bpline_para(Mat src, int row_block_num, int col_block_num, Mat &grid_points, float min, float max)
-
{
-
int grid_rows = row_block_num + BPLINE_BOARD_SIZE;
-
int grid_cols = col_block_num + BPLINE_BOARD_SIZE;
-
-
int grid_size = grid_rows*grid_cols;
-
grid_points.create(Size(
2*grid_size,
1), CV_32FC1);
-
-
-
float *grid_points_data = grid_points.ptr<
float>(
0);
-
srand((
unsigned
int)time(
NULL));
-
for (
int i =
0; i < grid_size; i++)
-
{
-
grid_points_data[i] = randf(min, max);
//x
-
-
-
int cnt =
100000000;
-
while(cnt--);
-
-
grid_points_data[i+grid_size] = randf(min, max);
//y
-
-
-
cnt =
100000000;
-
while(cnt--);
-
}
-
}
(3) 基于C++与CUDA实现的FFD变换代码
-
__global__ void Bspline_Ffd_kernel(uchar *srcimg, uchar *dstimg, int row_block_num, int col_block_num, float *grid_points, int row, int col)
-
{
-
int x = threadIdx.x + blockDim.x * blockIdx.x;
//col
-
int y = threadIdx.y + blockDim.y * blockIdx.y;
//row
-
-
-
if(x < col && y < row)
-
{
-
int grid_rows = row_block_num + BPLINE_BOARD_SIZE;
-
int grid_cols = col_block_num + BPLINE_BOARD_SIZE;
-
int grid_size = grid_rows*grid_cols;
-
float delta_x = col*
1.0/col_block_num;
-
float delta_y = row*
1.0/row_block_num;
-
float x_block = x / delta_x;
-
float y_block = y / delta_y;
-
-
int j =
floor(x_block);
-
int i =
floor(y_block);
-
float u = x_block - j;
-
float v = y_block - i;
-
-
-
float pX[
4], pY[
4];
-
pX[
0] = (
1 - u*u*u +
3*u*u -
3*u) /
6.0;
-
pX[
1] = (
4 +
3*u*u*u -
6*u*u) /
6.0;
-
pX[
2] = (
1 -
3*u*u*u +
3*u*u +
3*u) /
6.0;
-
pX[
3] = u*u*u /
6.0;
-
pY[
0] = (
1 - v*v*v +
3*v*v -
3*v) /
6.0;
-
pY[
1] = (
4 +
3*v*v*v -
6*v*v) /
6.0;
-
pY[
2] = (
1 -
3*v*v*v +
3*v*v +
3*v) /
6.0;
-
pY[
3] = v*v*v /
6.0;
-
-
-
float Tx =
0;
-
float Ty =
0;
-
for (
int m =
0; m <
4; m++)
//行
-
{
-
for (
int n =
0; n <
4; n++)
//列
-
{
-
int control_point_x = j + n;
-
int control_point_y = i + m;
-
-
-
float temp = pY[m] * pX[n];
-
-
Tx += temp*grid_points[control_point_y*grid_cols+control_point_x];
//x
-
Ty += temp*grid_points[control_point_y*grid_cols+control_point_x+grid_size];
//y
-
}
-
}
-
-
float src_x = x + Tx;
-
float src_y = y + Ty;
-
int x1 =
floor(src_x);
-
int y1 =
floor(src_y);
-
-
if (x1 <
1 || x1 >= col
-1 || y1 <
1 || y1 >= row
-1)
//越界
-
{
-
dstimg[y*col+x] =
0;
-
}
-
else
-
{
-
//dstimg[y*col+x] = srcimg[y1*col+x1]; //最邻近插值
-
int x2 = x1 +
1;
//双线性插值
-
int y2 = y1 +
1;
-
uchar pointa = srcimg[y1*col+x1];
-
uchar pointb = srcimg[y1*col+x2];
-
uchar pointc = srcimg[y2*col+x1];
-
uchar pointd = srcimg[y2*col+x2];
-
uchar gray = (uchar)((x2 - src_x)*(y2 - src_y)*pointa - (x1 - src_x)*(y2 - src_y)*pointb - (x2 - src_x)*(y1 - src_y)*pointc + (x1 - src_x)*(y1 - src_y)*pointd);
-
dstimg[y*col+x] = gray;
-
}
-
-
-
}
-
-
-
}
-
-
-
void Bspline_Ffd_cuda(Mat srcimg, Mat &dstimg, int row_block_num, int col_block_num, Mat grid_points)
-
{
-
-
dim3 Bpline_Block(16, 16);
//每个线程块有16*16个线程
-
int M = (srcimg.cols+Bpline_Block.x
-1)/Bpline_Block.x;
-
int N = (srcimg.rows+Bpline_Block.y
-1)/Bpline_Block.y;
-
dim3 Bpline_Grid(M, N);
-
-
-
int grid_rows = row_block_num + BPLINE_BOARD_SIZE;
-
int grid_cols = col_block_num + BPLINE_BOARD_SIZE;
-
int grid_size = grid_rows*grid_cols;
-
int img_size = srcimg.cols*srcimg.rows;
-
-
-
uchar *srcimg_cuda;
-
uchar *dstimg_cuda;
-
float *grid_points_cuda;
-
-
-
cudaMalloc((
void**)&srcimg_cuda, img_size);
-
cudaMalloc((
void**)&dstimg_cuda, img_size);
-
cudaMalloc((
void**)&grid_points_cuda,
2*grid_size*
sizeof(
float));
-
cudaMemcpy(srcimg_cuda, srcimg.data, img_size, cudaMemcpyHostToDevice);
-
cudaMemcpy(grid_points_cuda, grid_points.data,
2*grid_size*
sizeof(
float), cudaMemcpyHostToDevice);
-
-
-
Bspline_Ffd_kernel<< <Bpline_Grid, Bpline_Block >> >(srcimg_cuda, dstimg_cuda, row_block_num, col_block_num, grid_points_cuda, srcimg.rows, srcimg.cols);
-
-
-
Mat tmp(srcimg.size(), CV_8UC1);
-
cudaMemcpy(tmp.data, dstimg_cuda, img_size, cudaMemcpyDeviceToHost);
-
tmp.copyTo(dstimg);
-
-
-
cudaFree(srcimg_cuda);
-
cudaFree(dstimg_cuda);
-
cudaFree(grid_points_cuda);
-
}
(4) 目标函数代码
-
float F_fun_bpline(Mat S1, Mat Si, int row_block_num, int col_block_num, Mat grid_points)
-
{
-
double result;
-
Mat Si_tmp;
-
-
Bspline_Ffd_cuda(Si, Si_tmp, row_block_num, col_block_num, grid_points);
-
-
-
result = cal_cc_block(S1, Si_tmp,
5,
5);
//默认分为5*5块计算互相关
-
-
-
return result;
-
}
(5) 求取梯度代码
-
void cal_gradient(Mat S1, Mat Si, int row_block_num, int col_block_num, Mat grid_points, Mat &gradient)
-
{
-
float EPS =
1;
//1e-4f;
-
-
-
gradient.create(grid_points.size(), CV_32FC1);
-
-
float a1 = F_fun_bpline(S1, Si, row_block_num, col_block_num, grid_points);
-
-
-
Mat grid_p = grid_points.clone();
-
for(
int i =
0; i < grid_points.cols; i++)
-
{
-
grid_p.at<
float>(
0, i) += EPS;
-
float a2 = F_fun_bpline(S1, Si, row_block_num, col_block_num, grid_p);
-
grid_p.at<
float>(
0, i) -= EPS;
-
gradient.at<
float>(
0, i) = (a2 - a1) / EPS;
-
}
-
}
(6) 使用梯度来更新控制参数的代码
-
void update_grid_points(Mat &grid_points, Mat gradient, float alpha)
-
{
-
for(
int i =
0; i < grid_points.cols; i++)
-
{
-
grid_points.at<
float>(
0, i) = grid_points.at<
float>(
0, i) - gradient.at<
float>(
0, i)*alpha;
-
}
-
}
(7) 梯度下降法代码
-
int bpline_match(Mat S1, Mat Si, Mat &M, int row_block_num, int col_block_num, Mat &grid_points)
-
{
-
int max_iter =
5000;
//最多迭代次数
-
Mat gradient, pre_gradient;
-
Mat pre_grid_points;
-
double e =
0.000005;
//定义迭代精度
-
float ret1 =
0.0;
-
float ret2 =
0.0;
-
int cnt =
0;
-
float alpha =
8000000;
-
//求梯度
-
cal_gradient(S1, Si, row_block_num, col_block_num, grid_points, gradient);
//求梯度
-
int out_cnt =
0;
-
-
-
while (cnt < max_iter)
-
{
-
-
-
pre_grid_points = grid_points.clone();
-
update_grid_points(grid_points, gradient, alpha);
//更新输入参数
-
-
-
ret1 = F_fun_bpline(S1, Si, row_block_num, col_block_num, pre_grid_points);
//F_fun(S1, Si, S1_entropy, delta_x, delta_y, pre_grid_points);
-
ret2 = F_fun_bpline(S1, Si, row_block_num, col_block_num, grid_points);
//F_fun(S1, Si, S1_entropy, delta_x, delta_y, grid_points);
-
-
-
if (ret2 > ret1)
//如果当前轮迭代的目标函数值大于上一轮的函数值，则减小步长并重新计算梯度、重新更新参数
-
{
-
alpha *=
0.8;
-
grid_points = pre_grid_points.clone();
-
continue;
-
}
-
-
-
cout<<
"f="<<ret2<<
", alpha="<<alpha<<
endl;
-
-
-
if (
abs(ret2 - ret1) < e)
-
{
-
out_cnt++;
-
if(out_cnt >=
4)
//如果连续4次目标函数值不改变，则认为求到了最优解，停止迭代
-
{
-
Bspline_Ffd_cuda(Si, M, row_block_num, col_block_num, grid_points);
-
return
0;
-
}
-
}
-
else
-
{
-
out_cnt =
0;
-
}
-
gradient.copyTo(pre_gradient);
-
//求梯度
-
cal_gradient(S1, Si, row_block_num, col_block_num, grid_points, gradient);
//求梯度
-
-
//如果梯度相比上一次迭代有所下降，则增大步长
-
if(norm(gradient, NORM_L2) <= norm(pre_gradient, NORM_L2))
-
alpha *=
2;
-
-
cnt++;
-
}
-
-
-
return
-1;
-
-
-
}
(8) 测试代码
-
void ffd_match_test(void)
-
{
-
Mat img1 = imread(
"lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
-
Mat img2 = imread(
"lena_out.jpg", CV_LOAD_IMAGE_GRAYSCALE);
-
-
-
int row_block_num =
30;
-
int col_block_num =
30;
-
-
Mat grid_points;
-
init_bpline_para(img1, row_block_num, col_block_num, grid_points,
-0.001,
0.001);
-
-
-
Mat
out;
-
bpline_match(img1, img2,
out, row_block_num, col_block_num, grid_points);
-
-
-
imshow(
"img1", img1);
-
imshow(
"img2", img2);
-
imshow(
"out",
out);
-
imshow(
"img1-img2", abs(img1-img2));
-
imshow(
"img1-out", abs(img1-
out));
-
waitKey();
-
}
运行上述代码，对扭曲的Lena图像进行配准，结果如下图所示。可以看到，配准之后，配准图像没有那么扭曲了，也即与基准图像的位置更加对齐了。
基准图像
浮动图像
配准图像
基准图像与浮动图像的差值图
基准图像与配准图像的差值图
目标函数值的下降过程
微信公众号如下，欢迎关注：
转载：https://blog.csdn.net/shandianfengfan/article/details/113750401
查看评论