1 基于改进的三维区域生长的骨骼组织外轮廓提取传统的医学影像处理都是基于每层影像单独处理的.近年来,越来越多的专家****都提出从三维空间中整体处理医学影像信息.从三维空间中寻找感兴趣的组织结构相比单纯的从二维空间中查找需要更大的计算量和需要考虑更复杂的空间结构,但是也相应地额外获取了层与层之间的空间信息.尤其是有些医学影像数据集合中,某些组织器官如果单纯从每层中独立分离提取,它们就会形成断裂.有些细小的骨骼组织末端由于容积效应,使得它的像素的灰度值比普通骨骼组织的灰度值低很多,并且它在当前层的影像中不与其他骨骼组织有链接,这就使得在传统的二维区域生长或其他算法下,它将不会被标记为骨骼组织像素.三维区域生长[8, 9, 10, 11, 12, 13, 14, 15, 16]有效地解决了这个问题,它在纵向的生长过程中可以找到骨骼组织在纵向上的末端并将其标记.
基本的三维区域生长算法就是把图像分割中经典的二维区域生长算法扩展到三维空间中.将在二维空间中搜索使用的8-邻域在三维区域生长算法中扩展至26-邻域,26-邻域中心像素点周围上层影像和下层影像中的各9个像素和同层影像中的8个相邻像素点,构成了该中心像素点的9+9+8=26个邻域像素集合.
本文结合头部CTA影像特征从三维区域生长初始种子点确定和对生长判定条件这两个方面对三维区域生长方法进行改进.改进后的三维区域生长算法初始种子的确定采用了基于大量临床实验和CTA医学影像特点的先验知识,全自动地有效地确定初始种子点,不需要人为手工干扰.对于生长判定条件的改进是基于概率论的相关知识和生理学相关知识而进行的.
1.1 改进的初始种子点选取方法传统的区域生长算法的初始种子点的选取往往需要人工在图像中通过点击操作或者其他方法进行设定,种子点的选取是否准确在很大程度上影响算法的性能,这就要求相关的操作人员临床经验要很高并降低了算法的鲁棒性.本文提出改进的种子点选取方法,不需要人工交互,并且是根据不同的CTA影像数据集合可确定不同的骨骼像素的峰值.
图 1是一组实例中的不同组织的灰度分布直方图.图 1展示了一组典型的头部CTA数据集的主要3种类型的像素、血管组织像素、骨骼组织像素和其他组织像素统计直方图,其中横坐标是灰度值,纵坐标是该灰度值相对应的组织像素点数量.如图 1所示,肌肉组织、脂肪组织等其他组织的灰度值较低,主要集中于灰度范围在500~1 200之间,这些组织的灰度区间和骨骼血管组织的灰度区间交叉覆盖很少,基本可以忽略不计.CTA医学影像之所以不能单纯地通过灰度信息来分割,主要是由于血管组织的灰度范围1 250~1 750和骨骼组织的灰度范围1 200~2 800严重重叠,由图 1可以看出,骨骼组织像素的灰度范围基本将血管组织像素的灰度范围全部覆盖,因此不能单纯地通过灰度信息来区分血管组织和骨骼组织.但是,通过对大量临床头部CTA医学影像数据集的观察,发现如图 1黑色十字中心点所示的峰值点(该点为骨骼组织像素的数量的峰值),位于该点右侧的像素点全部都是骨骼像素点.全部的非骨骼像素点的灰度范围都位于该峰值左侧.因此,得到了快速全自动确定CTA头部数据集三维区域生长算法初始种子点的方法.
 |
