在前面的文章中，我們講到使用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*1low = (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圖像進行配準，結果如下圖所示。

基準圖像

浮動圖像

配準圖像

基準圖像與浮動圖像的差值圖

基準圖像與配準圖像的差值圖

目標函數值的下降過程

由上圖可知，LM算法的優化結果比梯度下降法更理想，其收斂速度更快，且優化結果更接近最優解。

本人微信公眾號如下，會不定時更新更精彩的內容，歡迎掃碼關注：

總結

以上是生活随笔為你收集整理的图像配准系列之基于FFD形变与LM算法的图像配准的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。