一种胸部X光图像中骨骼抑制的方法

摘要

本发明属于X线图像处理技术领域，具体涉及一种胸部X光图像中骨骼抑制的方法。本发明对双能剪影图像中的正常胸片建立肺部轮廓模型，再应用三阶B样条小波变换对图像做三尺度小波分解，将分解后的图像使用2-jet方法提取特征图像，再使用支持向量回归机建立骨骼模型，将骨骼模型应用到X光图像中，从而提取对应的预测骨骼图像，最后用X光图像与预测的骨骼图像做图像减法从而得到预测的软组织图像。因为抑制骨骼方法大都采用神经网络方法，而此方法易陷入局部最优，训练结果不太稳定，所以本专利采用支持向量回归机来进行骨骼预测，得到了更优的结果。

说明书

技术领域

本发明属于X线图像处理技术领域，具体涉及一种胸部X光图像中骨骼抑制的方法。 

背景技术

肺癌是当今对人类健康危害最大的恶性肿瘤之一。特别是近半个世纪以来，随着空气污染造成的环境不断恶化及吸烟人口数量的大量增加，各国肺癌的发病率和病死率都在急剧上升。由于肺属于人体内部脏器，多数肺癌在开始的时候只是在身体内悄悄生长，患者没有任何感觉。当患者因咳嗽、咯血及胸痛等临床症状就诊时，大多数已经处于中晚期，错过了治疗的最佳时机。所以，肺癌的早期诊断与治疗是提高肺癌患者生存率的关键。 

医学影像检查不仅能通过图像发现病灶的存在，定位病灶在全肺的位置，还能观察病灶的大小、形态及密度等物理特征。胸部疾病的影像诊断是以胸部x线摄影为主体的，但在胸部x线图像中，经常会由于病症位于图像中的肋骨或者锁骨结构遮挡区域而造成漏诊。 

双能量减影(Dual Energy Subtraction,DES)是在数字胸部x线摄影基础上发展的一种较新的成像技术。它利用骨与软组织对x线光子的能量衰减方式不同，以及不同原子量的物质的光电吸收效应差别，利用数字摄影将两种吸收效应的信息进行分离，选择性去除骨或软组织的衰减信息，进而获得纯粹的软组织像和骨组织像。因此应用双能量剪影技术可以得到三张图像，一张正常胸片，一张软组织图像和一张骨骼图像。 

DES的软组织减影图像对胸部小结节性病变，钙化结节性病变，肺纹理病变，特别是位于肺肋骨、锁骨重叠处结节性病变的清晰显示率明显高于普通胸部x光图像。这主要是减影区的肺部软组织图像除去了骨骼等重叠因素的影响，对于肋骨，锁骨相重叠的病灶显示清晰，这对于肋软骨钙化较多的患者更为有益。 

国内外的一些研究机构对胸部x光图像中骨骼抑制问题进行研究。一些国外学者利用正常的胸部x光图像来估计骨骼结构模型，产生骨骼抑制图像来图高病变的检测性能。Ahmed及Rasheed等人提出利用独立成分分析(Independent Component Analysis,ICA)方法来抑制后侧肋骨，进而提高结节的可见度。Jiann等人提出一种非参数的正常胸片中肋骨抑制算法。该算法由肺分割、肋骨分割、肋骨灰度估计及抑制等部分组成，提出利用实数编码遗传算法(Real-coded genetic algorithm,RCGA)来估计肋骨的灰度，进而实现肋骨的抑制。 

还有一些学者利用DES图像中的正常胸部图像及骨骼图像，构造回归模型，进而达到生成骨骼抑制图像的目的。Loog等人提出基于非线性回归的图像滤波算法来从x光图像中获得骨骼结构及软组织图像。该算法通使用k最近邻回归算法建立DES肺部图像与骨骼图像的回归模型。Kenji Suzuki等人提出基于多分辨率大规模训练人工神经元网络(Multiresolution massive training artificial neural network,MTANN)的图像处理技术来抑制x光图像中骨骼区域的对比度。该算法采用双能图像中的x光图像与骨结构图像作为训练数据，利用MTANN建立回归模型。在生成的骨骼抑制图像中，结节及血管等解剖结构清晰可见。 