| 图 1 头部影像灰度直方图Fig. 1 Histogram of head image in grayscale |
| 图选项 |
首先,通过头部分层算法,得到CTA头部影像数据集中躯干像素的集合,并统计全体躯干像素的灰度直方图.然后,根据该直方图计算出骨骼像素的临界峰值,记录该峰值的灰度数值.最后,扫描影像数据,将所有大于该峰值的像素点记录为三维区域生长算法的初始种子点.通过十几组临床数据的测试,确定了这种基于CTA影像数据灰度直方图特点的快速自动种子点选取方法的有效性和准确性.
1.2 改进的生长判定条件传统的区域生长算法判定条件往往不能充分地考虑局部信息，并且没有针对人体解剖学特性做出任何相应的响应.由于人体组织具有很强的连续性和封闭性,因此,本文针对这些特性对三维区域生长算法的判定条件做出相应的改进,使得改进后的三维区域生长算法可以更加良好地解决医学影像领域的问题.
制订针对于CTA影像特点的生长判定条件并引入了概率方法来辅助确定邻域点是否为骨骼组织像素点.通常的区域生长算法的判定条件一般有两种:一种是判定新加入点的灰度值与原区域的平均灰度是否在可接受范围内,数学形式如式(1);另一种是判定新加入点的灰度值和邀请它加入的原种子点的灰度绝对值之差是否在可接受范围内,数学表达式如式(2):
对于表达式(1)所示的第1种区域生长算法,存在的问题是某些灰度值很高的致密骨骼像素点有可能因为远离平均灰度值而不能被加入到骨骼组织区域中.对于式(2)所示的第2种区域生长算法,φ的取值十分困难,如果φ的取值过大很有可能会使得区域沿着一条灰度均匀下降的路线形成过生长.如果φ的取值过小,可能使得相同类别之间相邻像素的灰度变化过大而不能被加入到目标区域中.因此,对三维区域生长算法的判定条件进行了适合CTA头部影像数据的合理修改.
将判定条件设计成两部分:第1部分是根据新加入点的像素的灰度值来进行判定的；第2部分是根据该点周围像素点的归属来进行判定的.例如,有理由相信,如果一个新的需要判定的像素的周围26-邻域中有大量的骨骼像素点存在,则这个新的像素点很有可能就是骨骼像素点.其极端的例子就是,加入一个新的像素点具有很低的灰度值,但是它在某一个平面上的8-邻域全部都已经被判定为骨骼像素点,则该点应该属于灰度值较低的骨髓组织,该点应该被加入到骨骼组织区域中.因此,设计了新的加入判定规则:
式中:m为判定点的26-邻域内骨骼像素点的个数;n为当前已经被判定为骨骼区域的像素点的集合;λ为参数,用来调节区域生长限制调节的严格程度.
从式(3)中可以清楚地了解到,如果当前判定点的周围具有较多的已经判定为骨骼组织的像素点,则m的值越大,式(3)右侧的数值就会越小,因此,X0就更容易满足判定条件.同时,加入了另一个判定条件,即:如果在三维平面中的某一个平面下,该判定点的8-邻域内有7个以上像素点已经被判定为骨骼像素点,则该点就直接被判定为骨骼组织.因为,在同一平面内,8-邻域中的7个或以上的邻域点一定可以包围中心点及骨骼组织的中心点一定是高灰度的骨骼组织或低灰度的骨髓组织.
如图 2骨髓像素点示意图所示,矩形框区域中心像素点具有和周围软组织相近的较低的灰度值,如果按照原始的区域生长算法,它们将不会被加入到骨骼组织区域中.如果按照本文设计的改进的三维区域生长算法判定条件,它们将会被当作骨骼像素点而正确地加入到骨骼区域中.
 |
