CN101303767B - 基于块图像内容自适应分类的数字剪影图像配准方法 - Google Patents

Abstract

本发明属于医学图像处理技术领域，具体为一种块图像内容自适应分类的数字剪影图像(DSA图像)配准方法。数字减影血管造影术，是通过含有血管信息的活片和不含血管信息的蒙片之间的减法操作，得到一幅血管特征清晰的数字减影图像。本发明根据活片的血流特性，将原图像进行分块，并对不同内容的块图像采用不同的相似性测度：针对含血管 的块图像，采用互信息作为相似性测度，可有效地减小血管成分对配准的干扰；而对不含血管的块图像，采用模糊不变相位相关法作为相似性测度，能有效减小图像模糊对配准的干扰。本发明以高维互信息作为块图像分类的判据，对图像灰度变化有较强的鲁棒性。分析和实验表明本发明可有效地提高DSA图像整体配准的配准精度。

Description

基于块图像内容自适应分类的数字剪影图像配准方法

技术领域

本发明属于医学图像处理技术领域，具体为一种块图像内容自适应分类的数字剪影图像(DSA图像)配准方法。 

背景技术

数字减影血管造影术(digital subtraction angiography，DSA)，是血管可视化的关键技术。研究表明，血管的可视化有助于在临床诊断中用于研究病变血管的结构以及制定相关的治疗措施。在血管的可视化技术中，未注射血管造影剂之前的X光片称为蒙片，注入血管造影剂之后拍摄的X光片称为活片或者盈片。数字减影血管造影术就是希望通过蒙片和活片的减法操作，能够得到一幅血管特征清晰的数字减影图像。 

由于呼吸、心跳、吞咽、肌肉收缩、镜头震颤等影响，病人和镜头间会有相对运动.如果直接对蒙片和活片进行减法操作，这种相对运动会使得减影图像中产生伪影及可能的图像模糊。因此通常在进行减法操作之前，需要对蒙片和活片图像进行配准。研究表明：在临床应用中，DSA图像配准过程中主要碰到两种干扰，即：第一，活片中的血管存在，导致蒙片和活片图像对应位置的拓扑结构发生了变化，这样使得应用相似性测度对它们之间进行配准的鲁棒性下降；第二，病人和镜头间的相对运动会造成运动伪影和图像模糊。此外，由于DSA图像的变形是非刚性的，这些都给配准带来了极大困难和挑战。 

目前实用性较强的DSA配准方法^[3]是：首先，在活片中提取控制点，这些控制点一般按照一定密度取自图像的边缘点或者框架点。然后，根据某一相似性测度，在蒙片中对应位置的临域内，搜索及判定相关的控制点。这样蒙片和活片中的每对控制点都关联了原图像相同部分的某些区域，称该区域为子块图像或者块图像。根据活片图像的控制点在蒙片中搜寻相关控制点的过程，实际上就是对两幅对应块图像进行配准的过程。虽然DSA图像整体变形为非刚性的，但如果块图像的区域比较小，那么该块图像活片相对于蒙片的变形可近似为仿射变换，这意味着蒙片和活片之间的块图像配准可采用成熟的刚性配准技术进行配准。最后，拟合活片与蒙片图像在控制点位置之间对应关系的非线性变换，将该变换作用于蒙片，得到矫正后的蒙片，通过减法操作，得到清晰的减影图像。 

发明内容

本发明的目的在于提供一种块图像内容自适应分类的数字剪影图像配准方法，以提高配准结果的整体精度和鲁棒性。 

不同相似性测度对不同内容的DSA块图像的配准，具有不同的配准准确度和鲁棒性。如果根据块图像内容对待配块图进行分类，然后对不同类型的块图像采取不同的相似性测度，可以期待去获得一个整体精准度和鲁棒性较高的配准结果。 

如果对不同内容的块图像采取不同的相似性测度进行配准，那么首先要确定块图像的内容。我们知道发现该块图像是否有血管，主要的依据是人体的血流特性。文献[3]中的实验表明，活片块图像中的血管会随着心脏跳动呈现周期性强烈变化的灰度，且出现与心跳脉冲信号同步且类似的峰值，而普通的伪影则变化缓慢，此性质是本方法对块图像进行内容分类的主要依据。 

