路径规划鼻部环境建模方法及装置与流程

1.本发明涉及医用图像处理的技术领域，尤其涉及一种路径规划鼻部环境建模方法，以及路径规划鼻部环境建模装置。


背景技术：

2.手术过程需避开鼻部周围重要组织。但由于鼻部环境复杂多变，鼻部环境地图难以构建。
3.人体或动物体不同组织对x线的吸收系数不同，进行电子计算机断层扫描后得到的三维空间的灰度分布，即为ct图像。在ct图像空间中，通常x轴为左右方向，与x轴垂直、将图像切分为左右两部分的切面为矢状面；y轴为前后方向，与y轴垂直、将图像切分为前后两部分的切面为冠状面；z轴为上下方向，将图像切分为上下两部分、与z轴垂直的切面为水平面(横断面)。


技术实现要素：

4.为克服现有技术的缺陷，本发明要解决的技术问题是提供了一种路径规划鼻部环境建模方法，其能够实现鼻部环境地图的精准构建。
5.本发明的技术方案是：这种路径规划鼻部环境建模方法，其包括以下步骤：
6.(1)对ct图像进行滤波去噪预处理；
7.(2)利用轴向区域生长法，结合边缘检测提取分割目标沿坐标轴方向的边缘点对，计算边缘点对之间距离，以距离最大边缘点对之间的像素点作为种子点，沿轴向生长至目标边缘，实现重要组织器官的快速分割；
8.(3)提取重要组织器官的表面轮廓作为障碍物和地图边界，利用改进的marching cubes算法进行三维表面重建。
9.本发明通过对ct图像进行滤波去噪预处理，轴向区域生长法，结合边缘检测提取分割目标沿坐标轴方向的边缘点对，计算边缘点对之间距离，以距离最大边缘点对之间的像素点作为种子点，沿轴向生长至目标边缘，实现脑组织、眼眶内组织、下鼻甲和鼻中隔等重要组织器官的快速分割，进而提取重要组织的表面轮廓作为障碍物和地图边界，利用改进的marching cubes算法进行三维表面重建，实现鼻部环境地图的精准建模，因此能够实现鼻部环境地图的精准构建。
10.还提供了路径规划鼻部环境建模装置，其包括：
11.预处理模块，其配置来对ct图像进行滤波去噪预处理；
12.轴向区域生长模块，其配置来利用轴向区域生长法，结合边缘检测提取分割目标沿坐标轴方向的边缘点对，计算边缘点对之间距离，以距离最大边缘点对之间的像素点作为种子点，沿轴向生长至目标边缘，实现重要组织器官的快速分割；
13.重建模块，其配置来提取重要组织器官的表面轮廓作为障碍物和地图边界，利用改进的marching cubes算法进行三维表面重建。
附图说明
14.图1是根据本发明的路径规划鼻部环境建模方法的流程图。
15.图2示出了用于检测被骨质包围软组织的边缘点对的边缘检测子。
具体实施方式
16.如图1所示，这种路径规划鼻部环境建模方法，其包括以下步骤：
17.(1)对ct图像进行滤波去噪预处理；
18.(2)利用轴向区域生长法，结合边缘检测提取分割目标沿坐标轴方向的边缘点对，计算边缘点对之间距离，以距离最大边缘点对之间的像素点作为种子点，沿轴向生长至目标边缘，实现重要组织器官的快速分割；
19.(3)提取重要组织器官的表面轮廓作为障碍物和地图边界，利用改进的marching cubes算法进行三维表面重建。
20.本发明通过对ct图像进行滤波去噪预处理，轴向区域生长法，结合边缘检测提取分割目标沿坐标轴方向的边缘点对，计算边缘点对之间距离，以距离最大边缘点对之间的像素点作为种子点，沿轴向生长至目标边缘，实现脑组织、眼眶内组织、下鼻甲和鼻中隔等重要组织器官的快速分割，进而提取重要组织的表面轮廓作为障碍物和地图边界，利用改进的marching cubes算法进行三维表面重建，实现鼻部环境地图的精准建模，因此能够实现鼻部环境地图的精准构建。
21.优选地，所述步骤(1)中，通过旋转变换将倾斜的ct图像调正，以原点为旋转中心，ct图像中任一体素绕x轴、y轴、z轴分别旋转α、β、γ角度，其位置坐标的旋转变换矩阵为：
[0022][0023]
按照右手系定则，大拇指指向坐标轴正向时四指方向为旋转的正方向，旋转角度为正；倾斜校正以图像中心为旋转中心，将图像中心平移至原点，绕坐标轴旋转后再平移回原位置；假设ct图像中心位置坐标为：p
c
＝[x
c y
c z
c
]，某一体素位置坐标为：p＝[x y z]，经过旋转变换后其位置坐标为：
[0024]
p
′
＝(p
‑
p
c
)
·
r(α，β，γ) p
c
[0025]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2.2)
[0026]
校正前后ct图像特征发生变化：绕z轴旋转可校正左右倾斜，校正后颅脑在水平面或矢状面投影最长，面部皮肤边缘最前点y轴坐标最小；绕x轴旋转校正上下倾斜，校正后最前端皮肤边缘由下巴变为鼻尖或上唇处，面部皮肤边缘最前点的水平高度明显变大；绕y轴旋转校正两侧歪斜，校正后颅脑左右对称，骨质边缘最左侧点和最右侧点的水平高度近似相等；据此确定α、β、γ的值：遍历γ可能值
‑
90
°
～90
°
进行旋转变换并检测面部皮肤边缘最
前点，当其y轴坐标最小时对应角度即为γ取值；遍历α可能值
‑
90
°
～0
°
进行旋转变换并检测面部皮肤边缘最前点，当其z轴坐标显著增大时对应角度即为α取值；遍历β可能值
‑
90
°
～90
°
进行旋转变换并检测骨质边缘最左侧点和最右侧点，当其z轴坐标近似相等时对应角度即为β取值；然后按公式(2.1)、(2.2)更新ct图像所有体素点坐标，实现自动倾斜校正。
[0027]
优选地，所述步骤(1)中采用中值滤波法对ct图像进行去噪预处理。
[0028]
优选地，所述步骤(2)中，对于大小为m
×
n
×
h的三维ct图像f(x，y，z)，使用轴向区域生长法进行分割的步骤如下：
[0029]
(2.1)自动定位或者交互选择分割目标所在的感兴趣区域；
[0030]
(2.2)选择轴向边缘点对距离较大的坐标轴a1，沿另一轴a2将三维图像切分为一系列二维切片；
[0031]
(2.3)扫描二维切片每一行或列，检测和求解沿a1方向距离最大边缘点对，若某行边缘点对不存在则距离为0，若存在多对则取其中距离最大一对；
[0032]
(2.4)所有切片沿a1方向边缘点对最大距离构成一维序列，求解其大于阈值的最长连续子序列，子序列两端点为沿a2方向距离最大边缘点对，实现a2轴上的分割；
[0033]
(2.5)对子序列中切片，以距离最大边缘点对连线上点为种子点，分别沿a3正负向生长至目标边缘，实现a3轴上的分割，若未检测到边缘则与相邻种子点的边缘平齐；
[0034]
(2.6)求解沿a3轴方向距离最大边缘点对，以其连线上点为种子点，分别沿a1正负向生长至目标边缘，实现a1轴上的分割，若未检测到边缘点则与相邻种子点的边缘平齐；
[0035]
(2.7)分别沿坐标轴a1、a3切分图像，对每个切片上分割结果进行连通区域分析，去除离群点，完成分割。
[0036]
优选地，所述步骤(3)中，对于重要组织表面轮廓的体素，参照公式(2.5)计算其邻域内所有像素灰度值的中值v
med
，作为所需等值面值
[0037]
g(i，j，k)＝median({f(i u，j v，k w):
‑
n≤u，v，w≤n})
ꢀꢀ
(2.5)
[0038]
其中f(i u，j v，k w)为原像素邻域内像素值；
[0039]
对体素邻域内所有像素进行多项式曲线拟合以代替线性插值，多项式函数为：
[0040][0041]
其中w＝[w0，w1，
…
，w
m
]为多项式系数集合，x为体素x轴坐标，m为多项式阶数，为减小计算量取为3；拟合误差函数为：
[0042][0043]
其中v
n
为体素邻域内像素灰度值，对于三维ct图像可与其相邻的所有像素：n＝27，λ为正则化惩罚项的调节因子，取0.01。
[0044]
优选地，所述步骤(3)中，f(x，w)为w的线性函数，e(w)为凸函数，其最小值在一阶导数为0处取得，最小化误差函数e(w)求解出多项式系数w；令f(x，w)＝v
med
，求解该方程位于体素邻域内的解x
med
，为重建表面点的x轴坐标，y轴和z轴坐标的求解与此相同；最后连接所有相邻表面点得到三角面片，拼接所有面片获得重建表面。
[0045]
本领域普通技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，包括上述实施例方法的各步骤，而所述的存储介质可以是：rom/ram、磁碟、光盘、存储卡等。因此，与本发明的方法相对应的，本发明还同时包括一种路径规划鼻部环境建模装置。如图1所示，该装置包括：
[0046]
预处理模块，其配置来对ct图像进行滤波去噪预处理；
[0047]
轴向区域生长模块，其配置来利用轴向区域生长法，结合边缘检测提取分割目标沿坐标轴方向的边缘点对，计算边缘点对之间距离，以距离最大边缘点对之间的像素点作为种子点，沿轴向生长至目标边缘，实现重要组织器官的快速分割；
[0048]
重建模块，其配置来提取重要组织器官的表面轮廓作为障碍物和地图边界，利用改进的marching cubes算法进行三维表面重建。
[0049]
以下更详细地说明本发明。
[0050]
1、图像预处理
[0051]
1.1倾斜校正
[0052]
人体或动物体不同组织对x线的吸收系数不同，进行电子计算机断层扫描后得到的三维空间的灰度分布，即为ct图像。在ct图像空间中，通常x轴为左右方向，与x轴垂直、将图像切分为左右两部分的切面为矢状面；y轴为前后方向，与y轴垂直、将图像切分为前后两部分的切面为冠状面；z轴为上下方向，将图像切分为上下两部分、与z轴垂直的切面为水平面(横断面)。在进行ct扫描时，若患者体位不正会导致ct图像倾斜、左右不对称，不仅不利于医生观察，还会影响重要组织和病灶目标的定位和分割。在动物实验中需要先将动物麻醉，再进行ct扫描，更容易发生倾斜现象。旋转变换能够将倾斜的ct图像调正。
[0053]
以原点为旋转中心，ct图像中任一体素绕x轴、y轴、z轴分别旋转α、β、γ角度，其位置坐标的旋转变换矩阵为：
[0054][0055]
按照右手系定则，大拇指指向坐标轴正向时四指方向为旋转的正方向，旋转角度为正。倾斜校正需要以图像中心为旋转中心，因此将图像中心平移至原点，绕坐标轴旋转后再平移回原位置。假设ct图像中心位置坐标为：p
c
＝[x
c y
c z
c
]，某一体素位置坐标为：p＝[x y z]，经过旋转变换后其位置坐标为：
[0056]
p
′
＝(p
‑
p
c
)
·
r(α，β，γ) p
c
ꢀꢀꢀꢀꢀ
(2.2)
[0057]
结合先验知识可知校正前后ct图像特征发生变化：绕z轴旋转可校正左右倾斜，校正后颅脑在水平面或矢状面投影最长，面部皮肤边缘最前点y轴坐标最小；绕x轴旋转可校正上下倾斜(通常为向上倾斜)，校正后最前端皮肤边缘由下巴(或下唇)变为鼻尖(人体)或上唇(猴子)处，面部皮肤边缘最前点的水平高度(z轴坐标)明显变大；绕y轴旋转可校正两
侧歪斜，校正后颅脑左右对称，骨质边缘最左侧点和最右侧点的水平高度近似相等。据此确定α、β、γ的值：遍历γ可能值
‑
90
°
～90
°
进行旋转变换并检测面部皮肤边缘最前点，当其y轴坐标最小时对应角度即为γ取值；遍历α可能值
‑
90
°
～0
°
进行旋转变换并检测面部皮肤边缘最前点，当其z轴坐标显著增大时对应角度即为α取值；遍历β可能值
‑
90
°
～90
°
进行旋转变换并检测骨质边缘最左侧点和最右侧点，当其z轴坐标近似相等时对应角度即为β取值。然后按式2.1、2.2更新ct图像所有体素点坐标，实现自动倾斜校正。实际应用中也可交互调整旋转角度，实时观察校正后ct图像，直至左右对称、能够清晰观察到视神经管等重要组织和病灶区域。
[0058]
1.2滤波去噪
[0059]
受外界环境、成像设备质量等因素的影响，计算机断层扫描成像过程中会产生噪声，使得ct图像的信噪比下降，甚至会导致图像某些细节特征被噪声所掩盖，给后续的特征提取、边缘检测和器官分割带来困难，因此需要对ct图像进行滤波去噪。ct图像中噪声来源很多，有随机噪声、电子噪声和光学噪声等，主要表现为高斯噪声(白噪声)和脉冲噪声(椒盐噪声)
[75]
。常用的ct图像去噪方法如下：
[0060]
(1)邻域平均法
[0061]
使用邻域内所有像素的平均值替换图像中原像素的值。对于ct图像f，可以设定一个空间滤波器作为三维卷积核，用该卷积核在图像上滑动计算，即可更新所有像素的灰度值：
[0062][0063]
式中f(i，j，k)为原像素值，2n 1可取值为3、5等奇数，卷积核大小为(2n 1)3。
[0064]
(2)高斯滤波法
[0065]
对图像f进行加权平均，使用邻域内所有像素灰度值的加权平均替换原像素值，像素灰度值更新公式为：
[0066][0067]
其中k(u，v，w)为权重矩阵每个元素的值，该权重矩阵为高斯卷积核，可以由高斯函数求出。例如，二维高斯函数为权重矩阵中下标为(u，v)的元素值为：令σ＝1.5，n＝1，则可求得3
×
3的二维权重矩阵为：
[0068][0069]
(3)中值滤波法
[0070]
使用邻域内所有像素灰度值的中值替换原像素的值：
[0071]
g(i，j，k)＝median({f(i u，j v，k w):
‑
n≤u，v，w≤n})
ꢀꢀꢀ
(2.5)
[0072]
式中f(i u，j v，k w)为原像素邻域内像素值。
[0073]
(4)自适应维纳滤波法
[0074]
假设图像和噪声都是随机变量，找到未被噪声污染图像f的一个估计使得它们之间的均方误差最小，该误差度量为：
[0075][0076]
该误差函数的最小值在频率域中表示为：
[0077][0078]
其中，h(u，v)为退化函数，s
η
(u，v)为噪声功率谱，s
f
(u，v)为未退化图像的功率谱，对于白噪声，其功率谱密度在整个频域内是常数，上式可近似为：
[0079][0080]
式中，k是加到|h(u，v)|2所有项上的特定常数。
[0081]
(5)小波去噪法
[0082]
对大小为m
×
n的图像f(x，y)进行离散小波变换，得到低频和高频两部分：
[0083][0084][0085]
其中是二维尺度函数，为i＝{h，v，d}时对应的二维小波，h、v、d分别度量列方向、行方向、对角线方向的变化，j0为任意开始尺度，通常取0。噪声能量一般分布在所有小波系数上，图像信号能量主要集中在少数小波系数上，对小波系数进行非线性阈值处理，如软阈值处理，以去除噪声分量：
[0086][0087]
最后再利用处理后的小波系数进行离散小波反变换：
[0088][0089]
完成重构得到去噪后的图像。
[0090]
分别利用上述方法对ct图像进行去噪处理。不同方法各有优劣，结合算法原理和实验结果分析，邻域平均法由于平均效果会导致图像细节模糊；在抑制服从正态分布的噪声时，高斯滤波比较有效；中值滤波保持图像细节的同时，也能在一定程度上消除噪声；自适应维纳滤波法能够根据图像的局部方差自适应调整输出，局部方差越大滤波器的平滑作用就越强，不适合处理非平稳信号；小波变换去噪虽然能有效抑制图像噪声，但会产生振铃效应，不能很好地处理边缘轮廓信息，且算法实现比较复杂。综上所述，本文采用中值滤波法对ct图像进行去噪预处理。
[0091]
2、结合边缘检测的轴向区域生长
[0092]
提出一种结合边缘检测的轴向区域生长策略，提取图像中与坐标轴平行的边缘点对的最大距离作为分割目标的主要特征，以距离最大的轴向边缘点对连线线段上像素点作为种子点，沿该轴向的法线方向生长至分割目标边缘。
[0093]
定义分割目标的轴向边缘点对：若两个边缘点连线与坐标轴平行，且连线线段上像素点均位于分割目标内，则这两个点是分割目标的轴向边缘点对。在ct图像中，不同组织的ct值差异较大，如表2.1所示。
[0094]
表2.1
[0095][0096]
据此可构造与sobel算子等梯度算子类似的边缘检测子，用于检测目标边缘点对。例如，当软组织被骨质包围时，其边缘点对为骨质——软组织边缘和软组织——骨质边缘。令j(i)＝j
‑2·
j
‑1·
j0·
j1·
j2，当且仅当f(i)位于软组织和骨质边缘时，j(i)响应为1，否则为0。特别地，当超出图像边界范围时，令f(i)＝300。当两个边缘点i1、i2之间的体素为软组织时：
‑
300<f(i)<300，i1<i<i2，二者为软组织的边缘点对。
[0097]
对于大小为m
×
n
×
h的三维ct图像f(x，y，z)，使用轴向区域生长法进行分割的步骤如下：
[0098]
(1)自动定位或者交互选择分割目标所在的感兴趣区域；
[0099]
(2)选择轴向边缘点对距离较大的坐标轴a1，沿另一轴a2将三维图像切分为一系列二维切片；
[0100]
(3)扫描二维切片每一行(列)，检测和求解沿a1方向距离最大边缘点对，若某行边
缘点对不存在则距离为0，若存在多对则取其中距离最大一对；
[0101]
(4)所有切片沿a1方向边缘点对最大距离构成一维序列，求解其大于阈值的最长连续子序列(子序列两端点为沿a2方向距离最大边缘点对)实现a2轴上的分割；
[0102]
(5)对子序列中切片，以距离最大边缘点对连线上点为种子点，分别沿a3正负向生长至目标边缘，实现a3轴上的分割，若未检测到边缘则与相邻种子点的边缘平齐；
[0103]
(7)求解沿a3轴方向距离最大边缘点对，以其连线上点为种子点，分别沿a1正负向生长至目标边缘，实现a1轴上的分割，若未检测到边缘点则与相邻种子点的边缘平齐；
[0104]
(8)分别沿坐标轴a1、a3切分图像，对每个切片上分割结果进行连通区域分析，去除离群点，完成分割。
[0105]
该方法是一种实现简单、快速有效的图像分割方法。对于分割目标轮廓边缘的未封闭区域，以相邻轮廓边缘作为边缘，可能会带来分割误差，因此分割目标轮廓边缘封闭性越好则分割效果越好。使用该方法能够实现脑组织、眼眶、下鼻甲、鼻中隔等重要组织的快速分割，进而提取其表面轮廓用于经鼻视神经管减压术路径规划的鼻部环境建模。
[0106]
3、基于表面重建的鼻部环境建模
[0107]
用于路径规划环境地图建模的方法有很多，主要包括栅格地图法、可视图法、voronoi图法等。由于鼻部环境复杂多变，重要组织等障碍物形状不规则，无明确顶点，不适合采用可视图法等图方法。三维ct图像和重要组织分割结果在空间和亮度上均已离散栅格化，可采用栅格地图法构建鼻部环境模型。
[0108]
栅格地图法将环境空间标记为不可通行的障碍物栅格和可通行栅格，简单有效，栅格越多地图精度越高越容易找到最优路径，但大量栅格需要较大的存储空间，搜索效率可能会降低。ct图像数据量庞大，例如，尺寸为512
×
512
×
100的ct图像像素点数目超过两千万。因此，提取脑组织等重要组织的表面轮廓，代替其作为障碍物和地图边界，以减少地图栅格数量，提高搜索效率。为进一步提高地图精度，可基于表面重建技术构建鼻部环境模型。
[0109]
marching cubes算法是最常用的三维表面重建方法。首先对ct图像所有体素顶点进行判断和标记：获取体素顶点的灰度值，将大于预设阈值的顶点标记为1，否则标记为0。每个顶点有两种状态，8个顶点共有2^8＝256种情况。根据立方体的旋转对称性和互补对称性(等值面拓扑结构不变)，可将等值面分布的256种情况合并为15种。
[0110]
然后通过线性插值计算等值面与体素交点位置：
[0111][0112]
式中v
th
为等值面值，p1、p2为等值面与体素交线的顶点，f(p1)、f(p2)为其灰度值，p为所求交点。
[0113]
marching cubes算法线性插值时未考虑邻域内体素影响，且其等值面值通常需要提前设定。针对这些问题进行改进，自适应重建表面。对于重要组织表面轮廓的体素，参照式(2.5)计算其邻域内所有像素灰度值的中值v
med
，作为所需等值面值。对体素邻域内所有像素进行多项式曲线拟合以代替线性插值，多项式函数为：
[0114][0115]
式中w＝[w0，w1，
…
，w
m
]为多项式系数集合，x为体素x轴坐标，m为多项式阶数，为减小计算量可取为3。拟合误差函数为：
[0116][0117]
式中v
n
为体素邻域内像素灰度值，对于三维ct图像可取与其相邻的所有像素(包含自身)：n＝27，λ为正则化惩罚项的调节因子，可取接近0的微小值如0.01，防止多项式系数过大。
[0118]
易知f(x，w)为w的线性函数，e(w)为凸函数，其最小值在一阶导数为0处取得，最小化误差函数e(w)即可求解出多项式系数w。令f(x，w)＝v
med
，求解该方程位于体素邻域内的解x
med
，即为重建表面点的x轴坐标，y轴和z轴坐标的求解与此类似。最后连接所有相邻表面点得到三角面片，拼接所有面片获得重建表面。已知鼻部周围障碍物表面信息即可完成路径规划，无需再显示求解可通行区域，因此采用该方法对所有障碍物进行三维表面重建，得到所需鼻部环境模型。
[0119]
以上所述，仅是本发明的较佳实施例，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属本发明技术方案的保护范围。