发明内容
本发明的目的在于克服上述不足之处,从而提供一种基于量子行为粒子群算法的多分辨率医学图像配准方法,解决了基于互信息的目标函数存在许多局部极值从而给配准的优化过程带来的困难,大大地提高了配准精度和速度,达到了亚像素级,可以应用于临床诊断、放射治疗和图像引导的外科手术等领域;并且为提高配准过程的速度和精度以及鲁棒性,提出了使用多分辨率的方法来进行图像的配准。
按照本发明提供的技术方案,基于量子行为粒子群算法的多分辨率医学图像配准方法,特征是采用以下配准步骤:
(1)、图像的背景去除:首先求取一个阈值,之后采用种子填充算法去除参考图像和待配图像的背景部份:具体步骤如下:
(a)、求出图像中的最大和最小灰度z1和zk,令阈值初始值为
(b)、根据阈值Tk将图像分割成R1和R2两部分,分别求出两部分的平均灰度值z0和ZB:
式中,z(i,j)是图像上(i,j)点的灰度值,N(i,j)是(i,j)点的权重系数,这里取N(i,j)=1.0;
(c)、求出新的阈值
如果Tk=Tk+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 |
Δθ、Δtx和Δty分别代表旋转角的误差(单位为角度)、X轴平移和Y轴平移(单位为像素)的误差
具体实施方式:
为了更好的理解本发明的技术方案,以下对本发明的实施方式做进一步的介绍。
1、图像的背景去除
为了使图像免除噪声的干扰,需要将图像的背景去除;输入图像f,先求出图像中的最大和最小灰度,并令阈值初始为它们的平均值,然后根据阈值,将参考图像和待配图像分割成两部分,分别求出两部分的平均灰度值,由这两部分的平均灰度值求出新的阈值,最后使用种子填充算法去除参考图像和待配图像的背景部分。
具体步骤如下:
(1)、求出图像中的最大和最小灰度Z1和Zk,令阈值初始值为
(2)、根据阈值Tk将图像分割成R1和R2两部分,分别求出两部分的平均灰度值Z0和ZB:
式中,z(i,j)是图像上(i,j)点的灰度值,N(i,j)是(i,j)点的权重系数,这里取N(i,j)=1.0。
(3)、求出新的阈值
如果Tk=Tk+1,则结束,否则k=k+1,迭代执行上述步骤。
最后使用种子填充法从图像的左上角开始填充小于阈值的点,去除背景。
2、通过小波变换得到低分辨率图像
去除背景后的参考图像和待配图像进行小波变换,利用小波变换方法经一次或多次变换后得到两幅低分辨率的图像:为了降低图像的分辨率,需要将图像进行小波变换;图像通过它与低通滤波器的卷积形成一个平滑信号,与高通滤波器的卷积形成一个细节信号,从而把原图像分解为低分辨率的图像。实验中应用有平滑信号的那幅图像。
对于一维离散小波变换,在m级上的离散信号fm,分别通过它与h(低通)滤波器的卷积以形成一个平滑信号fm+1,与g(高通)滤波器的卷积以形成一个细节信号f′m+1,从而把它分解为m+1级,这个过程可以用下列公式表示,使用了Mallat提出的金字塔算法:
其中,fm+1是在分辨率级m+1上的平滑信号;f′m+1是细节信号。在fm上的离散点的总数等于fm+1和f′m+1的离散点的和,因此,fm+1和f′m+1两者必须在上式所描述的操作之后,在每隔一个数据点上采样。相同的过程进一步应用到fm+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)为二者的联合熵。由于互信息对重叠区域的变化比较敏感,采用如下两种归一化互信息的计算方式:
归一化互信息能更好地反映配准函数的变化。
在这里将前一步得到的低分辨率图像作为归一化互信息的两个变量,采用归一化互信息作为QPSO算法求解的目标函数。
3.2使用QPSO算法完成一次配准
经过处理的两幅低分辨率的图像,以归一化互信息作为目标函数,即使用量子行为粒子群算法求解得到低分辨率图像配准参数;即首先在解空间初始化一组粒子,计算粒子的目标函数值,即归一化互信息值,然后粒子通过追寻个体最优位置与全局最优位置经过一定的迭代次数后完成寻优过程;
QPSO算法与其他进化类算法相类似,具有进化和群体智能的特点。在QPSO算法中,每个候选解称为‘粒子’,若干个候选解就构成了群体。每个粒子没有重量和体积,通过目标函数确定它的适应值。每个粒子在解空间中运动,粒子通过追随自身的个体极值与群体的极值来动态的调整自己的位置信息。
QPSO算法的描述如下:
假设算法的搜索空间为D维,粒子群的规模为N,每个粒子包含下列信息:
xi=(xi1,xi2,…xiD):粒子的当前位置;
Pi=(Pi1,Pi2,…PiD):粒子i的当前最优位置,也可记为pbest;
Pg=(Pg1,Pg2,…PgD):粒子群的全局最优位置,也可记为gbest。
每个粒子都按照以下的进化公式来更新自己的位置信息:
pid(t)=φ·Pid(t)+(1-φ)·Pgd(t),φ=rand (9)
xid(t+1)=pid(t)±α·|mbestd(t)-xid(t)|·1n(1/u),u=rand (10)
其中,t是当前的迭代次数,mbest称为平均最优位置,它是所有粒子自身最优位置的中心点;pid为Pid与Pgd构成的超矩形中的一个随机点;参数α称为压缩-扩张因子,可以用来控制粒子的收敛速度,采用如下的取值方式,
α=(1.0-0.5)×(MAXITER-t)/MAXITER+0.5 (11)
其中,t是当前迭代次数,MAXITER算法的最大迭代次数。
QPSO算法完成图像配准的步骤:
步骤1:初始化算法参数,包括粒子数,问题维数,初始化空间及搜索空间,粒子的初始位置,初始最优值等;
步骤2:由式(8)计算群体的平均最优位置mbest(t);
步骤3:由式(9)计算随机位置pid(t+1);
步骤4:由式(10)计算粒子的新位置xid(t+1);
步骤5:由式(7)计算粒子新位置的适应度fitness(xi(t+1));
步骤6:更新粒子的当前最优位置,即:如果
fitness(xi(t+1))<fitness(pbesti(t)),则pbesti(t+1)=xi(t+1),否则,pbesti(t+1)=pbesti(t);
步骤7:更新群体的最优位置,即:如果fitness(pbesti(t+1))<fitness(gbest(t)),则gbest(t+1)=pbesti(t+1);
步骤8:循环步骤2~7,直至满足一定的结束条件,然后输出群体的全局最优位置gbest,即为配准参数。
3.3三线性PV算法实现插值
在配准过程中,由于浮动图像上的点通过空间变换后得到的点的坐标不一定是整数,需要通过插值方法来得到变换点的灰度值,浮动图f的样本a在某种空间变换下对应的参考图r的像素为b,通常b的空间坐标与任意一个实际的参考图像并不重叠。三线性PV插值算法不是通过邻居点确定b的灰度,而是按照周围8个像素与b点的空间距离分配权重,使周围8个像素点贡献于联合灰度分布统计,即
其中,r(i)是8个邻点的灰度,Wi是权重。
这种插值方法可以使互信息的计算更为精确,对于优化过程中的局部极值问题有所缓解
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°)。