图像配准算法之demons算法

demons算法是一种全局坐标变换模型的配准算法，该算法使用参考图像的梯度以及参考图像与浮动图像的灰度差值来计算每一个点的坐标偏移量，从而得到参考图像与浮动图像的整幅图的坐标偏移量，并使用坐标偏移量对浮动图像进行重采样和插值，重复迭代此过程直到迭代次数达到设定次数或参考图像与浮动图像相似度达到设定阈值为止。

demons算法的核心思想

demons算法是一种基于光流理论发展起来的配准方法，因此该算法与其它光流算法一样，具有亮度恒定和小运动的约束条件。假设t时刻图像I上点(x,y)的亮度值为I(x,y,t)，根据运动过程中亮度恒定有：

对上式等号右边泰勒展开得到：

因为前提条件是小运动，dx，dy，dt都很小，因此可以把余项ε忽略，于是有：

记：

那么有：

上式就是demons算法的核心思想，该式表明使用梯度信息和图像随时间变化的差值信息（也就是参考图像与浮动图像的帧间差）可以计算两帧图像之间点的偏移量。需要注意的是必须在小运动条件之下上式中的余项ε才可以忽略，因此当运动很大的时候把ε忽略掉会严重影响计算结果的精确度，由此可知大运动对demons算法配准结果的影响还是较大的。

—

原始demons算法原理

根据上式，该算法的迭代过程也可以看作是浮动图像中各个像素点逐步向参考图像中对应位置点扩散的过程，其中两图像对应点的灰度差值是扩散的外力，而参考图像对应点的梯度则是扩散的内力。

假设参考图像为S，浮动图像为M，则参考图像上点(x,y)的偏移量(U_x,U_y)计算如下式，其中S_x和S_y分别是参考图像上点(x,y)处x方向和y方向的梯度，▽f则是点(x,y)处参考图像与浮动图像的灰度差。

计算得到整幅图像的坐标偏移量之后，为了使偏移量在全局范围内是平滑连续的，每次迭代时还对整幅图的坐标偏移量进行高斯平滑，这样可以有效避免重采样之后图像出现的毛刺现象。

—

增加扩散速度系数α的demons算法

在原始demons算法的基础上，P. Cachier提出了增加扩散速度系数α来控制坐标偏移量（扩散速度）的大小，如下式所示，α越大偏移量越小，反之偏移量越大。在算法的迭代过程中为了提高配准精度，通常会随着迭代次数增加而逐渐增大α，即减小偏移量。

—

Active demons算法

随后H.Wang等人提出了Active demons算法，把浮动图像的梯度加入偏移量的计算中。原始算法中驱动扩散的内力为参考图像的梯度，这时增加了浮动图像的梯度作为新的内力，因此加快了迭代的收敛速度，其计算偏移量如下式所示，其中M_x和M_y分别是浮动图像上点(x,y)处x方向和y方向的梯度。

—

Inertial demons算法

后来Santos-Ribeiro A等人在Active demons算法的基础上提出了Inertial demons算法，即把上一层迭代计算得到的偏移量加入到当前层迭代的偏移量计算当中，进一步提升了收敛速度和配准精度。其计算如下式所示，其中k为当前的迭代次数，β为0~1的系数：

—

Inertial demons代码实现

1. 计算梯度代码

实际的实现过程中，我们使用Sobel梯度来代替差分梯度，由于Sobel梯度比差分梯度信息更丰富，因此可加快收敛速度：

void get_mat_gradient(Mat src, Mat &Fx, Mat &Fy)
{
  Mat src_tmp;
  Mat Fx_tmp;
  Mat Fy_tmp;
  //定义Sobel梯度算子
  Mat ker_x = (Mat_<float>(3, 3) << -1.0, 0.0, 1.0, -2.0, 0.0, 2.0, -1.0, 0.0, 1.0);
  Mat ker_y = (Mat_<float>(3, 3) << -1.0, -2.0, -1.0, 0.0, 0.0, 0.0, 1.0, 2.0, 1.0);


  src.convertTo(src_tmp, CV_32F);
  //使用Sobel算子作卷积运算，得到x、y方向的梯度
  filter2D(src_tmp, Fx_tmp, src_tmp.depth(), ker_x);
  filter2D(src_tmp, Fy_tmp, src_tmp.depth(), ker_y);


  Fx = Fx_tmp.clone();
  Fy = Fy_tmp.clone();
}

2. 计算坐标偏移代码

此处实现的是上述的Active demons算法，即同时使用参考图像、浮动图像来计算坐标偏移。

