在前面的文章中，我们讲到使用FFD形变作为坐标变换模型，使用梯度下降法作为优化算法来寻求FFD的最优控制参数：
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在上篇文章中已经讲过，此处不再重复贴出来：
求海塞矩阵代码：
-
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图像进行配准，结果如下图所示。
基准图像
浮动图像
配准图像
基准图像与浮动图像的差值图
基准图像与配准图像的差值图
目标函数值的下降过程
由上图可知，LM算法的优化结果比梯度下降法更理想，其收敛速度更快，且优化结果更接近最优解。
本人微信公众号如下，会不定时更新更精彩的内容，欢迎扫码关注：
转载：https://blog.csdn.net/shandianfengfan/article/details/113778069