在数字冠脉造影图像中分割血管数据的方法

摘要

本发明公开了一种在数字冠脉造影图像中分割血管数据的方法，其步骤为：①对冠脉造影图像序列进行频域增强，得到增强后的新图像序列f(x，y，t)；②在新图像序列f(x，y，t)中取出待分割单帧增强图像f(x，y)，对f(x，y)分别进行针对主动脉和细小血管的两个尺度下的圆盘结构低帽变换，得到灰度级图像g1(x，y)和g2(x，y)；③分别对g1(x，y)和g2(x，y)进行“相对阈值比较”提取灰度图象中的血管，得到二值图像B1(x，y)和B2(x，y)④合并两幅二值图，得到分割结果图。本发明方法充分利用了血管的运动、结构、尺度、形状、灰度五方面的信息，能对造影图像有效分割，操作简便，工作效率高。

说明书

技术领域

本发明属于医学成像技术领域，具体涉及一种在数字冠脉造影图像中分割血管数据的方法。

背景技术

数字血管造影技术在临床已应用20多年，是心脑血管疾病无创诊断与介入治疗手术导航的重要依据。冠脉造影术是使用最普遍的用于诊断冠状动脉疾病并评估狭窄程度的图像方法。造影图像处理中的一个关键任务就是进行图像分割，以使血管的生理特征能够更清楚的显示出来，结构分析、运动分析、三维可视化等后续操作，以及图像引导手术、治疗评估等应用研究都是以图像分割为基础的。然而由于X射线经过的组织厚度不同、血液中的造影剂浓度不均匀以及不同背景组织和目标混杂在一起，使得心脏冠状动脉造影图像的信噪比很低。另外，人体的组织结构和形状很复杂，而且人与人之间有相当大的差异。因此造影图像的血管分割是一项非常困难的任务。

C.Kirbas和F.Quek(“A review of vessel extraction techniques andalgorithms，”ACM Computing Surveys，vol.36 pp.81-121 June 2004.)提出，现有的的血管分割技术大致可以分为6大类：模式识别方法、基于模型的方法、基于跟踪的方法、人工智能方法、神经网络方法、管状物体检测法，而这些传统的针对造影图像的血管分割方法几乎都是利用单幅图像的灰度或者血管结构信息。由于冠脉造影图像之间的差异性很大并且单幅图像信噪比很低，使得这类方法的适应性存在很大问题并且不能提取出所有的血管，尤其是对比度低的部分。

发明内容

本发明的目的在于提供一种在数字冠脉造影图像中分割血管数据的方法，该方法能对造影图像有效分割，并且操作简便，工作效率高。

本发明提供的在数字冠脉造影图像中分割血管数据的方法，其步骤包括：

(1)对图像序列f₀(x，y，t)进行频域增强，得到新图像序列f(x，y，t)，其中图像序列f₀(x，y，t)的帧数为T，单幅图像大小为M×N，0≤x＜M，0≤y＜N，0≤t＜T，x，y，t均为整数；

(2)在新图像序列f(x，y，t)取出待分割单帧增强图像f(x，y)，对f(x，y)分别进行直径为D₁和D₂的圆盘结构低帽变换，得到g₁(x，y)和g₂(x，y)；其中，D₁＞D₂，D₁D₂均为正整数，D_MAX＜D₁≤(D_MAX+5)，D_MIN＜D₂≤(D_MIN+(D_MAX-D_MIN)/2)，D_MAX为最宽冠状动脉的直径，D_MIN为最窄冠状动脉的直径；

(3)按照下述步骤分别对g₁(x，y)和g₂(x，y)进行阈值比较，得到二值图像B₁(x，y)和B₂(x，y)：

(3.1)利用式(I)比较低帽变换图像g(x，y)和增强图像f(x，y)得到差值图d(x，y)，其中，g(x，y)为g₁(x，y)或g₂(x，y)；

(3.2)利用式(II)计算阈值T₁，T₁为d(x，y)图中的非零像素值的平均值

$>>>T>1>>=>>>>\Sigma >>d>>(>x>,>y>)>>>>0>\; >d>>(>x>,>y>)>>>>N>d>>>->->->>(>II>)>>>s>$

其中N_d为d(x，y)图中非零像素值的像素个数；

(3.3)若d(x，y)≥T₁，则该点为血管象素点，若d(x，y)＜T₁，则该点为背景点，得到二值图B₁(x，y)和B₂(x，y)；

(4)叠加二值图B₁(x，y)和B₂(x，y)，去除面积低于阈值T₂的块状噪声区域，T₂为根据图像大小和血管粗细选取的阈值，0＜T₂≤(D_MIN×D_MIN)，所得到的二值图，即为冠状动脉的提取结果。

