图像配准系列之基于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算法耗时。

本人微信公众号如下,会不定时更新更精彩的内容,欢迎扫码关注:

  • 8
    点赞
  • 17
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 7
    评论
评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

萌萌哒程序猴

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值