一种CT影像体表提取方法及系统与流程

本发明涉及医学影像处理技术领域，特别涉及一种CT影像体表提取方法及系统。

背景技术：

体表分割属于图像分割的范畴，常用的图像分割方法主要包括以下几类：基于区域的方法，基于边缘的方法，基于分类器、聚类、统计学的方法。依据待分割对象的不同特性，所采用或设计的分割方法也千差外别。针对体表的分割，在理想情况下使用边界跟踪或合理的阈值分割等方法就可以得到体表区域的边缘。但是临床的医学影像中，由于人体内外结构的异质性、CT成像过程的畸变或噪声的影响及其它不确定干扰因素，使得实际处理的医学影像体表往往边缘模糊、信号强度不均匀甚至噪声较大。另外，人体体表的结构比较复杂而且个体外部形态和生理差异较大，这也给医学影像的体表分割带来了很大的困难。

Manning Wang和Zhijian Song在Pattern Recognition Letters提出了采用全自动或手动确定阈值的分割方法提取体表结构。Giovanni Fattori等在IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING采用了基于Marching Cube分割的体表提取方法，主要选择体表的信号强度值的上限和下限阈值作为体表分割的依据。中国发明专利201310232262.6中基于灰度阈值的体表皮肤分割方法，首先计算皮肤的灰度阈值，对图像二值化处理；然后采用形态学方法处理二值图像，得到皮肤的像素级三维网格，接着还原表面像素区域的灰度值，最后利用以移动立方体算法对灰度图像三维重建，从而得到精确的皮肤三维网格。

上述基于阈值分割或边缘灰度分割方法的主要不足之处在于：(1)人体体表结构的复杂性和形态的差异，CT成像过程中的信号衰减、噪声、畸变等因素，都可以使体表成像的差异性较大，从而使精确阈值确定难度很大且很难获得满意的体表分割效果。(2)皮肤的分割结果如果存在孔洞或者包括皮肤以外的粘连组织结构时，如何对体表进行自动精细修复及剔除体表外的杂质，这是影响体表分割效果的重要因素，但是现有文献和发明专利中并没有考虑或陈述。

技术实现要素：

有鉴于此，有必要提供一种能够实现三维体表的快速、精确和鲁棒性分割提取的CT影像体表提取方法。

为实现上述目的，本申请采用下述技术方案：

一种CT影像体表提取方法，包括下述步骤：

获取三维CT影像；

对所述CT影像体表进行体表粗分割；

对所述CT影像体表进行体表细分割；

输出经细分割后的CT影像体表影像；

其中，体表粗分割包括下述步骤：

分别选取所述三维CT影像中体表的边缘和外部的两个种子点，分别记为第一种子点和第二种子点；

基于所述第一种子点和第二种子点采用孤立连接区域增长的方法区分体表内和体表外，获取体表的三维二值图像；

其中，体表细分割包括下述步骤：

采用边缘跟踪算法对所述三维二值图像的轴位影像分别逐层追踪最外边缘的轮廓，获得每一层切片的边缘轮廓坐标；

遍历所述三维二值图像所有层面的轴位影像进行逐层轮廓跟踪提取，获得所有层面的轮廓坐标集合。

优选地，所述三维CT影像为胸腹部三维CT影像。

优选地，其中，体表细分割中，边缘跟踪算法为8连通区域的轮廓跟踪算法。

优选地，所述8连通区域的轮廓跟踪算法，包括下述步骤：

步骤S110：从上到下，从左到右的顺序扫描所述三维二值图像的轴位影像，寻找最外围的第一个边界起始点，记为A0，其中，A0是具有最小行或列值的边缘点；

步骤S120：定义一个扫描方向变量dir,该变量用于记录上一步中沿着前一个边界点到当前边界点的移动方向，8连通区域时其初始化时取值为dir＝7；

步骤S130：按逆时针方向搜索当前像素的3×3邻域，其起始搜索方向设定如下：对8连通区域，若dir为奇数取(dir+7)取余8；若dir为偶数为(dir+6)取余8；

步骤S140：在3×3邻域中搜索到的第一个与当前像素值相同的像素便为新的边界点An，同时更新变量dir为新的方向值；

步骤S150：如果An等于第二个边界点A1且前一个边界点An-1等于第一个边界点A0，则停止搜索，结束跟踪，否则重复步骤S130继续搜索；

步骤S160：由边界点A0、A1、A2、…、An-2构成的边界便为要跟踪的边界。

优选地，体表细分割中，在完成获得每一层切片的边缘轮廓坐标后，还包括下述步骤：

对每一层切片的边缘轮廓进行各向异性扩散滤波处理。

另外，本发明还提供了一种CT影像体表提取系统，包括：

影像获取模块，用于获取三维CT影像；

体表粗分割模块，对所述CT影像体表进行体表粗分割；

体表细分割模块，对所述CT影像体表进行体表细分割；

