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

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

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

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

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

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

梯度下降法与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算法耗时。