综上所述，国内外对于胸部x光图像中骨骼抑制的研究才刚刚兴起。原有的大多数研究 是基于普通胸部x光图像中的骨骼提取、分割、模型建立及移去等。但是，由于普通x光图像中骨骼结构复杂，且边缘区域灰度的对比度较低等特点，仅利用普通x光图像实现骨骼的分割与抑制比较困难。而利用双能图像的骨骼回归模型可以从一个完全下颖的角度来获得骨骼图像，不要再对分割算法进行研究。由于在模型建立过程中，使用大量DES骨骼图像信息，因此获得的预测图像精度较高。 

发明内容

针对现有技术存在的不足，本发明提供一种稳定度高的胸部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) 

其中，$M(s,\theta )\left(\begin{array}{c}{x}_{\mathrm{jk}}\\ y_{\mathrm{jk}}\end{array}\right)=\left(\begin{array}{c}\left(s\cos \theta \right)x_{\mathrm{jk}}-\left(s\sin \theta \right)y_{\mathrm{jk}}\\ \left(s\sin \theta \right)x_{\mathrm{jk}}+\left(s\cos \theta \right)y_{\mathrm{jk}}\end{array}\right),$t=(t_x,t_y,...,t_x,t_y)^T。式中：W是一个点的加权对角矩阵，它由每一边界点的权值w_k构成。令R_kl表示第k个点和第l个点的距离，表示训练集中距离的变化，则第k个点的权值

如果令a_x=scosθ,a_y=ssinθ，则根据最小二乘法有四个线性等式 

$\left(\begin{array}{cccc}{X}_2& -Y_2& W& 0\\ Y_2& X_2& 0& W\\ Z& 0& X_2& Y_2\\ 0& Z& {-Y}_2& X_2\end{array}\right)\left(\begin{array}{c}{a}_x\\ a_y\\ t_x\\ t_y\end{array}\right)=\left(\begin{array}{c}{X}_1\\ Y_1\\ C_1\\ C_2\end{array}\right)---\left(4\right)$

其中，$X_i=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k{x}_{\mathrm{ik}},$$Y_i=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k{y}_{\mathrm{ik}},$$Z=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k(x_{2k}^2+y_{2k}^2),$$C_1=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k(x_{1k}{x}_{2k}+y_{1k}{y}_{2k}),$$C_2=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k(y_{1k}{x}_{2k}+x_{1k}{y}_{2k}),$$w=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k.$可以使用矩阵运算方法解出变量a_x,a_y,t_x,t_y。 

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

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

步骤2.1：提取特征图像；进行主成分分析对数据进行降维，对于图像I中点p附近的点p+△p的灰度可由Taylor展开获得，$I(p+\mathrm{\Delta p})\approx \underset{n=0}{\overset{N}{\Sigma }}\frac{1}{n!}(\Delta p_{i1},···,\Delta p_{\mathrm{in}}\frac{{\partial }^{n}I\left(p\right)}{{\partial i}_{1···}{\partial i}_n})---\left(2\right)$

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

$\left(\begin{array}{c}I(p+\mathrm{\Delta p})=\underset{n=0}{\overset{2}{\Sigma }}\frac{1}{n!}(\Delta p_{i1},···,\Delta p_{\mathrm{in}}\frac{{\partial }^{n}I\left(p\right)}{{\partial i}_1···{\partial i}_n})=I\left(p\right)+\Delta p_{i1}\frac{\partial I\left(p\right)}{{\partial i}_1}+\frac{1}{2!}\left(\Delta p_{i1}\Delta p_{i2}\frac{{\partial }^{2}I\left(p\right)}{{\partial i}_1{\partial i}_2}\right)\\ =I\left(p\right)+\Delta p_x\frac{\partial I\left(p\right)}{\partial x}+\frac{1}{2!}(\Delta p_x^2\frac{{\partial }^{2}I\left(p\right)}{{\partial x}^2}+\Delta p_x\Delta p_y\frac{{\partial }^{2}I\left(p\right)}{\partial x\partial y})\\ +\Delta p_y\frac{\partial I\left(p\right)}{\partial y}+\frac{1}{2!}(\Delta p_y^2\frac{{\partial }^{2}I\left(p\right)}{{\partial y}^2}+\Delta p_y\Delta p_x\frac{{\partial }^{2}I\left(p\right)}{\partial y\partial x})\end{array}\right)---\left(3\right)$

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

$A_{i_1,···i_n}(p,\sigma )={(I*G)}_{i_1,···i_n}(p,\sigma )=(I*G_{i_1,···i_n})(p,\sigma )---\left(5\right)$

