CN101604444A - 用图像质量评估作为相似性测度的血管减影图像配准方法 - Google Patents

Abstract

本发明属于医学图像处理技术领域，具体为一种用图像质量评估(SSIM)作为相似性侧度的血管减影图像配准方法。本发明方法包括特征点提取、特征点对应关系计算和非线性变换计算等，在特征点对应关系计算中，采用SSIM测度作为相似性测度，非线性变换中采用薄板样条函数进行图像配准。本发明所提出的方法可以有效的对血管减影图像进行配准，提高血管减影的效果。而且，本发明提出的方法的配准的准确率明显好于传统的配准算法，例如使用互信息作为相似性测度的配准方法。

Description

用图像质量评估作为相似性测度的血管减影图像配准方法

技术领域

本发明属于医学图像处理技术领域，具体涉及一种有效的血管减影图像配准方法。

背景技术

数字血管减影造影术(Digital Subtraction Angiography，DSA)，是血管可视化的关键技术。数字减影血管造影技术中，未注射血管造影剂之前的X光片称为蒙片(mask image)，注入血管造影剂之后拍摄的X光片称为活片(live image)。数字减影血管造影术就是希望通过蒙片和活片的减法操作，能够得到一幅血管特征清晰的数字减影图像^[1]，[2]。由于呼吸、心跳、吞咽、镜头震颤等影响，病人和镜头间会有相对运动。因此在进行减法操作之前，需要对蒙片和活片进行配准。

图像配准是寻找一种空间变换，使在两幅或更多幅图像的对应点达到空间上的一致。一般将需要配准的图像中一幅图像作为基准，称为参考图像，其他图像则称为待配准图像或者浮动图像。

现有的血管减影图像配准的主要有以下的方法。文献[3]中提出了一种基于小波的多分辨率的血管减影图像配准方法，这种算法使用小波来进行多分辨率处理，并且使用互信息作为相似性测度。文献[2]提出了一种自动为数字减影图像配准提出特征点的方法。然而，由于两幅数字减影图像的变换往往是严重非线性的，现有的方法并无法有效地进行图像配准。

非线性图像配准分为以下的步骤。1，需要在图像中提取出一些特征点，作为配准的参考点。2，通过在另一幅图像中寻找与这些特征点最接近的点，也就是相似性测度最高的点。得到了两幅图像之间特征点的对应关系。3，通过这些对应关系，可以得到两幅图像之间的非线性变换关系。最后，利用这个非线性变换关系，可以将两幅图像配准。

发明内容

本发明的目的在于提出一种配准效果好的血管减影图像配准方法。

本发明提出的血管减影图像配准方法是一种使用结构性图像质量评估因子(SSIM)作为相似性测度的图像配准方法。

结构性图像质量评估因子(SSIM)在文献[4]中提出的一种用来评价图像质量的方法。由于人眼的视觉对图像的结构信息比如说图像的纹理和图像的边缘比较敏感。因此使用图像的结构上的信息的差异可以很好的从人眼视觉的角度反映一幅图像相对参考图像的质量。两幅图像L，R之间的SSIM定义为：

$S(x,y)=l(L,R).c(L,R).s(L,R)$

$=\left(\frac{{2\mathrm{\μ}}_L{\mathrm{\μ}}_R+C_1}{{\mathrm{\μ}}_L^2+{\mathrm{\μ}}_R^2+C_1}\right)\·\left(\frac{{2\mathrm{\δ}}_L{\mathrm{\δ}}_R+C_2}{{\mathrm{\δ}}_L^2+{\mathrm{\δ}}_R^2+C_2}\right)\·\left(\frac{{2\mathrm{\δ}}_{\mathrm{LR}}+C_3}{{\mathrm{\δ}}_L{\mathrm{\δ}}_R+C_3}\right)---\left(1\right)$

其中μ_L和μ_R分别为图像L和R的局部均值，δ_L和δ_R分别为图像L和R的局部方差的均值，而δ_LR是图像L和R去除了均值后的互相关。C₁，C₂和C₃是为了防止均值，方差或者互相关过小的时候造成的数值上的不稳定而添加的常数项。

