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

摘要

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

主权项

一种胸部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的二阶Taylor展开为：$\begin{array}{c}I(p+\mathrm{\&Delta;}p)=\underset{n=0}{\overset{2}{\&Sigma;}}\frac{1}{n!}({\mathrm{\&Delta;p}}_{i1},\mathrm{...},{\mathrm{\&Delta;p}}_{in}\frac{{\&part;}^{n}I\left(p\right)}{\&part;i_1\mathrm{...}\&part;i_n})=I\left(p\right)+{\mathrm{\&Delta;p}}_{i1}\frac{\&part;I\left(p\right)}{\&part;i_1}+\frac{1}{2!}\left({\mathrm{\&Delta;p}}_{i1}{\mathrm{\&Delta;p}}_{i2}\frac{{\&part;}^{2}I\left(p\right)}{\&part;i_1\&part;i_2}\right)\\ =I\left(p\right)+{\mathrm{\&Delta;p}}_x\frac{\&part;I\left(p\right)}{\&part;x}+\frac{1}{2!}({\mathrm{\&Delta;p}}_x^2\frac{{\&part;}^{2}I\left(p\right)}{\&part;x^2}+{\mathrm{\&Delta;p}}_x{\mathrm{\&Delta;p}}_y\frac{{\&part;}^{2}I\left(p\right)}{\&part;x\&part;y})\\ +{\mathrm{\&Delta;p}}_y\frac{\&part;I\left(p\right)}{\&part;y}+\frac{1}{2!}({\mathrm{\&Delta;p}}_y^2\frac{{\&part;}^{2}I\left(p\right)}{\&part;y^2}+{\mathrm{\&Delta;p}}_y{\mathrm{\&Delta;p}}_x\frac{{\&part;}^{2}I\left(p\right)}{\&part;y\&part;x})\end{array}---\left(3\right)$对图像进行高斯平滑处理，对于图像I(p)进行高斯滤波后的图像为A(p,σ)＝(I*G)(p,σ)，滤波图像的n阶偏导数为则根据卷积性质可知：$A_{i_1,...i_n}(p,\mathrm{\&sigma;})={(I*G)}_{i_1,...i_n}(p,\mathrm{\&sigma;})=(I*G_{i_1,...i_n})(p,\mathrm{\&sigma;})---\left(5\right)$式中：角标i₁,...i_n表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；利用局部2‑jet特征图像来获取边界点的候选点及各个候选点的灰度代价，局部2‑jet特征为：J²[I](p,σ)＝{I,I_x,I_y,I_xx,I_xy,I_yy}              (7)其中：I表示原始图像，就是输入图像；Ix表示对原始图像求x方向的导数；Iy表示对原始图像求y方向的导数；Ixx表示对原始图像求x方向二次导数；Ixy表示对原始图像先求x方向的导数再求y方向的导数；Iyy表示对原始图像求y方向二次导数；步骤2.2：选取边界点的候选点：对于初始肺边界的每一个点，计算所有特征图像中初始肺边界的每一个点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度，选出相似程度最大的点，作为该边界点的候选点；相似程度为所有特征图像中初始肺边界的每一个点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度集合的马氏距离h_i；步骤2.3：使用动态规划进行肺区域分割；在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度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}_{1,1}& \mathrm{...}& h_{1,k}& \mathrm{...}& h_{1,m}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ h_{i,1}& \mathrm{...}& h_{i,k}& \mathrm{...}& h_{i,m}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ h_{n,1}& \mathrm{...}& h_{n,k}& \mathrm{...}& 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\mathrm{\&pi;}}{n_c}\left(k-1\right)\right)\\ sin\left(\frac{2\mathrm{\&pi;}}{n_c}\left(k-1\right)\right)\end{array}\right]---\left(11\right)$式中：为在特征图像上，以点p_i为圆心，r_c为半径的圆上的n_c个点的灰度，k＝1......n_c；分别为训练图像的第j个特征图像中第i个边界点的周围像素点灰度的均值及协方差，N为特征图像总数；搜索最优轮廓即找一条最佳路径，在矩阵C中每行选一个元素，沿着路径选择的时候，灰度与形状相似性代价的总和最小，即：$J\left(k_1\mathrm{.....}{k}_n\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{h}_i+\mathrm{\&gamma;}\underset{i=1}{\overset{n}{\&Sigma;}}{f}_i---\left(12\right)$其中，γ为灰度与形状相似性代价系数；调整γ值，使得这两种代价在边界搜索过程中发挥大致相同的作用；步骤3：使用B样条多尺度小波变换的特征图像提取；在模型建立过程中，提取有效描述不同尺度骨骼的特征图像是模型建立的基础；步骤3.1：B样条多尺度小波变换：B样条函数随着样条阶数的增加而快速收敛于高斯函数，其一阶导数可以逼近最优边缘检测算子；利用B样条小波进行多尺度边缘增强；利用函数：$H\left(e^{i\mathrm{\&omega;}}\right)=\&Sigma;h_k{e}^{-ik\mathrm{\&omega;}}=\frac{{\widehat{\mathrm{\&theta;}}}_m\left(2\mathrm{\&omega;}\right)}{{\widehat{\mathrm{\&theta;}}}_m\left(\mathrm{\&omega;}\right)}=\frac{e^{-i\mathrm{\&epsiv;}\mathrm{\&omega;}}{\left(\frac{\sin \mathrm{\&omega;}}{\mathrm{\&omega;}}\right)}^{m+1}}{e^{-i\frac{\mathrm{\&epsiv;}\mathrm{\&omega;}}{2}}{\left(\frac{\sin (\mathrm{\&omega;}/2)}{\mathrm{\&omega;}/2}\right)}^{m+1}}=e^{-i\frac{\mathrm{\&epsiv;}\mathrm{\&omega;}}{2}}{\left(\cos \frac{\mathrm{\&omega;}}{2}\right)}^{m+1}---\left(18\right)$$G\left(\mathrm{\&omega;}\right)=\frac{{\widehat{\mathrm{\&psi;}}}_m\left(2\mathrm{\&omega;}\right)}{{\widehat{\mathrm{\&theta;}}}_m\left(\mathrm{\&omega;}\right)}=\underset{k}{\&Sigma;}{g}_k{e}^{-ik\mathrm{\&omega;}}---\left(20\right)$求取B样条小波变换滤波器系数；式中：m表示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ⁿ，Rⁿ表示n维的实数空间；样本集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{|\&lt;w,x\&gt;+b-y|}{\sqrt{1+\left|\right|w|{|}^2}}---\left(23\right)$因为S集是ε‑线性近似，所以有|<w,x>+b‑f(x_i)|≤ε,i＝1,...,m；可以得到:$\frac{|\&lt;w,x\&gt;+b-y|}{\sqrt{1+\left|\right|w|{|}^2}}\&le;\frac{\mathrm{\&epsiv;}}{\sqrt{1+\left|\right|w|{|}^2}},i=1,...,l---\left(24\right)$于是有$d_i\&le;\frac{\mathrm{\&epsiv;}}{\sqrt{1+\left|\right|w|{|}^2}},i=1,...,l---\left(25\right)$(25)式表明是S中的点到超平面的距离的上界；ε‑不敏感损失函数的线性近似集S的最优近似超平面是通过最大化S中的点到超平面距离的上界而得到超平面；最优超平面是通过最大化得到的；因此，只要最小化||w||²就可以得到最优近似超平面；于是线性回归问题就化为：求下面的最优化问题：$min\frac{1}{2}\left|\right|w|{|}^2---\left(26\right)$约束为|＜w，x＞+b‑yi|≤ε，i＝1，...，m；另外，考虑拟合误差的情况，引入松弛因子ξ_i，当划分有误差时，ξ_i，都大于0，误差不存在取0；因而，标准ε不敏感支持向量回归机可以表示为：$\left\{\begin{array}{c}\underset{w,b,\mathrm{\&xi;}}{min}\frac{1}{2}\left|\right|w|{|}^2+C\underset{i=1}{\overset{l}{\&Sigma;}}({\mathrm{\&xi;}}_i+{\mathrm{\&xi;}}_i^{*})\\ \begin{array}{cc}s.t.& y_i-w\&CenterDot;x_i-b\&le;\mathrm{\&epsiv;}+{\mathrm{\&xi;}}_i\\ & w\&CenterDot;x_i+b-y_i\&le;\mathrm{\&epsiv;}+{\mathrm{\&xi;}}_i^{*}\\ & {\mathrm{\&xi;}}_i\&GreaterEqual;0\\ & {\mathrm{\&xi;}}_i^{*}\&GreaterEqual;0,i=1,2,...,l\end{array}\end{array}\right.---\left(27\right)$其中，C＞0为平衡因子；引入拉格朗日函数为：$\begin{array}{c}L(w,b,\mathrm{\&xi;},\mathrm{\&alpha;},{\mathrm{\&alpha;}}^{*},\mathrm{\&gamma;})=\frac{1}{2}\left|\right|w|{|}^2+C\underset{i=1}{\overset{n}{\&Sigma;}}({\mathrm{\&xi;}}_i+{\mathrm{\&xi;}}_i^{*})-\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&alpha;}}_i\&lsqb;{\mathrm{\&xi;}}_i+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)\&rsqb;\\ -\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&alpha;}}_i^{*}\&lsqb;{\mathrm{\&xi;}}_i+\mathrm{\&epsiv;}+y_i-f\left(x_i\right)\&rsqb;-\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&gamma;}}_i({\mathrm{\&xi;}}_i{\mathrm{\&gamma;}}_i+{\mathrm{\&xi;}}_i^{*}{\mathrm{\&gamma;}}_i^{*})\end{array}---\left(28\right)$其中，α_i，γi≥0，i＝1，...，n；函数L的极值应满足条件于是得到下面的式子：$w=\underset{i=1}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})x_i---\left(29\right)$$\underset{i=1}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})=0---\left(30\right)$C‑α_i‑γ_i＝0，i＝1，...，l                         (31)$C-{\mathrm{\&alpha;}}_i^{*}-{\mathrm{\&gamma;}}_i^{*}=0,i=1,...,l---\left(32\right)$将(29)‑(32)代入到(28)中，得到优化问题的对偶形式为：$\max \&lsqb;\frac{1}{2}\underset{i,j=1}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})({\mathrm{\&alpha;}}_j-{\mathrm{\&alpha;}}_j^{*})\&lt;x_i,x_j\&gt;+\underset{i=1}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})y_i-\underset{i=1}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i+{\mathrm{\&alpha;}}_i^{*})\mathrm{\&epsiv;}\&rsqb;---\left(33\right)$约束为：$\underset{i=1}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})=0,0\&le;{\mathrm{\&alpha;}}_i,{\mathrm{\&alpha;}}_i^{*}\&le;C,i=1,...,l---\left(34\right)$由于骨骼灰度预测为非线性回归，所以先使用一个非线性映射φ，把数据映射到一个高维特征空间，然后在高维特征空间进行线性回归建立模型；由于在上面的优化过程中，只考虑到高维特征空间中的内积运算，因此用一个核函数k(x,y)代替<φ(x),φ(y)>就可以实现非线性回归；于是，非线性回归的优化方程为最大化下面的函数：$\begin{array}{c}W(\mathrm{\&alpha;},{\mathrm{\&alpha;}}^{*})=-\frac{1}{2}\underset{i,j}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})({\mathrm{\&alpha;}}_j-{\mathrm{\&alpha;}}_j^{*})K(x_i,x_j)\\ +\underset{i,j}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})y_i-\underset{i,j}{\overset{n}{\&Sigma;}}({\mathrm{\&alpha;}}_i+{\mathrm{\&alpha;}}_i^{*})\mathrm{\&epsiv;}\end{array}---\left(35\right)$其约束条件为(32)，求解出α的值后，就可以得到f(x)：$f\left(x\right)=\underset{i=1}{\overset{N}{\&Sigma;}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x_i,x)+b---\left(36\right)$通常情况下，大部分α_i或的值将为零，不为零的α_i或所对应的样本称为支持向量，根据最优化条件，在鞍点有：$\left\{\begin{array}{cc}{\mathrm{\&alpha;}}_i\&lsqb;{\mathrm{\&xi;}}_i+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)\&rsqb;=0,& i=1,...l\\ {\mathrm{\&alpha;}}_i^{*}\&lsqb;{\mathrm{\&xi;}}_i^{*}+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)\&rsqb;=0,& i=1,...l\\ (C-{\mathrm{\&alpha;}}_i){\mathrm{\&xi;}}_i=0,& i=1,...l\\ (C-{\mathrm{\&alpha;}}_i^{*}){\mathrm{\&xi;}}_i^{*}=0,& i=1,...l\end{array}\right.---\left(37\right)$于是得到b的计算式如下：当α_j∈(0,C)当用任意一个支持向量可以计算出b的值，得到回归函数f(x)；建立预测模型后，可以根据肺部X光图像的灰度值分布来预测产生骨结构图像；步骤5：对正常胸片与预测的骨骼图像做图像减法来预测软组织图像。