基于亚像素的血管分割方法及其分割系统与流程

本发明涉及血管分割技术领域，尤其涉及一种基于亚像素的血管分割方法及其分割系统。

背景技术：

近年来，随着生物医学成像技术及医学影像设备的迅猛发展，医学图像的分辨率逐日提高，使得医学影像成为医学技术当中发展最为迅速的领域之一。因此，医学图像处理技术在医学科研和临床实践的作用、影响极其重要。医学图像分割是图像分割技术在医学领域的重要应用，是医学图像处理中重要的一部分，能够给医生提供重要的临床参考价值。在医学图像分割中，血管分割一直是研究的热点。

随着人类生活水平的提高，心脑血管疾病成为人类的头号死因，对血管疾病的早期定量诊断和风险评估对延长人类预期寿命，提高人类的身体抵抗力，起着非常关键的作用。x射线血管造影(x-rayangiogram,xra)是通过造影剂显示血管的一种方法。从经过ct扫描得到的造影图像中准确分割出血管树是进行可视化三维血管重建及血管疾病诊断和治疗的重要前提，也是进行辅助诊断和外科手术的重要保障。ct血管图像进行分割必须根据ct图像的特点寻求最佳的可行方法，以提高医生诊断识别的正确率。鉴于血管复杂的形态结构，以及成像环境中存在的噪音、造影剂衰减和分布不均匀等情况，导致图像中常存血管目标模糊和局部形态特异等各种问题，从而无法对目标对象(特别是血管细小末梢)进行准确分割。

因此，有必要提供一种新的基于亚像素的血管分割方法及其分割系统来解决上述问题。

技术实现要素：

本发明需要解决的技术问题是提供一种基于亚像素的血管分割方法及其分割系统，其可以对血管进行准确分割。

本发明提供了一种基于亚像素的血管分割方法，包括以下步骤：

s1、对ct扫描获取的dicom文件序列通过高位截取、反变换和调窗处理三个环节进行预处理，将所述dicom文件序列转换为计算机通用的bmp文件序列；

s2、对所述bmp文件序列使用canny算法进行粗边缘检测，并生成粗边缘检测结果；

s3、将所述粗边缘检测结果生成链表记录下来，并保留原始灰度；

s4、使用zernike矩对记录下来的所述粗边缘检测结果进行亚像素的边缘定位以获取边缘坐标；

s5、根据所述边缘坐标使用漫水填充算法将血管区域填充，并通过阈值分割获取三维血管结构；

s6、通过鼠标在所述三维血管结构中点击获取初始种子点，采用26邻域对所述三维血管结构进行区域生长分割，实现交互式分割出指定血管结构。

优选的，在步骤s1中，调窗处理环节将变换公式

和center＝c/rs-ri,width＝w/rs-rimin＝center-width,max＝center+width中的输入值划分为三类，划分标temp＝1024/(max-min)×gm

准关联图像的截距和斜率，分别为pixel_in≤min、pixel_in≥max和other，其中，c为窗位，w为窗宽，rs、ri分别为斜率和截距，gm为显示器的最大显示值，pixel_in为输入的图像数据原始值，pixel_out为输出的显示灰度值。

优选的，在步骤s2中，对于canny算法中的阈值选取，结合otsu算法，自适应的获取高低阈值，将canny算法得到的非极大值抑制图像n(x,y)中的像素划分为n1、n2和n3三个类别，根据三个类别的灰度值期望，定义类间方差σ²(k,t)，根据公式获得两个最佳阈值k^*和t^*，其中，n1代表非极大值抑制图像n(x,y)中的非边缘点，n2代表非极大值抑制图像n(x,y)中需要判断边缘点的点，n3代表非极大值抑制图像n(x,y)中的边缘点。

优选的，在步骤s4中，选取zernike矩的7x7模板进行亚像素边缘定位。

优选的，在步骤s6中，通过鼠标在指定血管中选取表现特征明显的点作为用来进行方向引导的初始种子点，将初始种子点与其邻域的体像素点之间的灰度差作为新种子点的判断依据一；将各相邻初始种子点之间的单位向量与初始种子点到其邻域点方向的单位向量的数量积作为新种子点判断依据二；根据判断依据一和判断依据二实现交互式分割指定血管结构。