由于SSIM对于非结构(比如图像的亮度)的失真不敏感，而对图像的结构失真相对敏感。其十分适合用于作为配准的相似性测度。因为在图像配准的时候希望相似性测度对图像亮度等非结构性信息不敏感。

本发明方法分为三步：特征点提取，特征点对应关系计算，非线性变换计算。

首先，进行特征点提取。由于血管减影的图像的大小一般都比较大，直接对每个像素进行配准从计算复杂度的角度来说并不现实。因此，只计算一些选定的特征点之间的对应关系。这些特征点采用手动选取，或者通过特征点选择器选取。这里，使用Harris检测器[5]来提取特征点。因为Harris检测器检测出的特征点相对比较稳定。而且，由于对所有的特征点进行配准的运算量比较大，仅仅对特征点的一部分进行配准。即在每个m×m像素的窗口中只采一个特征点。m的值越大，特征点的数量越小。通过实验，发现m的值可以取到20-100的整数，优选m为50～70。

然后，进行特征点对应关系计算。即对提取出的特征点在另一幅图像中寻找对应点。对于一个特征点，另一幅图像中的局部图像和这个特征点最相似的点被选择为这个特征点的对应点。也就是，对于每个特征点(x，y)，在另一幅图像中找到满足以下条件的对应点(x_p，y_p)：

$(x_p,y_p)=\arg {\max }_{x_p,y_p}S(I_l(x,y),I_m(x_p,y_p))---\left(2\right)$

其中，S(.，.)是一种相似性测度，在本发明中取SSIM测度，即如式(1)所示。I_m和I_l是在蒙片和活片上对应的点的局部图像区域。

由于蒙片和活片之间的形变一般不大，限制|x_p-x|和|y-y_p|分别不超过n。达到最大相似度的特征点的对应关系是通过穷举这个区域内的所有点的相似度，并且找出其中的相似度最大的点来达到的。通过实验，根据变形的程度n一般可以取10-50的整数，优选n＝25～35。

需要注意的是，在寻找特征点对应关系的时候，采用了SSIM测度作为相似性测度，这也是提出的配准方法与传统的配准算法的主要不同。

最后，根据得到的特征点的对应关系计算出两幅图片之间的非线性变换，并且通过这个非线性变换配准两幅图像。本发明使用薄板样条函数来实现配准[6]。得到非线性变换的过程如下：

根据得到的特征点的对应关系U＝{u_a：a＝1，2，...，n}和V＝{v_a：a＝1，2，...，n}。薄板样条函数可以表示成以下的形式：

$E_{\mathrm{TPS}\left(f\right)}=\underset{a=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|u_a-f\left(v_a\right)\left|\right|}^2+\mathrm{\λ}\∫\∫[\left(\frac{{\∂}^{2}f}{{\∂x}^2}\right)+2\left(\frac{{\∂}^{2}f}{\∂x\∂y}\right)+\left(\frac{{\∂}^{2}f}{{\∂y}^2}\right)]\mathrm{dxdy}---\left(3\right)$

其中f是在v_a和u_a之间的变换，λ是一个用来调节薄板样条函数平滑性的约束函数。λ取的比较大，变换的平滑性就越好。λ可以取到0.1-10。对于固定的λ，存在着函数f(v)可以表示成：

f(v)＝v.d+φ(v).w    (4)

其中v是需要变换的点集，d是一个3×3的仿射变换，w是一个n×3的非仿射变换矩阵。φ(v)是一个由薄板样条核生成的1×n向量。对于v中的每一个点，存在φ_a(v)，可以表示为φ_a(v)＝c‖v-v_a‖²log‖v-v_a‖，其中c是一个常量。

将(4)式代入(3)式，可以得到

E_TPS(d，w)＝‖U-Vd-φw‖²+λtrace(w^TΦw)         (5)

其中U和V分别是由u_a和v_a构成的矩阵。Φ是一个由φ(v_a)构成的n×n的矩阵。

然后通过QR分解可以分离出仿射和非仿射的变换：