| 图 2 骨髓像素点示意图Fig. 2 Schematic diagram of marrow pixel points |
| 图选项 |
通过改进的三维区域生长算法,已经可以在很大程度上得到CTA影像数据集中的大多数骨骼组织部分.然后,通过形态学的膨胀操作来增大得到的骨骼组织外轮廓,这样可以使得该轮廓处于真正的骨骼组织轮廓的外侧并且很接近于骨骼组织轮廓.
2 改进的主动轮廓模型算法结合医学影像的特征和功能分析的特殊需求对经典Snake模型[17, 18]进行了两个方面的改进：第一个方面是对Snake模型的初始轮廓的确定,不采用人工交互式方式确定初始轮廓,而是采用改进的三维区域生长算法的结果经过膨胀操作后直接作为改进的Snake模型的初始轮廓,这一改进既减少了算法对于人工交互的要求,也极大地提高了初始轮廓的准确性,这样可以提高Snake算法的稳定性和健壮性;第二个方面是在能量函数中引进了图像力这一参数项,这一项更加有效地利用了CTA医学影像的骨骼组织边缘灰度信息,可以更好地引导主动轮廓向骨骼组织边缘移动.
2.1 改进的初始轮廓确定初始轮廓的确定又在很大程度上影响Snake模型的性能,但经典Snake模型初始轮廓的确定往往比较困难并且不易实现自动化[19].本文通过改进的三维区域生长算法来确定Snake模型的初始轮廓.首先通过改进的三维区域生长算法对医学影像中的骨骼组织进行预分割处理,然后通过形态学操作使得它更加满足Snake模型的初始外轮廓要求,进而解决Snake模型的初始轮廓确定的自动化问题.
由于改进的三维区域生长算法得到的骨骼区域已经比较接近真实的骨骼区域分割结果,当它经过有限次数的膨胀操作之后,它将完全地将真正的骨骼区域完全包括在它的外边界之中,这样就能较好地满足Snake模型对初始轮廓位置的要求.
2.2 改进的能量函数Snake算法主要思想是通过定义一条初始的能量函数曲线,并将其初始位置定义在待分割物体的边界外围.通过求解该能量函数的极小值,能量函数曲线经过不断地演化,最终将收敛于目标区域的边界处,得到目标图像的轮廓.主动轮廓模型的定义为
该曲线是由一组控制点组成的,这些控制点首尾相接构成了Snake模型的轮廓线.x(s)和y(s)表示控制点在影像中的平面坐标位置,其中,s是描述边界的傅里叶变换形式的自变量.在Snake模型的控制点集合上定义模型的能量函数为
式中:Eint为轮廓线的内部能量;Eext为外部能量.
主动轮廓模型Snake的能量Esnake由反应Snake曲线弹性势能和弯曲势能的内部能量Eint和反应当前图像特征的外部能量Eext两部分组成.内部能量通常可以写为
式中:一阶微分项V′(s)表示对应能量曲线的张力;二阶微分项V″(s)表示对应能量曲线的曲率;α和β则是相应的势能系数,可以是随曲线变化的表达式,也可以是常量,通常这两个系数用常量表示.在实际使用主动轮廓模型时,经常需要确定轮廓上的多个离散的控制点,α和β的数值决定能量曲线在这些控制点的位置上可以伸展和弯曲的程度.Eint的第1项可以控制曲线的应力,即在控制点的总数量不变的情况下,减小该项意味着能量曲线将不断地收缩.例如,当β的取值较大时,会使二阶微分项的权重增加,这会使曲线更加趋于平滑,对于无其他条件干扰的封闭的轮廓,它将进化成为一个光滑的圆形.假设,Vi－1、Vi和Vi+1是主动轮廓曲线上相邻的3个点,曲线的一阶微分项和二阶微分项可以通过式(7)和式(8)近似地求解:
外部能量Eext是定义在整幅图像上的标量函数,用于反映图像中的某些特性,例如边缘特征等.基于边缘的主动轮廓模型是研究历史最长并且应用最广泛的主动轮廓模型,该模型采用图像的灰度梯度信息作为外部能量,可以使主动轮廓曲线在内部能量和外部能量的共同引导下,向灰度梯度较大的目标区域的边界移动.基于边缘的主动轮廓模型的外部能量Eext是用图像的灰度梯度模的负值丨Δ I丨构成的图像势能.它形成图像势能立场,驱使能量曲线向图像灰度梯度高的区域边缘移动.因此,外部能量可表述为
式中:Δ为梯度算子;I为矩阵;λ为外部能量的权值系数.因为,参数λ为正数,所以要想使Esnake(V)取得较小的数值,则Eext(V)也需要向数值小的方向移动,即轮廓线会向梯度数值较大的区域边缘移动.外部能量就这样驱动Snake模型的轮廓线向区域的边界移动.
本文在传统的主动轮廓模型中加入了影像能量项Eima(V(s))，影像能量项表示该像素点作为边缘的资质,如果该像素点的影像能量数值较大,则说明该点比较满足边缘条件.改进后的主动轮廓模型的能量函数为
式(10)中的Eint(V(s))和Eext(V(s))与传统的主动轮廓模型中的内部能量、外部能量完全一致,内部能量来表示曲线的弹性势能和弯曲势能,外部能量来驱使轮廓曲线向目标区域移动.影像能量的计算表达式为
式中：
其中:G为原始图像I与5×5的滤波器K的卷积;γ为调节系数,用来调控影像能量对改进的主动轮廓模型的影像的大小.滤波器K可以用来检测图像中的边缘,它具有各向同性的特性,因此,它可以无差别地检测图像中任何方向的边缘.它在灰度不变区域和灰度均匀变化区域均无响应,即在这些区域,G的值接近0.而在边缘区域,它将在灰度值较低的边缘外侧得到较高的数值,因此,当它和γ相乘的时候,将得到一个绝对值较大的负值,它将驱动轮廓曲线向这里移动.γ和α,β,λ相似,它们的数值需要通过实验确定.
3 算法实验与结果分析本实验采用CTA临床数据20组.它们所提供的CTA影像数据均符合DICOM3.0标准,影像数据分辨率为512像素×512像素,像素点间距为0.625 mm,影像层间距为0.900 mm,影像层数在400~700层之间.图 3给出了根据本文设计的头部骨骼组织提取算法得到的实验结果在三维可视化下的视图效果.
 |