式中：角标i₁,...i_n表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；对于一个给定尺度参数σ的n-jet特征可以定义为： 

$J^N\left[I\right](p,\sigma )=\left\{I_{i_1,...i_n}\right|n=0,...,N\},N\in Z---\left(6\right)$

由此可以看出，当N→∞时，在某一尺度空间中，图像I的局部特征由所有偏导数的组合构成。因此，可以采用高斯偏导数构造一个特征向量来描述局部特征。当偏导数阶数越高时，对图像的描述越准确。但是，如果使用过多的偏导数描述图像，会使得特征向量维数增加。 

利用局部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个边界点的形状相似性代价定义为： 

$f_i={(v_i-u_{v_i})}^T{\left({\mathrm{COV}}_{v_i}\right)}^{-1}(v_i-u_{v_i})---\left(8\right)$

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

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

$C=\left(\begin{array}{ccccc}{h}_{\mathrm{1,1}}& ···& h_{1,k}& ···& h_{1,m}\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ h_{i,1}& ···& h_{i,k}& ···& h_{i,m}\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ h_{n,1}& ···& h_{n,k}& ···& h_{n,m}\end{array}\right)---\left(9\right)$

$h_i=\underset{j=1}{\overset{N}{\Sigma }}{(g_i^j-u_{g_i^j})}^T{\left(s_{g_i^j}\right)}^{-1}(g_i^j-u_{g_i^j})---\left(10\right)$

$g_i^j=p_i+r_c\left(\begin{array}{c}\cos \left(\frac{2\pi }{n_c}(k-1)\right)\\ \sin \left(\frac{2\pi }{n_c}(k-1)\right)\end{array}\right)---\left(11\right)$

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

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

$J(k_1.....k_n)=\underset{i=1}{\overset{n}{\Sigma }}{h}_i+\gamma \underset{i=1}{\overset{n}{\Sigma }}{f}_i---\left(12\right)$

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

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

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

$N_m\left(x\right)=\frac{1}{(m-1)!}\underset{t=0}{\overset{m}{\Sigma }}{(-1)}^t\left(\begin{array}{c}m\\ t\end{array}\right){(x-t)}_{+}^{m-1},\forall n\in Z^{+}---\left(13\right)$

其中，$\left(\begin{array}{c}m\\ t\end{array}\right)=\frac{t!}{m!(t-m)!}.$由此可知，零阶B样条为： 

当m≥1时，m阶B样条函数N_m(x)可用卷积迭代关系表示： 

显然，N_m(x)是非负且紧支撑的，其支撑宽度为suppN_m=[0,m+1]，支撑中心为(m+1)/2，且N_m关于(m+1)/2中心对称。 

零价B样条的傅里叶变换为则m次B样条的傅里叶形式为： 

${\hat{N}}_m\left(\omega \right)=e^{-i\frac{m+1}{2}\omega }{\left(\frac{\sin (\omega /2)}{\omega /2}\right)}^{m+1}---\left(16\right)$

若将m次B样条基N_m(x)平移到N_m[x+(m+1)/2]，则对称中心移至坐标原点，傅里叶变换将不再出现上式中的线性相位。但是，当m为偶数时，将出现半整数节点。为避免出现半整数节点，当m为偶数时，将N_m(x)平移到N_m[x+(m+1)/2]，则对称中心移至1/2。将N_m(x)整数平移后的函数记为θ_m(x)，其傅里叶变换为： 

${\hat{\theta }}_m\left(\omega \right)=e^{-i\frac{\mathrm{ϵ\omega }}{2}}{\left(\frac{\sin (\omega /2)}{\omega /2}\right)}^{m+1}---\left(17\right)$

由多分辨率分析得到二尺度关系式可求得： 

$H\left(e^{\mathrm{i\omega }}\right)=\Sigma h_k{e}^{-\mathrm{ik\omega }}=\frac{{\hat{\theta }}_m\left(2\omega \right)}{{\hat{\theta }}_m\left(\omega \right)}=\frac{e^{-\mathrm{iϵ\omega }}{\left(\frac{\sin \omega }{\omega }\right)}^{m+1}}{e^{-i\frac{\mathrm{ϵ\omega }}{2}}{\left(\frac{\sin (\omega /2)}{\omega /2}\right)}^{m+1}}=e^{-i\frac{\mathrm{ϵ\omega }}{2}}{\left(\cos \frac{\omega }{2}\right)}^{m+1}---\left(18\right)$