$V=\left[Q_1{Q}_2\right]\left(\begin{array}{c}R\\ 0\end{array}\right)---\left(6\right)$

其中Q₁和Q₂是正交矩阵。R是一个上三角阵。最后，w和d的解可以表示为：

$w=Q_2{(Q_2^T\mathrm{\Φ}{Q}_2+{\mathrm{\λI}}_{N-3})}^{-1}{Q}_2^{T}U---\left(7\right)$

$d=R^{-1}{Q}_1^T(U-\mathrm{\Φw})---\left(8\right)$

由此得非线性变换函数f(v)，然后函数对图像进行非线性变换，即可以得到配准后的图像。

附图说明

图1为原图像(左)和经过变换后的图像(右)。

图2为利用Harris提取的特征点(左)和计算得到的特征点位移结果(右)。

图3第一行从左到右分别为未配准的图像、使用互信息配准后的图像和使用SSIM配准后的图像。第二行中是第一行中的图像和原有的图像相减后得到的图像经过归一化后的结果。

图4为配准前的减影结果(左)和使用本文中的方法配准后的减影结果(右)。

具体实施方式

为了测试本文中的配准方法的效果，将一幅活片图像手动的经过非线性变换，得到需要配准图像。首先，在活片图像中提取特征点(x_i，y_i)，其中i＝1，2，...，K，K是特征点的数量。然后将每个特征点在x和y方向上分别平移dx，dy个像素，其中dx，dy在-δ和δ范围内。δ在本实验取20。经过平移后的特征点的坐标为(x_i+dx_i，y_i+dy_i)。对于得到的特征点，通过薄板样条对活片图像进行变换，得到待配准的图像。原图像和变换后的图像分别如图1的左、右图所示。

使用本发明中的方法和使用互信息作为相似性测度的方法对两幅图像进行配准。首先进行特征点提取，在特征值点提取的时候，取m为60，利用Harris提取到特征值点的结果如图2左图所示。以SSIM作为相似性测度得出的特征点的位移结果如图2右图所示。然后，分别通过互信息和SSIM作为相似性测度对两幅图像进行配准。得到的结果如图3所示。图3的第一行从左到右分别为待配准的图像、通过互信息配准后的图像和通过本文中的方法配准后的图像；第二行中是第一行中的图像和原有的图像相减后得到的图像经过归一化后的结果。可以看到，未配准得图像的减影后的图像的伪影最多，而通过本发明方法配准后的结果明显比使用互信息配准后的图像的伪影少。

另外，通过计算配准后的图像和原有活片之间的均方误差来评价配准的性能。使用6组图片来进行实验。实验的结果如下表所示。可以看出，使用本发明中的方法配准的结果明显比使用SSIM方法配准后的方法均方误差小。

图像编号	1	2	3	4	5	6	7	8
									使用SSIM配准后的均方误差	2.3037	2.3037	2.1117	2.4365	3.1938	3.2432	2.7182	2.5329
使用互信息配准后的均方误差	3.9968	4.3434	3.0851	3.2851	4.8936	4.7780	4.2132	4.3354

表1.使用SSIM和互信息作为相似性测度的配准图像均方误差结果

然后，把本发明方法应用于实际的血管减影图像中。配准前和配准后的剪影结果如图4的左右图所示。因此，本文中的方法可以有效的配准血管减影图像。而且配准的结果比现有的方法性能好。

参考文献

[1]E.H.W.Meijering，W.J.Niessen，and M.A.Viegever，“Retrospective motioncorrection in digital subtraction angiography：a review，”IEEE Transactions on MedicalImaging，vol.18，pp.2-21，1999.

[2]Y.Bentoutou and N.Taleb，“Automatic extraction of control points for digitalsubtraction angiography image enhancement，”IEEE Transactions on Nuclear Science，vol.52，pp.238-246，2005.

[3]J.Yang，Y.Wang，S.Tang，S.Zhou，Y.Liu，and W.Chen，“Multiresolution ElasticRegistration of X-Ray Angiography Images Using Thin-Plate Spline，”IEEE Transactions onNuclear Science，vol.54，pp.152-166，2007.

