说明书自适应增强数字减影血管造影图像的方法
技术领域
本发明属于医学成像技术领域，具体涉及一种自适应增强数字减影血管造影图像的方法。
背景技术
数字减影血管造影DSA(Digital Subtraction Angiography)技术在临床已应用20多年，是心脑血管疾病无创诊断与介入治疗手术导航的重要依据。DSA图像处理中的一个关键任务就是进行图像增强，以使血管的生理特征能够更清楚的显示出来，方便结构分析、运动分析、三维可视化等后续操作，以及图像引导手术、肿瘤放射治疗、治疗评估等应用研究。由于X射线经过的组织厚度以及血液中的造影剂浓度不均匀，经过减影后的图像不能完全消除人体组织引起的噪声信号，背景和目标混杂在一起。另外，人体的组织结构和形状很复杂，而且人与人之间有相当大的差异。因此血管的增强是一项困难的任务。为了增强血管，抑制背景，目前普遍采用各向异性平滑对图像进行增强。Perona和Malik(参见P.Perona and J.Malik：‘Scale-space and edge detectionusing anisotropic diffusion’.IEEE Trans.Pattern Anal.MachineIntell.，12(7)，629-639(1990))引入基本的各向异性扩散方程为图像进行平滑处理。在图像噪声不大，并且选择较好参数的情况下，利用各向异性平滑能够取得不错的效果。但是当噪声较大时，该方法存在噪声处理缺陷。因此Catte(参见F.Catte，P.L.Lions，J.M.Morel，and T.Coll：‘Imageselective smoothing and edge detection by nonlinear diffusion’.SIAMJ.Numer.Anal.，29，182-193(1992))对其引入高斯平滑处理，使得各向异性平滑在应用时减少了噪声的影响。高斯平滑处理的采用解决了噪声的问题，但同时也引入了新的问题：(1)在平滑噪声的同时平滑了边缘。(2)每一步迭代都要进行一次高斯平滑，计算速度非常缓慢。(3)随着迭代的进行，噪声越来越少，在每一步迭代都需要重新估计高斯函数的方差。
另外传统的各向异性增强方法没有确定如何选取传导参数，大部分情况下都是通过实验来确定。通常情况下，传导参数过大会导致在不想执行平滑的地方产生平滑，传导参数过小则易导致在需要平滑的区域产生增强，而且过小的传导参数由于扩散作用不大而导致需要很多的迭代次数。不同的图像要想找到最佳传导参数，并不是一件容易的事情，因此有必要使得传导参数能够自适应地变化。因此寻找一种自适应的血管图像增强方法一直是人们研究的一个重要方向。
发明内容
本发明的目的在于提供一种在数字减影血管造影图像中增强血管数据的方法，该方法在整个迭代过程中避免了人为选择传导参数，即传导参数能够自适应地变化，从而能够自适应的对DSA图像进行增强，并降低了噪声的影响，同时增强了边缘，提高了图像增强的效果。
本发明提出的在数字减影血管造影图像中增强血管数据的方法按照下述步骤进行：
(1)设定迭代次数；
(2)将原始数字减影血管造影图像每个像素点(i，j)用小面模型进行拟和，并计算像素点(i，j)在拟合后的一阶偏导数和二阶偏导数
(3计算每个像素点(i，j)在拟合后的Hessian矩阵，用矩阵H(i，j)表示：
H ( i , j ) = I ^ ( i , j ) xx I ^ ( i , j ) xy I ^ ( i , j ) xy I ^ ( i , j ) yy ]]>
并计算矩阵H(i，j)的两个特征值λ(i，j)1，λ(i，j)2；
(4)针对每个像素点(i，j)选取传导参数T(i，j)的形式：
T ( i , j ) = ( λ ( i , j ) 1 2 + λ ( i , j ) 2 2 ) 1 2 ]]>
(5)针对每个像素点(i，j)选取传递函数的形式：
c ( | ▿ I ^ | ) ( i , j ) = 1 1 + ( | ▿ I ^ ( i , j ) | T ( i , j ) ) 2 ]]>
其中，
| ▿ I ^ | ( i , j ) = I ^ ( i , j ) x 2 + I ^ ( i , j ) y 2 ]]>
(6)依据各个像素点的传递函数采用各向异性扩散模型更新各个像素点灰度，从而更新原始图像灰度；
(7)重复步骤(2)至步骤(6)直到达到设定的迭代次数。
综合考虑计算效率和计算精度，所述步骤(2)选用5×5的核函数对原始数字减影血管造影图像进行拟和。
本发明利用小面模型对图像进行拟和来达到降噪的效果，并根据图像灰度的二阶导数计算Hessian矩阵并计算传导参数，这样在整个迭代过程中不需要人为选择传导参数，使得本发明自适应的对DSA图像进行增强，降低了噪声的影响，并增强了边缘，提高了增强的效果。
附图说明
图1为本发明流程图；
图2示出了本发明与传统各向异性扩散方法分别处理弱噪声DSA图像结果的对比图，其中
图2a是从原始的弱噪声DSA图像选取的一部分(像素为200×200)，
图2b是采用Catte方法得到的图像增强效果图(传导参数T＝30，50次迭代，高斯方差取为0.3)，
图2c是采用本发明迭代50次得到的图像增强效果图，
图2d是采用Catte方法得到的图像增强效果图(传导参数T＝30，150次迭代，高斯方差取为0.3)，
图2e是采用本发明迭代150次得到的图像增强效果图；
图3示出了本方法与传统各向异性扩散方法分别处理强噪声DSA图像结果的对比图，其中
图3a是从原始的强噪声DSA图像选取的一部分(像素为200×200)，
图3b是采用传统各向异性扩散Perona-Malik方法得到的图像增强效果图(传导参数T＝30，50次迭代)，
图3c和图3d是采用Catte方法得到的图像增强效果图(传导参数T＝30，50次迭代，高斯方差分别取为0.3和0.4)，
图3e是利用本发明迭代50次得到的图像增强效果图，
图3f是采用Catte方法得到的图像增强效果图(传导参数T＝30，150次迭代，高斯方差取为0.3)，
图3g是采用本发明迭代150次得到的图像增强效果图。
具体实施方式
本发明包括以下步骤：
(1)设定迭代次数，一般选择30到50次迭代。
(2)将原始DSA图像每个像素点用小面模型进行拟和。
将图像看作曲面，根据小面模型，原始图像中每一像素点的某邻域灰度曲面可以用离散正交多项式作曲面最佳拟合。综合考虑计算效率和计算精度，本发明选用5×5的核函数对原始图像进行拟和。例如，在横坐标方向设定区间R＝{-2 -1 0 1 2}，纵坐标方向设定区间C＝{-2 -1 0 1 2}，R和C为原始图像中某一像素点(i，j)具有元素数为5的对称邻域R×C内的坐标序数集，令r和c分别代表区间R和C中的值，小面模型的离散正交多项式的基底可以表示为：
1 , r , c , r 2 - 2 , rc , c 2 - 2 , r 3 - 17 5 r , ( r 2 - 2 ) c , r ( c 2 - 2 ) , c 3 - 17 5 c ]]>
因此，利用离散正交多项式构造的双变量立方函数f(r，c)(i，j)来估计对应邻域内点的灰度值为：
f ( r , c ) ( i , j ) = K 1 + K 2 r + K 3 c + K 4 ( r 2 - 2 ) + K 5 rc + K 6 ( c 2 - 2 ) + K 7 ( r 3 - 17 5 r ]]>
+ K 8 ( r 2 - 2 ) c + K 9 r ( c 2 - 2 ) + K 10 ( c 3 - 17 5 c ) ]]>
其中K1，K2，…Km…K10。为双变量离散多项式的系数，由下式表示：
K m = Σ ( r , c ) ∈ S g m ( r , c ) I ( r , c ) ( i , j ) Σ ( r , c ) ∈ S g m 2 ( r , c ) ]]>
S是定义在R×C二维对称邻域，I(r，c)(i，j)是在某一像素点(i，j)的对称邻域内位置(r，c)处((r，c)∈S)的图像灰度值，{g0(r，c)，g1(r，c)，…，g10(r，c))代表二维离散正交多项式的基底，其值分别为 1 , r , c , r 2 - 2 , rc , c 2 - 2 , r 3 - 17 5 r , ( r 2 - 2 ) c , r ( c 2 - 2 ) , c 3 - 17 5 c . ]]>
小面拟和后的每一个像素点的灰度值可以表示为在拟合曲面的中心点的函数值，在此例子中，即取r＝0，c＝0时，f(0，0)(i，j)＝K1-2*K4-2*K6。像素点(i，j)在拟合后的一阶偏导数和二阶偏导数可以用像素点(i，j)对称邻域内的灰度估计值f(r，c)(i，j)在r＝0，c＝0处的一阶和二阶偏导数分别表示为：
I ^ ( i , j ) x = ∂ f ( i , j ) ∂ r , I ^ ( i , j ) y = ∂ f ( i , j ) ∂ c ]]>
I ^ ( i , j ) xx = ∂ 2 f ( i , j ) ∂ r 2 , I ^ ( i , j ) xy = ∂ 2 f ( i , j ) ∂ r ∂ c , I ^ ( i , j ) yy = ∂ 2 f ( i , j ) ∂ c 2 ]]>
其中，
∂ f ( i , j ) ∂ r = K 2 - 17 5 K 7 , ∂ f ( i , j ) ∂ c = K 3 - 17 5 K 10 ]]>
∂ 2 f ( i , j ) ∂ r 2 = 2 * K 4 , ∂ 2 f ( i , j ) ∂ r ∂ c = K 5 , ∂ 2 f ( i , j ) ∂ c 2 2 * K 6 ]]>
(3计算每个像素点(i，j)在拟合后的Hessian矩阵，用矩阵H(i，j)表示：
H ( i , j ) = I ^ ( i , j ) xx I ^ ( i , j ) xy I ^ ( i , j ) xy I ^ ( i , j ) yy ]]>
并计算矩阵H(i，j)的两个特征值λ(i，j)1，λ(i，j)2；
(4)针对每个像素点(i，j)选取传导参数T(i，j)的形式：
T ( i , j ) = ( λ ( i , j ) 1 2 + λ ( i , j ) 2 2 ) 1 2 ]]>
(5)针对每个像素点(i，j)选取传递函数的形式：
c ( | ▿ I ^ | ) ( i , j ) = 1 1 + ( | ▿ I ^ ( i , j ) | T ( i , j ) ) 2 ]]>
其中，
| ▿ I ^ | ( i , j ) = I ^ ( i , j ) x 2 + I ^ ( i , j ) y 2 ]]>
其中，为像素点(i，j)梯度的模。
| ▿ I ^ | ( i , j ) = I ^ ( i , j ) x 2 + I ^ ( i , j ) y 2 ]]>
(6)依据各个像素点的传递函数采用各向异性扩散模型更新各个像素点灰度，从而更新原始图像灰度。
计算方法为：
I n + 1 ( i , j ) = I n ( i , j ) + Δt 4 d n ( i , j ) ]]>
上标n代表当前迭代次数，i和j分别代表原始图像中横纵坐标的取值，In(i，j)代表原始图像的像素点(i，j)在第n次迭代后的灰度值，Δt为采样间隔(一般取0.05-0.25)，dn(i，j)的表达形式为：
dn(i，j)＝cn(i，j-1)[In(i，j-1)-In(i，j)]+cn(i-1，j)[In(i-1，j)-In(i，j)]
+cn(i，j+1)[In(i，j+1)-In(i，j)]+cn(i+1，j)[In(i+1，j)-In(i，j)]
其中：
c ( i , j ) = c ( | ▿ I ^ | ) ( i , j ) ]]>
(7)重复步骤(2)至步骤(6)直到达到设定的迭代次数。
图2示出了本发明与传统各向异性扩散方法分别处理弱噪声DSA图像结果的对比图。
图2a是从原始的弱噪声DSA图像选取的一部分(像素为200×200)，图2b是采用用Catte算法得到的图像增强效果图(传递参数T＝30，50次迭代，高斯方差取为0.3)，图2d是采用Catte算法得到的图像增强效果图(传递参数T＝30，150次迭代，高斯方差取为0.3)，其结果是有的模糊了图像边缘，有的没有很好的抑制噪声。
图2c和图2f为本发明对原始图像分别迭代50次和150次的效果图，结果表明本发明对噪声抑制，图像边缘增强都起到非常好的作用。
图3示出了本发明与传统各向异性扩散方法分别处理强噪声DSA图像结果的对比图。
图3a是从原始的强噪声DSA图像选取的一部分(像素为200×200)，图3b是采用传统各向异性扩散Perona-Malik方法得到的图像增强效果图(传导参数T＝30，50次迭代)，图3c和图3d是采用Catte方法得到的图像增强效果图(传导参数T＝30，50次迭代，高斯方差分别取为0.3和0.4)，图3f是采用Catte方法得到的图像增强效果图(传导参数T＝30，150次迭代，高斯方差取为0.3)。从效果图可以看出，这些方法有的模糊了图像边缘，有的对噪声的抑制不强，都不能达到很好的图像增强效果。
图3e和图3g为本发明采用5×5区间进行小面模型拟合，分别对强噪声DSA图像迭代50次和150次的图像增强效果图，结果表明其有效降低了噪声的影响，增强了边缘，提高了增强效果。