图像配准系列之基于FFD形变与LM算法的图像配准

大家好，欢迎来到IT知识分享网。

在前面的文章中，我们讲到使用FFD形变作为坐标变换模型，使用梯度下降法作为优化算法来寻求FFD的最优控制参数：

图像配准系列之基于FFD形变与梯度下降法的图像配准

LM算法可以看作是梯度下降法与高斯-牛顿法的结合算法，它既具有梯度下降法的稳健性，又具有高斯-牛顿法的快速收敛性，而且相比来说更不容易陷入局部极值。有关LM算法的数学公式推导，我们也在前文中有详细说明：

最优化算法之牛顿法、高斯-牛顿法、LM算法

本文我们将使用LM算法作为优化算法来寻求FFD变换的最优控制参数。并与梯度下降法比较配准结果。

1. LM算法的计算过程

(1) 差分法求控制参数的梯度。

假设FFD变换模型有r+3行c+3列的控制点，每个控制点有x方向、y方向的两个控制参数，因此总共有N=2*(r+3)*(c+3)个控制参数需要优化，所有控制参数组成一个1行N列的向量X：

对于每个控制参数，其梯度的计算如下式，其中F为目标函数，EPS为一个较小的数，一般取1或者0.5即可。

所有控制点的梯度组成一个1行N列的梯度向量：

(2) 计算当前控制参数的目标函数值f0、矩阵gk、海塞矩阵H。

(3) 使用海塞矩阵计算矩阵G、矩阵h。

上式中，u为LM算法的控制步长，通常u取一个较小的初始值（比如0.001），在迭代优化过程中u的值根据情况而改变，详细如何改变在后续步骤说明。矩阵I为N*N的单位矩阵：

(4) 判断是否满足以下条件，如果满足则认为找到最优控制参数，因此停止循环迭代：

上式中，norm函数为L2范数，e通常取一个很小的值，比如10-12。比如X的L2范数可按下式计算：

(5) 更新X得到X’，然后计算X’的目标函数值f1，并计算步长因子ρ。

(6) 判断ρ是否大于0。

I. 如果ρ大于0则减小u的值：

则把X’赋值给X，并改变u的值。接着判断迭代次数是否达到设定的最大次数，如果达到则停止迭代，否则跳转到第(1)步执行。

II. 如果ρ小于等于0则增大u和v：

接着判断迭代次数是否达到设定的最大次数，如果达到则停止迭代，否则跳转到第(3)步执行。

上式中，在迭代开始之前通常把v初始化值为2。

2. C++实现代码

cal_gradient、F_fun_bpline、Bspline_Ffd_cuda、init_bpline_para在上篇文章中已经讲过，此处不再重复贴出来：

图像配准系列之基于FFD形变与梯度下降法的图像配准

求海塞矩阵代码：

Mat cal_hessian_matrix(Mat gradient) { Mat gradient_trans;   Mat hessian;   transpose(gradient, gradient_trans);   //转置   hessian = gradient_trans*gradient; return hessian; }

LM算法代码：

void bpline_match_LM(Mat S1, Mat Si, Mat &M, int row_block_num, int col_block_num, Mat &grid_points) { int iter = 150; //迭代次数   double e = 1e-12;  // 误差限 double u = 1e-3; double v = 2.0; Mat H; Mat G; Mat new_grid_points; Mat I = cv::Mat::eye(grid_points.cols, grid_points.cols, CV_32FC1); //单位矩阵 Mat h, h_T; Mat gradient, gradient_T; float f0, f1; Mat gk; double low = 1.0;    cal_gradient(S1, Si, row_block_num, col_block_num, grid_points, gradient); //求梯度 f0 = F_fun_bpline(S1, Si, row_block_num, col_block_num, grid_points); transpose(gradient, gradient_T); gk = gradient_T*f0; H = cal_hessian_matrix(gradient); //由雅克比矩阵计算海塞矩阵 for(int i = 0; i < iter; i++) { G = H + u*I; Mat G_T; invert(G, G_T, DECOMP_SVD); h = -G_T*gk;     if(norm(h, NORM_L2) < e*(norm(grid_points, NORM_L2)+e)) break; transpose(h, h_T); new_grid_points = grid_points + h_T; f1 = F_fun_bpline(S1, Si, row_block_num, col_block_num, new_grid_points); Mat L_0_h = 0.5*h_T*(u*h-gk); //1_N*N_1 = 1*1 low = (f0 - f1)/L_0_h.at<float>(0, 0); if(low > 0) { grid_points = new_grid_points.clone(); double tmp = 2*low-1;//5*low-1; tmp = 1 - tmp*tmp*tmp; double t = 0.; u = u*(t > tmp ? t : tmp); v = 2; cal_gradient(S1, Si, row_block_num, col_block_num, grid_points, gradient); //求梯度 f0 = F_fun_bpline(S1, Si, row_block_num, col_block_num, grid_points); transpose(gradient, gradient_T); gk = gradient_T*f0;       H = cal_hessian_matrix(gradient);    //由雅克比矩阵计算海塞矩阵 } else { u *= v; v *= 2; } cout<<"f1="<<f1<<", low="<<low<<", u="<<u<<endl; }   Bspline_Ffd_cuda(Si, M, row_block_num, col_block_num, grid_points); }

测试代码：

void ffd_match_LM_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_LM(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图像进行配准，结果如下图所示。

基准图像

浮动图像

配准图像

基准图像与浮动图像的差值图

基准图像与配准图像的差值图

目标函数值的下降过程

由上图可知，LM算法的优化结果比梯度下降法更理想，其收敛速度更快，且优化结果更接近最优解。