医院实地采集的活片图像，并不是只有一幅静态图像，而是连续拍摄的一段录影。在以往DSA配准的研究中，往往只取整个录影中的某一帧图像作为活片。这样意味着并未充分利用活片序列的全部信息，因为活片中的血管部分由于心跳会有周期性的灰度变化，而这些信息只有在时间序列上才能体现出来。根据心脏跳动的周期和X光机的采片率，可把连续的10-12帧图像作为一个单元，此单元则表征一次血流变化的周期性质。 

本发明提出的数字剪影图像的配准方法，其中块图像配准的步骤如下：根据血流特性，以高维互信息作为判据对块图像分为含血管块图像和不含血管块图像，然后针对含血管的块图像和不含血管的块图像采取不同的相似性测度：对于含血管的块图像，采用高维互信息作为相似性测度；对于不含血管的块图像，采用模糊不变相位相关法作为相似性测度。 

辨识块图像内容最简单的方法是根据所取图像整体灰度变化的周期性来分辨是否含有血管，但是由于除血管外仍然有不等量的背景信息存在，很难排除背景灰度变化对块图像整体灰度的影响。由于血管部分灰度变化比较大，那么一个单元内的两幅有血管图像的互信息值应该明显小于两幅无血管图像，基于这种思想，本文采取高维互信息的方法进行对是否有无血管的块图像进行划分。 

文献[4]中引入了高维互信息的定义，其定义互信息矩阵如下： 

$I_{\mathrm{matrix}}=\left[\begin{array}{cccc}{I}_{11}& I_{12}& \·\·\·& I_{1N}\\ I_{21}& \·\·\·& \·\·\·& I_{2N}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ I_{N1}& I_{N2}& \·\·\·& I_{\mathrm{NN}}\end{array}\right]---\left(1\right)$

式中I_ij，i，j＝1，2，...，N，表示第i幅图和第j幅图之间的归一化的互信息，N为一个单元内的图片数目。设λ_i(i＝1，2，...，N)为互信息矩阵的特征值。则对应的高维互信息I为： 

$I({\mathrm{\λ}}_1,{\mathrm{\λ}}_2,...,{\mathrm{\λ}}_N)=1+\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}\lg \left(\frac{{\mathrm{\λ}}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}\right)}{\mathrm{lgN}}---\left(2\right)$

根据其定义可见，高维互信息包含了一个单元内任意两幅块图像之间的信息，对图像间的整体灰度变化有较高的鲁棒性。有血管的图像单元，求得的高维互信息值比较小，反之较大。设定阈值，即可区分块图像中有无血管。 

具体步骤为： 

(1)在造影图像序列中，确定一幅图像为活片；以活片为基准，向前或者向后选取连续的N幅图为一个活片单元，N可取10到12之间； 

(2)在块图像配准过程中，取该单元其他图片中其相应对应位置的块图像，构成N幅块图像组成的块图像单元； 

(3)计算该块图像单元N幅块图像的高维互信息I； 

(4)设定阈值φ，如果I≥φ，则判定该块图像为无血管块图像，如果I＜φ，则判定该块图像为含血管块图像。 

我们的大量实验表明，阈值φ的经验值为0.4-0.5，可根据实际情况分析选取，最优值为0.46。 

在对块图像内容进行分类之后，对不同类型的块图像采用最适应的相似性测度，从而获得对图像整体配准的高鲁棒性。 

在块图像配准过程中，针对含血管块图像，互信息作为相似性测度可有效地减小血管成分对配准的干扰，优于其他相似性测度。而对血管成分很少的块图像，血管对配准的影响已经很小，为了抵抗模糊干扰，采用模糊不变相位相关法，能够取得更好的效果。 