式中，当m为奇数时，ε=0；当m为偶数时，ε=1。由欧拉公式和二项式展开可得低通滤波器为： 

当m为奇数时，

当m为偶数时，

选取B样条函数的一阶导数作为小波函数，即m+1次B样条函数在2^-1尺度上一阶微分函数为m次B样条小波函数： 

${\psi }_m\left(x\right)=2\frac{d}{\mathrm{dx}}{\theta }_{m+1}\left(2x\right)---\left(19\right)$

其傅里叶变换为且ψ_m(x)满足小波的容许性条件。由多分辨率分析中小波函数与尺度函数的二尺度关系可求得： 

$G\left(\omega \right)=\frac{\widehat{{\psi }_m}}{{\hat{\theta }}_m\left(\omega \right)}=\underset{k}{\Sigma }{g}_k{e}^{-\mathrm{ik\omega }}---\left(20\right)$

当m是奇数时，$G\left(\omega \right)={4\mathrm{ie}}^{-\frac{\mathrm{i\omega }}{2}}\sin \frac{\omega }{2};$当m是偶数时，$G\left(\omega \right)={4\mathrm{ie}}^{\frac{\mathrm{i\omega }}{2}}\sin \frac{\omega }{2}.$由此可知，时域有限脉冲响应为： 

当m=3时，各项滤波器系数为k=-2,-1,0,1,2，g_k={0,0,-2,2,0}。由此可见，三阶B样条小波的滤波器的有限长度为5。 

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

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

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

$d_i=\frac{|<w,x>+b-y|}{\sqrt{1+{\left|\right|w\left|\right|}^2}}---\left(23\right)$

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

$\frac{|<w,x>+b-y|}{\sqrt{1+{\left|\right|w\left|\right|}^2}}\le \frac{ϵ}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,...,l---\left(24\right)$

于是有 

$d_i\le \frac{ϵ}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,···,l---\left(25\right)$

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

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

求下面的最优化问题： 

$\min \frac{1}{2}{\left|\right|w\left|\right|}^2---\left(26\right)$

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

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

$\left(\begin{array}{c}\underset{w,b,\xi }{\min }\frac{1}{2}{\left|\right|w\left|\right|}^2+C\underset{i=1}{\overset{l}{\Sigma }}({\xi }_i+{\xi }_i^{*})\\ s.t.y_i-w·x_i-b\le ϵ+{\xi }_i\\ w·x_i+b-y_i\le ϵ+{\xi }_i^{*}\\ {\xi }_i\ge 0\\ {\xi }_i^{*}\ge 0,i=\mathrm{1,2},...,l\end{array}\right)---\left(27\right)$

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

$\left(\begin{array}{c}L(w,b,\xi ,\alpha ,{\alpha }^{*},\gamma )=\frac{1}{2}{\left|\right|w\left|\right|}^2+C\underset{i=1}{\overset{n}{\Sigma }}({\xi }_i+{\xi }_i^{*})-\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i[{\xi }_i+ϵ-y_i+f\left(x_i\right)]\\ -\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i^{*}[{\xi }_i+ϵ+y_i-f\left(x_i\right)]-\underset{i=1}{\overset{n}{\Sigma }}{\gamma }_i({\xi }_i{\gamma }_i+{\xi }_i^{*}{\gamma }_i^{*})\end{array}\right)---\left(28\right)$

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

$w=\underset{i=0}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})x_i---\left(29\right)$

$\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})=0---\left(30\right)$

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

$C-{\alpha }_i^{*}-{\gamma }_i^{*}=0,i=1,...,l---\left(32\right)$

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

$\max -\frac{1}{2}\underset{i,j=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})({\alpha }_j-{\alpha }_j^{*})<x_i,x_j>+\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})y_i-\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i+{\alpha }_i^{*})ϵ---\left(33\right)$

约束为： 

$\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})=\mathrm{0,0}\le {\alpha }_i,{\alpha }_i^{*}\le C,i=1,...,l---\left(34\right)$

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

$\left(\begin{array}{c}W(\alpha ,{\alpha }^{*})=-\frac{1}{2}\underset{i,j}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})({\alpha }_j-{\alpha }_j^{*})K(x_i,x_j)\\ +\underset{i,j}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})y_i-\underset{i,j}{\overset{n}{\Sigma }}({\alpha }_i+{\alpha }_i^{*})ϵ\end{array}\right)---\left(35\right)$

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

