demons算法是一种全局坐标变换模型的配准算法，该算法使用参考图像的梯度以及参考图像与浮动图像的灰度差值来计算每一个点的坐标偏移量，从而得到参考图像与浮动图像的整幅图的坐标偏移量，并使用坐标偏移量对浮动图像进行重采样和插值，重复迭代此过程直到迭代次数达到设定次数或参考图像与浮动图像相似度达到设定阈值为止。
01
demons算法的核心思想
demons算法是一种基于光流理论发展起来的配准方法，因此该算法与其它光流算法一样，具有亮度恒定和小运动的约束条件。假设t时刻图像I上点(x,y)的亮度值为I(x,y,t)，根据运动过程中亮度恒定有：
对上式等号右边泰勒展开得到：
因为前提条件是小运动，dx，dy，dt都很小，因此可以把余项ε忽略，于是有：
记：
那么有：
上式就是demons算法的核心思想，该式表明使用梯度信息和图像随时间变化的差值信息（也就是参考图像与浮动图像的帧间差）可以计算两帧图像之间点的偏移量。需要注意的是必须在小运动条件之下上式中的余项ε才可以忽略，因此当运动很大的时候把ε忽略掉会严重影响计算结果的精确度，由此可知大运动对demons算法配准结果的影响还是较大的。
02
原始demons算法原理
根据上式，该算法的迭代过程也可以看作是浮动图像中各个像素点逐步向参考图像中对应位置点扩散的过程，其中两图像对应点的灰度差值是扩散的外力，而参考图像对应点的梯度则是扩散的内力。
假设参考图像为S，浮动图像为M，则参考图像上点(x,y)的偏移量(Ux,Uy)计算如下式，其中Sx和Sy分别是参考图像上点(x,y)处x方向和y方向的梯度，▽f则是点(x,y)处参考图像与浮动图像的灰度差。
计算得到整幅图像的坐标偏移量之后，为了使偏移量在全局范围内是平滑连续的，每次迭代时还对整幅图的坐标偏移量进行高斯平滑，这样可以有效避免重采样之后图像出现的毛刺现象。
03
增加扩散速度系数α的demons算法
在原始demons算法的基础上，P. Cachier提出了增加扩散速度系数α来控制坐标偏移量（扩散速度）的大小，如下式所示，α越大偏移量越小，反之偏移量越大。在算法的迭代过程中为了提高配准精度，通常会随着迭代次数增加而逐渐增大α，即减小偏移量。
04
Active demons算法
随后H.Wang等人提出了Active demons算法，把浮动图像的梯度加入偏移量的计算中。原始算法中驱动扩散的内力为参考图像的梯度，这时增加了浮动图像的梯度作为新的内力，因此加快了迭代的收敛速度，其计算偏移量如下式所示，其中Mx和My分别是浮动图像上点(x,y)处x方向和y方向的梯度。
05
Inertial demons算法
后来Santos-Ribeiro A等人在Active demons算法的基础上提出了Inertial demons算法，即把上一层迭代计算得到的偏移量加入到当前层迭代的偏移量计算当中，进一步提升了收敛速度和配准精度。其计算如下式所示，其中k为当前的迭代次数，β为0~1的系数：
06
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();
}07
配准测试结果
1. Lena图像配准结果
参考图像
浮动图像
多次迭代过程中浮动图像的变化过程
多次迭代过程中浮动图像与参考图像的差值图变化过程
2. 脑部图像配准结果
参考图像
浮动图像
多次迭代过程中浮动图像的变化过程
多次迭代过程中浮动图像与参考图像的差值图变化过程
配准时α参数很重要，通常需要多试几个合适的参数，或者设定一个较小的初始值，随着迭代次数的增加，逐渐增大α值，这样配准效果会更好。高斯滤波参数通常取59*59窗口、10.0的sigma值。迭代次数取50~100次就好了。
本文就讲到这里，demons算法后面还发展了很多变种算法，比如微分同胚demons算法，这个我们后续再进一步研究。
欢迎关注“萌萌哒程序猴”微信公众号，接下来会不定时更新更加精彩的内容噢～
本文暂时没有评论，来添加一个吧(●'◡'●)