CN104574413A - 一种肺部ct图像的血管分叉点提取方法及*** - Google Patents

Abstract

本发明属于图像处理领域，提供了一种肺部CT图像的血管分叉点提取方法及***。该方法及***充分利用了血管分叉点在几何结构上的球型特征、以及张量投票的累加性，通过对增强后的肺CT图像上属于血管的像素点进行球张量投票，求解球张量显著性系数的局部最大值以获得备选血管分叉点，之后利用主成分分析算法去除血管的末端点，得到准确的肺部血管分叉点。该方法及***无需图像前期操作，对噪声具有较强的鲁棒性，为后期的图像配准和肺部运动估计提供重要依据，有利于后期的肺部CT图像分析，辅助临床诊断与治疗。

Description

一种肺部CT图像的血管分叉点提取方法及***

技术领域

本发明属于图像处理领域，尤其涉及一种肺部CT图像的血管分叉点提取方法及***。

背景技术

肺部CT图像是利用人体组织对X射线的吸收与透过率的不同扫描获得的重建图像。提取肺部CT图像的特征点具有重要应用意义，是视觉处理和图像处理的前提和基础。肺部CT图像的特征主要有角点，边缘点，拐点，边界和纹理等，在肺部CT图像中，血管分叉点因其特殊的几何结构，是肺运动估计、图像配准等计算机视觉领域中任务的基础，对图像分析中的自动标定和配准的精度产生直接影响。

现有技术提出的肺部CT图像的血管分叉点提取方法可以大致分为三类：第一类是基于血管中心线的分叉点提取方法，这类方法的一般流程为，先用血管分割算法提取出肺部图像的血管结构，再用血管细化算法进行骨架提取以得到单像素血管树形结构，而从血管树中提取到的节点作为血管分叉点。这种算法依赖于前期血管分割的精度以及细化的准确性，容易产生许多误检分叉点，或是对弱血管存在漏检情况；第二类是基于Hessian矩阵的分叉点提取方法，这类方法首先要计算图像二阶导数矩阵也即Hessian矩阵，继而求解矩阵特征值，用特征组合构建一个滤波器，以在分叉点处得到最大响应。这类方法对噪声敏感，同时难以区分血管分叉点与血管端点；第三类方法是基于机器学习算法的分叉点提取方法，这类方法主要是利用机器学习算法构建一个分类器，对从图像中提取到的特征(如特征向量直方图、灰度值局部均值等)进行训练以从候选点中划分出分叉点。这类方法由于使用了机器学习算法，在训练阶段会占用比较多时间，提取效率较低。

发明内容

本发明实施例的目的在于提供一种肺部CT图像的血管分叉点提取方法，旨在解决现有的血管交叉点提取方法对细化结果要求高以及对噪声敏感的问题。

本发明实施例是这样实现的，一种肺部CT图像的血管分叉点提取方法，所述方法包括以下步骤：

对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点；

每一个属于血管的像素点分别向投票窗口内的其它属于血管的像素点进行球张量投票；

根据球张量投票结果，对每一点累积的球张量票数进行分解，得到备选血管分叉点的集合；

利用主成分分析算法对每一备选血管分叉点进行分析，以去除血管的末端点，得到准确的血管分叉点。

本发明实施例的另一目的在于提供一种肺部CT图像的血管分叉点提取***，所述***包括：

血管像素点确定模块，用于对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点；

球张量投票模块，用于控制每一个属于血管的像素点分别向投票窗口内的其它属于血管的像素点进行球张量投票；

备选血管分叉点集合获取模块，用于根据球张量投票结果，对每一点累积的球张量票数进行分解，得到备选血管分叉点的集合；

血管分叉点提取模块，用于利用主成分分析算法对每一备选血管分叉点进行分析，以去除血管的末端点，得到准确的血管分叉点。

本发明提出的肺部CT图像的血管分叉点提取方法及***充分利用了血管分叉点在几何结构上的球型特征、以及张量投票的累加性，通过对增强后的肺CT图像上属于血管的像素点进行球张量投票，求解球张量显著性系数的局部最大值以获得备选血管分叉点，之后利用主成分分析算法去除血管的末端点，得到准确的肺部血管分叉点。该方法及***执行效率较高，无需图像前期处理操作，对噪声具有较强的鲁棒性，可以为后期的图像配准和肺部运动估计提供重要依据，有利于后期的肺部CT图像分析，辅助临床诊断与治疗。

附图说明

图1是本发明实施例提供的肺部CT图像的血管分叉点提取方法的流程图；

图2是本发明实施例中，确定属于血管的像素点的第一种详细流程图；

图3是本发明实施例中，确定属于血管的像素点的第二种详细流程图；

图4是本发明实施例中，对特征值和特征向量进行重构的详细流程图；

图5是本发明实施例中，确定准确的血管分叉点的详细流程图；

图6是本发明实施例提供的肺部CT图像的血管分叉点提取***的结构图；

