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

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

【技术实现步骤摘要】
一种胸部X光图像中骨骼抑制的方法
本专利技术属于X线图像处理
，具体涉及一种胸部X光图像中骨骼抑制的方法。
技术介绍
肺癌是当今对人类健康危害最大的恶性肿瘤之一。特别是近半个世纪以来，随着空气污染造成的环境不断恶化及吸烟人口数量的大量增加，各国肺癌的发病率和病死率都在急剧上升。由于肺属于人体内部脏器，多数肺癌在开始的时候只是在身体内悄悄生长，患者没有任何感觉。当患者因咳嗽、咯血及胸痛等临床症状就诊时，大多数已经处于中晚期，错过了治疗的最佳时机。所以，肺癌的早期诊断与治疗是提高肺癌患者生存率的关键。医学影像检查不仅能通过图像发现病灶的存在，定位病灶在全肺的位置，还能观察病灶的大小、形态及密度等物理特征。胸部疾病的影像诊断是以胸部x线摄影为主体的，但在胸部x线图像中，经常会由于病症位于图像中的肋骨或者锁骨结构遮挡区域而造成漏诊。双能量减影(DualEnergySubtraction,DES)是在数字胸部x线摄影基础上发展的一种较新的成像技术。它利用骨与软组织对x线光子的能量衰减方式不同，以及不同原子量的物质的光电吸收效应差别，利用数字摄影将两种吸收效应的信息进行分离，选择性去除骨或软组织的衰减信息，进而获得纯粹的软组织像和骨组织像。因此应用双能量剪影技术可以得到三张图像，一张正常胸片，一张软组织图像和一张骨骼图像。DES的软组织减影图像对胸部小结节性病变，钙化结节性病变，肺纹理病变，特别是位于肺肋骨、锁骨重叠处结节性病变的清晰显示率明显高于普通胸部x光图像。这主要是减影区的肺部软组织图像除去了骨骼等重叠因素的影响，对于肋骨，锁骨相重叠的病灶显示清晰，这对于肋软骨钙化较多的患者更为有益。国内外的一些研究机构对胸部x光图像中骨骼抑制问题进行研究。一些国外学者利用正常的胸部x光图像来估计骨骼结构模型，产生骨骼抑制图像来提高病变的检测性能。Ahmed及Rasheed等人提出利用独立成分分析(IndependentComponentAnalysis,ICA)方法来抑制后侧肋骨，进而提高结节的可见度。Jiann等人提出一种非参数的正常胸片中肋骨抑制算法。该算法由肺分割、肋骨分割、肋骨灰度估计及抑制等部分组成，提出利用实数编码遗传算法(Real-codedgeneticalgorithm,RCGA)来估计肋骨的灰度，进而实现肋骨的抑制。还有一些学者利用DES图像中的正常胸部图像及骨骼图像，构造回归模型，进而达到生成骨骼抑制图像的目的。Loog等人提出基于非线性回归的图像滤波算法来从x光图像中获得骨骼结构及软组织图像。该算法通使用k最近邻回归算法建立DES肺部图像与骨骼图像的回归模型。KenjiSuzuki等人提出基于多分辨率大规模训练人工神经元网络(Multiresolutionmassivetrainingartificialneuralnetwork,MTANN)的图像处理技术来抑制x光图像中骨骼区域的对比度。该算法采用双能图像中的x光图像与骨结构图像作为训练数据，利用MTANN建立回归模型。在生成的骨骼抑制图像中，结节及血管等解剖结构清晰可见。综上所述，国内外对于胸部x光图像中骨骼抑制的研究才刚刚兴起。原有的大多数研究是基于普通胸部x光图像中的骨骼提取、分割、模型建立及移去等。但是，由于普通x光图像中骨骼结构复杂，且边缘区域灰度的对比度较低等特点，仅利用普通x光图像实现骨骼的分割与抑制比较困难。而利用双能图像的骨骼回归模型可以从一个完全新颖的角度来获得骨骼图像，不要再对分割算法进行研究。由于在模型建立过程中，使用大量DES骨骼图像信息，因此获得的预测图像精度较高。
技术实现思路
针对现有技术存在的不足，本专利技术提供一种稳定度高的胸部X光图像中骨骼抑制的方法。本专利技术采取按如下步骤进行：步骤1：模型肺部初始轮廓位置的确定；步骤1.1：标记训练集中每张图像肺边界的边界点；步骤1.2：对齐训练形状向量；通过缩放、旋转和平移的操作使训练形状对齐，使它们尽可能对齐紧密，设xi是训练集中第i个形状中n个点的向量，xi＝(xi0,yi0,…xik,yik,…,xin-1,yin-1)T。其中，(xik,yik)是第i个形状中第k个点，给定两个相似形状xi和xj，选择旋转角度θ、缩放s、平移(tx,ty)，则M(s,θ)[x]代表旋转角度为θ和缩放比例为s的变换，使以下加权和最小化将xi映射为xj：Ej＝(xi-M(s,θ)[xj]-t)TW(xi-M(s,θ)[xj]-t)(1)其中，t＝(tx,ty,…,tx,ty)T。式中：W是一个点的加权对角矩阵，它由每一边界点的权值wk构成；令Rkl表示第k个点和第l个点的距离，表示训练集中距离的变化，则第k个点的权值如果令ax＝scosθ,ay＝ssinθ，则根据最小二乘法有四个线性等式其中，可以使用矩阵运算方法解出变量ax,ay,tx,ty。步骤1.3：建立肺部初始轮廓位置的模型；训练形状向量对齐后，利用主成分分析方法找出形状变化的统计信息，据此建立模型；步骤2：结合灰度信息和形状信息的肺实质分割：同时利用多张特征图像中边界点的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似；步骤2.1：提取特征图像；进行主成分分析对数据进行降维，对于图像I中点p附近的点p+Δp的灰度可由Taylor展开获得，由此可得，图像I的二阶Taylor展开为：对图像进行高斯平滑处理，对于图像I(p)进行高斯滤波后的图像为A(p,σ)＝(I*G)(p,σ)；滤波图像的n阶偏导数为则根据卷积性质可知：式中：角标i1,…in表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；对于一个给定尺度参数σ的n-jet特征可以定义为：...

