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

在前面的文章中，我们讲到使用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.333333;
      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图像进行配准，结果如下图所示。