CN102129686B - 一种基于体素级轮廓粗定位的亚体素表面检测方法 - Google Patents

Abstract

本发明公开了一种基于体素级轮廓粗定位的亚体素表面检测方法，对待处理的连续切片图像进行滤波去噪；生成标记模板与切片图像像素相对应；分配一个检测空间；按一定顺序将切片图像读入检测空间；对切片图像进行阈值分割并按八方向搜索策略进行轮廓提取，得到体素级轮廓点；逐点处理体素级轮廓点，输出亚体素表面轮廓点；判断切片图像是否全部处理完；将所有切片图像处理所得的亚体素表面轮廓点合并，得到亚体素精度的物体表面点云。本发明可从含有较多噪声和伪影的实际CT连续切片图像中提取完整的、高精度的物体表面点云，并通过标记模板大大减少了计算量。

Description

一种基于体素级轮廓粗定位的亚体素表面检测方法

技术领域

本发明属于CT图像表面重建技术领域，涉及基于CT图像体素级轮廓粗略定位进行高精度快速表面检测的一整套解决方案。

背景技术

点云是物体表面检测的核心数据，获取点云的方法有很多种，其中基于锥束CT技术的获取方式凭借其诸多优点得到了广泛应用。计算机断层成像技术(ComputedTomography，CT)是通过对物体不同角度的射线投影重建而获取被测物体断层图像信息的成像技术。锥束CT利用锥形束X射线源和面阵探测器采集被测物体的投影数据，与传统二维CT相比，锥束CT一次扫描即可重建出数百甚至上千个连续的断层图像，具有射线利用率高、切片连续、切片内和切片间空间分辨率相同、精度高等特点。断层图像通常也称为切片图像。

表面检测的两个基本要求是完整性和精度，检测速度也是实际应用中的一个关键性问题。根据其精度，目前实际应用的表面检测方法分为体素级和亚体素级两类，体素级检测结果完整性高，一般可达95％以上，但精度较低，一般在0.5个体素左右；亚体素级检测结果精度高，一般可达0.2～0.3个体素，但完整性不足，且时间复杂度较高。因此要想在保证完整性的前提下取得较高的提取精度，必须结合体素级检测方法的完整性优势和亚体素检测方法的高精度优势。

从目前的文献资料来看，体素级表面检测方法主要借助于经典边缘检测算子和图像分割算法。基于边缘检测算子的一类方法的检测精度受图像信噪比影响较大；在附加合理去噪算法的条件下，基于图像分割的一类方法一般可以取得完整性较好的体素级表面轮廓。高精度的亚体素级表面检测方法主要分为以下两类：一是基于Facet模型的表面检测方法，该类方法通过三次多项式拟合Facet模型，再根据梯度方向上的二阶导数零交叉点确定亚体素表面点位置。王凯、张定华、刘晶等人在《计算机辅助设计与图形学学报》(2007，19(9)：1100-1106)的文章“采用3D Gaussian Facet模型的亚体素表面检测”中扩展了原始Haralick Facet模型并建立了三维Gaussian Facet模型及其计算公式，该方法能有效降低邻近边缘干涉对检测结果的影响，能提取尺寸更小的结构边缘；二是基于矩的表面检测方法，该类方法通过建立理想边缘的参数化模型，并根据矩不变量架设理想边缘灰度分布和离散图像灰度分布之间的关系，由不变关系建立方程即可确定理想边缘的参数。刘亚威在其硕士学位论文《空间矩亚体素图像测量算法的研究》中给出了基于空间矩的2D亚体素边缘检测修正算法，3D校正算法推导过于复杂，尚无较好的修正方法。对于仿真CT图像数据，现有算法可以获得平均误差小于0.2个体素的完整表面轮廓点，但对于包含噪声和模糊效应的实际CT图像，亚体素表面检测方法会出现参数选择困难、轮廓提取不完整、伪边严重等问题。