本发明还提供了一种基于亚像素的血管分割系统，包括：

图像序列预处理模块，用于对ct扫描获取的dicom文件序列通过高位截取、反变换和调窗处理三个环节进行预处理，并将所述dicom文件序列转换为计算机通用的bmp文件序列；

亚像素边缘检测血管分割模块，用于对bmp文件序列使用canny算法和zemike矩进行处理以得到三维血管结构；以及

交互式分割指定血管结构模块，用于采用26领域对三维血管结构进行区域生长分割，实现交互式分割出指定血管结构。

优选的，所述亚像素边缘检测血管分割模块包括：

粗边缘检测模块，用于使用canny算法对对所述bmp文件序列进行粗边缘检测；

记录模块，用于将粗边缘检测的结果生成链表记录下来，并保留原始灰度；

边缘坐标获取模块，用于使用zernike矩对记录下来的粗边缘检测结果进行亚像素的边缘定位以获取边缘坐标；以及

三维血管结构形成模块，用于根据所述边缘坐标使用漫水填充算法将血管区域填充，并通过阈值分割获取三维血管结构。

与相关技术相比，本发明提供的基于亚像素的血管分割方法及其分割系统具有如下有益效果：

1).对dicom文件转换为bmp文件中的调窗处理算法进行改进：对常用的线性变换公式进行改进，将输入值和rs、ri关联，解决了在调窗处理后像素值可能发生的越界现象，有效的改善了显示效果；

2)提出canny_zernike矩法进行亚像素边缘检测分割精确的血管结构，对处理后的图像利用改进的canny算法进行粗边缘检测，然后使用zernike矩对粗边缘检测的结果进行亚像素级的边缘定位，可精确地分割血管的细小末梢；

3)提出自适应canny算子粗边缘检测算法，结合otsu算法(最大类间方差)的思想，自适应的获取高低阈值，解决了阈值选取难的问题，并且血管边缘更加清晰，细节更加完善；

4)能交互式的分割指定血管区域，可以将血管分割达到亚像素级精度。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图，其中：

图1为本发明基于亚像素的血管分割系统的结构框图；

图2图1中所示亚像素边缘检测血管分割模块的结构框图；

图3为本发明基于亚像素的血管分割系统的工作流程图；

图4为本发明基于亚像素的血管分割方法的流程图；

图5(a)为未经过图像序列预处理过的原始图像；

图5(b)为经过图像序列预处理过的bmp图像；

图6(a)为本发明理想亚像素边缘模型的立体示意图；

图6(b)为本发明理想亚像素边缘模型的平面示意图；

图6(c)为本发明理想亚像素边缘模型的平面旋转示意图；

图7(a)为整个血管结构图；

图7(b)对图6(a)所示整个血管结构经过本发明基于亚像素的血管分割方法及其分割系统分割出的指定血管结构图。

具体实施方式

下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

请参阅图1，本发明提供了一种基于亚像素的血管分割系统100，包括图像序列预处理模块1、亚像素边缘检测血管分割模块2以及交互式分割指定血管结构模块3。

所述图像序列预处理模块1用于对ct扫描获取的dicom文件序列通过高位截取、反变换和调窗处理三个环节进行预处理，并将所述dicom文件序列转换为计算机通用的bmp文件序列。通过对调窗处理算法进行改进，对常用的线性变换公式进行改进，将输入值和rs、ri关联，解决了在调窗处理后像素值可能发生的越界现象，有效的改善了显示效果，突出显示了图像的关键部分。

所述亚像素边缘检测血管分割模块2用于对bmp文件序列使用canny算法和zemike矩进行处理以得到三维血管结构，从而可精确地分割血管的细小末梢。另外，提出自适应canny算法粗边缘检测算法，结合了otsu算法的思想，自适应的获取高低阈值，解决了阈值选取难的问题，并且血管边缘更加清晰，细节更加完善。

结合图2所示，所述亚像素边缘检测血管分割模块2具体包括：

