一种利用序列数字减影血管造影图像分割血管数据的方法 【技术领域】
本发明属于医学成像技术领域，具体涉及一种在数字减影血管造影图像中分割血管数据的方法。
背景技术
数字减影血管造影(DSA：Digital Subtraction Angiography)技术在临床中已应用20多年，是心脑血管疾病无创诊断与介入治疗手术导航的重要依据。DSA图像处理中的一个关键任务就是进行图像分割，以使血管的生理特征能够更清楚地显示出来。结构分析，运动分析，三维可视化等后续操作，以及图像引导手术，肿瘤放射治疗，治疗评估等应用研究都是以图像分割为基础的。由于X射线经过的组织厚度以及血液中的造影剂浓度不均匀，经过减影后的图像不能完全消除人体组织引起的噪声信号，背景和目标混杂在一起。另外，人体的组织结构和形状很复杂，而且人与人之间有相当大的差异。因此血管的分割是一项困难的任务。
现有的血管分割技术大致可以分为6大类：模式识别方法、基于模型的方法、基于跟踪的方法、人工智能方法、神经网络方法、管状物体检测法，其中一些类型还可以细分(见“血管提取技术和算法的综述”，ACMComputing Surveys，36(2)：81-121，2004)。而阈值分割是最常见的直接检测区域的分割方法，其优点是简单，同时对于不同类的物体灰度值或其他特征值相差很大时，它能很有效的对图像进行分割，但它对图像中不存在明显灰度差异或灰度值范围有较大重叠的图像分割问题难以得到准确的结果。另外，由于它仅仅考虑图像的灰度信息而不考虑图像的空间信息，所以对噪声和灰度不均匀很敏感。
阈值分割方法分为全局阈值分割和局部阈值分割技术。早期用得较普遍的是全值阈值分割方法，但其并不适用于某些情况。以亮目标和暗背景的图像为例，目标的某些区域可能比背景的部分区域要暗，则任意一个确定的门限无法完全将目标和背景分离开来。也就是说，现有的全局分割技术很难在血管和背景灰度范围接近的局部地方取得好的分割效果，不适用于灰度直方图呈细窄的单峰分布的DSA图像，很难分得连续的血管段。
而另一方面，现有的局部阈值分割技术通常对图像进行分块后，对每个子块分别计算分割门限。由于未判断子块中是否存在血管，因此对不含血管的背景区域也将得到一个分割结果，这显然是不合理的。虽然直方图双峰法计算分割门限时，可利用直方图的双峰性在一定程度上对子块中是否存在血管作出判决，但其要求子块内所含血管点数和背景点数相当，这是十分苛刻很难满足的，且当子块较小，包含的象素数较少时，直方图双峰性的判断本身也存在困难。
Peng和Li(见“Knowledge-based Adaptive Thresholding Segmentation ofDigital Subtraction Angiography Images”，Image and Vision Computing25(2007)1236-1270)提出了基于繁忙度测度的局部阈值分割方法，在血管直径先验知识的指导下将图像分成适当大小的若干子块，采用有重叠的图像分块技术对各子块进行局部阈值分割，再融合各子块分割结果，充分考虑了DSA图像的灰度连续性，因而分割效果较好但是检测不到部分微细血管，也不能达到较好的去噪效果。
【发明内容】
本发明的目的在于提供一种利用序列数字减影血管造影图像分割血管数据的方法，该方法可以有效地提取出更完整的血管结构，并且操作简便，工作效率高。
本发明提供的利用序列数字减影血管造影图像分割血管数据的方法，其步骤包括：
(1)选取一个图像序列，记为I(x，y，t)，该图像序列I(x，y，t)包括开始注入造影剂到扩散到血管中的过程中的所有图像，并且该一系列图像为在时间上进行配准以后的数字减影血管造影图像；
(2)计算图像序列I(x，y，t)中的每一点在其时间域上的灰度标准差S(x，y)’S(x，y)是点(x，y)的时域统计方差；
(3)将图像中每一个点的灰度标准差S(x，y)映射到0～255的灰度空间，得到一幅血管树特征图像；
(4)对步骤(3)得到的血管树特征图像采用拉普拉斯滤波器进行锐化处理，增强边缘；
(5)对锐化后的血管树特征图像根据公式In+1(i,j)=In(i,j)+Δt4dn(i,j),]]>
求得n次迭代后的灰度值I(i，j)，迭代次数n初始值的取值范围为30～50，根据处理后图像血管区域被增强，背景区域被削弱的标准，不断调整n的取值，得到最佳的图像处理结果，作为各向异性平滑后的血管特征图；其中，In(i，j)是第n次迭代后像素点(i，j)处的灰度值，Δt是迭代的间隔时间，dn(i，j)采用下式计算：
dn(i，j)＝cn(i，j-1)[In(i，j-1)-In(i，j)]+cn(i-1，j)[In(i-1)-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(x,y)=1/1+(Ix2+Iy2/n)2]]>
Ix是点(i，j)处的x方向的边缘梯度值，Iy是点(i，j)处y方向的边缘梯度值；
(6)采用基于繁忙度测度的局部阈值分割方法对步骤(5)各向异性平滑后的血管特征图进行分割，得到二值图像；
(7)对步骤(6)分割后的得到的二值图像，采用标记连通域的方法，提取二值图像中面积最大的连通域，得到包含完整血管树结构的分割结果。
本发明针对现有全局阈值分割技术以及部分局部阈值分割技术对DSA图像处理的缺陷，根据DSA图像的特点，提出一种基于在时间上配准后的序列图像，利用统计特性方差的分割方法。本发明根据序列图像在时间上的变化，利用了统计量方差来描述图像上每一点的像素灰度值在时间上的变化信息，通过计算每一个像素点的方差值，然后映射到0～255的灰度空间，得到一幅描述每一点像素变化的特征图像，这样在时间上像素灰度值变化比较大的区域，即目标血管区域，方差值会比较大，反映在特征图上，目标被增强了；而在时间上像素灰度值变化比较小的区域，即不含血管的区域，方差值会比较小，反映在特征图上，背景被削弱了。本发明在得到了特征图以后，采用了基于血管直径的先验知识的重叠分块的局部阈值分割技术，对特征图进行二值分割。本发明还采用了最大连通域的标记技术，利用血管树在结构上的连通性，将血管树较完整得分割出来。
【附图说明】
图1为本方法的流程图；
图2为本发明实施例中某一序列图像，象素点的时间—灰度曲线，其中，图2a和图2b分别代表背景区域象素点和血管区域象素点；
图3为基于象素方差值映射后得到的特征图；
图4为对特征图进行拉普拉斯锐化后的结果图；
图5为对图4进行各向异性平滑后的结果图；
图6为在典型实施例中采用繁忙度测度的分割算法的结果图；
图7为对分割后的二值图像进行连通域标记后，得到的血管树；
图8为只采用序列中的最后一帧图片，进行与特征图相同的处理，首先拉普拉斯锐化，各向异性平滑，然后再用繁忙度算法分割以及连通域标记后的对比结果图。
【具体实施方式】
下面结合附图来说明本发明的典型实施例
本发明实例包括以下步骤：
(1)选取一个图像序列，记为I(x，y，t)，该图像序列I(x，y，t)包括开始注入造影剂到扩散到血管中的过程中的所有图像，并且该一系列图像为在时间上进行配准以后的数字减影血管造影图像；
(2)计算图像上的每一点在时间域上的灰度标准差值，计算公式为：S(x,y)=1TΣt=1T(I(x,y,t)-I‾(x,y))2,]]>S(x，y)其中是点(x，y)的时域统计标准差值，T是时间序列的长度，I(x，y，t)是第t帧图像中点(x，y)的灰度值，是点(x，y)在整个时间序列中的灰度平均值；
(3)计算出来的灰度标准差值范围大致在0～61.7之间，通过直线映射的方式将灰度标准差值变换到0～255的灰度空间中，形成了所要得到的血管树特征图，如图3所示；
(4)对特征图采用一个拉普拉斯的锐化滤波器进行处理，即用一个3×3的矩阵0101-41010]]>，与图像上的每一个点进行卷积相乘，将模块的中心与像素点对应，将乘积的结果赋给当前的像素点；
(5)对整幅图像采用各向异性平滑来消除噪声，根据公式In+1(i,j)=In(i,j)+Δt4dn(i,j),]]>，求得n次迭代后的邛度值I(i，j)，迭代次数n初始值的取值范围为30～50，根据处理后图像血管区域被增强，背景区域被削弱的标准，不断调整n的取值，得到最佳的图像处理结果，In(i，j)是第n次迭代后像素点(i，j)处的灰度值，Δt是迭代的间隔时间(一般取Δt＝0.05)，迭代次数n一般取30～50，dn(i，j)采用下式计算：
dn(i，j)＝cn(i，j-1)[In(i，j-1)-In(i，j)]+cn(i-1，j)[In(i-1)-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(x,y)=1/1+(Ix2+Iy2/n)2]]>
Ix是点(i，j)处的x方向的边缘梯度值，Iy是点(i，j)处y方向的边缘梯度
图像边缘梯度较大时，系数c(x，y)较小，扩散较弱，保持边缘；边缘梯度较小时，扩散系数大，扩散强度大，图像被模糊，从而实现各向异性的扩散能力；
(6)采用Peng和Li提出的基于繁忙度测度的局部阈值分割方法，对步骤(5)各向异性平滑后的血管特征图进行分割，得到二值图像。具体实施步骤如下：
(6.1)将DSA图像划分成N个互不重叠的a1×a2大小的子块E1、E2、……、EN，其中，d≤a1≤2d，d≤a2≤2d，d为血管直径的最大值；b1l、b2l、b3l、b4l为大小为b1×b2且包含子块E1在内的四个区域，其中，max(a1，a2)≤b1≤3max(a1，a2)，max(a1，a2)≤b2≤3max(a1，a2)；
(6.2)子区域b1l、b2l、b3l、b4l分别根据OTSU分割方法(最大类间房差)，求出每个区域的分割阈值Tλl，并对四个子区域进行分割；
在对区域bkl分割后所得的二值图对应的象素强度由fkl(i，j)来表示：
fkl(i,j)=0ifbkl(i,j)≥Tkl1ifbkl(i,j)<Tkl]]>
(6.3)采用繁忙度作为判断子区域分割是否合理的测度，来接受或者拒绝对子区域分割的假设门限。首先，利用如下公式：
busy(bkl)=Σi=1,j=1i=r,j=ctkl(i,j)]]>
然后，利用如下公式计算血管存在性标准的繁忙度门限Hσ
Hσ＝k·busymax+(1-k)·busymin
其中，busymax、busymin表示理想情况下b1×b2大小的区域繁忙度的最大或最小值，k是一个取值范围为[0，1]的系数；
如果busy(bkl)＞Hσ，则判定该区域bkl只包含背景，将二值区域各象素值fkl(i，j)置为背景所对应的逻辑值；否则，即判定该区域包含血管，保留二值区域的结果，即fkl(i，j)值保持不变；
(6.4)如果四个区域bkl均被判定为背景，其中k＝1、2、3和4，则将这四个区域的重叠子块El判为背景；如果四个区域中只有一个区域bλl被判定包含血管，其中1≤λ≤4，则将该区域的分割门限Tλl作为子块El的分割门限；如果四个区域中不只一个区域被判定包含血管，则对于针对不同区域bkl已求出的门限Tkl，选取其中最优的一个门限作为El的分割门限Tl；
(6.5)对子块El重复步骤(6.2)-(6.4)，完成整幅DSA图像的分割，得到分割后的二值图像
(7)对分割以后的图像，采用标记连通域的方法，提取二值图像中面积最大的连通域，得到包含较为完整血管树的分割结果。这里我们采用的是4-连通标记。
首先对分割得到的二值图像从左向右、从上向下进行扫描，如果当前像素值为0，就移到下一个扫描位置。如果当前像素值为1，就检查它左边和上边的2个近邻像素，如果他们全都为0，就给当前像素一个新的标记，如果上述两个近邻元素有一个像素值为1，就把该像素的标记赋给当前像素，如果他们的值都为1，且具有相同的标记，就将给标记赋给当前像素，如果他们的值为1但是不具有相同的标记，就将其中的1标记赋给当前像素并作个记号表明这2个标记等价。在扫描终结时将所有等价的标记对归入等价组，对各个组赋1个不同的标记，然后二次扫描图像，将每个标记用它所在等价组的标记代替。标记后得到的血管树的结果图如图7所示。
利用本方法可以得到较为完整的血管树的结构，完成将血管分割出来的工作，为随后的三维重建作基础。
通常，我们都是采用某一幅相对血管比较清晰的图像用来分割，但是由于造影剂随血液的流动，不同部位血管中的造影剂浓度不同，使得单帧中包含的血管树不连续，致使分割出来的血管不完整。本发明利用了造影剂随时间不断在血管中流动，血管区域的像素灰度值不断变化的信息，将完整的血管结构通过像素点的标准差信息反映出来，得到了一幅富含完整血管结构信息的特征图。对特征图进行拉普拉斯锐化增强边缘，用各向异性平滑消除噪声，然后用基于繁忙度测度的分割算法进行分割，对分割得到的二值图像进行连通域标记，最终将血管树分割出来。该方法相比对于单帧图像进行相同的处理，能够更好地分割出完整的血管树，取得了很好的效果。