| 图 3 头部骨骼三维可视化效果图Fig. 3 Head skeleton 3D visualization figure |
| 图选项 |
根据图 3的实验效果图可以看出,本文提出的头部骨骼组织提取算法可以准确有效地提取出CTA影像数据集中的骨骼组织部分,并且由于引入了概率理论和人体生理解剖学特征到三维区域生长算法中,使得骨骼外边缘像素的提取判断更加准确.
精准的骨骼组织提取结果为之后的三维重建工程提供了良好的表面信息和边界信息,进而使得经过三维重建之后的骨骼图像可以具有非常接近其真实生理形态下的光滑外表面.
本文设计的改进的Snake模型与GVF-Snake模型实验结果对比图如图 4所示.
 |
| 图 4 不同的Snake算法对比图Fig. 4 Comparison between different Snake algorithms |
| 图选项 |
在图 4中,将会看到,由于对骨骼组织边界提取不准确的传统的基于阈值的骨骼提取算法不能提供良好的边界信息,使得同样的三维重建技术重建出来的骨骼模型具有相对粗糙的表面和相对模糊的细节信息.
图 4给出了分别经过GVF-Snake算法处理和本文改进的Snake算法处理后的头部CTA影像.影像中作为实例被标记出来的两块骨骼结构属于下颌骨的一部分.从图 4(b)中可以看出,GVF-Snake模型,可以陷入到骨骼区域内部的深度凹陷区域,但是,它显然把其中被硬质骨骼组织包围的内部组织排除在外.它会去掉一些未被高灰度的硬质骨骼完全包围的低密度骨骼组织,与其对应的,本文所采用的改进的Snake算法则完全地保留住了这部分的低密度骨骼组织,但是,它也相应地多保留了一些被骨骼半包围的软组织等不属于骨骼组织的像素.
在实际的医学应用领域,医生更希望得到像图 4(a)中所示的分割提取结果.因为图 4(a)的中间位置左侧和右侧所示的下颌骨柱状骨骼组织的分割提取结果良好地保留了硬质骨骼内部软组织、骨髓组织等信息,并且可以更好地服务与后续的三维可视化工程中,使得该部分的骨骼提取结果在三维可视化中展现出非常良好的外表面.因为无论是临床医师还是研究****都更需要良好的外侧的观察视角,而不是从区域内部的视点出发去观察这个区域的细节.这就好像在观察人体骨架时,是从外侧观察它的外表面所具有的形态和特征,而不是从内部骨髓中去观察骨骼组织内部组成细节.
由于克服了基于传统阈值的骨骼提取算法的种种不良因素,本文提出的医学影像头部骨骼提取算法具有更加良好的操作性、稳定性和健壮性.通过成功地去除了人工种子点选取步骤,本算法有效地去除了由于医疗器械操作人员误操作而产生错误种子点进而导致错误的提取结果这一常见问题,因而大大提高了医疗器械的使用性能,得到了操作试用者的广泛好评.同时,对比示意图 5展示了:本文提出的骨骼组织分离提取算法由于采用了区域圆形度判定等预处理步骤,提前去除了一定数量的明显具有血管组织特征的血管组织,这就为之后的骨骼提取算法打下了良好的基础,使得本文提出的骨骼分割提取算法对血管组织的误提取数量比较少,反映在对比示意图中比较明显的一处就是,图 5(a)在上颌骨附近比图 5(b)具有更加混乱的组织提取结果.
 |