完整的基于模板匹配法的DSA配准流程主要分为三步：一、提取控制点；二、根据相似性测度，寻找相关控制点，这个过程实际上是块图像配准过程；三、根据控制点对拟合非线性变换，将该变换作用于整幅图像。本方法关注的重点是在提取活片的控制点之后，如何更精确和鲁棒地在蒙片中搜寻相关的控制点，进行块图像匹配，提高整体算法的鲁棒性。本方法结合边缘控制点提取和薄板样条非线性变换，构成一完整流程，对实际图像进行了大量实验，获得良好效果。 

附图说明

图1中，(a)模拟蒙片图像，(b)模拟活片图像，(c)矫正后的模拟蒙片图片，(d)未配准时的DSA图片，(e)本文方法的配准结果。 

图2中，(a)蒙片图像，(b)活片图像，(c)未配准前的DSA图像，(d)配准后的DSA 图像(上排为第一组实验结果，下排为第二组实验结果)。 

具体实施方式

一、模拟造影图像的配准实验 

考虑到真实DSA图像的配准结果难以进行参数化的评估，视觉上的观察比较也不非常可靠，这是多模态医学图像配准的一个通病和未解决的难题。为了对本发明的优越性进行验证，我们先对模拟造影图像进行配准实验，以便与其他方法进行数值比较。 

在模拟造影图像配准实验中，任意抽取一幅医学造影图像作为蒙片，对其进行特定非线性变换获得的图像作为活片。因为该非线性变换的相关参数已知，我们通过比较真实值和本文方法的配准结果之间的偏差，来测试本文算法的性能。具体实验步骤如下： 

(1)在DSA图像序列中抽取一幅N×N的图像，将其作为蒙片； 

(2)在蒙片中随机选取一组K个控制点，对每个控制点的坐标各自随机加上[-d，+d]的一个偏差值，获得一组新的控制点坐标； 

(3)根据两组控制点，拟合薄板样条变换，将该变换作用于蒙片生成一幅新的图片，并且通过Matlab函数进行一定程度的模糊化，作为本实验的活片； 

(4)获得实验图片以后，按照前面所述的流程对活片和蒙片进行配准。 

本实验中采用的图片大小为512×512，控制点偏移最大值d＝5，随机变换的控制点数目K＝25。配准过程中，每个控制点的搜索范围是以该点为中心的一个10×10的区域，与控制点对应的块图像大小定为60×60。 

为了证明本方法的优越性，我们与对块图像内容不分类而仅仅采用互信息的单相似性测度的匹配方案进行了比较。因为本项实验中我们已知两组控制点的真实偏差，故可以用控制点之间的均方根误差(RMSE)来比较配准精度^[3]。均方根误差(RMSE)的定义如式(3)： 

$\mathrm{RMSE}=$

$\sqrt{\frac{\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{(a_{11}{x}_i+a_{12}{y}_i+d_x+\underset{j=1}{\overset{K}{\mathrm{\Σ}}}{W}_{\mathrm{xj}}\mathrm{\φ}\left({\mathrm{\λ}}_{\mathrm{ij}}\right)-x_i^{\′})}^2+\underset{i=1}{\overset{K}{\mathrm{\Σ}}}{(a_{21}{x}_i+a_{22}{y}_i+d_y+\underset{j=1}{\overset{K}{\mathrm{\Σ}}}{W}_{\mathrm{yj}}\mathrm{\φ}\left({\mathrm{\λ}}_{\mathrm{ij}}\right)-y_i^{\′})}^2}{K}}---\left(3\right)$

其中，a₁₁，a₂₁，a₁₂，a₂₂，d_x，d_y为仿射变换参数，φ(λ)＝λ²logλ²，W为非线性加权系数，K为控制点集的势即控制点数目。(x，y)为蒙片控制点的坐标，(x′，y′)为活片控制点的坐标。 

实验结果如附图1所示。配准的RMSE结果如表1所示。可见，采用基于块图像内容的多相似性测度的配准方案，可有效地提高控制点的配准精度。 

表1 

	单相似性测度方  案	本文方法
			RMSE(pixel)	0.3951	0.2504

二、真实DSA图像的配准实验 

附图2给出了真实DSA图像的配准结果。 

