[发明专利]一种胸部X光图像中骨骼抑制的方法

权利要求书

1.一种胸部X光图像中骨骼抑制的方法，具体步骤如下：

步骤1：模型肺部初始轮廓位置的确定；

步骤1.1：标记训练集中每张图像肺边界的边界点；

步骤1.2：对齐训练形状向量；通过缩放、旋转和平移的操作使训练形状对齐，使它们尽可能对齐紧密，设x_i是训练集中第i个形状中n个点的向量，x_i=(x_i0,y_i0,...x_ik,y_ik,...,x_in-1,y_in-1)^T，其中，(x_ik,y_ik)是第i个形状中第k个点，给定两个相似形状x_i和x_j，选择旋转角度θ、缩放s、平移(t_x,t_y)，则M(s,θ)[x]代表旋转角度为θ和缩放比例为s的变换，使以下加权和最小化将x_i映射为x_j：

E_j=(x_i-M(s,θ)[x_j]-t)^TW(x_i-M(s,θ)[x_j]-t)               (1)

其中，t=(t_x,t_y,...,t_x,t_y)^T，式中:W是一个点的加权对角矩阵，它由每一边界点的权值w_k构成。令R_kl表示第k个点和第l个点的距离，表示训练集中距离的变化，则第k个点的权值

步骤1.3：建立肺部初始轮廓位置的模型；训练形状向量对齐后，利用主成分分析方法找出形状变化的统计信息，据此建立模型，

其特征在于：还包括如下步骤：

步骤2：结合灰度信息和形状信息的肺实质分割：同时利用多张特征图像中边界点的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似；

步骤2.1：提取特征图像；进行主成分分析对数据进行降维，对于图像I中点p附近的点p+△p的灰度可由Taylor展开获得，I(p+Δp)≈Σn=0N1n!(Δpi1,···,Δpin∂nI(p)∂i1···∂in)---(2)]]>

由此可得，图像I的二阶Taylor展开为：

I(p+Δp)=Σn=021n!(Δpi1,···,Δpin∂nI(p)∂i1···∂in)=I(p)+Δpi1∂I(p)∂i1+12!(Δpi1Δpi2∂2I(p)∂i1∂i2)=I(p)+Δpx∂I(p)∂x+12!(Δpx2∂2I(p)∂x2+ΔpxΔpy∂2I(p)∂x∂y)+Δpy∂I(p)∂y+12!(Δpy2∂2I(p)∂y2+ΔpyΔpx∂2I(p)∂y∂x)---(3)]]>

对图像进行高斯平滑处理，对于图像I(p)进行高斯滤波后的图像为A(p,σ)=(I*G)(p,σ)，滤波图像的n阶偏导数为则根据卷积性质可知：

Ai1,···in(p,σ)=(I*G)i1,···in(p,σ)=(I*Gi1,···in)(p,σ)---(5)]]>

式中：角标i₁,...i_n表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；

利用局部2-jet特征图像来获取边界点的候选点及各个候选点的灰度代价，局部2-jet特征为：

J²[I](p,σ)={I,I_x,I_y,I_xx,I_xy,I_yy}            (7)

步骤2.2：选取边界点的候选点：对于初始肺边界的每一个点，计算所有特征图像中该点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度，选出相似程度最大的点，作为该边界点的候选点。相似程度为所有特征图像中该点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度集合的马氏距离h_i：

步骤2.3：使用动态规划进行肺区域分割；在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度h_i在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度h_i；图像中第i个边界点的形状相似性代价定义为：

fi=(vi-uvi)T(COVvi)-1(vi-uvi)---(8)]]>

式中，v_i=p_i+1-p_i，表示第i个边界点的形状特征，分别表示所有训练图像中第i个边界点形状特征的均值与协方差，

对于测试图像中第i个边界点p_i，在指定的搜索区域内存在m个具有较小灰度相似性代价候选点，则n个边界点就会产生一个n×m的灰度代价矩阵：

C=h1,1···h1,k···h1,m·········hi,1···hi,k···hi,m·········hn,1···hn,k···hn,m---(9)]]>

hi=Σj=1N(gij-ugij)T(sgij)-1(gij-ugij)---(10)]]>

gij=pi+rccos(2πnc(k-1))sin(2πnc(k-1))---(11)]]>

式中：为在特征图像上，以点p_i为圆心，r_c为半径的圆上的n_c个点的灰度，k=1......n_c。分别为训练图像的第j个特征图像中第i个边界点的周围像素点灰度的均值及协方差，N为特征图像总数，

搜索最优轮廓即找一条最佳路径，在矩阵C中每行选一个元素，沿着路径选择的时候，灰度与形状相似性代价的总和最小，即：

J(k1.....kn)=Σi=1nhi+γΣi=1nfi---(12)]]>

其中，γ为灰度与形状相似性代价系数。调整γ值，使得这两种代价在边界搜索过程中发挥大致相同的作用；

步骤3：使用B样条多尺度小波变换的特征图像提取；在模型建立过程中，提取有效描述不同尺度骨骼的特征图像是模型建立的基础；

步骤3.1：B样条多尺度小波变换：B样条函数随着样条阶数的增加而快速收敛于高斯函数，其一阶导数可以逼近最优边缘检测算子；利用B样条小波进行多尺度边缘增强；利用函数：

H(eiω)=Σhke-ikω=θ^m(2ω)θ^m(ω)=e-iϵω(sinωω)m+1e-iϵω2(sin(ω/2)ω/2)m+1=e-iϵω2(cosω2)m+1---(18)]]>