采用传统亚体素表面检测方法对锥束CT重建的三维图像(即连续切片图像的层叠)进行轮廓提取时，需要对空间内每个体素进行拟合计算，时间复杂度非常高。大量的含噪背景区域对表面轮廓提取并没有贡献，并且导致不必要的计算，因此实际应用中需要进行方法加速。亚体素表面检测的加速方法一般有基于可分滤波器的加速方案和基于感兴趣区域(Region of Interest，ROI)的加速方案。ROI法是高精度表面检测中一种非常实用的降低计算量的方法，它通过缩小图像中待处理区域的范围，仅处理感兴趣区域内的体素，以达到缩减图像总处理时间的目的。王凯、张定华、张顺利等人在《中国图形图像学报》(2009，14(2)：329-333)的文章“对基于facet模型的表面检测的加速技术研究”利用递归的可分滤波器结合分段包围盒的加速方案将原算法的时间复杂度从O(N²·(2M+1)³)降低到O(N²)，其中N为切片图像尺寸，M为拟合模板尺寸。

发明内容

为了克服现有技术将亚体素表面检测方法应用于实际锥束CT三维图像时出现的参数选择困难、轮廓提取不完整、伪边严重和时间复杂度高等不足，本发明提供一种基于体素级轮廓粗定位的亚体素表面检测方法，在确保提取的表面轮廓完整性的前提下，利用3D Gaussian Facet模型进一步精确定位表面点，达到快速高精度的提取被检测物体完整表面轮廓点的目的。

本发明解决其技术问题所采用的技术方案包括以下步骤：

(1)对待处理的连续切片图像进行滤波去噪；

(2)设待处理的连续切片图像大小均为M×N像素，生成一个M×N大小的数组，数组元素为结构体，结构体成员包括计算标识、表面体素标识、基于3D Gaussian Facet模型的表面检测方法中的拟合系数、梯度方向、梯度模值G和二阶导数过零点距离ρ，该数组称为标记模板，与切片图像像素相对应；

(3)设基于3D Gaussian Facet模型的表面检测方法的拟合邻域大小为H×H×H，分配一个可存储H幅切片图像的内存空间，该空间称为检测空间，并设W＝int(H/2)；

(4)设当前待处理的切片图像为第K幅，若K＝1，将第1幅切片图像读入检测空间，并将该切片图像复制W份作为与其相邻的下面W层切片图像，然后将与该切片图像相邻的上面W层切片图像读入检测空间；否则，读入第(K+W)幅切片图像替换检测空间中的第(K-W-1)幅切片图像，若(K+W)大于顶层切片图像序号，则复制顶层切片图像作为需读入的切片图像；

(5)对第K幅切片图像进行OTSU阈值分割，得到其二值图像，并按八方向搜索策略对该二值图像进行轮廓提取，得到体素级轮廓点；

(6)以上步所得体素级轮廓点为粗定位基准，逐点处理体素级轮廓点，输出亚体素表面轮廓点；

(7)判断：若切片图像未全部处理完，则K＝K+1，转第(4)步；否则，继续下一步；

(8)将所有切片图像处理所得的亚体素表面轮廓点合并，得到亚体素精度的物体表面点云。

在上述方法的第(6)步中，逐点处理体素级轮廓点的具体步骤如下：

1)对标记模板数组元素结构体成员全部置“否”或0；

2)设当前待处理的体素级轮廓点为J，其在切片图像中的位置为(x，y)，判断：若标记模板(x，y)位置的计算标识为“否”，则在检测空间中按基于3D Gaussian Facet模型的表面检测方法，计算J的拟合系数、梯度方向、梯度模值G、二阶导数过零点距离ρ，并将结果存入标记模板(x，y)位置的数组元素，相应的计算标识置为“是”；

3)判断：若标记模板(x，y)位置的表面体素标识为“否”，则在J的|ρ|＜0.5时，将其表面体素标识置为“是”，转第7)步；否则，设置邻域边长L＝3；

4)在检测空间中按基于3D Gaussian Facet模型的表面检测方法，分别计算J的L×L邻域内，在标记模板中计算标识为“否”的体素的拟合系数、梯度方向、梯度模值G、二阶导数过零点距离ρ，并将结果分别存入标记模板中各体素所对应的数组元素，相应的计算标识置为“是”；

5)将J的L×L邻域内所有体素按其梯度模值G从大到小排序，并查找第1个|ρ|＜0.5且表面体素标识为“否”的体素，若有，则将其相应表面体素标识置为“是”，转第7)步；若无，则判断当前邻域边长，若L＝3，则设置L＝5并转第4)步，否则继续下一步；

6)将标记模板(x，y)位置的表面体素标识置为“是”，并置ρ＝0；