影像输出模块，输出经细分割后的CT影像体表影像；

其中，体表粗分割模块包括：

分别选取所述三维CT影像中体表的边缘和外部的两个种子点，分别记为第一种子点和第二种子点；

基于所述第一种子点和第二种子点采用孤立连接区域增长的方法区分体表内和体表外，获取体表的三维二值图像；

其中，体表细分割模块包括：

采用边缘跟踪算法对所述三维二值图像的轴位影像分别逐层追踪最外边缘的轮廓，获得每一层切片的边缘轮廓坐标；

遍历所述三维二值图像所有层面的轴位影像进行逐层轮廓跟踪提取，获得所有层面的轮廓坐标集合。

优选地，所述三维CT影像为胸腹部三维CT影像。

优选地，其中，体表细分割模块中，边缘跟踪算法为8连通区域的轮廓跟踪算法。

优选地，还包括滤波处理模块，用于对每一层切片的边缘轮廓进行各向异性扩散滤波处理

本发明采用上述技术方案，其有益效果在于：

本发明提出了一种CT影像体表提取方法及系统，采用粗-精结合的体表分割方法，粗分割采用基于孤立连接的区域增长区分体表内和体表外，从而大致获得体表结构，细分割过程针对粗分割三维体表的轴位轮廓进行修复和细节保持平滑处理，从而实现了三维体表的快速、精确和鲁棒性分割提取。

附图说明

图1为本发明提供的CT影像体表提取方法的步骤流程图。

图2为本申请实施例提供的输入CT影像。

图3为本申请实施例提供的种子点选择和孤立连接区域增长示意图。

图4(a)为二值化轴位影像，图4(b)为边缘增长结果。

图5为所述8连通区域的轮廓跟踪算法步骤流程图。

图6为二维轮廓各向异性扩散滤波处理示意图。

图7为体表分割结果示意图。

图8为本申请提供的一种CT影像体表提取系统示意图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

请参阅图1，为本发明提供的CT影像体表提取方法100的步骤流程图，包括下述步骤：

步骤S10：获取三维CT影像；

优选地，所述三维CT影像为胸腹部三维CT影像。可以理解，本申请中给出的三维CT影像还可以是其他部位。请参阅图2，为本申请实施例提供的输入CT影像。

步骤S20：对所述CT影像体表进行体表粗分割；

其中，体表粗分割包括下述步骤：

分别选取所述三维CT影像中体表的边缘和外部的两个种子点，分别记为第一种子点和第二种子点；

基于所述第一种子点和第二种子点采用孤立连接区域增长的方法区分体表内和体表外，获取体表的三维二值图像；

可以理解，本申请提供的孤立连接区域增长是一种改进型的区域增长法，通过分别设置一个在待分割区域(前景区域)的种子点和一个非待分割区域(背景区域)的种子点，并设定一个区域增长的最低阈值，通过不断区域增长待分割区域的种子点而不包含非待分割区域(背景区域)的种子点增长区域，从而获得理想的分割结果。在体表分割中，由于体表及体表内区域与体表外区域的灰度值相差很大，通过在体表及体表内选择第一个种子点而在体表外选择第二个种子点，孤立连接区域增长的分割结果可以尽可能的包含体表及体表内区域而不包含体表外区域。请参阅图3，为本申请实施例提供的种子点选择和孤立连接区域增长示意图。

本申请采用孤立连接区域增长做体表提取的粗分割主要有以下有点：(i)尽可能的在分割结果中包含体表及体表内区域，避免体表提取的部分区域缺失。(ii)很好的去除体表外区域的干扰，滤除了体表外的噪声或杂质。需要注意的是，体表内部的结构是相当复杂的，但是粗分割的目标并不是针对体表内部，而是在于区分体表外和体表内，使体表及体表内形成一个分割的整体。

步骤S30：对所述CT影像体表进行体表细分割；

其中，体表细分割包括下述步骤：

采用边缘跟踪算法对所述三维二值图像的轴位影像分别逐层追踪最外边缘的轮廓，获得每一层切片的边缘轮廓坐标；

遍历所述三维二值图像所有层面的轴位影像进行逐层轮廓跟踪提取，获得所有层面的轮廓坐标集合。

可以理解，边缘跟踪是指通过顺序找出目标区域边缘点从而实现感兴趣区域的分割。体表粗分割后的影像是三维二值图像，体表细分割过程就是从三维二值图像中提取最外围的边缘。由于粗分割结果中包含的体表及体表内区域为1，而体表外的背景区域为0，通过边缘追踪获得粗分割结果最外围的轮廓点集合，即可以获得体表的信息。请参阅图4(a)为二值化轴位影像，图4(b)为边缘增长结果(白色轮廓线所示，采用黑色增强显示效果)。

本发明通过对上述三维二值体数据的轴位影像分别逐层追踪最外边缘的轮廓，获得每一层切片的边缘轮廓坐标，然后将所有层面的轮廓坐标集合起来即获得了三维体表的曲面坐标。