G(ω)=ψm^(2ω)θ^m(ω)=Σkgke-ikω---(20)]]>

求取B样条小波变换滤波器系数：

式中：h_k为低通高斯滤波器系数，g_k为高通高斯滤波器系数，利用所得h_k和g_k对图像进行快速小波分解，原始图像的行分别与h_k和g_k进行卷积，再对其列进行下采样，可以得到两幅子图像，然后沿着列的方向对两幅子图像做滤波并进行下采样，就可以得到四幅1/4大小的输出子图像，对角线细节图像，垂直细节图像，水平细节图像和近似图像；

步骤3.2：高斯滤波局部2-jet特征图像提取，由于高斯偏导数构造的特征向量能用来描述图像的局部特征，根据公式(7)式中，提取不同尺度小波变换后的细节图像的2-jet特征(I,I_x,I_y,I_xx,I_xy,I_yy)图来建立回归模型；

步骤4：使用支持向量回归进行骨骼图像的预测，通过求得回归函数f(x)=ω·x+b，建立骨骼模型；设样本集为(x_i,y_i)，i=1,2,...,l，x_i∈Rⁿ，样本集S是ε-不敏感损失函数的线性近似。当线性回归函数f(x)=ω·x+b拟合(x_i,y_i)时，假设所有训练数据的拟合误差精度为ε，即：

|y_i-f(x_i)|≤εi=1,2,...,l                    (22)

让d_i表示点(x_i,y_i)∈S到超平面f(x)的距离：

di=|<w,x>+b-y|1+||w||2---(23)]]>

因为S集是ε-线性近似，所以有|<w,x>+b-f(x_i)|≤ε,i=1,...,m。可以得到:

|<w,x>+b-y|1+||w||2≤ϵ1+||w||2,i=1,...,l---(24)]]>

于是有

di≤ϵ1+||w||2,i=1,···,l---(25)]]>

(23)式表明是S中的点到超平面的距离的上界。ε-不敏感损失函数的线性近似集S的最优近似超平面是通过最大化S中的点到超平面距离的上界而得到超平面；

最优超平面是通过最大化得到的(即最小化)。因此，只要最小化||w||²就可以得到最优近似超平面。于是线性回归问题就化为：

求下面的最优化问题：

min12||w||2---(26)]]>

约束为|<w,x>+b-y_i|≤ε,i=1,...,m。另外，考虑拟合误差的情况，引入松弛因子当划分有误差时，都大于0，误差不存在取0；

因而，标准ε不敏感支持向量回归机可以表示为：

minw,b,ξ12||w||2+CΣi=1l(ξi+ξi*)s.t.yi-w·xi-b≤ϵ+ξiw·xi+b-yi≤ϵ+ξi*ξi≥0ξi*≥0,i=1,2,...,l---(27)]]>

其中，C>0为平衡因子。引入拉格朗日函数为：

L(w,b,ξ,α,α*,γ)=12||w||2+CΣi=1n(ξi+ξi*)-Σi=1nαi[ξi+ϵ-yi+f(xi)]-Σi=1nαi*[ξi+ϵ+yi-f(xi)]-Σi=1nγi(ξiγi+ξi*γi*)---(28)]]>

其中，α_i，γ_i≥0，i=1,...,n。函数L的极值应满足条件于是得到下面的式子：

w=Σi=0n(αi-αi*)xi---(29)]]>

Σi=1n(αi-αi*)=0---(30)]]>

C-α_i-γ_i=0,i=1,...,l                   (31)

C-αi*-γi*=0,i=1,...,l---(32)]]>

将(29)-(32)代入到(28)中，得到优化问题的对偶形式为：

max-12Σi,j=1n(αi-αi*)(αj-αj*)<xi,xj>+Σi=1n(αi-αi*)yi-Σi=1n(αi+αi*)ϵ---(33)]]>

约束为：

Σi=1n(αi-αi*)=0,0≤αi,αi*≤C,i=1,...,l---(34)]]>

由于骨骼灰度预测为非线性回归，所以先使用一个非线性映射φ，把数据映射到一个高维特征空间，然后在高维特征空间进行线性回归建立模型。由于在上面的优化过程中，只考虑到高维特征空间中的内积运算，因此用一个核函数k(x,y)代替<φ(x),φ(y)>就可以实现非线性回归。于是，非线性回归的优化方程为最大化下面的函数：

W(α,α*)=-12Σi,jn(αi-αi*)(αj-αj*)K(xi,xj)+Σi,jn(αi-αi*)yi-Σi,jn(αi+αi*)ϵ---(35)]]>

其约束条件为(32)，求解出α的值后，就可以得到f(x)：

f(x)=Σi=1N(αi-αi*)K(xi,x)+b---(36)]]>

通常情况下，大部分α_i或的值将为零，不为零的α_i或所对应的样本称为支持向量，根据最优化条件，在鞍点有：

αi[ξi+ϵ-yi+f(xi)]=0,i=1,...lai*[ξi*+ϵ-yi+f(xi)]=0,i=1,...l(C-αi)ξi=0,i=1,...l(C-αi*)ξi*=0,i=1,...l---(37)]]>

于是得到b的计算式如下：

b=yi-ϵ-Σi=1l(αi-αi*)K(xj,xi),]]>当α_j∈(0,C)

b=yi+ϵ-Σi=1l(αi-αi*)K(xj,xi),]]>当

用任意一个支持向量可以计算出b的值，得到回归函数f(x)；建立预测模型后，可以根据肺部x光图像的灰度值分布来预测产生骨结构图像；

步骤5：对正常胸片与预测的骨骼图像做图像减法来预测的软组织图像。