针对传统的数字冠状动脉造影图像分割方法几乎都是采用单幅图像的灰度信息，所利用信息有限的缺陷，本发明首先利用冠状动脉图像序列中血管和背景的不同运动变化，采用时域傅立叶高通变换来减弱低频背景噪声并增强血管。同时，还考虑到血管的管状结构、直径信息，采用符合血管结构的多尺度形态学滤波进一步增强待分割单帧图像。最后，根据血管与背景在不同图像之间的灰度差异一致性，分割出二值图像。本发明方法充分利用了血管的运动、结构、尺度、形状、灰度五方面的信息，对不同图像之间的较大差异适应性好，冠状动脉分割结果完整，对于细小血管或者血管末尾等对比度弱的部分也有很好表现。

附图说明

图1为本发明方法的流程图；

图2为左冠状动脉(LCA)树的造影图像序列(48*512*512大小，8位)；

图3为图2经过步骤(1)(β＝1/5)后的增强图像序列；

图4(a)为从原始图像序列(图2)中取出的第22帧图像，图4(b)为从增强图像序列(图3)中取出的第22帧图像；

图5为以图4(b)为例的动脉提取过程。图5(a)示出直径D₁＝35的圆盘结构低帽(Bottom-Hat)变换变换图；图5(b)示出直径D₁＝5的圆盘结构低帽(Bottom-Hat)变换变换图；图5(c)，图5(d)分别为图5(a)，图5(b)“相对阈值比较”结果图；

图6所示的反映了在本发明的典型实施例中使用上述方法与单幅图像分割方法的一个比较结果，图6(a)为原始序列图中第22帧图像，6(b)表示本方法的分割结果；图6(c)表示用S.Ehio基于形态学方法的分割结果；图6(d)表示用全局阈值分割方法的结果；图6(a)中的白矩形框内的暗状长条形区域用以标注用于注入造影剂的导管。

具体实施方式

下面结合附图和实例对本发明作进一步详细的说明。

如图1所示，本发明方法的步骤包括：

(1)对图像序列f₀(x，y，t)进行频域增强，该图像序列f₀(x，y，t)的帧数为T，单幅图像大小为M×N，0≤x＜M，0≤y＜N，0≤t＜T，x，y，t均为整数。

对图像序列进行频域增强可以采用多种方法，现举例说明其中一种方法的具体步骤：

(1.1)对图像序列f₀(x，y，t)进行M×N次沿时间轴t的离散傅立叶变换，得到频域信号F₀(x，y，k)，0≤k＜T，离散傅立叶变换公式如式(1)：其中j为虚数单位

$>>>F>0>>>(>x>,>y>,>k>)>>=>>1>T>over>>\Sigma >>t>=>0>>>T>->1>>>>f>0>>>(>x>,>y>,>t>)>>>e>>->j>2>\pi kt>/>T>\; >>s>$

            0≤k＜T，0≤x＜M，0≤y＜N    (1)

(1.2)对频域信号F₀(x，y，k)采用式(2)进行M×N次沿k轴的高通滤波，高通滤波函数为H(x，y，k)＝(1-e^-βt)，其中β是预先给定的滤波常数，用以抑制某些频率的背景噪声，经过多次实验，β值的优选范围为1/8--1/5。

            F(x，y，k)＝F₀(x，y，k)(1-e^-βm)

            0≤k＜T，0≤x＜M，0≤y＜N，(2)

(1.3)利用离散傅立叶反变换(式(3)，其中j为虚数单位)作用于滤波后的频域信号，获得新的图像序列：

$>>f>>(>x>,>y>,>t>)>>=over>>\Sigma >>k>=>0>>>T>->1>>ver>>F>‾>>>(>x>,>y>,>k>)>>>e>>j>2>\pi k>/>T>\; >>s>(3)$

            0≤t＜T，0≤x＜M，0≤y＜N

(2)在新图像序列f(x，y，t)取出待分割单帧增强图像f(x，y)，对f(x，y)分别进行直径为D₁和D₂的圆盘结构低帽(Bottom-Hat)变换，得到g₁(x，y)和g₂(x，y)；其中，D₁＞D₂，D₁ D₂均为正整数，D_MAX＜D₁≤(D_MAX+5)，D_MIN＜D₂≤(D_MIN+(D_MAX-D_MIN)/2)，D_MAX为最宽冠状动脉的直径，D_MIN为最窄冠状动脉的直径。

形态学低帽(Bottom-Hat)变换是一种常见的用于提取图像暗结构的形态学算子，该算子用g表示，定义为g＝I-(I·b)，其中I是输入图像，b是结构元素，f·b代表闭操作。

(3)分别对g₁(x，y)和g₂(x，y)进行“相对阈值比较”，得到二值图像B₁(x，y)和B₂(x，y)。过程如下：

(3.1)利用式(4)比较低帽变换图像g(x，y)和增强图像f(x，y)得到差值图d(x，y)，其中，g(x，y)为g₁(x，y)或g₂(x，y)；

