CN101216939B - 基于量子行为粒子群算法的多分辨率医学图像配准方法 - Google Patents

Abstract

本发明涉及一种基于量子行为粒子群算法的多分辨率医学图像配准方法，特征是：将参考图像和待配图像的背景去除后；经过处理的两幅低分辨率的图像，以归一化互信息作为目标函数，使用量子行为粒子群算法求解得到低分辨率图像配准参数；利用三线性PV算法实现插值；计算互信息时，对出界点进行修正处理；经过量子行为粒子群算法求得的配准参数作为初始值，将低分辨率图像的配准参数作为Powell方法的输入，从而完成图像配准。本发明解决了基于互信息的目标函数存在许多局部极值，大大地提高了配准精度和速度，达到了亚像素级；在临床诊断的图像判别、放射治疗的图像定位和外科手术的图像引导等领域有广泛的应用。

Description

基于量子行为粒子群算法的多分辨率医学图像配准方法

技术领域

本发明涉及一种基于量子行为粒子群算法的多分辨率医学图像配准方法，具体地说是可以用来解决使用互信息的目标函数存在的诸多局部极值问题，在临床诊断的图像判别、放射治疗的图像定位和外科手术的图像引导等领域有广泛的应用。

背景技术

在已有技术中，医学图像的配准技术是90年代发展起来的医学图像处理的一个重要分支，是医学图像处理的一项基本任务，对于临床诊断和治疗有重要意义，受到了医学界和工程界的重视。医学图像配准是指对于一幅医学图像寻求一种(或一系列)空间变换，使它与另一幅医学图像上的对应点达到空间上的一致，这种一致是指人体上的同一解剖点在两张匹配图像上有相同的空间位置。配准的结果应使两幅图像上所有解剖点，或至少是所有具有诊断意义的点及感兴趣的点都达到匹配。医学图像配准主要有两大类方法，基于灰度的方法和基于特征的方法。基于灰度的配准方法直接利用图像的灰度数据进行配准，从而避免了因分割而带来的误差，具有精度较高、鲁棒性强、不需要预处理而能实现自动配准，在基于灰度的配准方法中，基于互信息的方法包括互信息和归一化互信息方法已经被广泛使用并具有最高的精度。

基于互信息的医学图像配准中使用较多的优化算法有Powell法、单纯形法、遗传算法、模拟退火算法、粒子群优化算法等。这些优化算法各有优点，但也存在不足之处。例如Powell法与遗传算法都是无需求导数的直接优化法，但遗传算法的收敛速度较慢，而Powell法的优化速度虽然很快，但容易陷入局部最优，而粒子群优化算法虽然收敛快，但由于不是全局优化算法，因此也容易陷入局部最优解。在实际应用中，经常将多种优化算法混合使用，即开始时使用粗略的快速算法，然后使用精确的慢速算法。

发明内容

本发明的目的在于克服上述不足之处，从而提供一种基于量子行为粒子群算法的多分辨率医学图像配准方法，解决了基于互信息的目标函数存在许多局部极值从而给配准的优化过程带来的困难，大大地提高了配准精度和速度，达到了亚像素级，可以应用于临床诊断、放射治疗和图像引导的外科手术等领域；并且为提高配准过程的速度和精度以及鲁棒性，提出了使用多分辨率的方法来进行图像的配准。

按照本发明提供的技术方案，基于量子行为粒子群算法的多分辨率医学图像配准方法，特征是采用以下配准步骤：

(1)、图像的背景去除：首先求取一个阈值，之后采用种子填充算法去除参考图像和待配图像的背景部份：具体步骤如下：

(a)、求出图像中的最大和最小灰度z₁和z_k，令阈值初始值为

$T_0=\frac{Z_1+Z_k}{2}$

(b)、根据阈值T_k将图像分割成R₁和R₂两部分，分别求出两部分的平均灰度值z₀和Z_B：