粗边缘检测模块20，用于使用canny算法对所述bmp文件序列进行粗边缘检测；

记录模块21，用于将粗边缘检测的结果生成链表记录下来，并保留原始灰度；

边缘坐标获取模块22，用于使用zernike矩对记录下来的粗边缘检测结果进行亚像素的边缘定位以获取边缘坐标；以及

三维血管结构形成模块23，用于根据所述边缘坐标使用漫水填充算法将血管区域填充，并通过阈值分割获取三维血管结构。

所述交互式分割指定血管结构模块3用于采用26领域对三维血管结构进行区域生长分割，实现交互式分割出指定血管结构。

结合图3所示，本发明提供的所述基于亚像素的血管分割系统100先通过所述图像序列预处理模块1对图像序列进行预处理，然后通过所述亚像素边缘检测血管分割模块2进行血管边缘检测，具体为像素级canny算法粗定位边缘和亚像素zernike矩精确定位边缘，从而获得三维血管结构，最后通过所述交互式分割指定血管结构模块3基于交互式的区域生长获取指定血管结构。

结合图4所示，本发明还提供了一种根据上述基于亚像素的血管分割系统100的分割方法，该分割方法包括以下步骤：

s1、对ct扫描获取的dicom文件序列通过高位截取、反变换和调窗处理三个环节进行预处理，将所述dicom文件序列转换为计算机通用的bmp文件序列；

s2、对所述bmp文件序列使用canny算法进行粗边缘检测，并生成粗边缘检测结果；

s3、将所述粗边缘检测结果生成链表记录下来，并保留原始灰度；

s4、使用zernike矩对记录下来的所述粗边缘检测结果进行亚像素的边缘定位以获取边缘坐标；

s5、根据所述边缘坐标使用漫水填充算法将血管区域填充，并通过阈值分割获取三维血管结构；

s6、通过鼠标在所述三维血管结构中点击获取初始种子点，采用26邻域对所述三维血管结构进行区域生长分割，实现交互式分割出指定血管结构。

在步骤s1中，将dicom文件序列转换为bmp文件序列主要包括如下步骤：

(1)读取dicom文件序列，并保存图像显示所必须的图像数据信息(图像帧数、图像宽度、图像高度、分配位数(bitsallocated)、存储位数(bitsstored)、最高位数(highbit)以及图像窗位(windowscenter)、窗宽(windowwidth)、截距(rescaleintercept)、斜率(rescaleslope)等)和参数信息；

(2)图像数据处理，包括高位截取、反变换和调窗处理三个环节；

(3)填充bmp文件信息，将处理后的图像数据和参数按照bmp文件格式依次填充并显示。

所谓调窗处理(c/w)，即窗宽窗位调整，原理上是根据预知的窗宽和窗位值，获得需要显示的窗口大小(窗宽w)和中心位置(窗位c)，从而将窗口内的值转换成显示时的最亮和最暗范围内的值，将高于窗口上限的灰度值设置为显示器灰度最高值，低于窗口下限的灰度值设置为灰度最低值，通过动态改变窗宽窗位值移动窗口来显示不同的灰度段。

在调窗处理后，像素值可能发生越界现象(即处理后的像素值超出了实际的像素值范围)，本发明针对对此种情况，对常用的线性变换公式进行改进，将输入值和rs、ri关联，有效的改善原有算法的显示效果，公式如下：

和

center＝c/rs-ri,width＝w/rs-ri

min＝center-width,max＝center+width

temp＝1024/(max-min)×gm

调窗处理环节将上述变换公式中的输入值划分为三类，划分标准关联图像的截距和斜率，分别为pixel_in≤min、pixel_in≥max和other，其中，c为窗位，w为窗宽，rs、ri分别为斜率和截距，gm为显示器的最大显示值，pixel_in为输入的图像数据原始值，pixel_out为输出的显示灰度值。如图5(a)和图5(b)所示，通过截取、反变换和调窗处理三个环节进行预处理的图像与未经过预处理的图像的对比，可见，经过预处理后，有效改善了显示效果，且像素值不发生越界，突出显示了图像的关键部分。