【技术保护点】
一种胸部X光图像中骨骼抑制的方法，具体步骤如下：步骤1：模型肺部初始轮廓位置的确定；步骤1.1：标记训练集中每张图像肺边界的边界点；步骤1.2：对齐训练形状向量；通过缩放、旋转和平移的操作使训练形状对齐，使它们尽可能对齐紧密，设xi是训练集中第i个形状中n个点的向量，xi=(xi0,yi0,...xik,yik,...,xin?1,yin?1)T，其中，(xik,yik)是第i个形状中第k个点，给定两个相似形状xi和xj，选择旋转角度θ、缩放s、平移(tx,ty)，则M(s,θ)[x]代表旋转角度为θ和缩放比例为s的变换，使以下加权和最小化将xi映射为xj：Ej=(xi?M(s,θ)[xj]?t)TW(xi?M(s,θ)[xj]?t)???????????????(1)其中，t=(tx,ty,...,tx,ty)T，式中:W是一个点的加权对角矩阵，它由每一边界点的权值wk构成。令Rkl表示第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)]]>式中：角标i1,...in表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；利用局部2?jet特征图像来获取边界点的候选点及各个候选点的灰度代价，局部2?jet特征为：J2[I](p,σ)={I,Ix,Iy,Ixx,Ixy,Iyy}????????????(7)步骤2.2：选取边界点的候选点：对于初始肺边界的每一个点，计算所有特征图像中该点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度，选出相似程度最大的点，作为该边界点的候选点。相似程度为所有特征图像中该点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度...

