在之前的文章中，我们分别使用了梯度下降发与LM算法来优化FFD形变的控制参数，达到图像配准的目的：
在本文中，我们改为使用粒子群算法来来优化FFD形变的控制参数（相似度衡量指标不变）。粒子群算法的原理，我们在之前的文章也有讲过：
梯度下降法与LM算法都是单线程的寻找最优解，而粒子群算法则不一样，其有多个解（也即多个粒子）同时进行优化，每轮迭代时都在多个粒子中记录最接近最优解的粒子，粒子群算法相当于一种多线程寻求最优解的算法。因此该算法更不容易陷入局部极值。
同样假设FFD变换模型有r+3行c+3列的控制点，每个控制点有x方向、y方向的两个控制参数，因此总共有N=2*(r+3)*(c+3)个控制参数需要优化，也就是说，每个粒子的数据维度为N。
下面直接上代码。
1. 一些全局变量：
-
const
int row_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;
-
-
-
const
int NUM =
300;
//粒子数
-
const
float c1 =
1.8;
//粒子群参数1
-
const
float c2 =
1.8;
//粒子群参数2
-
//控制参数被初始化为-1到1之间的随机数
-
float xmin =
-1;
-
float xmax =
1;
-
//粒子群的速度范围被钳制在为-20到20之间，这是经验值，合适的钳制范围可以加快收敛速度
-
const
float vmin =
-10;
-
const
float vmax =
10;
-
-
-
-
-
//定义粒子群，粒子个数为NUM，每个粒子为一个结构体
-
struct particle
-
{
-
float x[DATA_SIZE];
//当前粒子包含的控制参数
-
-
-
float bestx[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));
-
-
-
int cntt =
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);
-
Mat img2 = imread(
"lena_out.jpg", CV_LOAD_IMAGE_GRAYSCALE);
-
-
-
Mat grid_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算法耗时。
本人微信公众号如下，会不定时更新更精彩的内容，欢迎扫码关注：
转载：https://blog.csdn.net/shandianfengfan/article/details/113812923
查看评论