在步骤s2中，对于canny算法中的阈值选取，结合otsu算法，自适应的获取高低阈值，将canny算法得到的非极大值抑制图像n(x,y)中的像素划分为n1，n2和n3三个类别，其中n1中包含灰度值为{m1,m2,...,mk}的像素，代表着原图中的非边缘点；n2中包含灰度值为{mk+1,mk+2,...,mt}，代表着原图中需要判断边缘点的点；n3中包含灰度值为{mt+1,mt+2,...,ml}，代表着原图中的边缘点。另外，设图像中总的像素数为n，灰度值为mi对应的像素个数为ni，那么其概率为则整个区间的灰度值期望为在n1，n2，n3三个类内的灰度值期望分别为：

令

则可以定义类间方差为：

σ²(k,t)＝p1(k)·(e1(k)-e)²+p2(k,t)·(e2(k,t)-e)²+p3(t)·(e3(t)-e)²

并且以下关系成立：

p1(k)·e1(k)+p2(k,t)·e2(k,t)+p3(t)·e3(t)＝e和p1(k)+p2(k,t)+p3(t)＝1

两个最佳阈值k^*和t^*是使得σ²(k,t)最大的值，用下式寻找最佳阈值：

该过程由选择第一个k值开始(该值是1，因为在0灰度值处寻找阈值没有意义，增量值为整数)。接着，t的所有制在大于k和小于l的范围内增加。然后，将k增大到其下一个值，t的所有值再次在大于k的所有值范围内增加。重复该过程，直到k＝l-2为止。该处理的结果是一个二维阵列σ²(k,t)，最后一步是在该阵列中寻找最大值。对应于最大值的k值和t值就是最佳阈值k^*和t^*。如果存在几个最大值，则对应于k和t的值被平均以得到最终的阈值。

在步骤s4中，选取zernike矩的7x7模板进行亚像素边缘定位。图像f(x,y)的n阶m次zernike矩^[ii]定义为：

若图像旋转角，则旋转之前的zernike矩z′nm和旋转之后的zernike矩znm之间有如下关系：

从式中可以看出，图像旋转前后只是相角发生了变化，而zernike矩的模不变，这个就是zernike矩的旋转不变性。利用旋转后的zernike矩可以很容易计算出边缘参数，从而实现对比边缘的亚像素定位。

基于zernike矩的边缘检测算法原理是通过建立zernike矩与理想的亚像素模型的4个边缘参数的关系，分别求解矩而得到图像模型的4个边缘参数，再将参数与预设的阈值进行比较判断，进而精确定位边缘点。理想亚像素边缘模型如图6(a)～(c)所示，其中，k为灰度阶跃高度，h为背景灰度，l(l∈[-1,1])为中心到边缘的距离，为l与x轴之间的夹角。单位圆内虚线是理想阶跃边缘，直线两侧的灰度分别为目标(h+k)和背景(k)。图6(c)是图6(b)顺时针旋转的结果。

根据图6(a)～(c)中的模型，计算旋转角度后不同阶次的zernike矩表示为：

由上述公式可以得到边缘的3个参数，分别表示为：

im[z′11]是关于y的奇函数，因此故

在离散数字图像情况下，zernike矩的计算可以采用模板与图像灰度的卷积。本文考虑到模板效应，若采用单位圆进行采样，在n×n像素区域内，当模板在图像上移动进行卷积时，模板覆盖的是模板中心周围n²个像素，此时单位半径变为了因此需要把单位圆上计算出来的垂直距离l放大倍，故像素点的亚像素坐标为：

虽然zernike矩模板n值越大，亚像素边缘定位精度越高，但是也增大了计算量，因此本发明优选的选取zernike矩算子的7×7模板进行亚像素边缘定位。

在对整个图像利用改进的canny算法进行自适应阈值边缘检测，并将粗边缘检测的结果生成链表记录下来，保留原始灰度后，得到的新图像，利用zernike矩对保存在数组中可能的边缘点图像重新定位图像边缘，步骤如下：

(1)计算7×7模板系数；