[4]Z.Wang，A.C.Bovik，H.R.Sheikh，and E.P.Simoncelli，“Image quality assessment：From error visibility to structural similarity，”IEEE transactions on image processing，vol.13，pp.600-612，2004.

[5]C.Harris and M.Stephens，“A combined corner and edge detector，”Alvey visionconference，1988，p.50.

[6]F.L.Bookstein，“Principal warps：Thin-plate splines and the decomposition ofdeformations，”IEEE Transactions on pattern analysis and machine intelligence，vol.11，pp.567-585，1989.

Claims (1)

1.一种用图像质量评估作为相似性侧度的的血管减影图像配准方法，其特征在于具体步骤如下：

首先，进行特征点提取；特征点采用手动选取，或者通过特征点选择器选取；在每个m×m像素的窗口中只采一个特征点，m为20-100的整数；

然后，进行特征点对应关系计算；即对提取出的特征点在另一幅图像中寻找对应点；对于一个特征点，另一幅图像中的局部图像和这个特征点最相似的点被选择为这个特征点的对应点；也就是，对于每个特征点(x，y)，在另一幅图像中找到满足以下条件的对应点(x_p，y_p)：

$(x_p,y_p)={\arg \max }_{x_p,y_p}S(I_l(x,y),I_m(x_p,y_p))---\left(1\right)$

其中，S(.，.)是SSIM测度，I_m和I_l是在蒙片和活片上对应的点的局部图像区域；

限制|x_p-x|和|y-y_p|分别不超过n；达到最大相似度的特征点的对应关系是通过穷举这个区域内的所有点的相似度，并且找出其中的相似度最大的点来达到的；n取10-50的整数；

最后，根据得到的特征点的对应关系计算出两幅图片之间的非线性变换，并且通过这个非线性变换配准两幅图像；其中，使用薄板样条函数来实现配准；得到非线性变换的过程如下：

根据得到的特征点的对应关系U＝{u_a：a＝1，2，...，n}和V＝{v_a：a＝1，2，...，n}；薄板样条函数表示成以下的形式：

$E_{\mathrm{TPS}\left(f\right)}=\underset{a=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|u_a-f\left(v_a\right)\left|\right|}^2+\mathrm{\λ}\∫\∫[\left(\frac{{\∂}^{2}f}{{\∂x}^2}\right)+2(\frac{{\∂}^{2}f}{\∂x\∂y})+(\frac{{\∂}^{2}f}{{\∂y}^2}\left)\right]\mathrm{dxdy}---\left(2\right)$

其中f是在v_a和u_a之间的变换，λ是一个用来调节薄板样条函数平滑性的约束函数；λ取值为0.1-10；对于固定的λ，存在着函数f(v)可以表示成：

f(v)＝v.d+φ(v).w                    (3)

其中v是需要变换的点集，d是一个3×3的仿射变换，w是一个n×3的非仿射变换矩阵；φ(v)是一个由薄板样条核生成的1×n向量；对于v中的每一个点，存在φ_a(v)，表示为φ_a(v)＝c||v-v_a||²log||v-v_a||，其中c是一个常量；

将(3)式代入(2)式，可以得到

E_TPS(d，w)＝‖U-Vd-φw||²+λtrace(w^TΦw)         (4)

其中U和V分别是由u_a和v_a构成的矩阵；Φ是一个由φ(v_a)构成的n×n的矩阵；

然后通过QR分解分离出仿射和非仿射的变换：

$V=\left[Q_1{Q}_2\right]\left(\begin{array}{c}R\\ 0\end{array}\right)---\left(5\right)$

其中Q₁和Q₂是正交矩阵；R是一个上三角阵；最后，w和d的解表示为：

$w=Q_2{(Q_2^T\mathrm{\Φ}{Q}_2+\mathrm{\λ}{I}_{N-3})}^{-1}{Q}_2^{T}U---\left(6\right)$

$d=R^{-1}{Q}_1^T(U-\mathrm{\Φw})---\left(7\right)$

由此得非线性变换函数f(v)，然后函数对图像进行非线性变换，即得到配准后的图像。

Images