7)按基于3D Gaussian Facet模型的亚体素表面轮廓点计算方法，计算当前表面体素标识为“是”的体素的亚体素表面轮廓点；

8)判断：若体素级轮廓点未全部处理完，转第2)步处理下一个体素级轮廓点；否则，结束。

本发明的有益效果是：本发明方法结合了体素级表面检测方法和亚体素表面检测方法的优势，可从含有较多噪声和伪影的实际CT连续切片图像中提取完整的、高精度的物体表面点云，并通过标记模板大大减少了计算量，整体可获得3倍以上的加速比。

下面结合附图和实施例对本发明进一步说明。

附图说明

图1为本发明方法进行亚体素表面检测的实现流程图。

具体实施方式

对一组锥束CT重建的100层连续切片图像，在Intel Core II 2.33GHz处理器、2G内存的计算机上，应用本发明方法进行亚体素表面检测，执行以下步骤如下：

(1)采用非局部平均算法对待处理的连续切片图像进行滤波去噪；

(2)待处理的连续切片图像大小均为1024×1024像素，生成一个1024×1024大小的数组，数组元素为结构体，结构体成员包括计算标识、表面体素标识、基于3D GaussianFacet模型的表面检测方法中的拟合系数、梯度方向、梯度模值G、二阶导数过零点距离ρ，该数组称为标记模板，与切片图像像素相对应；

(3)设基于3D Gaussian Facet模型的表面检测方法的拟合邻域大小为5×5×5，分配一个可存储5幅切片图像的内存空间，该空间称为检测空间，并设W＝int(5/2)＝2；

(4)设当前待处理的切片图像为第K幅，若K＝1，将第1幅切片图像读入检测空间，并将该切片图像复制W份作为与其相邻的下面W层切片图像，然后将与该切片图像相邻的上面W层切片图像读入检测空间；否则，读入第(K+W)幅切片图像替换检测空间中的第(K-W-1)幅切片图像，若(K+W)大于顶层切片图像序号，则复制顶层切片图像作为需读入的切片图像；

(5)对第K幅切片图像进行OTSU阈值分割，得到其二值图像，并按八方向搜索策略对该二值图像进行轮廓提取，得到体素级轮廓点；

(6)以上步所得体素级轮廓点为粗定位基准，按一定规则逐点处理，输出亚体素表面轮廓点，具体步骤为：

1)对标记模板数组元素结构体成员全部置“否”或0；

2)设当前待处理的体素级轮廓点为J，其在切片图像中的位置为(x，y)，判断：若标记模板(x，y)位置的计算标识为“否”，则在检测空间中按基于3D Gaussian Facet模型的表面检测方法，计算J的拟合系数、梯度方向、梯度模值G、二阶导数过零点距离ρ，并将结果存入标记模板(x，y)位置的数组元素，相应的计算标识置为“是”；

3)判断：若标记模板(x，y)位置的表面体素标识为“否”，则在J的|ρ|＜0.5时，将其表面体素标识置为“是”，转第7)步；否则，设置邻域边长L＝3；

4)在检测空间中按基于3D Gaussian Facet模型的表面检测方法，分别计算J的L×L邻域内，在标记模板中计算标识为“否”的体素的拟合系数、梯度方向、梯度模值G、二阶导数过零点距离ρ，并将结果分别存入标记模板中各体素所对应的数组元素，相应的计算标识置为“是”；

5)将J的L×L邻域内所有体素按其梯度模值G从大到小排序，并查找第1个|ρ|＜0.5且表面体素标识为“否”的体素，若有，则将其相应表面体素标识置为“是”，转第7)步；若无，则判断当前邻域边长，若L＝3，则设置L＝5并转第4)步，否则继续下一步；

6)将标记模板(x，y)位置的表面体素标识置为“是”，并置ρ＝0；

7)按基于3D Gaussian Facet模型的亚体素表面轮廓点计算方法，计算当前表面体素标识为“是”的体素的亚体素表面轮廓点；

8)判断：若体素级轮廓点未全部处理完，转第2)步处理下一个体素级轮廓点；否则，结束。

(7)判断：若切片图像未全部处理完，则K＝K+1，转第(4)步；否则，继续下一步；

(8)将所有切片图像处理所得的亚体素表面轮廓点合并，得到亚体素精度的物体表面点云。