(2)将每个像素与模板系数求卷积得到zernike矩的z00,z11和z20；

(3)根据上述公式计算出旋转角

(4)根据上述公式以及计算出参数l和k；

(5)确定亚像素边缘的判断依据，即k≥kt∩l≤lt(kt,lt作为判断阈值，kt为k矩阵取均值所得，lt为l矩阵取均值的倍所得)，若像素点满足该条件，则该像素点是边缘点，然后利用上述公式计算出亚像素坐标。最后，根据获取的边缘坐标使用漫水填充算法将血管区域填充，然后通过阈值分割获取三维血管结构。

在步骤s6中，关于初始种子点的选取：用户通过鼠标在三维血管结构中点击获取若干个“关键点”，一般来说在相对趋于直线的血管区域上指定少量关键点，在血管弯曲区域指定较多关键点，以提高分割的准确性。这些“关键点”作为初始种子点，种子点在三维情况下称为体素。同时建立一个空队列，将种子点压入其中，则队列被初始化为a1,a2,...an，其中a1,a2,...an表示选取的初始种子点或种子区域。对于一个断层大小为m×n像素的三维图像，且初始种子点ai到ai+1在z方向上共h层切片，其体数据表示为v＝{f(i,j,h)|i＝1,2,...,m；j＝1,2,...,n；h＝1,2,...h}。

通过鼠标在指定血管中选取表现特征明显的点作为用来进行方向引导的初始种子点，将初始种子点与其邻域的体像素点之间的灰度差作为新种子点的判断依据一；将各相邻初始种子点之间的单位向量与初始种子点到其邻域点方向的单位向量的数量积作为新种子点判断依据二；根据判断依据一和判断依据二实现交互式分割指定血管结构。

迭代过程：从种子点出发，研究其邻域内尚未经过处理的体像素点t.t可表示为其中n(x)表示体素点x的邻域。对于三维图像，其邻域类型分为6邻域、18邻域和26邻域。

本发明优选的采用26邻域对数据进行处理，种子点和其邻域之间的相似性判断依据为|g(p)-g(q)|<τ，其中，g(p)表示种子点p处的灰度值，g(q)表示种子点某一邻域点q处的灰度值，τ为选定的阈值。

另外，种子点ai的空间坐标为(xi,yi,zi)，种子点ai+1的空间坐标为(xi+1,yi+1,zi+1)，将种子点ai到ai+1方向的单位向量定义为向量种子点ai到ai+1在方向上h层切片中任意种子点到其邻域内尚未经过处理的体像素点t方向的单位向量定义为向量将两向量的数量积作为获取指定区域的判断标准，其中δ为选定的阈值。即判断依据为如下

当检测点满足上述判断依据公式时，将该点作为新种子点归入队列中。随着迭代的进行，当没有更多的新种子点产生时，则区域生长结束，当前队列中所有的像素所构成的区域就是分割结果。

本发明使用扫描的股前外侧皮瓣的血管ct数据，得到指定血管的三维结构图，如图7(a)和(b)所示，图7(a)为整个血管结构，图7(b)为指定的血管结构，可见，本发明可以精确的分割出指定的血管结构。

与相关技术相比，本发明提供的基于亚像素的血管分割方法及其分割系统具有如下有益效果：

1).对dicom文件转换为bmp文件中的调窗处理算法进行改进：对常用的线性变换公式进行改进，将输入值和rs、ri关联，解决了在调窗处理后像素值可能发生的越界现象，有效的改善了显示效果；

2)提出canny_zernike矩法进行亚像素边缘检测分割精确的血管结构，对处理后的图像利用改进的canny算法进行粗边缘检测，然后使用zernike矩对粗边缘检测的结果进行亚像素级的边缘定位，可精确地分割血管的细小末梢；

3)提出自适应canny算子粗边缘检测算法，结合otsu算法(最大类间方差)的思想，自适应的获取高低阈值，解决了阈值选取难的问题，并且血管边缘更加清晰，细节更加完善；

4)能交互式的分割指定血管区域，可以将血管分割达到亚像素级精度。

以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领域，均同理包括在本发明的专利保护范围内。