$Z_0=\frac{\underset{z(i,j)<T_k}{\mathrm{\&Sigma;}}z(i,j)\&times;N(i,j)}{\underset{z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

$Z_B=\frac{\underset{z(i,j)>T_k}{\mathrm{\&Sigma;}}z(i,j)\&times;N(i,j)}{\underset{z(i,j)>T_k}{\mathrm{\&Sigma;}}N(i,j)}$

式中，z(i，j)是图像上(i，j)点的灰度值，N(i，j)是(i，j)点的权重系数，这里取N(i，j)＝1.0；

(c)、求出新的阈值

$T_{k+1}=\frac{Z_0+Z_B}{2}$

如果T_k＝T_k+1，则结束，否则k＝k+1，迭代执行上述步骤；

最后使用种子填充法从图像的左上角开始填充小于阈值的点，去除背景；

(2)、通过小波变换得到低分辨率图像：小波变换得到两幅低分辨率的图像：去除背景后的参考图像和待配图像进行小波变换，利用小波变换方法经一次或多次变换后得到两幅低分辨率的图像；

(3)、利用量子行为粒子群算法求取配准参数：经过处理的两幅低分辨率图像，以归一化互信息作为目标函数，使用量子行为粒子群算法，即使用量子行为粒子群算法求解得到两幅低分辨率图像的配准参数；

(4)、利用Powell方法求取配准参数：经过求得的配准参数作为初始值，以高分辨率图像作为对象，使用归一化互信息作为目标函数，将低分辨率图像的配准参数作为Powell方法的输入，利用共轭方向并以此作为搜索方向，在一定迭代次数后输出最终的配准参数，从而完成图像配准；

所述使用量子行为粒子群算法求取配准参数：

(1)首先在解空间初始化一组粒子，计算粒子的目标函数值，即归一化互信息值，然后粒子通过量子行为粒子群算法的搜索，经过一定的迭代次数后完成寻优过程；

(2)迭代求解过程中，利用三线性PV算法实现插值；

(3)在计算互信息时，对出界点进行修正处理。

所述的利用三线性PV算法实现插值：在求解过程中，由于浮动图像上的点通过空间变换后得到的点坐标不是整数，需要通过插值来得到变换点的灰度值；三线性PV插值不会引入新的灰度值，浮动图像中的一点的灰度对联合直方图是由参考图像中的点的周围最近邻的8个点取与三线性插值相同的权重加权而得到。

所述的出界点进行修正处理：计算互信息时，对出界点进行修正处理，出界点是当浮动图中的样本点经过一定的空间变换后的对应点落在参考图之外的点，互信息的计算必须考虑出界点；处理出界点时使出界点的灰度等于距其最近的边界像素点的灰度。

本发明与已有技术相比具有以下优点：

本发明误配率低，配准速度快，配准精度高。因此，本发明解决了基于互信息的目标函数存在许多局部极值从而给配准的优化过程带来的困难，大大地提高了配准精度和速度，达到了亚像素级。应用于临床诊断、放射治疗和图像引导的外科手术，为他们提供了有效的辅助手段。

本发明的优点可以从如下的实验中得到验证：

选取人头部MRI图像1幅(256像素×256像素)，如图3所示：把图像去除背景，如图4所示：去除背景后的图像顺时针旋转30°，再向下平移、向右各平移5个像素，如图5所示：将去除背景的原图像作为参考图像，变换后的图像作为浮动图像，分别用4种算法，即Powell方法，PSO算法，QPSO算法，QPSO算法与Powell方法结合的方法，分别进行配准，每种算法分别运行10次，得到的配准结果：

1、误配率低：从表1的统计数据可以看出，QPSO参与的方法能够保证配准的结果，具有较强的鲁棒性；

表1各种算法的误配率

算法	误配率(％)
		Powell	50
PSO	20
		QPSO	0
QPSO+Powell	0

2、配准速度快：如图2所示：给出了PSO算法与QPSO算法求解配准参数的目标函数值的收敛曲线，从附图2的曲线可以看出，与PSO算法相比，QPSO算法能够更快的收敛到稳定解，也就表明，QPSO算法具有更快的配准速度；而实验结果表明，QPSO-POWELL结合算法比QPSO配准所花费的计算时间更少。

3、配准精度高：表2中记录了4种方法的配准参数与实际值的均方误差，从中可以看出，QPSO算法参与的配准实验得到的参数的值与实际值最为接近，具有更高的求解精度。

表2配准后的参数

算法	Δθ	Δt<sub>x</sub>	Δt<sub>y</sub>
				Powell	0.1564	0.4837	0.2486
PSO	0.0325	0.1133	0.3755
				QPSO	0.0183	0.0115	0.1732
QPSO+Powell	0.0174	0.0167	0.1815

Δθ、Δt_x和Δt_y分别代表旋转角的误差(单位为角度)、X轴平移和Y轴平移(单位为像素)的误差

附图说明

图1为本发明的工艺流程简图。

图2为已有技术的PSO算法与本发明的QPSO算法求解配准参数的目标函数值的收敛曲线比较图。

图3为参考图像。

图4为本发明去除背景后的参考图像。

图5为本发明使用的浮动图像。

图6为本发明二级小波分解后的参考图像。

图7为本发明二级小波分解后的浮动图像。

具体实施方式：

为了更好的理解本发明的技术方案，以下对本发明的实施方式做进一步的介绍。

1、图像的背景去除

为了使图像免除噪声的干扰，需要将图像的背景去除；输入图像f，先求出图像中的最大和最小灰度，并令阈值初始为它们的平均值，然后根据阈值，将参考图像和待配图像分割成两部分，分别求出两部分的平均灰度值，由这两部分的平均灰度值求出新的阈值，最后使用种子填充算法去除参考图像和待配图像的背景部分。

具体步骤如下：

(1)、求出图像中的最大和最小灰度Z₁和Z_k，令阈值初始值为

$T_0=\frac{Z_1+Z_k}{2}---\left(1\right)$

(2)、根据阈值T_k将图像分割成R₁和R₂两部分，分别求出两部分的平均灰度值Z₀和Z_B：

$Z_0=\frac{\underset{z(i,j)<T_k}{\mathrm{\&Sigma;}}z(i,j)\&times;N(i,j)}{\underset{z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}---\left(2\right)$

$Z_B=\frac{\underset{z(i,j)>T_k}{\mathrm{\&Sigma;}}z(i,j)\&times;N(i,j)}{\underset{z(i,j)>T_k}{\mathrm{\&Sigma;}}N(i,j)}---\left(3\right)$

式中，z(i，j)是图像上(i，j)点的灰度值，N(i，j)是(i，j)点的权重系数，这里取N(i，j)＝1.0。

(3)、求出新的阈值

$T_{k+1}=\frac{Z_0+Z_B}{2}---\left(4\right)$

如果T_k＝T_k+1，则结束，否则k＝k+1，迭代执行上述步骤。

最后使用种子填充法从图像的左上角开始填充小于阈值的点，去除背景。

2、通过小波变换得到低分辨率图像

去除背景后的参考图像和待配图像进行小波变换，利用小波变换方法经一次或多次变换后得到两幅低分辨率的图像：为了降低图像的分辨率，需要将图像进行小波变换；图像通过它与低通滤波器的卷积形成一个平滑信号，与高通滤波器的卷积形成一个细节信号，从而把原图像分解为低分辨率的图像。实验中应用有平滑信号的那幅图像。

对于一维离散小波变换，在m级上的离散信号f_m，分别通过它与h(低通)滤波器的卷积以形成一个平滑信号f_m+1，与g(高通)滤波器的卷积以形成一个细节信号f′_m+1，从而把它分解为m+1级，这个过程可以用下列公式表示，使用了Mallat提出的金字塔算法：

$f_{m+1}\left(n\right)=\underset{k}{\mathrm{\&Sigma;}}h(2n-k)f_m\left(k\right)$

$f_{m+1}^{\&prime;}\left(n\right)=\underset{k}{\mathrm{\&Sigma;}}g(2n-k)f_m\left(k\right)$

其中，f_m+1是在分辨率级m+1上的平滑信号；f′_m+1是细节信号。在f_m上的离散点的总数等于f_m+1和f′_m+1的离散点的和，因此，f_m+1和f′_m+1两者必须在上式所描述的操作之后，在每隔一个数据点上采样。相同的过程进一步应用到f_m+1上，并在下一级分辨率创建细节的和平滑的信号，直到达到所期望的分辨率级别。

二维离散小波变换，它是对信号先在x方向上进行一维小波变换。然后再对2个子带在y方向上进行一维小波变换得到4个子带。对两次一维小波变换都是平滑信号的那幅图像再进行一次二维离散小波变换，得到原图像的二级小波分解图像。实验中用小波变换都是平滑信号的那幅图像来进行实验。

3、基于归一化互信息作为目标函数，使用QPSO算法求解一组配准参数

3.1归一化互信息的定义

互信息用来描述两个随机变量间的统计相关性，是一个变量包含另一个变量的信息量的多少的度量。它可用熵来描述：

I(A，B)＝H(A)+H(B)-H(A，B)    (5)

其中H(A)和H(B)分别为图像A和B的熵，H(A，B)为二者的联合熵。由于互信息对重叠区域的变化比较敏感，采用如下两种归一化互信息的计算方式：

$I(A,B)=\frac{H\left(A\right)+H\left(B\right)}{H(A,B)}---\left(6\right)$

$\mathrm{ECC}(A,B)=\frac{2I(A,B)}{H\left(A\right)+H\left(B\right)}---\left(7\right)$

归一化互信息能更好地反映配准函数的变化。

在这里将前一步得到的低分辨率图像作为归一化互信息的两个变量，采用归一化互信息作为QPSO算法求解的目标函数。

3.2使用QPSO算法完成一次配准

经过处理的两幅低分辨率的图像，以归一化互信息作为目标函数，即使用量子行为粒子群算法求解得到低分辨率图像配准参数；即首先在解空间初始化一组粒子，计算粒子的目标函数值，即归一化互信息值，然后粒子通过追寻个体最优位置与全局最优位置经过一定的迭代次数后完成寻优过程；

QPSO算法与其他进化类算法相类似，具有进化和群体智能的特点。在QPSO算法中，每个候选解称为‘粒子’，若干个候选解就构成了群体。每个粒子没有重量和体积，通过目标函数确定它的适应值。每个粒子在解空间中运动，粒子通过追随自身的个体极值与群体的极值来动态的调整自己的位置信息。

QPSO算法的描述如下：

假设算法的搜索空间为D维，粒子群的规模为N，每个粒子包含下列信息：

x_i＝(x_i1，x_i2，…x_iD)：粒子的当前位置；

P_i＝(P_i1，P_i2，…P_iD)：粒子i的当前最优位置，也可记为pbest；

P_g＝(P_g1，P_g2，…P_gD)：粒子群的全局最优位置，也可记为gbest。

每个粒子都按照以下的进化公式来更新自己的位置信息：

$\mathrm{mbest}\left(t\right)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{P}_i\left(t\right)=(\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{P}_{i1,}\left(t\right)...,\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{P}_{\mathrm{iD}}\left(t\right))---\left(8\right)$

p_id(t)＝φ·P_id(t)+(1-φ)·P_gd(t)，φ＝rand    (9)

x_id(t+1)＝p_id(t)±α·|mbest_d(t)-x_id(t)|·1n(1/u)，u＝rand    (10)

其中，t是当前的迭代次数，mbest称为平均最优位置，它是所有粒子自身最优位置的中心点；p_id为P_id与P_gd构成的超矩形中的一个随机点；参数α称为压缩-扩张因子，可以用来控制粒子的收敛速度，采用如下的取值方式，

α＝(1.0-0.5)×(MAXITER-t)/MAXITER+0.5    (11)

其中，t是当前迭代次数，MAXITER算法的最大迭代次数。

QPSO算法完成图像配准的步骤：

步骤1：初始化算法参数，包括粒子数，问题维数，初始化空间及搜索空间，粒子的初始位置，初始最优值等；

步骤2：由式(8)计算群体的平均最优位置mbest(t)；

步骤3：由式(9)计算随机位置p_id(t+1)；

步骤4：由式(10)计算粒子的新位置x_id(t+1)；

步骤5：由式(7)计算粒子新位置的适应度fitness(x_i(t+1))；

步骤6：更新粒子的当前最优位置，即：如果

fitness(x_i(t+1))＜fitness(pbest_i(t))，则pbest_i(t+1)＝x_i(t+1)，否则，pbest_i(t+1)＝pbest_i(t)；

步骤7：更新群体的最优位置，即：如果fitness(pbest_i(t+1))＜fitness(gbest(t))，则gbest(t+1)＝pbest_i(t+1)；

步骤8：循环步骤2～7，直至满足一定的结束条件，然后输出群体的全局最优位置gbest，即为配准参数。

3.3三线性PV算法实现插值

在配准过程中，由于浮动图像上的点通过空间变换后得到的点的坐标不一定是整数，需要通过插值方法来得到变换点的灰度值，浮动图f的样本a在某种空间变换下对应的参考图r的像素为b，通常b的空间坐标与任意一个实际的参考图像并不重叠。三线性PV插值算法不是通过邻居点确定b的灰度，而是按照周围8个像素与b点的空间距离分配权重，使周围8个像素点贡献于联合灰度分布统计，即

且

i＝1，2，…，8        (12)

其中，r(i)是8个邻点的灰度，W_i是权重。

这种插值方法可以使互信息的计算更为精确，对于优化过程中的局部极值问题有所缓解

3.4处理出界点

计算互信息时，对出界点进行修正处理，出界点是当浮动图中的某样本点经过一定的空间变换后的对应点落在参考图之外的点，互信息的计算必须考虑出界点；处理出界点的策略是使出界点的灰度等于距其最近的边界像素点的灰度，这样相当于扩大了参考图的背景，同时保持优化过程中的样本数不变，计算的互信息更为准确。

4、使用Powell方法求解精确的配准参数

经过量子行为粒子群算法す求得的配准参数作为初始值，以高分辨率图像作为对象，使用归一化互信息作为目标函数，将将低分辨率图像的配准参数作为Powell方法的输入，利用共轭方向并以此作为搜索方向，在一定迭代次数后输出最终的配准参数，从而完成图像配准。

即以第三步中求解结果gbest作为Powell的初始点，高分辨率的图像作为研究对象，即归一化互信息的参数为两幅高分辨率的图像，求解过程如下：步骤1：以gbest的三个参数作为变换t1(x，y，θ)的初始变量，取3个坐标轴(Z轴代表角度)单位向量p1，p2，p3作为初始方向，取门限极小值ε，沿着3个方向做三次一维搜索得到新的变换t2(x，y，θ)。

步骤2：求出新方向p4＝t2-t1，从t2沿着p4做一次一维搜索得到新的变换t3，如果||t3-t1||＜ε则停止，否则转步骤3。

步骤3：进行方向替换p1＝p2，p2＝p3，p3＝p4，t1＝t3，t1沿着p1，p2，p3方向做一维搜索得到新的变换t2，转步骤2。

其中一维搜索使用黄金分割法。通过以上步骤得到的更替方向都具有相互共轭的性质，保证了Powell算法寻优的准确性。

在本发明的实例验证中，采用的图像为图3，去除背景后图像为图4，配准时将图4作为参考图像。图4顺时针旋转30°，再向下、向右各平移5个像素，得到图5，图5作为浮动图像。小波变化后的低分辨率图像为图6和图7，经QPSO算法配准后的结果为(3.8743，3.6724，29.8278°)，经Powell方法求解的结果为(4.9853，5.01326，29.8625°)。

Claims (4)

1.一种基于量子行为粒子群算法的多分辨率医学图像配准方法，其特征是采用以下配准步骤：

(1)、图像的背景去除：首先求取一个阈值，之后采用种子填充算法去除参考图像和待配图像的背景部份，具体步骤如下：

(a)、求出图像中的最大和最小灰度z₁和z_k，令阈值初始值为

$T_0=\frac{Z_1+Z_k}{2}$

(b)、根据阈值T_k将图像分割成R₁和R₂两部分，分别求出两部分的平均灰度值z₀和Z_B：

$Z_0=\frac{\underset{z(i,j)<T_k}{\mathrm{\&Sigma;}}z(i,j)\&times;N(i,j)}{\underset{z(i,j)<T_k}{\mathrm{\&Sigma;}}N(i,j)}$

$Z_B=\frac{\underset{z(i,j)>T_k}{\mathrm{\&Sigma;}}z(i,j)\&times;N(i,j)}{\underset{z(i,j)>T_k}{\mathrm{\&Sigma;}}N(i,j)}$

式中，z(i，j)是图像上(i，j)点的灰度值，N(i，j)是(i，j)点的权重系数，这里取N(i，j)＝1.0；

(c)、求出新的阈值

$T_{k+1}=\frac{Z_0+Z_B}{2}$

如果T_k＝T_k+1，则结束，否则k＝k+1，迭代执行上述步骤；

最后使用种子填充法从图像的左上角开始填充小于阈值的点，去除背景；

(2)、通过小波变换得到低分辨率图像：小波变换得到两幅低分辨率的图像，去除背景后的参考图像和待配图像进行小波变换，利用小波变换方法经一次或多次变换后得到两幅低分辨率的图像；

(3)、利用量子行为粒子群算法求取配准参数：经过处理的两幅低分辨率图像，以归一化互信息作为目标函数，使用量子行为粒子群算法，即使用量子行为粒子群算法求解得到两幅低分辨率图像的配准参数；

(4)、利用Powell方法求取配准参数：经过求得的配准参数作为初始值，以高分辨率图像作为对象，使用归一化互信息作为目标函数，将低分辨率图像的配准参数作为Powell方法的输入，利用共轭方向并以此作为搜索方向，在一定迭代次数后输出最终的配准参数，从而完成图像配准。

2.根据权利要求1所述的基于量子行为粒子群算法的多分辨率医学图像配准方法，其特征在于所述使用量子行为粒子群算法求取配准参数：

(1)首先在解空间初始化一组粒子，计算粒子的目标函数值，即归一化互信息值，然后粒子通过量子行为粒子群算法的搜索，经过一定的迭代次数后完成寻优过程；

(2)迭代求解过程中，利用三线性PV算法实现插值；

(3)在计算互信息时，对出界点进行修正处理。

3.根据权利要求2所述的基于量子行为粒子群算法的多分辨率医学图像配准方法，其特征在于所述的利用三线性PV算法实现插值：在求解过程中，由于浮动图像上的点通过空间变换后得到的点坐标不是整数，需要通过插值来得到变换点的灰度值；三线性PV插值不会引入新的灰度值，浮动图像中的一点的灰度对联合直方图是由参考图像中的点的周围最近邻的8个点取与三线性插值相同的权重加权而得到。

4.根据权利要求2所述的基于量子行为粒子群算法的多分辨率医学图像配准方法，其特征在于所述的对出界点进行修正处理：计算互信息时，对出界点进行修正处理，出界点是当浮动图中的样本点经过一定的空间变换后的对应点落在参考图之外的点，互信息的计算必须考虑出界点；处理出界点时使出界点的灰度等于距其最近的边界像素点的灰度。

Images