图7是图6中血管像素点确定模块的第一种结构图；

图8是图6中血管像素点确定模块的第二种结构图；

图9是图8中重构子模块的结构图；

图10是图6中血管分叉点提取模块的结构图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提出的肺部CT图像的血管分叉点提取方法充分利用了血管分叉点在几何结构上的球型特征、以及张量投票的累加性，通过对增强后的肺CT图像上属于血管的像素点进行球张量投票，求解球张量显著性系数的局部最大值以获得备选血管分叉点，之后利用主成分分析算法去除血管的末端点，得到准确的肺部血管分叉点。

图1示出了本发明实施例提供的肺部CT图像的血管分叉点提取方法的流程，包括以下步骤：

S1：对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点。

作为第一种优选实现方式，如图2所示，步骤S1进一步包括以下步骤：

S11：利用多尺度高斯函数对肺CT图像进行平滑。

假设G(x,y,z；σ)为尺度为σ的三维高斯函数，则肺CT图像I(x,y,z)在尺度σ下的平滑结果I_σ(x,y,z)表示为：其中， $G(x,y,z;\mathrm{\σ})=\frac{1}{{\left(\sqrt{2\mathrm{\π}}{\mathrm{\σ}}^2\right)}^3}{e}^{-\frac{x^2+y^2+z^2}{2{\mathrm{\σ}}^2}}.$

S12：在每一尺度下，根据平滑结果计算肺CT图像中每一点的Hessian矩阵。

假设尺度σ下，肺CT图像中点(x,y,z)处的Hessian矩阵为H_σ(x,y,z)，则其表示为：

$H_{\mathrm{\σ}}(x,y,z)=\left[\begin{array}{ccc}\frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{{\∂x}^2}& \frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{\∂x\∂y}& \frac{\∂I_{\mathrm{\σ}}^2(x,y,z)}{\∂x\∂z}\\ \frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{\∂x\∂y}& \frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{{\∂y}^2}& \frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{\∂y\∂z}\\ \frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{\∂x\∂z}& \frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{\∂y\∂z}& \frac{{\∂I}_{\mathrm{\σ}}^2(x,y,z)}{{\∂z}^2}\end{array}\right]$

S13：对每一点的Hessian矩阵进行特征值分解，得到三个特征值以及与三个特征值分别一一对应的特征向量。

本发明实施例中，对于Hessian矩阵H_σ(x,y,z)进行分解后得到的三个特征值记为λ₁、λ₂、λ₃，且满足|λ₁|≤|λ₂|≤|λ₃|；与三个特征值分别一一对应的三个特征向量记为

Hessian矩阵的特征值和特征向量可以描述管状结构的几何特征。具体来说，对于属于管状结构上的点，其沿着血管走向的特征向量对应的特征值为三个特征值中较小的一个；而沿与血管走向垂直的切面方向的其余两个特征向量张成一个平面，且其余两个特征向量对应的特征值大小相近，并为三个特征值中较大的两个，即满足|λ₃|≈|λ₂|》|λ₁|≈0。

S14：根据每一点的特征值和特征向量，估计对应点在每一尺度下属于管状结构的可能性。

假设尺度σ下，图像中点(x,y,z)对应的三个特征值满足|λ₁|≤|λ₂|≤|λ₃|，点(x,y,z)属于管状结构的可能性为V_s(σ)，则满足：