(4)

(3.2)利用式(5)计算阈值T₁，T₁为d(x，y)图中的非零像素值的平均值

$>>>T>1>>=>>>>\Sigma >>d>>(>x>,>y>)>>>>0>\; >d>>(>x>,>y>)>>>>N>d>>>>s>(5)$

其中N_d为d(x，y)图中非零像素值的像素个数。

(3.3)若d(x，y)≥T₁，则该点为血管象素点，若d(x，y)＜T₁，则该点为背景点，得到二值图B₁(x，y)和B₂(x，y)。

(4)二值图合并：叠加二值图B₁(x，y)和B₂(x，y)，去除面积低于阈值T₂的块状噪声区域，T₂为根据图像大小和血管粗细选取的阈值，0＜T₂≤(D_MIN×D_MIN)，所得到的二值图即为冠状动脉的提取结果。

实例：

该实例图像如图2所示，图像中含有暗目标和亮背景，暗目标表示动脉血管，下面详细说明本实例的步骤：

(1)对图2所示的图像序列f₀(x，y，t)(0≤x＜512，0≤y＜512，0≤t＜48)进行频域增强，这个序列包含48幅从心脏收缩期到心脏扩张期再到心脏收缩期的图像。其过程为：

(1.1)对图像序列f₀(x，y，t)进行512×512次沿时间轴t的离散傅立叶变换，得到48×512×512点数的频域信号F₀(x，y，k)，0≤k＜48，离散傅立叶变换公式如下：

$>>>F>0>>>(>x>,>y>,>k>)>>=>>1>48>over>>\Sigma >>t>=>0>>47>>>f>0>>>(>x>,>y>,>t>)>>>e>>->j>2>\pi kt>/>48>\; >>s>$

               0≤k＜48，0≤x＜512，0≤y＜512   (1)

(1.2)对频域信号F₀(x，y，k)采用(2)式进行512×512次沿k轴的高通滤波，高通滤波函数为 $>>H>>(>x>,>y>,>k>)>>=>>(>1>->>e>>->>1>5>>t>\; >)>>>s>$

$>ver>>F>‾>>>(>x>,>y>,>k>)>>=>>F>0>>>(>x>,>y>,>k>)>>>(>1>->>e>>->>1>5>>m>\; >)>>>s>$

(2)

(1.3)对滤波后的频域信号F(x，y，k)进行512×512次离散傅立叶反变换(式(3))，获得如图4所示的新时域图像序列：

$>>f>>(>x>,>y>,>t>)>>=over>>\Sigma >>k>=>0>>47>ver>>F>‾>>>(>x>,>y>,>k>)>>>e>>->j>2>\pi k>/>48>\; >>s>$

               0≤t＜48，0≤x＜M，0≤y＜N    (3)

(2)从增强图像序列中，取出第22帧待分割单帧图像f(x，y)，对f(x，y)分别进行直径D₁＝35，D₂＝5的圆盘结构低帽(Bottom-Hat)变换，得到g₁(x，y)(图5(a))和g₂(x，y)图5(b)。

(3)分别对g₁(x，y)和g₂(x，y)进行“相对阈值比较”得到二值图像B₁(x，y)(图5(c))和B₂(x，y)(图5(d))。

(3.1)利用式(4)比较低帽变换图像g(x，y)和增强图像f(x，y)得到差值图d(x，y)，其中，g(x，y)为g₁(x，y)或g₂(x，y)；

(4)

(3.2)利用式(5)计算阈值T₁，T₁为d(x，y)图中的非零像素值的平均值

$>>>T>1>>=>>>>\Sigma >>d>>(>x>,>y>)>>>>0>\; >d>>(>x>,>y>)>>>>N>d>>>>s>(5)$

其中N_d为d(x，y)图中非零像素值的像素个数。

(3.3)若d(x，y)≥T₁，则该点为血管象素点，若d(x，y)＜T₁，则该点为背景点，得到二值图B₁(x，y)和B₂(x，y)。

(4)二值图合并：叠加二值图B₁(x，y)和B₂(x，y)，去除面积低于阈值T₂＝20的块状噪声区域，所得到的二值图(图6(b))即为冠状动脉的提取结果。

由图6(b)看出，采用本发明得到的分割结果(图6(b))分割完整，残留噪声少，图6(a)中白色矩形框内的静态造影剂导管部分绝大部分被去除；但仅根据单幅图像的灰度信息和结构信息方法得到的分割结果(图6(c)，图6(d))，血管结构不完整：对比度弱的小血管被去掉、分割引入很大噪声、静态导管存在于分割结果中。

采用本发明方法，有效的提高了图像的对比度，能够从低信噪比造影图像中自动地提取整支动脉树，对于细小血管部分或者血管末尾部分也有很好表现。