对于这组切片图像，表1给出了随机抽取的2幅切片图像分别采用3D GaussianFacet(3GF)表面检测方法和本发明方法处理后在表面轮廓完整性方面的比较，可见本发明方法获得更好的完整性；表2给出了3GF方法与本发明方法在提取精度方面的比较，可见本发明的提取精度明显较高；表3给出了3GF方法(可分滤波器+ROI加速)和本发明方法在计算速度上的比较，可见本发明方法获得了可观的加速比。实际处理结果验证了本发明方法的可行性和有效性。

表1亚体素表面检测完整性比较

	3GF方法	本发明方法
			断边	3	0
伪轮廓	12	0

表2亚体素表面检测精度比较

表3亚体素表面检测速度比较

Claims (1)

1.一种基于体素级轮廓粗定位的亚体素表面检测方法，其特征在于包括以下步骤：

(1)对待处理的连续切片图像进行滤波去噪；

(2)设待处理的连续切片图像大小均为M×N像素，生成一个M×N大小的数组，数组元素为结构体，结构体成员包括计算标识、表面体素标识、基于3D GaussianFacet模型的表面检测方法中的拟合系数、梯度方向、梯度模值G和二阶导数过零点距离ρ，该数组称为标记模板，与切片图像像素相对应；

(3)设基于3D Gaussian Facet模型的表面检测方法的拟合邻域大小为H×H×H，分配一个可存储H幅切片图像的内存空间，该空间称为检测空间，并设W＝int(H/2)；

(4)设当前待处理的切片图像为第K幅，若K＝1，将第1幅切片图像读入检测空间，并将该切片图像复制W份作为与其相邻的下面W层切片图像，然后将与该切片图像相邻的上面W层切片图像读入检测空间；否则，读入第(K+W)幅切片图像替换检测空间中的第(K-W-1)幅切片图像，若(K+W)大于顶层切片图像序号，则复制顶层切片图像作为需读入的切片图像；

(5)对第K幅切片图像进行OTSU阈值分割，得到其二值图像，并按八方向搜索策略对该二值图像进行轮廓提取，得到体素级轮廓点；

(6)以上步所得体素级轮廓点为粗定位基准，逐点处理体素级轮廓点，输出亚体素表面轮廓点，具体步骤如下：

1)对标记模板数组元素结构体成员全部置“否”；

2)设当前待处理的体素级轮廓点为J，其在切片图像中的位置为(x，y)，判断：若标记模板(x，y)位置的计算标识为“否”，则在检测空间中按基于3D GaussianFacet模型的表面检测方法，计算J的拟合系数、梯度方向、梯度模值G、二阶导数过零点距离ρ，并将结果存入标记模板(x，y)位置的数组元素，相应的计算标识置为“是”；

3)判断：若标记模板(x，y)位置的表面体素标识为“否”，则在J的|ρ|＜0.5时，将其表面体素标识置为“是”，转第7)步；在J的|ρ|≥0.5时，设置邻域边长L＝3；

4)在检测空间中按基于3D Gaussian Facet模型的表面检测方法，分别计算J的L×L邻域内，在标记模板中计算标识为“否”的体素的拟合系数、梯度方向、梯度模值G、二阶导数过零点距离ρ，并将结果分别存入标记模板中各体素所对应的数组元素，相应的计算标识置为“是”；

5)将J的L×L邻域内所有体素按其梯度模值G从大到小排序，并查找第1个|ρ|＜0.5且表面体素标识为“否”的体素，若有，则将其相应表面体素标识置为“是”，转第7)步；若无，则判断当前邻域边长，若L＝3，则设置L＝5并转第4)步，否则继续下一步；

6)将标记模板(x，y)位置的表面体素标识置为“是”，并置ρ＝0；

7)按基于3D Gaussian Facet模型的亚体素表面轮廓点计算方法，计算当前表面体素标识为“是”的体素的亚体素表面轮廓点；

8)判断：若体素级轮廓点未全部处理完，转第2)步处理下一个体素级轮廓点；否则，结束；

(7)判断：若切片图像未全部处理完，则K＝K+1，转第(4)步；否则，继续下一步；

(8)将所有切片图像处理所得的亚体素表面轮廓点合并，得到亚体素精度的物体表面点云。