void demons_one_wang(Mat S, Mat M, Mat Sx, Mat Sy, float alpha, int win_size, float sigma, Mat &Tx, Mat &Ty)
{
  Mat diff = M - S;   //计算差值图
  float a1, a2, ax, ay;
  Mat Tx_tmp(S.size(), CV_32FC1);
  Mat Ty_tmp(S.size(), CV_32FC1);


  Mat Mx, My;
  get_mat_gradient(M, Mx, My);    //求浮动图像的梯度


  for (int i = 0; i < S.rows; i++)
  {
    for (int j = 0; j < S.cols; j++)
    {
      a1 = pow(Sx.at<float>(i, j), 2) + pow(Sy.at<float>(i, j), 2) + pow(alpha, 2)*pow(diff.at<float>(i, j), 2);  //分母
      a2 = pow(Mx.at<float>(i, j), 2) + pow(My.at<float>(i, j), 2) + pow(alpha, 2)*pow(diff.at<float>(i, j), 2);
      if ((a1 > -0.0000001 && a1 < 0.0000001) || (a2 > -0.0000001 && a2 < 0.0000001))
      {
        Tx_tmp.at<float>(i, j) = 0.0;
        Ty_tmp.at<float>(i, j) = 0.0;
      }
      else
      {
        ax = Sx.at<float>(i, j) / a1 + Mx.at<float>(i, j) / a2;   //分子
        Tx_tmp.at<float>(i, j) = (-diff.at<float>(i, j)*ax);


        ay = Sy.at<float>(i, j) / a1 + My.at<float>(i, j) / a2;   //分子
        Ty_tmp.at<float>(i, j) = (-diff.at<float>(i, j)*ay);
      }


    }
  }


  Tx_tmp = Tx_tmp * 10;   //实际测试时，发现对计算得到的坐标偏移扩大一定倍数，可加快收敛速度
  Ty_tmp = Ty_tmp * 10;
  //对坐标偏移进行高斯平滑，减小毛刺
  GaussianBlur(Tx_tmp, Tx, Size(win_size, win_size), sigma, sigma);   
  GaussianBlur(Ty_tmp, Ty, Size(win_size, win_size), sigma, sigma);


}

3. 像素重采样代码

计算得到整张浮动图像的坐标偏移之后，再使用坐标偏移对浮动图像进行重采样。由于坐标偏移是浮点型数据，重采样时需要插值，关于重采样与插值算法原理，可参考以下文章。我们调用Opencv的remap函数可以方便地实现重采样。

常见图像插值算法的原理与C++实现

void movepixels_2d2(Mat src, Mat &dst, Mat Tx, Mat Ty, int interpolation)
{
  Mat src_tmp(src.size(), CV_32FC1);
  Mat dst_tmp(src.size(), CV_32FC1, 0.0);
  Mat Tx_map(src.size(), CV_32FC1, 0.0);
  Mat Ty_map(src.size(), CV_32FC1, 0.0);


  src.convertTo(src_tmp, CV_32F);


  for (int i = 0; i < src.rows; i++)
  {
    for (int j = 0; j < src.cols; j++)
    {
      //根据坐标偏移计算坐标映射
      float x = static_cast<float>(j + Tx.at<float>(i, j));
      float y = static_cast<float>(i + Ty.at<float>(i, j));
      
      //判断坐标映射是否超出范围
      if (x >= 0 && x < src.cols && y >= 0 && y < src.rows)
      {
        Tx_map.at<float>(i, j) = x;
        Ty_map.at<float>(i, j) = y;
      }
      else
      {
        Tx_map.at<float>(i, j) = 0;
        Ty_map.at<float>(i, j) = 0;
      }
    }
  }
  //像素重采样
  remap(src_tmp, dst_tmp, Tx_map, Ty_map, interpolation);    
  dst_tmp.copyTo(dst);
}

4. 迭代过程代码

整个配准过程是一个迭代过程，也即使用上一轮迭代的配准结果作为当前轮迭代的输入，在多次迭代过程中，浮动图像逐渐接近参考图像。此处，我们实现的是Inertial demons算法，也即把上一层迭代计算得到的偏移量加入到当前层迭代的偏移量计算当中。

void Inertial_demons(Mat S, Mat M, Mat &D, float alpha, int win_size, float sigma, int num)
{
  Mat Sx, Sy;
  Mat S_tmp, M_tmp;


  get_mat_gradient(S, Sx, Sy);    //求参考图像的梯度


  Mat Tx_tmp(S.size(), CV_32FC1, 0.0);
  Mat Ty_tmp(S.size(), CV_32FC1, 0.0);
  Mat Tx(S.size(), CV_32FC1, 0.0);
  Mat Ty(S.size(), CV_32FC1, 0.0);


  S.convertTo(S_tmp, CV_32F);   //图像转换为浮点型
  M.convertTo(M_tmp, CV_32F);


  for (int i = 0; i < num; i++)   //迭代过程
  {
    Tx_tmp = Tx.clone();  //保存上一轮迭代结果
    Ty_tmp = Ty.clone();
    demons_one_wang(S_tmp, M_tmp, Sx, Sy, alpha, win_size, sigma, Tx, Ty);
    Tx = Tx + Tx_tmp*0.75;  //把上一层迭代计算得到的偏移量加入到当前层迭代的偏移量计算当中
    Ty = Ty + Ty_tmp*0.75;
    movepixels_2d2(M_tmp, M_tmp, Tx, Ty, INTER_CUBIC);  //像素重采样
  }
  //将浮点型转换为uchar型输出
  M_tmp.convertTo(D, CV_8U);
}

5. 测试代码

运行以下测试代码，分别配准Lena图像和脑部图像：

void demons_test(void)
{
  
  Mat img1 = imread("brain3.png", CV_LOAD_IMAGE_GRAYSCALE);
  Mat img2 = imread("brain4.png", CV_LOAD_IMAGE_GRAYSCALE);
  //Mat img1 = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
  //Mat img2 = imread("lena_ffd.jpg", CV_LOAD_IMAGE_GRAYSCALE)
  
  Mat out;
  Inertial_demons(img1, img2, out, 10.0, 59, 10.0, 100);


  imshow("img1", img1);
  imshow("img2", img2);
  imshow("out", out);
  imshow("img2-img1", abs(img2 - img1));
  imshow("out-img1", abs(out - img1));
  waitKey();
}

—

配准测试结果

1. Lena图像配准结果