图像配准算法与三维重建（图像配准系列之基于FFD形变与粒子群算法的图像配准）

在之前的文章中，我们分别使用了梯度下降发与LM算法来优化FFD形变的控制参数，达到图像配准的目的：

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

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

在本文中，我们改为使用粒子群算法来来优化FFD形变的控制参数（相似度衡量指标不变）。粒子群算法的原理，我们在之前的文章也有讲过：

粒子群（PSO）算法的理解与应用

梯度下降法与LM算法都是单线程的寻找最优解，而粒子群算法则不一样，其有多个解（也即多个粒子）同时进行优化，每轮迭代时都在多个粒子中记录最接近最优解的粒子，粒子群算法相当于一种多线程寻求最优解的算法。因此该算法更不容易陷入局部极值。

同样假设FFD变换模型有r 3行c 3列的控制点，每个控制点有x方向、y方向的两个控制参数，因此总共有N=2*(r 3)*(c 3)个控制参数需要优化，也就是说，每个粒子的数据维度为N。

下面直接上代码。

1. 一些全局变量：

constintrow_block_num_pso=30;//FFD网格的行数 const int col_block_num_pso = 30; //FFD网格的列数 //FFD网格的控制参数个数 const int DATA_SIZE = (row_block_num_pso BPLINE_BOARD_SIZE)*(col_block_num_pso BPLINE_BOARD_SIZE)*2; constintNUM=300;//粒子数 constfloatc1=1.8; //粒子群参数1 constfloatc2=1.8;//粒子群参数2 //控制参数被初始化为-1到1之间的随机数 float xmin = -1; float xmax = 1; //粒子群的速度范围被钳制在为-20到20之间，这是经验值，合适的钳制范围可以加快收敛速度 constfloatvmin=-10; constfloatvmax=10; //定义粒子群，粒子个数为NUM，每个粒子为一个结构体 struct particle { floatx[DATA_SIZE];//当前粒子包含的控制参数 floatbestx[DATA_SIZE];//当前粒子的历史最优控制参数 float f; //前粒子包含的控制参数对应的目标函数值 float bestf; //当前粒子的历史最优控制参数对应的目标函数值 }swarm[NUM];

2. 延时函数代码。

这里的延时函数，是在连续获取随机数时，增加一定的延时间隔，或许能增加随机数的随机性（个人经验，有待考证）。

void delay_for(long int cnt) { while(cnt--); }

3. 粒子群优化代码。

void PSO(Mat S1, Mat Si, Mat &M, Mat &grid_points) { grid_points.create(1, DATA_SIZE, CV_32FC1); float *grid_points_p = grid_points.ptr<float>(0); for (int i = 0; i < DATA_SIZE; i )//初始化全局最优 { grid_points_p[i] = randf(xmin, xmax); } float gbestf_pre = 0; float gbestf = 100000000.0; Mat para_x(1, DATA_SIZE, CV_32FC1); //初始化粒子群 for (int i = 0; i < NUM; i ) { particle *p1 = &swarm[i]; for (int j = 0; j < DATA_SIZE; j ) { p1->x[j] = randf(xmin, xmax); p1->bestx[j] = randf(xmin, xmax); } memcpy((float *)para_x.data, p1->x, DATA_SIZE*sizeof(float)); p1->f = F_fun_bpline(S1, Si, row_block_num_pso, col_block_num_pso, para_x); p1->bestf = 100000000.0; } float *V = (float *)calloc(DATA_SIZE, sizeof(float)); const int cnt = 5000; float w = 0.0025/(cnt-1); srand((unsigned)time(NULL)); intcntt=0; for (int t = 0; t < cnt; t ) { for (int i = 0; i < NUM; i ) { particle*p1=&swarm[i]; for (int j = 0; j < DATA_SIZE; j ) //进化方程 { float d1 = randf(0, 1); delay_for(100000); float d2 = randf(0, 1); V[j] = w*(cnt-1-t)*V[j] c1*d1*(p1->bestx[j] - p1->x[j]) c2*d2*(grid_points.ptr<float>(0)[j] - p1->x[j]); V[j] = (V[j] < vmin) ? vmin : ((V[j] > vmax) ? vmax : V[j]); p1->x[j] = p1->x[j] V[j]; } memcpy((float *)para_x.data, p1->x, DATA_SIZE*sizeof(float)); p1->f = F_fun_bpline(S1, Si, row_block_num_pso, col_block_num_pso, para_x); if (p1->f < p1->bestf) //改变该粒子的历史最优 { for (int j = 0; j < DATA_SIZE; j ) { p1->bestx[j] = p1->x[j]; } p1->bestf = p1->f; } if (p1->bestf < gbestf) //改变所有例子的全局最优 { for (int j = 0; j < DATA_SIZE; j ) { grid_points.ptr<float>(0)[j] = p1->bestx[j]; } for (int j = 0; j < DATA_SIZE; j ) //把当前全局最优的粒子随机放到另一位置 { p1->x[j] = randf(xmin, xmax); } gbestf_pre = gbestf; gbestf = p1->bestf; printf("t=%d,gbestf=%lf\n",t,gbestf); } } if(abs(gbestf_pre-gbestf)<1e-6) { cntt ; if(cntt>=1000) { break; } } else { cntt = 0; } } free(V); Bspline_Ffd_cuda(Si,M,row_block_num_pso,col_block_num_pso,grid_points); write_data_file("gradient_list.m",d_list); }

4. 测试代码：

void ffd_match_pso_test(void) { Mat img1 = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE); Matimg2=imread("lena_out.jpg",CV_LOAD_IMAGE_GRAYSCALE); Matgrid_points; Mat out; PSO(img1, img2, out, grid_points); imshow("img1", img1); imshow("img2", img2); imshow("out", out); imshow("img1-img2", abs(img1-img2)); imshow("img1-out", abs(img1-out)); waitKey(); }

5. 运行结果：

运行上述代码，同样对扭曲的Lena图像进行配准，结果如下图所示。

参考图像

浮动图像

配准图像

参考图像与浮动图像的差值图

参考图像与配准图像的差值图

全局最优目标函数值的降低过程

至此，我们分别使用了三种优化算法来优化FFD形变的控制参数：梯度下降发、LM算法、粒子群算法，梯度下降法比较稳定，但容易陷入局部极值，LM算法兼具稳定与不容易陷入局部极值的特性，不过当参数量很大时LM算法计算海塞矩阵或者矩阵的逆很是耗时，相比来说粒子群算法的不容易陷入局部极值特性更好，而且也没有LM算法耗时。