优选地，边缘跟踪算法为8连通区域的轮廓跟踪算法。

请参阅图5，所述8连通区域的轮廓跟踪算法，包括下述步骤：

步骤S110：从上到下，从左到右的顺序扫描所述三维二值图像的轴位影像，寻找最外围的第一个边界起始点，记为A0，其中，A0是具有最小行或列值的边缘点；

步骤S120：定义一个扫描方向变量dir,该变量用于记录上一步中沿着前一个边界点到当前边界点的移动方向，8连通区域时其初始化时取值为dir＝7；

步骤S130：按逆时针方向搜索当前像素的3×3邻域，其起始搜索方向设定如下：对8连通区域，若dir为奇数取(dir+7)取余8；若dir为偶数为(dir+6)取余8；

步骤S140：在3×3邻域中搜索到的第一个与当前像素值相同的像素便为新的边界点An，同时更新变量dir为新的方向值；

步骤S150：如果An等于第二个边界点A1且前一个边界点An-1等于第一个边界点A0，则停止搜索，结束跟踪，否则重复步骤S130继续搜索；

步骤S160：由边界点A0、A1、A2、…、An-2构成的边界便为要跟踪的边界。

上述图像轮廓跟踪算法可以获取任一层面轴位二值化影像的体表轮廓，对三维体数据的所有层面轴位二值化影像进行逐层轮廓跟踪提取，则可以获得所有层面的轮廓坐标集合。

优选地，体表细分割中，在完成获得每一层切片的边缘轮廓坐标后，还包括下述步骤：

对每一层切片的边缘轮廓进行各向异性扩散滤波处理。

可以理解，在形态上，二值化影像的轮廓跟踪结果普遍存在阶梯状或陡峭的变化，这与实际体表的平滑形态是有差别的。主要原因是上述边缘跟踪过程是基于二值化影像的轮廓进行，容易受到影像数字离散化、噪声等因素的影像。因此，对所获得体表的平滑滤波处理显得尤为关键，需要强调的是，体表的平滑处理不能改变体表轮廓的原有变化形态，即首先需要保持轮廓的形态细节特征不变，同时平滑滤波处理噪声及离散化引起的不规则突变。

为此，本发明采用了对二维轮廓进行各向异性扩散滤波处理。在数字图像处理中，各向异性扩散滤波主要是用于进行图像平滑，各向异性扩散在平滑图像时可以尽可能保留图像边缘。请参阅图6为二维轮廓各向异性扩散滤波处理。

具体地，本发明并不是采用各向异性扩散滤波对二维或三维图像进行处理，而是使用各向异性扩散滤波对二维轮廓进行各向异性扩散滤波处理。这个处理过程，实质上是针对二维轮廓点集合的X坐标和Y坐标分别进行一维的各向异性扩散滤波处理。假定二维轮廓边缘离散点的坐标为

{(x_i，y_i：i＝0，...，n-1)}，那么各向异性扩散滤波处理就是分别针对轮廓离散点集合的{(x_i：i＝0，...，n-1)}和{(y：i＝0，...，n-1)}一维向量分别进行各向异性扩散滤波处理。

步骤S40：输出经细分割后的CT影像体表影像。

请参阅图7，为体表分割结果示意图。

请参阅图8，本申请还提供了一种CT影像体表提取系统，包括：

影像获取模块110，用于获取三维CT影像；

体表粗分割模块120，对所述CT影像体表进行体表粗分割；

体表细分割模块130，对所述CT影像体表进行体表细分割；

影像输出模块140，输出经细分割后的CT影像体表影像；

其中，体表粗分割模块120包括：

分别选取所述三维CT影像中体表的边缘和外部的两个种子点，分别记为第一种子点和第二种子点；

基于所述第一种子点和第二种子点采用孤立连接区域增长的方法区分体表内和体表外，获取体表的三维二值图像；

其中，体表细分割模块130包括：

采用边缘跟踪算法对所述三维二值图像的轴位影像分别逐层追踪最外边缘的轮廓，获得每一层切片的边缘轮廓坐标；

遍历所述三维二值图像所有层面的轴位影像进行逐层轮廓跟踪提取，获得所有层面的轮廓坐标集合。

还包括滤波处理模块150，用于对每一层切片的边缘轮廓进行各向异性扩散滤波处理。

本申请提供的CT影像体表提取系统，其详细的实现方式在上述已有详细描述，这里不再赘述。

本发明提出了一种CT影像体表提取方法及系统，采用粗-精结合的体表分割方法，粗分割采用基于孤立连接的区域增长区分体表内和体表外，从而大致获得体表结构，细分割过程针对粗分割三维体表的轴位轮廓进行修复和细节保持平滑处理，从而实现了三维体表的快速、精确和鲁棒性分割提取。

虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。