实验参数如下：第一组实验图片大小为512×512，控制点个数为64，块图像大小定为50×50，每个控制点的搜索范围是以该点为中心的一个10×10的区域，第二组实验图片大小为360×360，控制点个数为81，其他参数同第一组实验。(虽然第一组图片较大，但图片局部变形较小，故采用较少的控制点也可取得较好的效果，并减小了计算量)。 

实验结果表明，经过配准之后，伪影得到明显消除，同时血管部分也更加清晰，图像更加平滑。视觉观察中并未发现有信息缺失，配准后的DSA图像中血管比较清晰明显，运动伪影得到了较好的消除，图像噪声也得到了较好的抑制。 

参考文献： 

[1]Erik H.W.Meijering，Retrospective motion correction in digital subtraction angiography：areview，IEEE Transactions on Medical Imaging，January 1999，Vol.18，No.1：2-21 

[2]Ville Ojansivu and Janne Heikkil

，Image registration using blur-invariant phase correlation，IEEE Signal Processing Letters，July 2007，Vol.14，No.7：449-452 

[3]Youcef Bentoutou and N.Taleb，Automatic extraction of control points for digital subtractionangiography image enhancement，IEEE Transactions on Nuclear Science，February 2005，Vol.52，No.1：238-246 

[4]Bo Wang and Yi Shen，A method on calculating high-dimensional mutual information and itsapplication to registration of multiple ultrasound images，Ultrasonics(2006)，doi：10.1016/j.ultras.2006.07.012 

[5]F.L.Bookstein，Principal warps：thin-plate splines and the decomposition of deformations，IEEE Transactions on Pattern Anal.Mach.Intell.，1989，Vol.11，No.6：567-585 

[6]Zhiguo Cao，et al，DSA image registration based on multiscale gabor filters and mutualInformation，2005，Proceedings of the 2005 IEEE，International Conference on InformationAcquisition：105-110。 

Claims (2)

1.一种基于块图像内容自适应分类的数字剪影图像配准方法，包括提取控制点、图像块配准和根据控制点对拟合非线性变换，将该项变换作用于整幅图像三个步骤，其特征在于其中图像块配准的步骤如下：根据血流特性，以高维互信息作为判据对块图像分为含血管块图像和不含血管块图像，然后针对含血管的块图像和不含血管的块图像采取不同的相似性测度：对于含血管的块图像，采用高维互信息作为相似性测度；对于不含血管的块图像，采用模糊不变相位相关法作为相似性测度。

2.根据权利要求1所述的基于块图像内容自适应分类的数字剪影图像配准方法，其特征在于所述以高维互信息作为判据，将块图像分为含血管块图像和不含血管块图像的步骤如下：

(1)在造影图像序列中，确定一幅图像为活片；以活片为基准，向前或者向后选取连续的N幅图为一个活片单元，N取10到12之间；

(2)在块图像配准过程中，取该单元其他图片中其相应对应位置的块图像，构成N幅块图像组成的块图像单元；

(3)计算该块图像单元N幅块图像的高维互信息I；

(4)设定阈值φ，如果I≥φ，则判定该块图像为无血管块图像，如果I＜φ，则判定该块图像为含血管块图像；

其中，高维互信息I的计算式为：

设λ_i(i＝1，2，...，N)为互信息矩阵的特征值，则对应的高维互信息为：

$I({\mathrm{\λ}}_1,{\mathrm{\λ}}_2,...,{\mathrm{\λ}}_N)=1+\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\mathrm{\λ}}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}\lg \left(\frac{{\mathrm{\λ}}_i}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}\right)}{\mathrm{lgN}}---\left(2\right)$

这里λ_i，i＝1，2，...，N，下述互信息矩阵的特征值：

$I_{\mathrm{matrix}}=\left[\begin{array}{cccc}{I}_{11}& I_{12}& ...& I_{1N}\\ I_{21}& ...& ...& I_{2N}\\ ...& ...& ...& ...\\ I_{N1}& I_{N2}& ...& I_{\mathrm{NN}}\end{array}\right]---\left(1\right)$

式中I_ij，i，j＝1，2，...，N，表示第i幅图和第j幅图之间的归一化的互信息，N为一个单元内的图片数目。

Images