$f\left(x\right)=\underset{i=1}{\overset{N}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(x_i,x)+b---\left(36\right)$

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

$\left(\begin{array}{cc}{\alpha }_i[{\xi }_i+ϵ-y_i+f\left(x_i\right)]=0& ,i=1,...l\\ a_i^{*}[{\xi }_i^{*}+ϵ-y_i+f\left(x_i\right)]=0& ,i=1,...l\\ (C-{\alpha }_i){\xi }_i=0& ,i=1,...l\\ (C-{\alpha }_i^{*}){\xi }_i^{*}=0& ,i=1,...l\end{array}\right)---\left(37\right)$

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

$b=y_i-ϵ-\underset{i=1}{\overset{l}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(x_j,x_i),$当α_j∈(0,C) 

$b=y_i+ϵ-\underset{i=1}{\overset{l}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(x_j,x_i),$当∈(0,C) 

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

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

所述步骤1.2：对齐训练形状向量的具体步骤如下： 

所述的将分块的数据映射到特征空间的具体方法如下： 

步骤1.2.1：旋转、缩放和平移每个肺区域形状，使其与训练集中的第一个形状对齐； 

步骤1.2.2：根据对齐形状，计算平均形状； 

步骤1.2.3：旋转、缩放和平移平均形状使其与第一个形状对齐； 

步骤1.2.权：重新将每个形状与当前平均形状对齐； 

步骤1.2.5：如果过程收敛或者到指定循环次数，退出；否则转到步骤。 

所述步骤3中对胸部x光图像做三尺度小波分解。 

本发明的优点为：本发明先对双能剪影图像中的正常胸片建立肺部轮廓模型，再应用三阶B样条小波变换对图像做三尺度小波分解，将分解后的图像使用2-jet方法提取特征图像， 再使用支持向量回归机建立骨骼模型，将骨骼模型应用到X光图像中，从而提取对应的预测骨骼图像，最后用X光图像与预测的骨骼图像做图像减法从而得到预测的软组织图像。因为抑制骨骼方法大都采用神经网络方法，而此方法易陷入局部最优，训练结果不太稳定，所以本专利采用支持向量回归机来进行骨骼预测，得到了更优的结果。 

附图说明

图1是本发明的总体流程图； 

图2是本发明的肺区分割流程图； 

图3是本发明中高斯滤波局部2-jet提取特征图像示意图； 

图4是本发明中最优近似超平面的示意图； 

图5是本发明中支持向量回归示意图。 

具体实施方式

如图1至图5所示，本发明采取按如下步骤进行： 

步骤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) 

其中，$M(s,\theta )\left(\begin{array}{c}{x}_{\mathrm{jk}}\\ y_{\mathrm{jk}}\end{array}\right)=\left(\begin{array}{c}\left(s\cos \theta \right)x_{\mathrm{jk}}-\left(s\sin \theta \right)y_{\mathrm{jk}}\\ \left(s\sin \theta \right)x_{\mathrm{jk}}+\left(s\cos \theta \right)y_{\mathrm{jk}}\end{array}\right),$t=(t_x,t_y,...,t_x,t_y)^T式中：W是一个点的加权对角矩阵，它由每一边界点的权值w_k构成。令R_kl表示第k个点和第l个点的距离，表示训练集中距离的变化，则第k个点的权值

如果令a_x=scosθ,a_y=ssinθ，则根据最小二乘法有四个线性等式 

$\left(\begin{array}{cccc}{X}_2& -Y_2& W& 0\\ Y_2& X_2& 0& W\\ Z& 0& X_2& Y_2\\ 0& Z& {-Y}_2& X_2\end{array}\right)\left(\begin{array}{c}{a}_x\\ a_y\\ t_x\\ t_y\end{array}\right)=\left(\begin{array}{c}{X}_1\\ Y_1\\ C_1\\ C_2\end{array}\right)---\left(4\right)$

其中，$X_i=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k{x}_{\mathrm{ik}},$$Y_i=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k{y}_{\mathrm{ik}},$$Z=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k(x_{2k}^2+y_{2k}^2),$$C_1=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k(x_{1k}{x}_{2k}+y_{1k}{y}_{2k}),$$C_2=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k(y_{1k}{x}_{2k}+x_{1k}{y}_{2k}),$$w=\underset{k=0}{\overset{n-1}{\Sigma }}{w}_k.$可以使用矩阵运算方法解出变量a_x,a_y,t_x,t_y。 

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

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

步骤2.1：提取特征图像；进行主成分分析对数据进行降维，对于图像I中点p附近的点p+△p的灰度可由Taylor展开获得，$I(p+\mathrm{\Delta p})\approx \underset{n=0}{\overset{N}{\Sigma }}\frac{1}{n!}(\Delta p_{i1},···,\Delta p_{\mathrm{in}}\frac{{\partial }^{n}I\left(p\right)}{{\partial i}_{1···}{\partial i}_n})$

                                                                       (2) 

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

$\left(\begin{array}{c}I(p+\mathrm{\Delta p})=\underset{n=0}{\overset{2}{\Sigma }}\frac{1}{n!}(\Delta p_{i1},···,\Delta p_{\mathrm{in}}\frac{{\partial }^{n}I\left(p\right)}{{\partial i}_1···{\partial i}_n})=I\left(p\right)+\Delta p_{i1}\frac{\partial I\left(p\right)}{{\partial i}_1}+\frac{1}{2!}\left(\Delta p_{i1}\Delta p_{i2}\frac{{\partial }^{2}I\left(p\right)}{{\partial i}_1{\partial i}_2}\right)\\ =I\left(p\right)+\Delta p_x\frac{\partial I\left(p\right)}{\partial x}+\frac{1}{2!}(\Delta p_x^2\frac{{\partial }^{2}I\left(p\right)}{{\partial x}^2}+\Delta p_x\Delta p_y\frac{{\partial }^{2}I\left(p\right)}{\partial x\partial y})\\ +\Delta p_y\frac{\partial I\left(p\right)}{\partial y}+\frac{1}{2!}(\Delta p_y^2\frac{{\partial }^{2}I\left(p\right)}{{\partial y}^2}+\Delta p_y\Delta p_x\frac{{\partial }^{2}I\left(p\right)}{\partial y\partial x})\end{array}\right)---\left(3\right)$

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

$A_{i_1,···i_n}(p,\sigma )={(I*G)}_{i_1,···i_n}(p,\sigma )=(I*G_{i_1,···i_n})(p,\sigma )---\left(5\right)$

式中：角标i₁,...i_n表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；对于一个给定尺度参数σ的n-jet特征可以定义为： 

$J^N\left[I\right](p,\sigma )=\left\{I_{i_1,...i_n}\right|n=0,...,N\},N\in Z---\left(6\right)$

由此可以看出，当N→∞时，在某一尺度空间中，图像I的局部特征由所有偏导数的组合构成。因此，可以采用高斯偏导数构造一个特征向量来描述局部特征。当偏导数阶数越高时，对图像的描述越准确。但是，如果使用过多的偏导数描述图像，会使得特征向量维数增加。 

利用局部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个边界点的形状相似性代价定义为： 

$f_i={(v_i-u_{v_i})}^T{\left({\mathrm{COV}}_{v_i}\right)}^{-1}(v_i-u_{v_i})---\left(8\right)$

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

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

$C=\left(\begin{array}{ccccc}{h}_{\mathrm{1,1}}& ···& h_{1,k}& ···& h_{1,m}\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ h_{i,1}& ···& h_{i,k}& ···& h_{i,m}\\ ·& & ·& & ·\\ ·& & ·& & ·\\ ·& & ·& & ·\\ h_{n,1}& ···& h_{n,k}& ···& h_{n,m}\end{array}\right)---\left(9\right)$

$h_i=\underset{j=1}{\overset{N}{\Sigma }}{(g_i^j-u_{g_i^j})}^T{\left(s_{g_i^j}\right)}^{-1}(g_i^j-u_{g_i^j})---\left(10\right)$

$g_i^j=p_i+r_c\left(\begin{array}{c}\cos \left(\frac{2\pi }{n_c}(k-1)\right)\\ \sin \left(\frac{2\pi }{n_c}(k-1)\right)\end{array}\right)---\left(11\right)$

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

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

$J(k_1.....k_n)=\underset{i=1}{\overset{n}{\Sigma }}{h}_i+\gamma \underset{i=1}{\overset{n}{\Sigma }}{f}_i---\left(12\right)$

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

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

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

$N_m\left(x\right)=\frac{1}{(m-1)!}\underset{t=0}{\overset{m}{\Sigma }}{(-1)}^t\left(\begin{array}{c}m\\ t\end{array}\right){(x-t)}_{+}^{m-1},\forall n\in Z^{+}---\left(13\right)$

其中，$\left(\begin{array}{c}m\\ t\end{array}\right)=\frac{t!}{m!(t-m)!}.$由此可知，零阶B样条为： 

当m≥1时，m阶B样条函数N_m(x)可用卷积迭代关系表示： 

显然，N_m(x)是非负且紧支撑的，其支撑宽度为suppN_m=[0,m+1]，支撑中心为(m+1)/2，且N_m关于(m+1)/2中心对称。 

零阶B样条的傅里叶变换为则m次B样条的傅里叶形式为： 

${\hat{N}}_m\left(\omega \right)=e^{-i\frac{m+1}{2}\omega }{\left(\frac{\sin (\omega /2)}{\omega /2}\right)}^{m+1}---\left(16\right)$

若将m次B样条基N_m(x)平移到N_m[x+(m+1)/2]，则对称中心移至坐标原点，傅里叶变换将不再出现上式中的线性相位。但是，当m为偶数时，将出现半整数节点。为避免出现半整数节点，当m为偶数时，将N_m(x)平移到N_m[x+(m+1)/2]，则对称中心移至1/2。将N_m(x)整数 平移后的函数记为θ_m(x)，其傅里叶变换为： 

${\hat{\theta }}_m\left(\omega \right)=e^{-i\frac{\mathrm{ϵ\omega }}{2}}{\left(\frac{\sin (\omega /2)}{\omega /2}\right)}^{m+1}---\left(17\right)$

由多分辨率分析得到二尺度关系式可求得： 

$H\left(e^{\mathrm{i\omega }}\right)=\Sigma h_k{e}^{-\mathrm{ik\omega }}=\frac{{\hat{\theta }}_m\left(2\omega \right)}{{\hat{\theta }}_m\left(\omega \right)}=\frac{e^{-\mathrm{iϵ\omega }}{\left(\frac{\sin \omega }{\omega }\right)}^{m+1}}{e^{-i\frac{\mathrm{ϵ\omega }}{2}}{\left(\frac{\sin (\omega /2)}{\omega /2}\right)}^{m+1}}=e^{-i\frac{\mathrm{ϵ\omega }}{2}}{\left(\cos \frac{\omega }{2}\right)}^{m+1}---\left(18\right)$

式中，当m为奇数时，ε=0；当m为偶数时，ε=1。由欧拉公式和二项式展开可得低通滤波器为： 

当m为奇数时，

当m为偶数时，

选取B样条函数的一阶导数作为小波函数，即m+1次B样条函数在2^-1尺度上一阶微分函数为m次B样条小波函数： 

${\psi }_m\left(x\right)=2\frac{d}{\mathrm{dx}}{\theta }_{m+1}\left(2x\right)---\left(19\right)$

其傅里叶变换为且ψ_m(x)满足小波的容许性条件。由多分辨率分析中小波函数与尺度函数的二尺度关系可求得： 

$G\left(\omega \right)=\frac{\widehat{{\psi }_m}}{{\hat{\theta }}_m\left(\omega \right)}=\underset{k}{\Sigma }{g}_k{e}^{-\mathrm{ik\omega }}---\left(20\right)$

当m是奇数时，$G\left(\omega \right)={4\mathrm{ie}}^{-\frac{\mathrm{i\omega }}{2}}\sin \frac{\omega }{2};$当m是偶数时，$G\left(\omega \right)={4\mathrm{ie}}^{\frac{\mathrm{i\omega }}{2}}\sin \frac{\omega }{2}.$由此可知，时域有限脉冲响应为： 

当m=3时，各项滤波器系数为k=-2,-1,0,1,2，g_k={0,0,-2,2,0}。由此可见，三阶B样条小波的滤波器的有限长度为5。 

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

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

本实施例对胸部x光图像做三尺度小波分解，获得12张小波分解图像，其中3张近似图像，3张水平细节图像，3张垂直细节图像，3张对角线细节图像。由于在胸部x光图像中，很少有呈现对角线分布的骨骼结构。因此，不提取对角线细节图像的特征。 

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

$d_i=\frac{|<w,x>+b-y|}{\sqrt{1+{\left|\right|w\left|\right|}^2}}---\left(23\right)$

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

$\frac{|<w,x>+b-y|}{\sqrt{1+{\left|\right|w\left|\right|}^2}}\le \frac{ϵ}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,...,l---\left(24\right)$

于是有 

$d_i\le \frac{ϵ}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,···,l---\left(25\right)$

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

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

求下面的最优化问题： 

$\min \frac{1}{2}{\left|\right|w\left|\right|}^2---\left(26\right)$

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

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

$\left(\begin{array}{c}\underset{w,b,\xi }{\min }\frac{1}{2}{\left|\right|w\left|\right|}^2+C\underset{i=1}{\overset{l}{\Sigma }}({\xi }_i+{\xi }_i^{*})\\ s.t.y_i-w·x_i-b\le ϵ+{\xi }_i\\ w·x_i+b-y_i\le ϵ+{\xi }_i^{*}\\ {\xi }_i\ge 0\\ {\xi }_i^{*}\ge 0,i=\mathrm{1,2},...,l\end{array}\right)---\left(27\right)$

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

$\left(\begin{array}{c}L(w,b,\xi ,\alpha ,{\alpha }^{*},\gamma )=\frac{1}{2}{\left|\right|w\left|\right|}^2+C\underset{i=1}{\overset{n}{\Sigma }}({\xi }_i+{\xi }_i^{*})-\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i[{\xi }_i+ϵ-y_i+f\left(x_i\right)]\\ -\underset{i=1}{\overset{n}{\Sigma }}{\alpha }_i^{*}[{\xi }_i+ϵ+y_i-f\left(x_i\right)]-\underset{i=1}{\overset{n}{\Sigma }}{\gamma }_i({\xi }_i{\gamma }_i+{\xi }_i^{*}{\gamma }_i^{*})\end{array}\right)---\left(28\right)$

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

$w=\underset{i=0}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})x_i---\left(29\right)$