【技术特征摘要】
1.一种胸部X光图像中骨骼抑制的方法，具体步骤如下：步骤1：模型肺部初始轮廓位置的确定；步骤1.1：标记训练集中每张图像肺边界的边界点；步骤1.2：对齐训练形状向量；通过缩放、旋转和平移的操作使训练形状对齐，使它们尽可能对齐紧密，设xi是训练集中第i个形状中n个点的向量，xi＝(xi0,yi0,...xik,yik,...,xin-1,yin-1)T，其中，(xik,yik)是第i个形状中第k个点，给定两个相似形状xi和xj，选择旋转角度θ、缩放s、平移(tx,ty)，则M(s,θ)[x]代表旋转角度为θ和缩放比例为s的变换，使以下加权和最小化将xi映射为xj：Ej＝(xi-M(s,θ)[xj]-t)TW(xi-M(s,θ)[xj]-t)(1)其中，t＝(tx,ty,...,tx,ty)T，式中：W是一个点的加权对角矩阵，它由每一边界点的权值wk构成；令Rkl表示第k个点和第l个点的距离，表示训练集中距离的变化，则第k个点的权值步骤1.3：建立肺部初始轮廓位置的模型；训练形状向量对齐后，利用主成分分析方法找出形状变化的统计信息，据此建立模型；其特征在于：还包括如下步骤：步骤2：结合灰度信息和形状信息的肺实质分割：同时利用多张特征图像中边界点的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似；步骤2.1：提取特征图像；进行主成分分析对数据进行降维，对于图像I中点p附近的点p+△p的灰度可由Taylor展开获得，由此可得，图像I的二阶Taylor展开为：对图像进行高斯平滑处理，对于图像I(p)进行高斯滤波后的图像为A(p,σ)＝(I*G)(p,σ)，滤波图像的n阶偏导数为则根据卷积性质可知：式中：角标i1,...in表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；利用局部2-jet特征图像来获取边界点的候选点及各个候选点的灰度代价，局部2-jet特征为：J2[I](p,σ)＝{I,Ix,Iy,Ixx,Ixy,Iyy}(7)其中：I表示原始图像，就是输入图像；Ix表示对原始图像求x方向的导数；Iy表示对原始图像求y方向的导数；Ixx表示对原始图像求x方向二次导数；Ixy表示对原始图像先求x方向的导数再求y方向的导数；Iyy表示对原始图像求y方向二次导数；步骤2.2：选取边界点的候选点：对于初始肺边界的每一个点，计算所有特征图像中初始肺边界的每一个点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度，选出相似程度最大的点，作为该边界点的候选点；相似程度为所有特征图像中初始肺边界的每一个点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度集合的马氏距离hi；步骤2.3：使用动态规划进行肺区域分割；在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度hi；图像中第i个边界点的形状相似性代价定义为：式中，vi＝pi+1-pi，表示第i个边界点的形状特征，分别表示所有训练图像中第i个边界点形状特征的均值与协方差；对于测试图像中第i个边界点pi，在指定的搜索区域内存在m个具有较小灰度相似性代价候选点，则n个边界点就会产生一个n×m的灰度代价矩阵：式中：为在特征图像上，以点pi为圆心，rc为半径的圆上的nc个点的灰度，k＝1......nc；分别为训练图像的第j个特征图像中第i个边界点的周围像素点灰度的均值及协方差，N为特征图像总数；搜索最优轮廓即找一条最佳路径，在矩阵C中每行选一个元素，沿着路径选择的时候，灰度与形状相似性代价的总和最小，即：其中，γ为灰度与形状相似性代价系数；调整γ值，使得这两种代价在边界搜索过程中发挥大致相同的作用；步骤3：使用B样条多尺度小波变换的特征图像提取；在模型建立过程中，提取有效描述不同尺度骨骼的特征图像是模型建立的基础；步骤3.1：B样条多尺度小波变换：B样条函数随着样条阶数的增加而快速收敛于高斯函数，其一阶导数可以逼近最优边缘检测算子；利用B样条小波进行多尺度边缘增强；利用函数：求取B样条小波变换滤波器系数；式中：m表示B样条小波的阶数；hk为低通高斯滤波器系数，gk为高通高斯滤波器系数，利用所得hk和gk对图像进行快速小波分解，原始图像的行分别与hk和gk进行卷积，再对其列进行下采样，可以得到两幅子图像，然后沿着列的方向对两幅子图像做滤波并进行下采样，就可以得到四幅1/4大小的输出子图像，对角线细节图像，垂直细节图像，水平细节图像和近似图像；步骤3.2：高斯滤波局...

【专利技术属性】
技术研发人员：郭薇，张国栋，丛林，王柳，
申请(专利权)人：沈阳航空航天大学，
类型：发明
国别省市：