$V_s\left(\mathrm{\σ}\right)=\left\{\begin{array}{cc}0& {\mathrm{\λ}}_2>0\mathrm{or}{\mathrm{\λ}}_3>0\\ (1-e^{-\frac{R_A^2}{2{\mathrm{\α}}^2}})\·e^{-\frac{R_B^2}{2{\mathrm{\β}}^2}}\·(1-e^{-\frac{S^2}{2{\mathrm{\γ}}^2}})\·e^{-\frac{2{\mathrm{Coeff}}^2}{\left|{\mathrm{\λ}}_2\right|{\mathrm{\λ}}_3^2}}& \mathrm{otherwise}\end{array}\right.$

其中， $R_A=\frac{\left|{\mathrm{\λ}}_2\right|}{\left|{\mathrm{\λ}}_3\right|},$ $R_B=\frac{\left|{\mathrm{\λ}}_1\right|}{\sqrt{\left|{\mathrm{\λ}}_2{\mathrm{\λ}}_3\right|}},$ $S=\sqrt{{\mathrm{\λ}}_1^2+{\mathrm{\λ}}_2^2+{\mathrm{\λ}}_3^2},$ coeff是一个常数，α是常数，一般可取0.5，β是常数，一般可取0.5，γ是设定的常数。

S15：取每一点在不同尺度下的可能性的最大值作为对应点属于管状结构的可能性的最终值。

假设该可能性的最终值记为V，则有：其中，σ_min,σ_max分别是最小尺度和最大尺度。

S16：利用弥散函数对肺CT图像中每一可能性大于0的点的强度进行更新，直到更新次数达到最大迭代次数为止，得到增强后的血管强度。

本发明实施例中，利用VED算法对图中管状结构的可能性进行弥散，弥散函数表示为：其中，V_t是弥散后的管状结构强度，t是弥散时间，▽.是散度算子，D是弥散张量，且满足：

${\mathrm{\λ}}_1^{\′}=-+(\mathrm{\ω}-1)V^{\frac{1}{L}}$

${\mathrm{\λ}}_2^{\′}={\mathrm{\λ}}_3^{\′}=1+(\mathrm{\ϵ}-1)V^{\frac{1}{L}}$

其中，为Hessian矩阵H_σ(x,y,z)进行分解后得到的三个特征值，ω为一参数，用以表明各向异性弥散的强度，可取ω＝5，ε为一参数，用以保证弥散张量D是一个正定矩阵，可取ε＝0.01，L为一参数，用以控制弥散函数对血管影响的敏感性，可取L＝2。

S17：判定增强后的血管强度达到设定阈值的点为属于血管的像素点。

假设设定阈值为Th，若增强后的结果中，点(x,y,z)的增强后的血管强度V(x,y,z)满足V(x,y,z)≥Th，则判定点(x,y,z)为属于血管的像素点。其中，设定阈值Th可取增强后的肺CT图像中各点增强后的血管强度的平均值。

作为第二种优选实现方式，如图3所示，与图2所示不同，此时，在通过步骤S15得到可能性的最终值之后，在步骤S16之前，还包括以下步骤：

S18：以可能性大于0的点的特征值最小方向为法线方向，对邻域内其它可能性大于0的点进行棒张量投票，并根据投票结果对每一可能性大于0的点的特征值和特征向量进行重构，以确定每一可能性大于0的点的管状结构的走向方向。

进一步地，如图4所示，步骤S18又可包括以下步骤：

S181：以每一可能性大于0的点为投票点、以对应投票点的特征值最小的方向为法线方向，对邻域内其它可能性大于0的点为票数接收点进行棒张量投票。

假设可能性大于0的点(x,y,z)对应的三个特征值满足|λ₁|≤|λ₂|≤|λ₃|，特征向量为棒张量为S，板张量为P，球张量为B，则点(x,y,z)在对应尺度下的Hessian矩阵H可分解为棒张量、板张量和球张量之和，即有：H＝(λ₃-λ₂)S+(λ₂-λ₁)P+λ₁B，其中， (λ₃-λ₂)表示曲面显示性。

在本发明实施例中，假设可能性大于0的点(x,y,z)为棒张量投票点，以点(x,y,z)的特征值最小的方向为法线方向，对邻域内其他可能性大于0的点R进行投票，R为票数接收点，则点(x,y,z)向点R投出的票数为包含方向和强度的棒张量Stick(l,θ,π)，且满足：

$\mathrm{Stick}(l,\mathrm{\θ},\mathrm{\π})=({\mathrm{\λ}}_3-{\mathrm{\λ}}_2)\mathrm{DF}(s,k,\mathrm{\σ})\left[\begin{array}{c}-\sin \left(2\mathrm{\θ}\right)\\ \cos \left(2\mathrm{\θ}\right)\\ 0\end{array}\right]\left[\begin{array}{ccc}-\sin \left(2\mathrm{\θ}\right)& \cos \left(2\mathrm{\θ}\right)& 0\end{array}\right]$

其中， $\mathrm{DF}(s,k,\mathrm{\σ})=e^{-\frac{s^2+c{k}^2}{{\mathrm{\σ}}^2}}$ 为显著性衰减函数，且 $s=\frac{\mathrm{\θl}}{\sin \mathrm{\θ}},$ $k=\frac{2\sin \mathrm{\θ}}{l},$ θ为点(x,y,z)与点R的连线l与张成的平面夹角，张成的平面的法线方向为s为连线l的弧长，σ指定了投票的尺度范围，决定了投票窗口的大小，c为尺度范围σ的函数，用于制约曲率的退化程度，且满足：

本发明实施例中，通过以最小特征值对应特征向量为法线方向，又加上一个棒张量显著性(λ₃-λ₂)作为权重，进行棒张量投票。投票结束后，图像中的每一可能性大于0的都会获得周围邻域内其他点的投票累加。

S182：根据投票结果对每一可能性大于0的点的特征值和特征向量进行重构，以确定每一可能性大于0的点的管状结构的走向方向。

本发明实施例中，对票数接收点R处接收到的票数Stick(l,θ,π)进行累加，累加的过程包括张量大小和方向的累加，记T_R′(x,y,z)为接收点接收到的累加张量，对其进行特征分解：

其中|λ'₃|≤|λ'₂|≤|λ'₁|为T_R′(x,y,z)的特征值，为投票结束后累加张量T_R′(x,y,z)的特征向量，该三个新特征向量分别对应特征值最小、次小、最大的特征值，此时所得的特征向量的方向即是对原图扩散方向的纠正方向。此时，步骤S16是根据重构的特征向量代替原来的特征向量利用弥散函数对血管结构的血管强度进行弥散的。

相对于第一种实现方式，第二种实现方式利用重构的张量方向可以较好地减少血管强度沿血管切面的弥散，同时增强沿血管方向的弥散，达到抑制噪声，增强血管特征的效果。

S2：每一个属于血管的像素点分别向投票窗口内的其它属于血管的像素点进行球张量投票。

在本发明实施例中，假设属于血管的像素点(x,y,z)为球张量投票点，对投票窗口内的点Q进行投票，Q为票数接收点。选定连接点(x,y,z)和点Q的光滑球面，记球面弧长为s，则点(x,y,z)向点Q投出的票数为棒张量Stick(l,θ,π)在两个方向旋转后得到的球张量B，且满足：

$B(x,y,z)={\mathrm{\λ}}_1{\mathrm{\λ}}_2{\mathrm{\λ}}_3{I}_{\mathrm{\σ}}(x,y,z){\∫}_0^{2\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}{S}_{\mathrm{\θ\φ\ψ}}\mathrm{d\φd\ψ}{|}_{\mathrm{\θ}=0}$

其中，I_σ(x,y,z)表示在尺度σ下所述肺CT图像的平滑结果，S_θφψ表示棒张量Stick(l,θ,π)旋转角度φ和角度ψ后的结果，角度φ是棒张量Stick(l,θ,π)在xy平面内、以z轴为旋转中心的旋转角度，角度ψ是棒张量Stick(l,θ,π)在xz平面内、以y轴为旋转中心的旋转角度，棒张量Stick(l,θ,π)满足：

$\mathrm{Stick}(l,\mathrm{\θ},\mathrm{\π})=\mathrm{DF}(s,k,\mathrm{\σ})\left[\begin{array}{c}-\sin \left(2\mathrm{\θ}\right)\\ \cos \left(2\mathrm{\θ}\right)\\ 0\end{array}\right]\left[\begin{array}{ccc}-\sin \left(2\mathrm{\θ}\right)& \cos \left(2\mathrm{\θ}\right)& 0\end{array}\right]$

其中， $\mathrm{DF}(s,k,\mathrm{\σ})=e^{-\frac{s^2+c{k}^2}{{\mathrm{\σ}}^2}}$ 为显著性衰减函数，且 $s=\frac{\mathrm{\θl}}{\sin \mathrm{\θ}},$ $k=\frac{2\sin \mathrm{\θ}}{l},$ θ为点(x,y,z)与点Q的连线l与或张成的平面夹角，l是点(x,y,z)与点Q的连线，s为l的弧长，σ即为前述的尺度，其指定了投票的尺度范围，决定了投票窗口的大小，c为尺度范围σ的函数，用于制约曲率的退化程度，且满足： $c=\frac{-16\log \left(0.1\right)\×(\mathrm{\σ}-1)}{{\mathrm{\π}}^2}.$

本发明实施例中，投票窗口可以是一边长为ws的正方体，且边长ws满足：λ₁λ₂λ₃作为票数的权重，是因为球状几何结果的三个特征值都较大，三个特征值的乘积也较大。

像素点Q得到的球张量票数B_Q是所述像素点Q所在的投票窗口内所有属于血管的像素点向像素点Q投出的球张量之和，即：

S3：根据球张量投票结果，对每一点累积的球张量票数进行分解，得到备选血管分叉点的集合。

本发明实施例中，步骤S3又可以包括以下步骤：在球张量投票结束后，对每一属于血管的像素点(x,y,z)得到的球张量投票结果B进行特征值分解：

其中，λ₁″、λ₂″、λ₃″分别为特征值分解后的三个特征值，且有|λ₁″≤|λ₂″|≤|λ₃″|，分别为与所述三个特征值分别一一对应的三个特征向量；

若点(x,y,z)的最小特征值|λ₁″|为投票窗口内的最大值，则确定点(x,y,z)为备选血管分叉点。

S4：利用主成分分析(Principle Component Analysis，PCA)算法对每一备选血管分叉点进行分析，以去除血管的末端点，得到准确的血管分叉点。

由于在利用步骤S3得到的备选血管分叉点中，存在血管的末端点，利用PCA算法对备选血管分叉点进行筛选，以获取更加准确的血管分叉点。此时，进一步地，如图5所示，步骤S4又可包括以下步骤：

S41：以备选血管分叉点为中心，取一投票窗口内所有属于血管的像素点的坐标，构成一集合。

S42：计算该集合内的坐标均值，以及协方差矩阵。

假设集合X为X＝{(x_i,y_i,z_i)|i＝1,...,N}，其中，N是投票窗口内属于血管的像素点的坐标个数，坐标均值为(v_x,v_y,v_z)，则有 若协方差矩阵为C，则有：

$C=\frac{1}{N}\left[\begin{array}{ccc}\underset{i}{\mathrm{\Σ}}(x_i-v_x)(x_i-v_x)& \underset{i}{\mathrm{\Σ}}(x_i-v_x)(y_i-v_y)& \underset{i}{\mathrm{\Σ}}(x_i-v_x)(z_i-v_z)\\ \underset{i}{\mathrm{\Σ}}(y_i-v_y)(x_i-v_x)& \underset{i}{\mathrm{\Σ}}(y_i-v_y)(y_i-v_y)& \underset{i}{\mathrm{\Σ}}(y_i-v_y)(z_i-v_z)\\ \underset{i}{\mathrm{\Σ}}(z_i-v_z)(x_i-v_x)& \underset{i}{\mathrm{\Σ}}(z_i-v_z)(y_i-v_y)& \underset{i}{\mathrm{\Σ}}(z_i-v_z)(z_i-v_z)\end{array}\right]$

S43：对协方差矩阵C进行特征值分解：

其中，λ_1,C、λ_2,C、λ_3,C分别为特征值分解后的三个特征值，且满足λ_1,C|≤|λ_2,C|≤|λ_3,C|，分别为与三个特征值分别一一对应的三个特征向量，并取最大特征值|λ_3,C|对应的特征向量作为主方向。

S44：以主方向为法线、过备选血管分叉点作一平面，统计处于该平面一侧的属于血管的像素点个数、以及处于该平面另一侧的属于血管的像素点个数，根据处于该平面两侧的属于血管的像素点个数的比值，区分血管的末端点和血管分叉点。

本发明实施例中，假设以主方向为法线、过备选血管分叉点(x,y,z)作一平面，处于该平面一侧的属于血管的像素点个数为N₁，处于该平面另一侧的属于血管的像素点个数为N₂，则当或时，判定血管分叉点(x,y,z)为血管的末端点而不是血管分叉点。

图6示出了本发明实施例提供的肺部CT图像的血管分叉点提取***的结构，为了便于说明，仅示出了与本发明实施例相关的部分，该***可以是内置于其它各类图像变换***中的硬件单元、软件单元或软硬件单元的结合。

具体地，本发明实施例提供的肺部CT图像的血管分叉点提取***包括：血管像素点确定模块1，用于对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点；球张量投票模块2，用于控制每一个属于血管的像素点分别向投票窗口内的其它属于血管的像素点进行球张量投票；备选血管分叉点集合获取模块3，用于根据球张量投票结果，对每一点累积的球张量票数进行分解，得到备选血管分叉点的集合；血管分叉点提取模块4，用于利用主成分分析算法对每一备选血管分叉点进行分析，以去除血管的末端点，得到准确的血管分叉点。其中，球张量投票模块2的详细执行流程对应步骤S2所述，备选血管分叉点集合获取模块3的详细执行流程对应步骤S3所述，不赘述。

进一步地，如图7所示，血管像素点确定模块1的第一种结构可包括：平滑子模块11，用于利用多尺度高斯函数对肺CT图像进行平滑；第一计算子模块12，用于在每一尺度下，根据平滑结果计算肺CT图像中每一点的Hessian矩阵；第二计算子模块13，用于对每一点的Hessian矩阵进行特征值分解，得到三个特征值以及与三个特征值分别一一对应的特征向量；估算子模块14，用于根据每一点的特征值和特征向量，估计对应点在每一尺度下属于管状结构的可能性；取值子模块15，用于取每一点在不同尺度下的可能性的最大值作为对应点属于管状结构的可能性的最终值；弥散子模块16，用于利用弥散函数对肺CT图像中每一可能性大于0的点的强度进行更新，直到更新次数达到最大迭代次数为止，得到增强后的血管强度；第一判定子模块17，用于判定增强后的血管强度达到设定阈值的点为属于血管的像素点。其中，平滑子模块11、第一计算子模块12、第二计算子模块13、估算子模块14、取值子模块15、弥散子模块16、第一判定子模块17分别的详细执行流程对应步骤S11至步骤S17所述，不赘述。

血管像素点确定模块1的第二种结构如图8所示。与图7所示同步，第二种结构还可包括：重构子模块18，用于以取值子模块15得到的可能性大于0的点的特征值最小方向为法线方向，对邻域内其它可能性大于0的点进行棒张量投票，并根据投票结果对每一可能性大于0的点的特征值和特征向量进行重构，以确定每一可能性大于0的点的管状结构的走向方向。此时，弥散子模块16是根据重构子模块18重构的特征向量，利用弥散函数对肺CT图像中每一可能性大于0的点的强度进行更新的。

更进一步地，如图9所示，重构子模块18可包括：棒张量投票子模块181，用于以每一可能性大于0的点为投票点、以对应投票点的特征值最小的方向为法线方向，对邻域内其它可能性大于0的点为票数接收点进行棒张量投票；特征值和特征向量重构子模块182，用于根据投票结果对每一可能性大于0的点的特征值和特征向量进行重构，以确定每一可能性大于0的点的管状结构的走向方向。其中，棒张量投票子模块181和特征值和特征向量重构子模块182的详细执行流程对应步骤S181和步骤S182所述，不赘述。

进一步地，如图10所示，血管分叉点提取模块4可包括：集合获取子模块41，用于以备选血管分叉点为中心，取一投票窗口内所有属于血管的像素点的坐标，构成一集合；第三计算子模块42，用于计算集合内的坐标均值(v_x,v_y,v_z)，以及协方差矩阵C；主方向获取子模块43，用于对协方差矩阵C进行特征值分解：其中，λ_1,C、λ_2,C、λ_3,C分别为特征值分解后的三个特征值，且满足|λ_1,C|≤|λ_2,C|≤|λ_3,C|，分别为与所述三个特征值分别一一对应的三个特征向量，并取最大特征值|λ_3,C|对应的特征向量作为主方向；第二判定子模块44，用于以主方向为法线、过备选血管分叉点作一平面，统计处于该平面一侧的属于血管的像素点个数、以及处于该平面另一侧的属于血管的像素点个数，根据处于该平面两侧的属于血管的像素点个数的比值，区分血管的末端点和血管分叉点。其中，第三计算子模块42的详细执行流程对应步骤S42所述，第二判定子模块44的详细执行流程对应步骤S44所述，不赘述。

本发明实施例提出的肺部CT图像的血管分叉点提取方法及***充分利用了血管分叉点在几何结构上的球型特征、以及张量投票的累加性，通过对增强后的肺CT图像上属于血管的像素点进行球张量投票，求解球张量显著性系数的局部最大值以获得备选血管分叉点，之后利用主成分分析算法去除血管的末端点，得到准确的肺部血管分叉点。该方法及***执行效率较高，无需图像前期处理操作，对噪声具有较强的鲁棒性，可以为后期的图像配准和肺部运动估计提供重要依据，有利于后期的肺部CT图像分析，辅助临床诊断与治疗。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤是可以通过程序来控制相关的硬件完成，所述的程序可以在存储于一计算机可读取存储介质中，所述的存储介质，如ROM/RAM、磁盘、光盘等。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种肺部CT图像的血管分叉点提取方法，其特征在于，所述方法包括以下步骤：

对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点；

每一个属于血管的像素点分别向投票窗口内的其它属于血管的像素点进行球张量投票；

根据球张量投票结果，对每一点累积的球张量票数进行分解，得到备选血管分叉点的集合；

利用主成分分析算法对每一备选血管分叉点进行分析，以去除血管的末端点，得到准确的血管分叉点。

2.如权利要求1所述的肺部CT图像的血管分叉点提取方法，其特征在于，所述对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点的步骤又包括以下步骤：

利用多尺度高斯函数对肺CT图像进行平滑，若尺度为σ的三维高斯函数为G(x,y,z；σ)，所述肺CT图像为I(x,y,z)，所述肺CT图像I(x,y,z)在所述尺度σ下的平滑结果为I_σ(x,y,z)，则所述利用多尺度高斯函数对肺CT图像进行平滑的表示为： $I_{\mathrm{\σ}}(x,y,z)=I(x,y,z)\⊗G(x,y,z;\mathrm{\σ}),$ 其中， $G(x,y,z;\mathrm{\σ})=\frac{1}{{\left(\sqrt{3\mathrm{\π}}\mathrm{\σ}\right)}^3}{e}^{-\frac{x^2+y^2+z^2}{2{\mathrm{\σ}}^2}};$

在每一尺度下，根据平滑结果计算所述肺CT图像中每一点的Hessian矩阵；

对每一点的所述Hessian矩阵进行特征值分解，得到三个特征值以及与三个特征值分别一一对应的特征向量；

根据每一点的所述特征值和特征向量，估计对应点在每一尺度下属于管状结构的可能性；

取每一点在不同尺度下的可能性的最大值作为对应点属于管状结构的可能性的最终值；

利用弥散函数对所述肺CT图像中每一可能性大于0的点的强度进行更新，直到更新次数达到最大迭代次数为止，得到增强后的血管强度；

判定所述增强后的血管强度达到设定阈值的点为所述属于血管的像素点。

3.如权利要求1所述的肺部CT图像的血管分叉点提取方法，其特征在于，所述对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点的步骤又包括以下步骤：

利用多尺度高斯函数对肺CT图像进行平滑，若尺度为σ的三维高斯函数为G(x,y,z；σ)，所述肺CT图像为I(x,y,z)，所述肺CT图像I(x,y,z)在所述尺度σ下的平滑结果为I_σ(x,y,z)，则所述利用多尺度高斯函数对肺CT图像进行平滑的表示为： $I_{\mathrm{\σ}}(x,y,z)=I(x,y,z)\⊗G(x,y,z;\mathrm{\σ}),$ 其中， $G(x,y,z;\mathrm{\σ})=\frac{1}{{\left(\sqrt{3\mathrm{\π}}\mathrm{\σ}\right)}^3}{e}^{-\frac{x^2+y^2+z^2}{2{\mathrm{\σ}}^2}};$

在每一尺度下，根据平滑结果计算所述肺CT图像中每一点的Hessian矩阵；

对每一点的所述Hessian矩阵进行特征值分解，得到三个特征值以及与三个特征值分别一一对应的特征向量；

根据每一点的所述特征值和特征向量，估计对应点在每一尺度下属于管状结构的可能性；

取每一点在不同尺度下的可能性的最大值作为对应点属于管状结构的可能性的最终值；

以所述可能性大于0的点的特征值最小方向为法线方向，对邻域内其它可能性大于0的点进行棒张量投票，并根据投票结果对每一可能性大于0的点的特征值和特征向量进行重构，以确定所述每一可能性大于0的点的管状结构的走向方向；

根据重构的所述特征向量，利用弥散函数对所述肺CT图像中每一可能性大于0的点的强度进行更新，直到更新次数达到最大迭代次数为止，得到增强后的血管强度；

判定所述增强后的血管强度达到设定阈值的点为所述属于血管的像素点。

4.如权利要求3所述的肺部CT图像的血管分叉点提取方法，其特征在于，所述以所述可能性大于0的点的特征值最小方向为法线方向，对邻域内其它可能性大于0的点进行棒张量投票，并根据投票结果对每一可能性大于0的点的特征值和特征向量进行重构，以确定所述每一可能性大于0的点的管状结构的走向方向的步骤又包括以下步骤：

以每一可能性大于0的点为投票点、以对应投票点的特征值最小的方向为法线方向，对邻域内其它可能性大于0的点为票数接收点进行棒张量投票；

根据投票结果对所述每一可能性大于0的点的特征值和特征向量进行重构，以确定所述每一可能性大于0的点的管状结构的走向方向。

5.如权利要求1所述的肺部CT图像的血管分叉点提取方法，其特征在于，若所述属于血管的像素点(x,y,z)对应的三个特征值为λ₁、λ₂、λ₃，且满足|λ₁|≤|λ₂|≤|λ₃|，对应的特征向量为点(x,y,z)向所述投票窗口内其它属于血管的像素点Q投出的票数为棒张量Stick(l,θ,π)在两个方向旋转后得到的球张量B(x,y,z)，且满足：

$B(x,y,z)={\mathrm{\λ}}_1{\mathrm{\λ}}_2{\mathrm{\λ}}_3{I}_{\mathrm{\σ}}(x,y,z){\∫}_0^{2\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}{S}_{\mathrm{\θ\φ\ψ}}\mathrm{d\φd\ψ}{|}_{\mathrm{\θ}=0}$

其中，I_σ(x,y,z)表示在尺度σ下所述肺CT图像的平滑结果，S_θφψ表示棒张量Stick(l,θ,π)旋转角度φ和角度ψ后的结果，角度φ是棒张量Stick(l,θ,π)在xy平面内、以z轴为旋转中心的旋转角度，角度ψ是棒张量Stick(l,θ,π)在xz平面内、以y轴为旋转中心的旋转角度，所述棒张量Stick(l,θ,π)满足：

$\mathrm{Stick}(l,\mathrm{\θ},\mathrm{\π})=\mathrm{DF}(s,k,\mathrm{\σ})\left[\begin{array}{c}-\sin \left(2\mathrm{\θ}\right)\\ \cos \left(2\mathrm{\θ}\right)\\ 0\end{array}\right]\left[\begin{array}{ccc}-\sin \left(2\mathrm{\θ}\right)& \cos \left(2\mathrm{\θ}\right)& 0\end{array}\right]$

其中， $\mathrm{DF}(s,k,\mathrm{\σ})=e^{-\frac{s^2+{\mathrm{ck}}^2}{{\mathrm{\σ}}^2}}$ 为显著性衰减函数，且 $s=\frac{\mathrm{\θl}}{\sin \mathrm{\θ}},k=\frac{2\sin \mathrm{\θ}}{l},$ θ为点(x,y,z)与点Q的连线l与张成的平面夹角，l是点(x,y,z)与点Q的连线，s为l的弧长，σ指定了投票的尺度范围，决定了投票窗口的大小，c为尺度范围σ的函数，用于制约曲率的退化程度，且满足：所述投票窗口是一边长为ws的正方体，且边长ws满足：

所述像素点Q得到的球张量票数B_Q是所述像素点Q所在的投票窗口内所有属于血管的像素点向所述像素点Q投出的球张量之和，即：

6.如权利要求1所述的肺部CT图像的血管分叉点提取方法，其特征在于，所述根据球张量投票结果，对每一点累积的球张量票数进行分解，得到备选血管分叉点的集合的步骤又包括以下步骤：

对所述每一属于血管的像素点(x,y,z)得到的球张量投票结果B进行特征值分解：

其中，λ₁″、λ₂″、λ₃″分别为特征值分解后的三个特征值，且有|λ₁″|≤|λ₂″|≤|λ₃″|，分别为与所述三个特征值分别一一对应的三个特征向量；

若所述点(x,y,z)的最小特征值|λ₁″|为投票窗口内的最大值，则确定所述点(x,y,z)为所述备选血管分叉点。

7.如权利要求1所述的肺部CT图像的血管分叉点提取方法，其特征在于，所述利用主成分分析算法对每一备选血管分叉点进行分析，以去除血管的末端点，得到准确的血管分叉点的步骤又包括以下步骤：

以所述备选血管分叉点为中心，取一投票窗口内所有属于血管的像素点的坐标，构成一集合；

计算所述集合内的坐标均值(v_x,v_y,v_z)，以及协方差矩阵C；

对所述协方差矩阵C进行特征值分解：

其中，λ_1,C、λ_2,C、λ_3,C分别为特征值分解后的三个特征值，且满足|λ_1,C|≤|λ_2,C|≤|λ_3,C|，分别为与所述三个特征值分别一一对应的三个特征向量，并取最大特征值|λ_3,C|对应的特征向量作为主方向；

以所述主方向为法线、过所述备选血管分叉点作一平面，统计处于所述平面一侧的属于血管的像素点个数、以及处于所述平面另一侧的属于血管的像素点个数，根据处于所述平面两侧的属于血管的像素点个数的比值，区分血管的末端点和血管分叉点。

8.一种肺部CT图像的血管分叉点提取***，其特征在于，所述***包括：

血管像素点确定模块，用于对肺CT图像中的支气管和血管进行增强，并确定增强后的像素点中、属于血管的像素点；

球张量投票模块，用于控制每一个属于血管的像素点分别向投票窗口内的其它属于血管的像素点进行球张量投票；

备选血管分叉点集合获取模块，用于根据球张量投票结果，对每一点累积的球张量票数进行分解，得到备选血管分叉点的集合；

血管分叉点提取模块，用于利用主成分分析算法对每一备选血管分叉点进行分析，以去除血管的末端点，得到准确的血管分叉点。

9.如权利要求8所述的肺部CT图像的血管分叉点提取***，其特征在于，所述血管像素点确定模块包括：

平滑子模块，用于利用多尺度高斯函数对肺CT图像进行平滑，若尺度为σ的三维高斯函数为G(x,y,z；σ)，所述肺CT图像为I(x,y,z)，所述肺CT图像I(x,y,z)在所述尺度σ下的平滑结果为I_σ(x,y,z)，则所述利用多尺度高斯函数对肺CT图像进行平滑的表示为： $I_{\mathrm{\σ}}(x,y,z)=I(x,y,z)\⊗G(x,y,z;\mathrm{\σ}),$ 其中， $G(x,y,z;\mathrm{\σ})=\frac{1}{{\left(\sqrt{3\mathrm{\π}}\mathrm{\σ}\right)}^3}{e}^{-\frac{x^2+y^2+z^2}{2{\mathrm{\σ}}^2}};$

第一计算子模块，用于在每一尺度下，根据平滑结果计算所述肺CT图像中每一点的Hessian矩阵；

第二计算子模块，用于对每一点的所述Hessian矩阵进行特征值分解，得到三个特征值以及与三个特征值分别一一对应的特征向量；

估算子模块，用于根据每一点的所述特征值和特征向量，估计对应点在每一尺度下属于管状结构的可能性；

取值子模块，用于取每一点在不同尺度下的可能性的最大值作为对应点属于管状结构的可能性的最终值；

弥散子模块，用于利用弥散函数对所述肺CT图像中每一可能性大于0的点的强度进行更新，直到更新次数达到最大迭代次数为止，得到增强后的血管强度；

第一判定子模块，用于判定所述增强后的血管强度达到设定阈值的点为所述属于血管的像素点。

10.如权利要求8所述的肺部CT图像的血管分叉点提取***，其特征在于，所述血管分叉点提取模块包括：

集合获取子模块，用于以所述备选血管分叉点为中心，取一投票窗口内所有属于血管的像素点的坐标，构成一集合；

第三计算子模块，用于计算所述集合内的坐标均值(v_x,v_y,v_z)，以及协方差矩阵C；

主方向获取子模块，用于对所述协方差矩阵C进行特征值分解：

其中，λ_1,C、λ_2,C、λ_3,C分别为特征值分解后的三个特征值，且满足|λ_1,C|≤|λ_2,C|≤|λ_3,C|，分别为与所述三个特征值分别一一对应的三个特征向量，并取最大特征值|λ_3,C|对应的特征向量作为主方向；

第二判定子模块，用于以所述主方向为法线、过所述备选血管分叉点作一平面，统计处于所述平面一侧的属于血管的像素点个数、以及处于所述平面另一侧的属于血管的像素点个数，根据处于所述平面两侧的属于血管的像素点个数的比值，区分血管的末端点和血管分叉点。7 -->