$\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})=0---\left(30\right)$

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

$C-{\alpha }_i^{*}-{\gamma }_i^{*}=0,i=1,...,l---\left(32\right)$

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

$\max -\frac{1}{2}\underset{i,j=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})({\alpha }_j-{\alpha }_j^{*})<x_i,x_j>+\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})y_i-\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i+{\alpha }_i^{*})ϵ---\left(33\right)$

约束为： 

$\underset{i=1}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})=\mathrm{0,0}\le {\alpha }_i,{\alpha }_i^{*}\le C,i=1,...,l---\left(34\right)$

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

$\left(\begin{array}{c}W(\alpha ,{\alpha }^{*})=-\frac{1}{2}\underset{i,j}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})({\alpha }_j-{\alpha }_j^{*})K(x_i,x_j)\\ +\underset{i,j}{\overset{n}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})y_i-\underset{i,j}{\overset{n}{\Sigma }}({\alpha }_i+{\alpha }_i^{*})ϵ\end{array}\right)---\left(35\right)$

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

$f\left(x\right)=\underset{i=1}{\overset{N}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(x_i,x)+b---\left(36\right)$

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

$\left(\begin{array}{cc}{\alpha }_i[{\xi }_i+ϵ-y_i+f\left(x_i\right)]=0& ,i=1,...l\\ a_i^{*}[{\xi }_i^{*}+ϵ-y_i+f\left(x_i\right)]=0& ,i=1,...l\\ (C-{\alpha }_i){\xi }_i=0& ,i=1,...l\\ (C-{\alpha }_i^{*}){\xi }_i^{*}=0& ,i=1,...l\end{array}\right)---\left(37\right)$

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

$b=y_i-ϵ-\underset{i=1}{\overset{l}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(x_j,x_i),$当αj∈(0,C) 

$b=y_i+ϵ-\underset{i=1}{\overset{l}{\Sigma }}({\alpha }_i-{\alpha }_i^{*})K(x_j,x_i),$当∈(0,C) 

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

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

所述步骤1.2：对齐训练形状向量的具体步骤如下： 

所述的将分块的数据映射到特征空间的具体方法如下： 

步骤1.2.1：旋转、缩放和平移每个肺区域形状，使其与训练集中的第一个形状对齐； 

步骤1.2.2：根据对齐形状，计算平均形状； 

步骤1.2.3：旋转、缩放和平移平均形状使其与第一个形状对齐； 

步骤1.2.权：重新将每个形状与当前平均形状对齐； 

步骤1.2.5：如果过程收敛或者到指定循环次数，退出；否则转到步骤。 

所述步骤3中对胸部x光图像做三尺度小波分解。