| 图 5 头部骨骼分割算法对比示意图Fig. 5 Schematic diagram of comparison between different head bone segmentation algorithms |
| 图选项 |
4 结 论本文结合计算机断层血管造影影像的图像特点,提出基于改进的三维区域生长和改进的主动轮廓模型相结合的头部骨骼提取算法.
1) 结合CTA影像的灰度特点,提出了三维区域生长的快速的骨骼区域种子点提取方法,并对三维区域生长算法的判定条件进行了适合CTA头部影像数据的合理修改.
2) 对于经典的Snake模型,本文结合医学影像的特征和功能对其进行了两个方面的改进:一方面是对Snake模型的初始轮廓确定的改进,这一改进减少了算法对于人工交互的要求;另一方面是在能量函数中引进了图像力这一参数项,使得该模型可以更好地解决当前的问题.
3) 最后给出实验结果并和传统的基于阈值和区域生长的骨骼组织分割算法的结果进行对比分析,证实本文提出的骨骼组织分割提取算法效果良好.
参考文献
| [1] | Witkin A P.Scale-space filtering[C]//Proceeding International Joint Conference Artificial Intelligence.San Francisco, CA:Morgan Kaufmann Publishers Inc., 1983, 2:1019-1022. |
| [2] | DeanBidgood W.Understanding and using DICOM, the data interchange standard for biomedical imaging[J].Journal of the American Medical Informatics Association, 2007, 4(3):199-202. |
| Click to display the text | |
| [3] | Freeman E, Robson E, Bates B, et al.Head first design patterns[M].Cambridge:O'Reilly Media, 2005:35-40. |
| [4] | Behrens T, Rohr K, Stiehl H S.Robust segmentation of tubular structures in 3-D medical images by parametric object detection and tracking[J].IEEE Transactions on Medical Imaging, 2012, 21(8):924-933. |
| Click to display the text | |
| [5] | Rinkel G J E, Djibuti M, Algra A, et al.Prevalence and risk of rupture of intracranial aneurysms:A systematic review[J].Stroke, 2008, 29(1):251-256. |
| Click to display the text | |
| [6] | Wink O, Niessen W J, Viergever M A.Fast delineation and visualization of vessels in 3-D angiographic images[J].IEEE Transactions on Medical Imaging, 2000, 19(4):337-346. |
| Click to display the text | |
| [7] | Udupa J K, Samarasekara S.Fuzzy connectedness and object definition:Theory, algorithms, and applications in image segmentation[J].IEEE Transactions on Graphical Models and Image Processing, 1996, 58(3):246-261. |
| Click to display the text | |
| [8] | Chan F T, Vese L.Active contours without edges[J].IEEE Transactions on Image Processing, 2001, 10(2):266-277. |
| Click to display the text | |
| [9] | Kass M, Witkin A.Snakes:Active contour models[J].International Journal of Computer Vision, 2007, 1(4):259-268. |
| [10] | Terzopoulous D, Fleischer K.Deformable models[J].The Visual Computer, 2008, 79(1):12-49. |
| [11] | Li P H, Zhang T W.Review on active contour model (snake model)[J].Journal of Software, 2000, 11(6):751-757. |
| Click to display the text | |
| [12] | Li J, Yang X, Shi P.http://html.rhhz.net/BJHKHTDXXBZRB/A%20fast%20level%20set%20approach%20to%20image%20segmentation%20based%20on%20Mumford-Shah%20model[J].Journal of Computer, 2002, 25(11):1175-1183. |
| Click to display the text | |
| [13] | Cremers D, Kohlberer T, Schnorr C.Shape statistics in kernel space for variation image segmentation[J].Pattern Recognition, 2003, 36(9):1929-1943. |
| Click to display the text | |
| [14] | Stephen M P.Image processing and analysis at IPAG[J].IEEE Transaction on Medical Imaging, 2003, 22(12):1505-1518. |
| Click to display the text | |
| [15] | Pardo X M, Radeva P, Cabello D.Discriminated snakes for 3D reconstruction of anatomical organs[J].Medical Image Analysis, 2003, 7(3):293-310. |
| Click to display the text | |
| [16] | Montagnat J, Delingette H, Ayache N.A review of deformable surfaces:Topology, geometry and deformation[J].Image and Vision Computing, 2001, 19(14):1023-1040. |
| Click to display the text | |
| [17] | 李天庆, 张毅, 刘志, 等.Snake模型综述[J].计算机工程与应用, 2004, 40(11):43-45. Li T Q, Zhang Y, Liu Z, et al.A review on Snake models[J].Computer Engineering, 2004, 40(11):43-45(in Chinese). |
| Cited By in Cnki (118) | |
| [18] | 倪雅樱, 胡志忠.基于Snake模型的医学图像分割[J].上海生物医学工程, 2007, 28(3):165-170. Ni Y Y, Hu Z Z.Medical image segmentation based on snake model[J].Shanghai Journal of Biomedical Engineering, 2007, 28(3):165-170(in Chinese). |
| Cited By in Cnki (5) | |
| [19] | 王蒙, 戴亚平, 王庆林.一种新的FAST-Snake跟踪方法[J].自动化学报, 2014, 40(6):1108-1115. Wang M, Dai Y P, Wang Q L.A novel FAST-Snake object tracking approach[J].Acta Automatica Sinica, 2014, 40(6):1108-1115(in Chinese). |
| Cited By in Cnki |
