CN103824281A

CN103824281A - Bone inhibition method in chest X-ray image

Info

Publication number: CN103824281A
Application number: CN201410007747.XA
Authority: CN
Inventors: 郭薇; 张国栋; 丛林; 王柳
Original assignee: Shenyang Aerospace University
Current assignee: Shenyang Aerospace University
Priority date: 2014-01-07
Filing date: 2014-01-07
Publication date: 2014-05-28
Anticipated expiration: 2034-01-07

Abstract

The invention belongs to the technical field of X-ray image processing, and in particular relates to a bone suppression method in chest X-ray images. The present invention establishes a lung contour model for the normal chest radiograph in the dual-energy silhouette image, and then uses the third-order B-spline wavelet transform to decompose the image by three-scale wavelet, uses the 2-jet method to extract the characteristic image from the decomposed image, and then A support vector regression machine is used to establish a bone model, and the bone model is applied to the X-ray image to extract the corresponding predicted bone image. Finally, the X-ray image is subtracted from the predicted bone image to obtain the predicted soft tissue image. Because most of the bone suppression methods use the neural network method, and this method is easy to fall into a local optimum, and the training results are not stable, so this patent uses a support vector regression machine for bone prediction, and better results are obtained.

Description

A Method for Skeletal Suppression in Chest X-ray Images

技术领域 technical field

本发明属于X线图像处理技术领域，具体涉及一种胸部X光图像中骨骼抑制的方法。 The invention belongs to the technical field of X-ray image processing, and in particular relates to a bone suppression method in chest X-ray images. the

背景技术 Background technique

肺癌是当今对人类健康危害最大的恶性肿瘤之一。特别是近半个世纪以来，随着空气污染造成的环境不断恶化及吸烟人口数量的大量增加，各国肺癌的发病率和病死率都在急剧上升。由于肺属于人体内部脏器，多数肺癌在开始的时候只是在身体内悄悄生长，患者没有任何感觉。当患者因咳嗽、咯血及胸痛等临床症状就诊时，大多数已经处于中晚期，错过了治疗的最佳时机。所以，肺癌的早期诊断与治疗是提高肺癌患者生存率的关键。 Lung cancer is one of the most harmful malignant tumors to human health today. Especially in the past half century, with the continuous deterioration of the environment caused by air pollution and the large increase in the number of smokers, the incidence and mortality of lung cancer in various countries have been rising sharply. Since the lung is an internal organ of the human body, most lung cancers just grow quietly in the body at the beginning, and the patient does not feel anything. When patients go to the doctor due to clinical symptoms such as cough, hemoptysis, and chest pain, most of them are already in the middle and late stage, and the best time for treatment has been missed. Therefore, early diagnosis and treatment of lung cancer is the key to improve the survival rate of lung cancer patients. the

医学影像检查不仅能通过图像发现病灶的存在，定位病灶在全肺的位置，还能观察病灶的大小、形态及密度等物理特征。胸部疾病的影像诊断是以胸部x线摄影为主体的，但在胸部x线图像中，经常会由于病症位于图像中的肋骨或者锁骨结构遮挡区域而造成漏诊。 Medical imaging examination can not only detect the existence of lesions through images, locate the location of lesions in the whole lung, but also observe physical characteristics such as size, shape and density of lesions. The imaging diagnosis of chest diseases is mainly based on chest X-ray photography, but in chest X-ray images, the diagnosis is often missed because the disease is located in the area covered by the rib or clavicle structure in the image. the

双能量减影(Dual Energy Subtraction,DES)是在数字胸部x线摄影基础上发展的一种较新的成像技术。它利用骨与软组织对x线光子的能量衰减方式不同，以及不同原子量的物质的光电吸收效应差别，利用数字摄影将两种吸收效应的信息进行分离，选择性去除骨或软组织的衰减信息，进而获得纯粹的软组织像和骨组织像。因此应用双能量剪影技术可以得到三张图像，一张正常胸片，一张软组织图像和一张骨骼图像。 Dual Energy Subtraction (DES) is a relatively new imaging technology developed on the basis of digital chest X-ray photography. It takes advantage of the different energy attenuation methods of x-ray photons between bone and soft tissue, and the difference in the photoelectric absorption effect of substances with different atomic weights. It uses digital photography to separate the information of the two absorption effects, and selectively removes the attenuation information of bone or soft tissue. Obtain pure soft tissue images and bone tissue images. Therefore, three images can be obtained by applying dual-energy silhouette technology, one normal chest radiograph, one soft tissue image and one bone image. the

DES的软组织减影图像对胸部小结节性病变，钙化结节性病变，肺纹理病变，特别是位于肺肋骨、锁骨重叠处结节性病变的清晰显示率明显高于普通胸部x光图像。这主要是减影区的肺部软组织图像除去了骨骼等重叠因素的影响，对于肋骨，锁骨相重叠的病灶显示清晰，这对于肋软骨钙化较多的患者更为有益。 The soft tissue subtraction images of DES can clearly display small nodular lesions, calcified nodular lesions, and pulmonary marking lesions, especially the nodular lesions at the overlap of lung ribs and clavicles, which are significantly higher than ordinary chest X-ray images. This is mainly because the lung soft tissue image in the subtraction area removes the influence of overlapping factors such as bones, and the overlapping lesions of ribs and clavicles are clearly displayed, which is more beneficial for patients with more costal cartilage calcification. the

国内外的一些研究机构对胸部x光图像中骨骼抑制问题进行研究。一些国外学者利用正常的胸部x光图像来估计骨骼结构模型，产生骨骼抑制图像来图高病变的检测性能。Ahmed及Rasheed等人提出利用独立成分分析(Independent Component Analysis,ICA)方法来抑制后侧肋骨，进而提高结节的可见度。Jiann等人提出一种非参数的正常胸片中肋骨抑制算法。该算法由肺分割、肋骨分割、肋骨灰度估计及抑制等部分组成，提出利用实数编码遗传算法(Real-coded genetic algorithm,RCGA)来估计肋骨的灰度，进而实现肋骨的抑制。 Some research institutions at home and abroad have studied the problem of skeletal suppression in chest X-ray images. Some foreign scholars use normal chest X-ray images to estimate the bone structure model, and generate bone-suppressed images to improve the detection performance of lesions. Ahmed and Rasheed et al proposed to use the independent component analysis (Independent Component Analysis, ICA) method to suppress the posterior ribs, thereby improving the visibility of nodules. Jiann et al proposed a non-parametric rib suppression algorithm in normal chest radiographs. The algorithm is composed of lung segmentation, rib segmentation, rib gray level estimation and suppression, etc. It is proposed to use real-coded genetic algorithm (RCGA) to estimate the rib gray level, and then realize the rib suppression. the

还有一些学者利用DES图像中的正常胸部图像及骨骼图像，构造回归模型，进而达到生成骨骼抑制图像的目的。Loog等人提出基于非线性回归的图像滤波算法来从x光图像中获得骨骼结构及软组织图像。该算法通使用k最近邻回归算法建立DES肺部图像与骨骼图像的回归模型。Kenji Suzuki等人提出基于多分辨率大规模训练人工神经元网络(Multiresolution massive training artificial neural network,MTANN)的图像处理技术来抑制x光图像中骨骼区域的对比度。该算法采用双能图像中的x光图像与骨结构图像作为训练数据，利用MTANN建立回归模型。在生成的骨骼抑制图像中，结节及血管等解剖结构清晰可见。 Some scholars use normal chest images and bone images in DES images to construct regression models, and then achieve the purpose of generating bone suppression images. Loog et al proposed an image filtering algorithm based on nonlinear regression to obtain bone structure and soft tissue images from X-ray images. The algorithm uses the k-nearest neighbor regression algorithm to establish a regression model of DES lung images and bone images. Kenji Suzuki et al. proposed an image processing technology based on multiresolution massive training artificial neural network (MTANN) to suppress the contrast of bone regions in X-ray images. The algorithm uses X-ray images and bone structure images in dual-energy images as training data, and uses MTANN to establish a regression model. In the generated bone-suppressed images, anatomical structures such as nodules and blood vessels are clearly visible. the

综上所述，国内外对于胸部x光图像中骨骼抑制的研究才刚刚兴起。原有的大多数研究是基于普通胸部x光图像中的骨骼提取、分割、模型建立及移去等。但是，由于普通x光图像中骨骼结构复杂，且边缘区域灰度的对比度较低等特点，仅利用普通x光图像实现骨骼的分割与抑制比较困难。而利用双能图像的骨骼回归模型可以从一个完全下颖的角度来获得骨骼图像，不要再对分割算法进行研究。由于在模型建立过程中，使用大量DES骨骼图像信息，因此获得的预测图像精度较高。 To sum up, the research on skeletal suppression in chest X-ray images at home and abroad is just emerging. Most of the original research is based on bone extraction, segmentation, model building and removal in ordinary chest X-ray images. However, due to the complex structure of bones in ordinary X-ray images and the low gray contrast of edge regions, it is difficult to achieve bone segmentation and suppression using only ordinary X-ray images. However, the skeletal regression model using dual-energy images can obtain skeletal images from a completely novel perspective, so there is no need to study the segmentation algorithm. Since a large amount of DES bone image information is used in the process of model building, the accuracy of the predicted image obtained is high. the

发明内容 Contents of the invention

针对现有技术存在的不足，本发明提供一种稳定度高的胸部X光图像中骨骼抑制的方法。 Aiming at the deficiencies in the prior art, the present invention provides a method for bone suppression in chest X-ray images with high stability. the

本发明采取按如下步骤进行： The present invention takes to carry out as follows:

步骤1：模型肺部初始轮廓位置的确定； Step 1: Determination of the initial contour position of the model lung;

步骤1.1：标记训练集中每张图像肺边界的边界点； Step 1.1: Mark the boundary points of the lung boundaries of each image in the training set;

步骤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： Step 1.2: Align the training shape vectors; align the training shapes by scaling, rotating and translating operations so that they are aligned as closely as possible, let x _i be the vector of n points in the i-th shape in the training set, x _i =(x _i0 ,y _i0 ,...x _ik ,y _ik ,...,x _in-1 ,y _in-1 ) ^T . Among them, (xi _ik , y _ik ) is the kth point in the i-th shape, given two similar shapes x _i and x _j , choose the rotation angle θ, scale s, and translate (t _x , t _y ), then M(s,θ)[x] represents a transformation with rotation θ and scaling s that minimizes the following weighted sum mapping x _i to x _j :

E_j=(x_i-M(s,θ)[x_j]-t)^TW(x_i-M(s,θ)[x_j]-t) (1) E _j =( _xi -M(s,θ)[x _j ]-t) ^T W( _xi -M(s,θ)[x _j ]-t) (1)

其中， $M (s, θ) [\begin{matrix} x_{jk} \\ y_{jk} \end{matrix}] = (\begin{matrix} (s \cos θ) x_{jk} - (s \sin θ) y_{jk} \\ (s \sin θ) x_{jk} + (s \cos θ) y_{jk} \end{matrix}),$ t=(t_x,t_y,...,t_x,t_y)^T。式中：W是一个点的加权对角矩阵，它由每一边界点的权值w_k构成。令R_kl表示第k个点和第l个点的距离，

表示训练集中距离的变化，则第k个点的权值

in,

m (the s, θ) [\begin{matrix} x_{jk} \\ {the y}_{jk} \end{matrix}] = (\begin{matrix} (the s \cos θ) x_{jk} - (the s \sin θ) {the y}_{jk} \\ (the s \sin θ) x_{jk} + (the s \cos θ) {the y}_{jk} \end{matrix}),

t=(t _x ,t _y ,...,t _x ,t _y ) ^T . In the formula: W is a weighted diagonal matrix of a point, which is composed of the weight w _k of each boundary point. Let R _kl represent the distance between the kth point and the lth point,

Indicates the change of the distance in the training set, then the weight of the kth point

如果令a_x=scosθ,a_y=ssinθ，则根据最小二乘法有四个线性等式 If a _x =scosθ, a _y =ssinθ, then there are four linear equations according to the least square method

$(\begin{matrix} {X x}_{22} & - - {Y Y}_{22} & W W & 00 \\ {Y Y}_{22} & {X x}_{22} & 00 & W W \\ Z Z & 00 & {X x}_{22} & {Y Y}_{22} \\ 00 & Z Z & {- - Y Y}_{22} & {X x}_{22} \end{matrix}) (\begin{matrix} {a a}_{x x} \\ {a a}_{y the y} \\ {t t}_{x x} \\ {t t}_{y the y} \end{matrix}) = = (\begin{matrix} {X x}_{11} \\ {Y Y}_{11} \\ {C C}_{11} \\ {C C}_{22} \end{matrix}) - - - - - - ((44))$

其中， $X_{i} = Σ_{k = 0}^{n - 1} w_{k} x_{ik},$ $Y_{i} = Σ_{k = 0}^{n - 1} w_{k} y_{ik},$ $Z = Σ_{k = 0}^{n - 1} w_{k} (x_{2 k}^{2} + y_{2 k}^{2}),$ $C_{1} = Σ_{k = 0}^{n - 1} w_{k} (x_{1 k} x_{2 k} + y_{1 k} y_{2 k}),$ $C_{2} = Σ_{k = 0}^{n - 1} w_{k} (y_{1 k} x_{2 k} + x_{1 k} y_{2 k}),$ $w = Σ_{k = 0}^{n - 1} w_{k} .$ 可以使用矩阵运算方法解出变量a_x,a_y,t_x,t_y。 in, $x_{i} = Σ_{k = 0}^{no - 1} w_{k} x_{ik},$ $Y_{i} = Σ_{k = 0}^{no - 1} w_{k} {the y}_{ik},$ $Z = Σ_{k = 0}^{no - 1} w_{k} (x_{2 k}^{2} + {the y}_{2 k}^{2}),$ $C_{1} = Σ_{k = 0}^{no - 1} w_{k} (x_{1 k} x_{2 k} + {the y}_{1 k} {the y}_{2 k}),$ $C_{2} = Σ_{k = 0}^{no - 1} w_{k} ({the y}_{1 k} x_{2 k} + x_{1 k} {the y}_{2 k}),$ $w = Σ_{k = 0}^{no - 1} w_{k} .$ The variables a _x , a _y , t _x , t _y can be solved using the matrix operation method.

步骤1.3：建立肺部初始轮廓位置的模型；训练形状向量对齐后，利用主成分分析方法找出形状变化的统计信息，据此建立模型； Step 1.3: Establish the model of the initial contour position of the lung; after the training shape vectors are aligned, use the principal component analysis method to find out the statistical information of the shape change, and build the model accordingly;

步骤2：结合灰度信息和形状信息的肺实质分割：同时利用多张特征图像中边界点的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似； Step 2: Segmentation of lung parenchyma combining grayscale information and shape information: using the grayscale and shape information of boundary points in multiple feature images at the same time, so that the searched boundary grayscale and shape information are similar to the training images;

步骤2.1：提取特征图像；进行主成分分析对数据进行降维，对于图像I中点p附近的点p+△p的灰度可由Taylor展开获得， $I (p + Δp) \approx Σ_{n = 0}^{N} \frac{1}{n!} (Δ p_{i 1}, \cdot \cdot \cdot, Δ p_{in} \frac{{&PartialD;}^{n} I (p)}{{&PartialD; i}_{1 \cdot \cdot \cdot} {&PartialD; i}_{n}}) - - - (2)$ Step 2.1: Extract the feature image; perform principal component analysis to reduce the dimensionality of the data, and the gray level of the point p+△p near the point p in the image I can be obtained by Taylor expansion, $I (p + Δp) \approx Σ_{no = 0}^{N} \frac{1}{no!} (Δ p_{i 1}, &Center Dot; &Center Dot; &Center Dot;, Δ p_{in} \frac{{&PartialD;}^{no} I (p)}{{&PartialD; i}_{1 \cdot \cdot \cdot} {&PartialD; i}_{no}}) - - - (2)$

由此可得，图像I的二阶Taylor展开为： It can be obtained that the second-order Taylor expansion of image I is:

$\begin{matrix} I I ((p p + + Δp Δp)) = = {Σ Σ}_{n no = = 00}^{22} \frac{11}{n no!!} ((Δ Δ {p p}_{i i 11},, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;,, Δ Δ {p p}_{in in} \frac{{&PartialD; &PartialD;}^{n no} I I ((p p))}{{&PartialD; &PartialD; i i}_{11} \cdot \cdot \cdot \cdot \cdot \cdot {&PartialD; &PartialD; i i}_{n no}})) = = I I ((p p)) + + Δ Δ {p p}_{i i 11} \frac{&PartialD; &PartialD; I I ((p p))}{{&PartialD; &PartialD; i i}_{11}} + + \frac{11}{22!!} ((Δ Δ {p p}_{i i 11} Δ Δ {p p}_{i i 22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; i i}_{11} {&PartialD; &PartialD; i i}_{22}})) \\ = = I I ((p p)) + + Δ Δ {p p}_{x x} \frac{&PartialD; &PartialD; I I ((p p))}{&PartialD; &PartialD; x x} + + \frac{11}{22!!} ((Δ Δ {p p}_{x x}^{22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; x x}^{22}} + + Δ Δ {p p}_{x x} Δ Δ {p p}_{y the y} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{&PartialD; &PartialD; x x &PartialD; &PartialD; y the y})) \\ + + Δ Δ {p p}_{y the y} \frac{&PartialD; &PartialD; I I ((p p))}{&PartialD; &PartialD; y the y} + + \frac{11}{22!!} ((Δ Δ {p p}_{y the y}^{22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; y the y}^{22}} + + Δ Δ {p p}_{y the y} Δ Δ {p p}_{x x} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{&PartialD; &PartialD; y the y &PartialD; &PartialD; x x})) \end{matrix} - - - - - - ((33))$

对图像进行高斯平滑处理，对于图像I(p)进行高斯滤波后的图像为A(p,σ)=(I*G)(p,σ)。滤波图像的n阶偏导数为则根据卷积性质可知： Gaussian smoothing is performed on the image, and the image after Gaussian filtering on the image I(p) is A(p,σ)=(I*G)(p,σ). The nth order partial derivative of the filtered image is According to the convolution properties, we can know:

${A A}_{{i i}_{11},, \cdot \cdot \cdot &Center Dot; \cdot \cdot {i i}_{n no}} ((p p,, σ σ)) = = {((I I * * G G))}_{{i i}_{11},, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; {i i}_{n no}} ((p p,, σ σ)) = = ((I I * * {G G}_{{i i}_{11},, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; {i i}_{n no}})) ((p p,, σ σ)) - - - - - - ((55))$

式中：角标i₁,...i_n表示数据分别对变量求偏导数，把高斯偏导数看作是一个滤波器，高斯偏导数响应的组合构成一个完整的特征描述；对于一个给定尺度参数σ的n-jet特征可以定义为： In the formula: the subscripts i ₁ ,...i _n indicate that the data calculates the partial derivatives of the variables respectively, and the Gaussian partial derivatives are regarded as a filter, and the combination of the Gaussian partial derivative responses constitutes a complete feature description; for a given The n-jet feature of the scale parameter σ can be defined as:

${J J}^{N N} [[I I]] ((p p,, σ σ)) = = {{{I I}_{{i i}_{11},, . . . . . . {i i}_{n no}} | | n no = = 00,, . . . . . .,, N N}},, N N &Element; &Element; Z Z - - - - - - ((66))$

由此可以看出，当N→∞时，在某一尺度空间中，图像I的局部特征由所有偏导数的组合构成。因此，可以采用高斯偏导数构造一个特征向量来描述局部特征。当偏导数阶数越高时，对图像的描述越准确。但是，如果使用过多的偏导数描述图像，会使得特征向量维数增加。 It can be seen from this that when N → ∞, in a certain scale space, the local features of the image I are composed of the combination of all partial derivatives. Therefore, Gaussian partial derivatives can be used to construct a feature vector to describe local features. When the partial derivative order is higher, the description of the image is more accurate. However, if too many partial derivatives are used to describe the image, the dimension of the feature vector will increase. the

利用局部2-jet特征图像来获取边界点的候选点及各个候选点的灰度代价，局部2-jet特征为： Use the local 2-jet feature image to obtain the candidate points of the boundary points and the gray cost of each candidate point. The local 2-jet features are:

J²[I](p,σ)={I,I_x,I_y,I_xx,I_xy,I_yy} (7) J ² [I](p,σ)={I,I _x ,I _y ,I _xx ,I _xy ,I _yy } (7)

步骤2.2：选取边界点的候选点：对于初始肺边界的每一个点，计算所有特征图像中该点搜索区域内所有像素点的灰度与训练特征图像中相应点灰度的相似程度，选出相似程度最大的点，作为该边界点的候选点。相似程度为所有特征图像中该点周围像素点灰度到训练样本特征图像中相应点周围像素点灰度集合的马氏距离h_i： Step 2.2: Select the candidate point of the boundary point: For each point of the initial lung boundary, calculate the similarity between the gray levels of all pixels in the search area of this point in all feature images and the gray level of the corresponding point in the training feature image, and select The point with the greatest similarity is used as the candidate point of the boundary point. The degree of similarity is the Mahalanobis distance h _i from the gray level of the pixels around the point in all feature images to the set of gray levels of pixels around the corresponding point in the feature image of the training sample:

步骤2.3：使用动态规划进行肺区域分割；在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度h_i在边界点搜索区域内，像素点的灰度相似性代价为该点周围像素点灰度与训练图像中相应边界点的周围像素点灰度的相似程度h_i；图像中第i个边界点的形状相似性代价定义为： Step 2.3: Use dynamic programming to segment the lung area; in the boundary point search area, the gray similarity cost of the pixel point is the similarity between the gray level of the surrounding pixels of this point and the surrounding pixel gray level of the corresponding boundary point in the training image h _i In the boundary point search area, the gray similarity cost of the pixel point is the similarity h _i of the gray level of the surrounding pixels of this point and the surrounding pixel gray level of the corresponding boundary point in the training image; the i-th boundary in the image The shape similarity cost of a point is defined as:

${f f}_{i i} = = {(({v v}_{i i} - - {u u}_{{v v}_{i i}}))}^{T T} {(({COV COV}_{{v v}_{i i}}))}^{- - 11} (({v v}_{i i} - - {u u}_{{v v}_{i i}})) - - - - - - ((88))$

式中，v_i=p_i+1-p_i，表示第i个边界点的形状特征，

分别表示所有训练图像中第i个边界点形状特征的均值与协方差。 In the formula, v _i =p _i+1 -p _i represents the shape feature of the i-th boundary point,

Represent the mean and covariance of the i-th boundary point shape features in all training images, respectively.

对于测试图像中第i个边界点p_i，在指定的搜索区域内存在m个具有较小灰度相似性代价候选点，则n个边界点就会产生一个n×m的灰度代价矩阵： For the i-th boundary point p _i in the test image, there are m candidate points with a smaller gray-scale similarity cost in the specified search area, then n boundary points will generate an n×m gray-scale cost matrix:

$C C = = [\begin{matrix} {h h}_{1,1 1,1} & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & {h h}_{11,, k k} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {h h}_{11,, m m} \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ {h h}_{i i,, 11} & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & {h h}_{i i,, k k} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {h h}_{i i,, m m} \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ {h h}_{n no,, 11} & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & {h h}_{n no,, k k} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {h h}_{n no,, m m} \end{matrix}] - - - - - - ((99))$

${h h}_{i i} = = {Σ Σ}_{j j = = 11}^{N N} {(({g g}_{i i}^{j j} - - {u u}_{{g g}_{i i}^{j j}}))}^{T T} {(({s the s}_{{g g}_{i i}^{j j}}))}^{- - 11} (({g g}_{i i}^{j j} - - {u u}_{{g g}_{i i}^{j j}})) - - - - - - ((1010))$

${g g}_{i i}^{j j} = = {p p}_{i i} + + {r r}_{c c} [\begin{matrix} cos cos ((\frac{22 π π}{{n no}_{c c}} ((k k - - 11)))) \\ sin sin ((\frac{22 π π}{{n no}_{c c}} ((k k - - 11)))) \end{matrix}] - - - - - - ((1111))$

式中：

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

is the gray level of n _c points on a circle with point p _i as the center and r _c as the radius on the feature image, k=1...n _c . are the mean and covariance of the gray values of the surrounding pixels of the i-th boundary point in the j-th feature image of the training image, respectively. N is the total number of feature images,

搜索最优轮廓即找一条最佳路径，在矩阵C中每行选一个元素，沿着路径选择的时候，灰度与形状相似性代价的总和最小，即： To search for the optimal contour is to find an optimal path, select one element in each row in the matrix C, and when selecting along the path, the sum of the gray level and shape similarity costs is the smallest, namely:

$J J (({k k}_{11} . . . . . . . . . . {k k}_{n no})) = = {Σ Σ}_{i i = = 11}^{n no} {h h}_{i i} + + γ γ {Σ Σ}_{i i = = 11}^{n no} {f f}_{i i} - - - - - - ((1212))$

其中，γ为灰度与形状相似性代价系数。调整γ值，使得这两种代价在边界搜索过程中发挥大致相同的作用； Among them, γ is the gray level and shape similarity cost coefficient. Adjust the value of γ so that these two costs play roughly the same role in the boundary search process;

步骤3：使用B样条多尺度小波变换的特征图像提取；在模型建立过程中，提取有效描述不同尺度骨骼的特征图像是模型建立的基础； Step 3: Feature image extraction using B-spline multi-scale wavelet transform; in the process of model building, extracting feature images that effectively describe bones of different scales is the basis of model building;

步骤3.1：B样条多尺度小波变换：B样条函数随着样条阶数的增加而快速收敛于高斯函数，其一阶导数可以逼近最优边缘检测算子；利用B样条小波进行多尺度边缘增强；m阶B样条基函数N(m)可定义为： Step 3.1: B-spline multi-scale wavelet transform: The B-spline function quickly converges to a Gaussian function with the increase of the spline order, and its first-order derivative can approach the optimal edge detection operator; use B-spline wavelet for multi-scale Scale edge enhancement; m-order B-spline basis function N(m) can be defined as:

${N N}_{m m} ((x x)) = = \frac{11}{((m m - - 11))!!} {Σ Σ}_{t t = = 00}^{m m} {((- - 11))}^{t t} (\begin{matrix} m m \\ t t \end{matrix}) {((x x - - t t))}_{+ +}^{m m - - 11},, &ForAll; &ForAll; n no &Element; &Element; {Z Z}^{+ +} - - - - - - ((1313))$

其中， $(\begin{matrix} m \\ t \end{matrix}) = \frac{t!}{m! (t - m)!} .$ 由此可知，零阶B样条为： in, $(\begin{matrix} m \\ t \end{matrix}) = \frac{t!}{m! (t - m)!} .$ It can be seen that the zero-order B-spline is:

当m≥1时，m阶B样条函数N_m(x)可用卷积迭代关系表示： When m≥1, the m-order B-spline function N _m (x) can be expressed by the convolution iteration relationship:

显然，N_m(x)是非负且紧支撑的，其支撑宽度为suppN_m=[0,m+1]，支撑中心为(m+1)/2，且N_m关于(m+1)/2中心对称。 Obviously, N _m (x) is non-negative and tightly supported, its support width is suppN _m =[0,m+1], the support center is (m+1)/2, and N _m is about (m+1)/ 2 center symmetrical.

零价B样条的傅里叶变换为则m次B样条的傅里叶形式为： The Fourier transform of the zero-valent B-spline is Then the Fourier form of B-spline of degree m is:

${\overset{^^}{N N}}_{m m} ((ω ω)) = = {e e}^{- - i i \frac{m m + + 11}{22} ω ω} {((\frac{sin sin ((ω ω / / 22))}{ω ω / / 22}))}^{m m + + 11} - - - - - - ((1616))$

若将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)，其傅里叶变换为： If the m-time B-spline basis N _m (x) is translated to N _m [x+(m+1)/2], the center of symmetry will move to the coordinate origin, and the Fourier transform will no longer appear in the linear phase in the above formula . However, when m is even, half-integer nodes will appear. To avoid half-integer nodes, when m is an even number, translate N _m (x) to N _m [x+(m+1)/2], then the center of symmetry will move to 1/2. The function after integer translation of N _m (x) is recorded as θ _m (x), and its Fourier transform is:

${\overset{^^}{θ θ}}_{m m} ((ω ω)) = = {e e}^{- - i i \frac{ϵω ϵω}{22}} {((\frac{sin sin ((ω ω / / 22))}{ω ω / / 22}))}^{m m + + 11} - - - - - - ((1717))$

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

可求得： The two-scale relationship obtained by multiresolution analysis

can be obtained:

$H h (({e e}^{iω iω})) = = Σ Σ {h h}_{k k} {e e}^{- - ikω ikω} = = \frac{{\overset{^^}{θ θ}}_{m m} ((22 ω ω))}{{\overset{^^}{θ θ}}_{m m} ((ω ω))} = = \frac{{e e}^{- - iϵω iϵω} {((\frac{sin sin ω ω}{ω ω}))}^{m m + + 11}}{{e e}^{- - i i \frac{ϵω ϵω}{22}} {((\frac{sin sin ((ω ω / / 22))}{ω ω / / 22}))}^{m m + + 11}} = = {e e}^{- - i i \frac{ϵω ϵω}{22}} {((cos cos \frac{ω ω}{22}))}^{m m + + 11} - - - - - - ((1818))$

式中，当m为奇数时，ε=0；当m为偶数时，ε=1。由欧拉公式和二项式展开可得低通滤波器为： In the formula, when m is an odd number, ε=0; when m is an even number, ε=1. From Euler's formula and binomial expansion, the low-pass filter can be obtained as:

当m为奇数时，

When m is odd,

当m为偶数时，

When m is an even number,

选取B样条函数的一阶导数作为小波函数，即m+1次B样条函数在2^-1尺度上一阶微分函数为m次B样条小波函数： The first-order derivative of the B-spline function is selected as the wavelet function, that is, the first-order differential function of the m+1 B-spline function on the ^2-1 scale is the m-order B-spline wavelet function:

${ψ ψ}_{m m} ((x x)) = = 22 \frac{d d}{dx dx} {θ θ}_{m m + + 11} ((22 x x)) - - - - - - ((1919))$

其傅里叶变换为

且ψ_m(x)满足小波的容许性条件。由多分辨率分析中小波函数与尺度函数的二尺度关系

可求得： Its Fourier transform is

And ψ _m (x) satisfies the admissibility condition of wavelet. Two-scale relationship between wavelet function and scaling function in multi-resolution analysis

can be obtained:

$G G ((ω ω)) = = \frac{\overset{^^}{{ψ ψ}_{m m}}}{{\overset{^^}{θ θ}}_{m m} ((ω ω))} = = \underset{k k}{Σ Σ} {g g}_{k k} {e e}^{- - ikω ikω} - - - - - - ((2020))$

当m是奇数时， $G (ω) = {4 ie}^{- \frac{iω}{2}} \sin \frac{ω}{2};$ 当m是偶数时， $G (ω) = {4 ie}^{\frac{iω}{2}} \sin \frac{ω}{2} .$ 由此可知，时域有限脉冲响应为： When m is odd, $G (ω) = {4 ie}^{- \frac{iω}{2}} \sin \frac{ω}{2};$ When m is even, $G (ω) = {4 ie}^{\frac{iω}{2}} \sin \frac{ω}{2} .$ It can be seen that the time-domain finite impulse response is:

当m=3时，各项滤波器系数为k=-2,-1,0,1,2，g_k={0,0,-2,2,0}。由此可见，三阶B样条小波的滤波器的有限长度为5。 When m=3, each filter coefficient is k=-2,-1,0,1,2, g _k ={0,0,-2,2,0}. It can be seen that the finite length of the filter of the third-order B-spline wavelet is 5.

求取B样条小波变换滤波器系数： Calculate B-spline wavelet transform filter coefficients:

式中：h_k为低通滤波器系数，g_k为高通滤波器系数，利用所得h_k和g_k对图像进行快速小波分解，原始图像的行分别与h_k和g_k进行卷积，再对其列进行下采样，可以得到两幅子图像，然后沿着列的方向对两幅子图像做滤波并进行下采样，就可以得到四幅1/4大小的输出子图像，对角线细节图像，垂直细节图像，水平细节图像和近似图像； In the formula: h _k is the coefficient of the low-pass filter, g _k is the coefficient of the high-pass filter, and the obtained h _k and g _k are used to perform fast wavelet decomposition on the image, and the rows of the original image are convolved with h _k and g _k respectively, and then By downsampling its column, two sub-images can be obtained, and then the two sub-images are filtered and down-sampled along the direction of the column, and four output sub-images of 1/4 size can be obtained, and the diagonal detail image , the vertical detail image, the horizontal detail image and the approximate image;

步骤3.2：高斯滤波局部2-jet特征图像提取，由于高斯偏导数构造的特征向量能用来描述图像的局部特征，提取不同尺度小波变换后的细节图像的2-jet特征(I,I_x,I_y,I_xx,I_xy,I_yy)图来建立回归模型。根据公式(7)式中，高斯滤波的尺度为σ=2、4； Step 3.2: Gaussian filter local 2-jet feature image extraction, since the feature vector constructed by Gaussian partial derivative can be used to describe the local features of the image, extract the 2-jet features (I, I _x , I _y , I _xx , I _xy , I _yy ) graph to build a regression model. According to the formula (7), the scale of the Gaussian filter is σ=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)时，假设所有训练数据的拟合误差精度为ε，即： Step right: use support vector regression to predict the bone image, and establish the bone model by obtaining the regression function f(x)=ω·x+b: set the sample set as ( _xi , y _i ), i=1,2 ,...,l, _xi ∈ R ⁿ , the sample set S is a linear approximation of the ε-insensitive loss function. When the linear regression function f(x)=ω·x+b fits ( _xi , y _i ), it is assumed that the fitting error accuracy of all training data is ε, namely:

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

让d_i表示点(x_i,y_i)∈S到超平面f(x)的距离： Let d _i denote the distance of a point ( _xi , y _i ) ∈ S to the hyperplane f(x):

${d d}_{i i} = = \frac{| | ＜ < w w,, x x > > + + b b - - y the y | |}{\sqrt{11 + + {| | | | w w | | | |}^{22}}} - - - - - - ((23 twenty three))$

因为S集是ε-线性近似，所以有|<w,x>+b-f(xi)|≤ε,i=1,...,m。可以得到: Since the S set is an ε-linear approximation, there are |<w,x>+b-f(xi)|≤ε,i=1,...,m. Can get:

$\frac{| | ＜ < w w,, x x > > + + b b - - y the y | |}{\sqrt{11 + + {| | | | w w | | | |}^{22}}} ＜ < \frac{ϵ ϵ}{\sqrt{11 + + {| | | | w w | | | |}^{22}}},, i i = = 11,, . . . . . .,, l l - - - - - - ((24 twenty four))$

于是有 So there are

${d d}_{i i} ＜ < \frac{ϵ ϵ}{\sqrt{11 + + {| | | | w w | | | |}^{22}}},, i i = = 11,, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;,, l l - - - - - - ((2525))$

(23)式表明

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

is the upper bound on the distance of points in S to the hyperplane. The optimal approximation hyperplane of the linear approximation set S of the ε-insensitive loss function is obtained by maximizing the upper bound of the distance between the points in S and the hyperplane;

最优超平面是通过最大化

得到的(即最小化

)。因此，只要最小化||w||²就可以得到最优近似超平面。于是线性回归问题就化为： The optimal hyperplane is obtained by maximizing

obtained (i.e. minimized

). Therefore, the optimal approximate hyperplane can be obtained by simply minimizing ||w|| ² . So the linear regression problem becomes:

求下面的最优化问题： Find the following optimization problem:

$min min \frac{11}{22} {| | | | w w | | | |}^{22} - - - - - - ((2626))$

约束为|<w,x>+b-y_i|≤ε,i=1,...,m。另外，考虑拟合误差的情况，引入松弛因子

当划分有误差时，

都大于0，误差不存在取0； The constraints are |<w,x>+by _i |≤ε,i=1,...,m. In addition, considering the fitting error, the relaxation factor is introduced

When there is an error in the division,

Both are greater than 0, and if the error does not exist, take 0;

因而，标准ε不敏感支持向量回归机可以表示为： Therefore, the standard ε-insensitive support vector regression machine can be expressed as:

$\{\begin{matrix} \underset{w w,, b b,, ξ ξ}{min min} \frac{11}{22} {| | | | w w | | | |}^{22} + + C C {Σ Σ}_{i i = = 11}^{l l} (({ξ ξ}_{i i} + + {ξ ξ}_{i i}^{* *})) \\ s the s . . t t . . {y the y}_{i i} - - w w \cdot \cdot {x x}_{i i} - - b b ＜ < ϵ ϵ + + {ξ ξ}_{i i} \\ w w \cdot \cdot {x x}_{i i} + + b b - - {y the y}_{i i} ＜ < ϵ ϵ + + {ξ ξ}_{i i}^{* *} \\ {ξ ξ}_{i i} &GreaterEqual; &Greater Equal; 00 \\ {ξ ξ}_{i i}^{* *} &GreaterEqual; &Greater Equal; 00,, i i = = 1,2 1,2,, . . . . . .,, l l \end{matrix} - - - - - - ((2727))$

其中，C>0为平衡因子。引入拉格朗日函数为： Among them, C>0 is the balance factor. Introduce the Lagrangian function as:

$\begin{matrix} L L ((w w,, b b,, ξ ξ,, α α,, {α α}^{* *},, γ γ)) = = \frac{11}{22} {| | | | w w | | | |}^{22} + + C C {Σ Σ}_{i i = = 11}^{n no} (({ξ ξ}_{i i} + + {ξ ξ}_{i i}^{* *})) - - {Σ Σ}_{i i = = 11}^{n no} {α α}_{i i} [[{ξ ξ}_{i i} + + ϵ ϵ - - {y the y}_{i i} + + f f (({x x}_{i i}))]] \\ - - {Σ Σ}_{i i = = 11}^{n no} {α α}_{i i}^{* *} [[{ξ ξ}_{i i} + + ϵ ϵ + + {y the y}_{i i} - - f f (({x x}_{i i}))]] - - {Σ Σ}_{i i = = 11}^{n no} {γ γ}_{i i} (({ξ ξ}_{i i} {γ γ}_{i i} + + {ξ ξ}_{i i}^{* *} {γ γ}_{i i}^{* *})) \end{matrix} - - - - - - ((2828))$

其中，α_i，

γi≥0，i=1,...,n。函数L的极值应满足条件

于是得到下面的式子： Among them, α _i ,

γi≥0, i=1,...,n. The extremum of the function L should satisfy the condition

Then the following formula is obtained:

$w w = = {Σ Σ}_{i i = = 00}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) {x x}_{i i} - - - - - - ((2929))$

${Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) = = 00 - - - - - - ((3030))$

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

$C C - - {α α}_{i i}^{* *} - - {γ γ}_{i i}^{* *} = = 00,, i i = = 11,, . . . . . .,, l l - - - - - - ((3232))$

将(29)-(32)代入到(28)中，得到优化问题的对偶形式为： Substituting (29)-(32) into (28), the dual form of the optimization problem is obtained as:

$max max - - \frac{11}{22} {Σ Σ}_{i i,, j j = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) (({α α}_{j j} - - {α α}_{j j}^{* *})) ＜ < {x x}_{i i},, {x x}_{j j} > > + + {Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) {y the y}_{i i} - - {Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} + + {α α}_{i i}^{* *})) ϵ ϵ - - - - - - ((3333))$

约束为： Constraints are:

${Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) = = 0,0 0,0 ＜ < {α α}_{i i},, {α α}_{i i}^{* *} ＜ < C C,, i i = = 11,, . . . . . .,, l l - - - - - - ((3434))$

由于骨骼灰度预测为非线性回归，所以先使用一个非线性映射φ，把数据映射到一个高维特征空间，然后在高维特征空间进行线性回归建立模型。由于在上面的优化过程中，只考虑到高维特征空间中的内积运算，因此用一个核函数k(x,y)代替<φ(x),φ(y)>就可以实现非线性回归。于是，非线性回归的优化方程为最大化下面的函数： Since the bone grayscale prediction is a nonlinear regression, a nonlinear mapping φ is first used to map the data to a high-dimensional feature space, and then a linear regression is performed in the high-dimensional feature space to build a model. Since only the inner product operation in the high-dimensional feature space is considered in the above optimization process, nonlinear regression can be realized by replacing <φ(x),φ(y)> with a kernel function k(x,y) . Therefore, the optimization equation for nonlinear regression is to maximize the following function:

$\begin{matrix} W W ((α α,, {α α}^{* *})) = = - - \frac{11}{22} {Σ Σ}_{i i,, j j}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) (({α α}_{j j} - - {α α}_{j j}^{* *})) K K (({x x}_{i i},, {x x}_{j j})) \\ + + {Σ Σ}_{i i,, j j}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) {y the y}_{i i} - - {Σ Σ}_{i i,, j j}^{n no} (({α α}_{i i} + + {α α}_{i i}^{* *})) ϵ ϵ \end{matrix} - - - - - - ((3535))$

其约束条件为(32)。求解出α的值后，就可以得到f(x)： Its constraints are (32). After solving the value of α, f(x) can be obtained:

$f f ((x x)) = = {Σ Σ}_{i i = = 11}^{N N} (({α α}_{i i} - - {α α}_{i i}^{* *})) K K (({x x}_{i i},, x x)) + + b b - - - - - - ((3636))$

通常情况下，大部分α_i或

的值将为零，不为零的α_i或所对应的样本称为支持向量。根据最优化条件(KKT条件)，在鞍点有： Typically, most α _i or

The value of will be zero, non-zero α _i or The corresponding samples are called support vectors. According to the optimization condition (KKT condition), at the saddle point:

$\{\begin{matrix} {α α}_{i i} [[{ξ ξ}_{i i} + + ϵ ϵ - - {y the y}_{i i} + + f f (({x x}_{i i}))]] = = 00 & ,, i i = = 11,, . . . . . . l l \\ {a a}_{i i}^{* *} [[{ξ ξ}_{i i}^{* *} + + ϵ ϵ - - {y the y}_{i i} + + f f (({x x}_{i i}))]] = = 00 & ,, i i = = 11,, . . . . . . l l \\ ((C C - - {α α}_{i i})) {ξ ξ}_{i i} = = 00 & ,, i i = = 11,, . . . . . . l l \\ ((C C - - {α α}_{i i}^{* *})) {ξ ξ}_{i i}^{* *} = = 00 & ,, i i = = 11,, . . . . . . l l \end{matrix} - - - - - - ((3737))$

于是得到b的计算式如下： So the calculation formula of b is as follows:

$b = y_{i} - ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),$ 当α_j∈(0,C) $b = {the y}_{i} - ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),$ When α _j ∈ (0,C)

$b = y_{i} + ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),$ 当

∈(0,C)

b = {the y}_{i} + ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),

when

∈(0,C)

用任意一个支持向量可以计算出b的值，得到回归函数f(x)；建立预测模型后，可以根据肺部x光图像的灰度值分布来预测产生骨结构图像； The value of b can be calculated with any support vector, and the regression function f(x) can be obtained; after the prediction model is established, the bone structure image can be predicted according to the gray value distribution of the lung X-ray image;

步骤5：对正常胸片与预测的骨骼图像做图像减法来预测的软组织图像。 Step 5: Perform image subtraction on the normal chest radiograph and the predicted bone image to obtain the predicted soft tissue image. the

所述步骤1.2：对齐训练形状向量的具体步骤如下： The step 1.2: The specific steps of aligning the training shape vectors are as follows:

所述的将分块的数据映射到特征空间的具体方法如下： The specific method of mapping the block data to the feature space is as follows:

步骤1.2.1：旋转、缩放和平移每个肺区域形状，使其与训练集中的第一个形状对齐； Step 1.2.1: Rotate, scale and translate each lung region shape to align with the first shape in the training set;

步骤1.2.2：根据对齐形状，计算平均形状； Step 1.2.2: Calculate the average shape according to the aligned shape;

步骤1.2.3：旋转、缩放和平移平均形状使其与第一个形状对齐； Step 1.2.3: Rotate, scale and translate the average shape to align with the first shape;

步骤1.2.权：重新将每个形状与当前平均形状对齐； Step 1.2. Weight: Re-align each shape with the current average shape;

步骤1.2.5：如果过程收敛或者到指定循环次数，退出；否则转到步骤。 Step 1.2.5: If the process converges or reaches the specified number of cycles, exit; otherwise go to step. the

所述步骤3中对胸部x光图像做三尺度小波分解。 In the step 3, a three-scale wavelet decomposition is performed on the chest x-ray image. the

本发明的优点为：本发明先对双能剪影图像中的正常胸片建立肺部轮廓模型，再应用三阶B样条小波变换对图像做三尺度小波分解，将分解后的图像使用2-jet方法提取特征图像，再使用支持向量回归机建立骨骼模型，将骨骼模型应用到X光图像中，从而提取对应的预测骨骼图像，最后用X光图像与预测的骨骼图像做图像减法从而得到预测的软组织图像。因为抑制骨骼方法大都采用神经网络方法，而此方法易陷入局部最优，训练结果不太稳定，所以本专利采用支持向量回归机来进行骨骼预测，得到了更优的结果。 The advantages of the present invention are: the present invention first establishes the lung contour model on the normal chest radiograph in the dual-energy silhouette image, and then applies the third-order B-spline wavelet transform to perform three-scale wavelet decomposition on the image, and uses 2- The jet method extracts the feature image, and then uses the support vector regression machine to establish the bone model, applies the bone model to the X-ray image, thereby extracting the corresponding predicted bone image, and finally uses the X-ray image and the predicted bone image to perform image subtraction to obtain the prediction soft tissue images. Because most of the bone suppression methods use the neural network method, and this method is easy to fall into a local optimum, and the training results are not stable, so this patent uses a support vector regression machine for bone prediction, and better results are obtained. the

附图说明 Description of drawings

图1是本发明的总体流程图； Fig. 1 is overall flow chart of the present invention;

图2是本发明的肺区分割流程图； Fig. 2 is the segmentation flow chart of lung area of the present invention;

图3是本发明中高斯滤波局部2-jet提取特征图像示意图； Fig. 3 is a schematic diagram of feature image extracted by Gaussian filtering local 2-jet in the present invention;

图4是本发明中最优近似超平面的示意图； Fig. 4 is the schematic diagram of optimal approximation hyperplane among the present invention;

图5是本发明中支持向量回归示意图。 Fig. 5 is a schematic diagram of support vector regression in the present invention. the

具体实施方式Detailed ways

如图1至图5所示，本发明采取按如下步骤进行： As shown in Figures 1 to 5, the present invention takes the following steps to carry out:

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

in,

m (the s, θ) [\begin{matrix} x_{jk} \\ {the y}_{jk} \end{matrix}] = (\begin{matrix} (the s \cos θ) x_{jk} - (the s \sin θ) {the y}_{jk} \\ (the s \sin θ) x_{jk} + (the s \cos θ) {the y}_{jk} \end{matrix}),

t=(t _x ,t _y ,...,t _x ,t _y ) ^T where: W is a weighted diagonal matrix of a point, which is composed of the weight w _k of each boundary point. Let R _kl represent the distance between the kth point and the lth point, Indicates the change of the distance in the training set, then the weight of the kth point

步骤2：结合灰度信息和形状信息的肺实质分割：同时利用多张特征图像中边界点的灰度与形状信息，使得搜索到的边界灰度、形状信息与训练图像相似； Step 2: Lung parenchyma segmentation combining grayscale information and shape information: using the grayscale and shape information of boundary points in multiple feature images at the same time, so that the searched boundary grayscale and shape information are similar to the training images;

步骤2.1：提取特征图像；进行主成分分析对数据进行降维，对于图像I中点p附近的点p+△p的灰度可由Taylor展开获得， $I (p + Δp) \approx Σ_{n = 0}^{N} \frac{1}{n!} (Δ p_{i 1}, \cdot \cdot \cdot, Δ p_{in} \frac{{&PartialD;}^{n} I (p)}{{&PartialD; i}_{1 \cdot \cdot \cdot} {&PartialD; i}_{n}})$ Step 2.1: Extract the feature image; perform principal component analysis to reduce the dimensionality of the data, and the gray level of the point p+△p near the point p in the image I can be obtained by Taylor expansion, $I (p + Δp) \approx Σ_{no = 0}^{N} \frac{1}{no!} (Δ p_{i 1}, \cdot \cdot &Center Dot;, Δ p_{in} \frac{{&PartialD;}^{no} I (p)}{{&PartialD; i}_{1 &Center Dot; &Center Dot; &Center Dot;} {&PartialD; i}_{no}})$

(2) (2)

$\begin{matrix} I I ((p p + + Δp Δp)) = = {Σ Σ}_{n no = = 00}^{22} \frac{11}{n no!!} ((Δ Δ {p p}_{i i 11},, \cdot \cdot \cdot \cdot \cdot &Center Dot;,, Δ Δ {p p}_{in in} \frac{{&PartialD; &PartialD;}^{n no} I I ((p p))}{{&PartialD; &PartialD; i i}_{11} \cdot \cdot \cdot \cdot \cdot \cdot {&PartialD; &PartialD; i i}_{n no}})) = = I I ((p p)) + + Δ Δ {p p}_{i i 11} \frac{&PartialD; &PartialD; I I ((p p))}{{&PartialD; &PartialD; i i}_{11}} + + \frac{11}{22!!} ((Δ Δ {p p}_{i i 11} Δ Δ {p p}_{i i 22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; i i}_{11} {&PartialD; &PartialD; i i}_{22}})) \\ = = I I ((p p)) + + Δ Δ {p p}_{x x} \frac{&PartialD; &PartialD; I I ((p p))}{&PartialD; &PartialD; x x} + + \frac{11}{22!!} ((Δ Δ {p p}_{x x}^{22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; x x}^{22}} + + Δ Δ {p p}_{x x} Δ Δ {p p}_{y the y} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{&PartialD; &PartialD; x x &PartialD; &PartialD; y the y})) \\ + + Δ Δ {p p}_{y the y} \frac{&PartialD; &PartialD; I I ((p p))}{&PartialD; &PartialD; y the y} + + \frac{11}{22!!} ((Δ Δ {p p}_{y the y}^{22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; y the y}^{22}} + + Δ Δ {p p}_{y the y} Δ Δ {p p}_{x x} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{&PartialD; &PartialD; y the y &PartialD; &PartialD; x x})) \end{matrix} - - - - - - ((33))$

对图像进行高斯平滑处理，对于图像I(p)进行高斯滤波后的图像为A(p,σ)=(I*G)(p,σ)。滤波图像的n阶偏导数为

则根据卷积性质可知： Gaussian smoothing is performed on the image, and the image after Gaussian filtering on the image I(p) is A(p,σ)=(I*G)(p,σ). The nth order partial derivative of the filtered image is

According to the convolution properties, we can know:

${A A}_{{i i}_{11},, \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; {i i}_{n no}} ((p p,, σ σ)) = = {((I I * * G G))}_{{i i}_{11},, \cdot &Center Dot; \cdot \cdot \cdot \cdot {i i}_{n no}} ((p p,, σ σ)) = = ((I I * * {G G}_{{i i}_{11},, \cdot \cdot \cdot &Center Dot; \cdot \cdot {i i}_{n no}})) ((p p,, σ σ)) - - - - - - ((55))$

式中，v_i=p_i+1-p_i，表示第i个边界点的形状特征，

$C C = = [\begin{matrix} {h h}_{1,1 1,1} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {h h}_{11,, k k} & \cdot &Center Dot; \cdot \cdot \cdot \cdot & {h h}_{11,, m m} \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \\ {h h}_{i i,, 11} & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & {h h}_{i i,, k k} & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & {h h}_{i i,, m m} \\ \cdot \cdot & \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ {h h}_{n no,, 11} & \cdot \cdot \cdot \cdot \cdot \cdot & {h h}_{n no,, k k} & \cdot \cdot \cdot \cdot \cdot \cdot & {h h}_{n no,, m m} \end{matrix}] - - - - - - ((99))$

式中：为在特征图像上，以点p_i为圆心，r_c为半径的圆上的n_c个点的灰度，k=1......n_c。

分别为训练图像的第j个特征图像中第i个边界点的周围像素点灰度的均值及协方差。N为特征图像总数， In the formula: is the gray level of n _c points on a circle with point p _i as the center and r _c as the radius on the feature image, k=1...n _c .

are the mean and covariance of the gray values of the surrounding pixels of the i-th boundary point in the j-th feature image of the training image, respectively. N is the total number of feature images,

零阶B样条的傅里叶变换为则m次B样条的傅里叶形式为： The Fourier transform of the zero-order B-spline is Then the Fourier form of B-spline of degree m is:

由多分辨率分析得到二尺度关系式可求得： The two-scale relationship obtained by multiresolution analysis can be obtained:

当m为奇数时，

When m is odd,

当m为偶数时， When m is an even number,

其傅里叶变换为

可求得： Its Fourier transform is

can be obtained:

当m=3时，各项滤波器系数为k=-2,-1,0,1,2，

g_k={0,0,-2,2,0}。由此可见，三阶B样条小波的滤波器的有限长度为5。 When m=3, each filter coefficient is k=-2,-1,0,1,2,

g _k ={0,0,-2,2,0}. It can be seen that the finite length of the filter of the third-order B-spline wavelet is 5.

求取3阶B样条小波变换滤波器系数： Calculate the third-order B-spline wavelet transform filter coefficients:

本实施例对胸部x光图像做三尺度小波分解，获得12张小波分解图像，其中3张近似图像，3张水平细节图像，3张垂直细节图像，3张对角线细节图像。由于在胸部x光图像中，很少有呈现对角线分布的骨骼结构。因此，不提取对角线细节图像的特征。 In this embodiment, three-scale wavelet decomposition is performed on chest X-ray images to obtain 12 wavelet-decomposed images, including 3 approximate images, 3 horizontal detail images, 3 vertical detail images, and 3 diagonal detail images. Because in the chest x-ray image, there are few skeletal structures that show a diagonal distribution. Therefore, features of diagonal detail images are not extracted. the

步骤权：使用支持向量回归进行骨骼图像的预测，通过求得回归函数f(x)=ω·x+b，建立骨骼模型；设样本集为(x_i,y_i)，i=1,2,...,l，x_i∈Rⁿ，样本集S是ε-不敏感损失函数的线性近似。当线性回归函数f(x)=ω·x+b拟合(x_i,y_i)时，假设所有训练数据的拟合误差精度为ε，即： Step right: use support vector regression to predict the bone image, and establish the bone model by obtaining the regression function f(x)=ω·x+b; set the sample set as ( _xi , y _i ), i=1,2 ,...,l, _xi ∈ R ⁿ , the sample set S is a linear approximation of the ε-insensitive loss function. When the linear regression function f(x)=ω·x+b fits ( _xi , y _i ), it is assumed that the fitting error accuracy of all training data is ε, namely:

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

因为S集是ε-线性近似，所以有|<w,x>+b-f(x_i)|≤ε,i=1,...,m。可以得到: Since the set S is an ε-linear approximation, there is |<w,x>+bf(x _i )|≤ε,i=1,...,m. can get:

于是有 So there are

${d d}_{i i} ＜ < \frac{ϵ ϵ}{\sqrt{11 + + {| | | | w w | | | |}^{22}}},, i i = = 11,, \cdot \cdot \cdot \cdot \cdot \cdot,, l l - - - - - - ((2525))$

(23)式表明

最优超平面是通过最大化

得到的(即最小化

obtained (i.e. minimized

求下面的最优化问题： Find the following optimization problem:

$min min \frac{11}{22} {| | | | w w | | | |}^{22} - - - - - - ((2626))$

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

都大于0，误差不存在取0； The constraints are |<w,x>+by _i |≤ε,i=1,...,m. In addition, considering the fitting error, the relaxation factor is introduced When there is an error in the division,

Both are greater than 0, and if the error does not exist, take 0;

$\{\begin{matrix} \underset{w w,, b b,, ξ ξ}{min min} \frac{11}{22} {| | | | w w | | | |}^{22} + + C C {Σ Σ}_{i i = = 11}^{l l} (({ξ ξ}_{i i} + + {ξ ξ}_{i i}^{* *})) \\ s the s . . t t . . {y the y}_{i i} - - w w \cdot &Center Dot; {x x}_{i i} - - b b ＜ < ϵ ϵ + + {ξ ξ}_{i i} \\ w w \cdot &Center Dot; {x x}_{i i} + + b b - - {y the y}_{i i} ＜ < ϵ ϵ + + {ξ ξ}_{i i}^{* *} \\ {ξ ξ}_{i i} &GreaterEqual; &Greater Equal; 00 \\ {ξ ξ}_{i i}^{* *} &GreaterEqual; &Greater Equal; 00,, i i = = 1,2 1,2,, . . . . . .,, l l \end{matrix} - - - - - - ((2727))$

其中，α_i，

γi≥0，i=1,...,n。函数L的极值应满足条件

于是得到下面的式子： Among them, α _i ,

γi≥0, i=1,...,n. The extremum of the function L should satisfy the condition

Then the following formula is obtained:

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

约束为： Constraints are:

通常情况下，大部分α_i或

$b = y_{i} - ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),$ 当αj∈(0,C) $b = {the y}_{i} - ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),$ When αj∈(0,C)

$b = y_{i} + ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),$ 当∈(0,C) $b = {the y}_{i} + ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),$ when ∈(0,C)

Claims

1. A method for bone suppression in a chest X-ray image, the specific steps are as follows:

Step 1: Determination of the initial contour position of the model lung;

Step 1.1: mark the boundary points of the lung boundary of each image in the training set;

Step 1.2: Align the training shape vectors; align the training shapes by scaling, rotating and translating operations so that they are aligned as closely as possible, let x _i be the vector of n points in the i-th shape in the training set, x _i =(x _i0 ,y _i0 ,...x _ik ,y _ik ,...,x _in-1 ,y _in-1 ) ^T , where (x _ik ,y _ik ) is the kth point in the i-th shape, Given two similar shapes x _i and x _j , choose rotation angle θ, scaling s, and translation (t _x , t _y ), then M(s,θ)[x] represents the rotation angle θ and scaling ratio s transform, mapping _xi to _xj such that the following weighted sum is minimized:

E _j =( _xi -M(s,θ)[x _j ]-t) ^T W( _xi -M(s,θ)[x _j ]-t) (1)

in,

t=(t _x ,t _y ,...,t _x ,t _y ) ^T , where: W is a weighted diagonal matrix of a point, which is composed of the weight w _k of each boundary point. Let R _kl represent the distance between the kth point and the lth point, Indicates the change of the distance in the training set, then the weight of the kth point

Step 1.3: Establish the model of the initial contour position of the lung; after the training shape vectors are aligned, use the principal component analysis method to find out the statistical information of the shape change, and build the model accordingly.

It is characterized in that: it also includes the following steps:

Step 2: Segmentation of lung parenchyma combining grayscale information and shape information: using the grayscale and shape information of boundary points in multiple feature images at the same time, so that the searched boundary grayscale and shape information are similar to the training images;

Step 2.1: Extract the feature image; perform principal component analysis to reduce the dimensionality of the data, and the gray level of the point p+△p near the point p in the image I can be obtained by Taylor expansion,

I (p + Δp) \approx Σ_{no = 0}^{N} \frac{1}{no!} (Δ p_{i 1}, &Center Dot; &Center Dot; &Center Dot;, Δ p_{in} \frac{{&PartialD;}^{no} I (p)}{{&PartialD; i}_{1 \cdot \cdot \cdot} {&PartialD; i}_{no}}) - - - (2)

It can be obtained that the second-order Taylor expansion of image I is:

\begin{matrix} I I ((p p + + Δp Δp)) = = {Σ Σ}_{n no = = 00}^{22} \frac{11}{n no!!} ((Δ Δ {p p}_{i i 11},, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;,, Δ Δ {p p}_{in in} \frac{{&PartialD; &PartialD;}^{n no} I I ((p p))}{{&PartialD; &PartialD; i i}_{11} \cdot \cdot \cdot &Center Dot; \cdot \cdot {&PartialD; &PartialD; i i}_{n no}})) = = I I ((p p)) + + Δ Δ {p p}_{i i 11} \frac{&PartialD; &PartialD; I I ((p p))}{{&PartialD; &PartialD; i i}_{11}} + + \frac{11}{22!!} ((Δ Δ {p p}_{i i 11} Δ Δ {p p}_{i i 22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; i i}_{11} {&PartialD; &PartialD; i i}_{22}})) \\ = = I I ((p p)) + + Δ Δ {p p}_{x x} \frac{&PartialD; &PartialD; I I ((p p))}{&PartialD; &PartialD; x x} + + \frac{11}{22!!} ((Δ Δ {p p}_{x x}^{22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; x x}^{22}} + + Δ Δ {p p}_{x x} Δ Δ {p p}_{y the y} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{&PartialD; &PartialD; x x &PartialD; &PartialD; y the y})) \\ + + Δ Δ {p p}_{y the y} \frac{&PartialD; &PartialD; I I ((p p))}{&PartialD; &PartialD; y the y} + + \frac{11}{22!!} ((Δ Δ {p p}_{y the y}^{22} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{{&PartialD; &PartialD; y the y}^{22}} + + Δ Δ {p p}_{y the y} Δ Δ {p p}_{x x} \frac{{&PartialD; &PartialD;}^{22} I I ((p p))}{&PartialD; &PartialD; y the y &PartialD; &PartialD; x x})) \end{matrix} - - - - - - ((33))

Gaussian smoothing is performed on the image, and the image after Gaussian filtering for the image I(p) is A(p,σ)=(I*G)(p,σ), and the nth order partial derivative of the filtered image is

According to the convolution properties, we can know:

{A A}_{{i i}_{11},, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; {i i}_{n no}} ((p p,, σ σ)) = = {((I I * * G G))}_{{i i}_{11},, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; {i i}_{n no}} ((p p,, σ σ)) = = ((I I * * {G G}_{{i i}_{11},, \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; {i i}_{n no}})) ((p p,, σ σ)) - - - - - - ((55))

In the formula: the subscripts i ₁ ,...i _n represent the partial derivatives of the data to the variables respectively, and the Gaussian partial derivatives are regarded as a filter, and the combination of the Gaussian partial derivative responses constitutes a complete feature description;

Use the local 2-jet feature image to obtain the candidate points of the boundary points and the gray cost of each candidate point. The local 2-jet feature is:

J ² [I](p,σ)={I,I _x ,I _y ,I _xx ,I _xy ,I _yy } (7)

Step 2.2: Select the candidate point of the boundary point: For each point of the initial lung boundary, calculate the similarity between the gray levels of all pixels in the search area of this point in all feature images and the gray level of the corresponding point in the training feature image, and select The point with the greatest similarity is used as the candidate point of the boundary point. The degree of similarity is the Mahalanobis distance h _i from the gray level of the pixels around the point in all feature images to the set of gray levels of pixels around the corresponding point in the feature image of the training sample:

Step 2.3: Use dynamic programming to segment the lung area; in the boundary point search area, the gray similarity cost of the pixel point is the similarity between the gray level of the surrounding pixels of this point and the surrounding pixel gray level of the corresponding boundary point in the training image h _i In the boundary point search area, the gray similarity cost of the pixel point is the similarity h _i of the gray level of the surrounding pixels of this point and the surrounding pixel gray level of the corresponding boundary point in the training image; the i-th boundary in the image The shape similarity cost of a point is defined as:

{f f}_{i i} = = {(({v v}_{i i} - - {u u}_{{v v}_{i i}}))}^{T T} {(({COV COV}_{{v v}_{i i}}))}^{- - 11} (({v v}_{i i} - - {u u}_{{v v}_{i i}})) - - - - - - ((88))

In the formula, v _i =p _i+1 -p _i represents the shape feature of the i-th boundary point,

Represent the mean and covariance of the i-th boundary point shape features in all training images, respectively,

For the i-th boundary point p _i in the test image, there are m candidate points with a smaller gray-scale similarity cost in the specified search area, then n boundary points will generate an n×m gray-scale cost matrix:

C C = = [\begin{matrix} {h h}_{1,1 1,1} & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & {h h}_{11,, k k} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {h h}_{11,, m m} \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ {h h}_{i i,, 11} & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & {h h}_{i i,, k k} & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & {h h}_{i i,, m m} \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ {h h}_{n no,, 11} & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & {h h}_{n no,, k k} & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {h h}_{n no,, m m} \end{matrix}] - - - - - - ((99))

{h h}_{i i} = = {Σ Σ}_{j j = = 11}^{N N} {(({g g}_{i i}^{j j} - - {u u}_{{g g}_{i i}^{j j}}))}^{T T} {(({s the s}_{{g g}_{i i}^{j j}}))}^{- - 11} (({g g}_{i i}^{j j} - - {u u}_{{g g}_{i i}^{j j}})) - - - - - - ((1010))

{g g}_{i i}^{j j} = = {p p}_{i i} + + {r r}_{c c} [\begin{matrix} cos cos ((\frac{22 π π}{{n no}_{c c}} ((k k - - 11)))) \\ sin sin ((\frac{22 π π}{{n no}_{c c}} ((k k - - 11)))) \end{matrix}] - - - - - - ((1111))

In the formula:

is the gray level of n _c points on a circle with point p _i as the center and r _c as the radius on the feature image, k=1...n _c .

are the mean value and covariance of the surrounding pixel gray levels of the i-th boundary point in the j-th feature image of the training image, N is the total number of feature images,

Searching for the optimal contour is to find an optimal path, select one element in each row in the matrix C, and when selecting along the path, the sum of the gray level and shape similarity costs is the smallest, namely:

J J (({k k}_{11} . . . . . . . . . . {k k}_{n no})) = = {Σ Σ}_{i i = = 11}^{n no} {h h}_{i i} + + γ γ {Σ Σ}_{i i = = 11}^{n no} {f f}_{i i} - - - - - - ((1212))

Among them, γ is the gray level and shape similarity cost coefficient. Adjust the value of γ so that these two costs play roughly the same role in the boundary search process;

Step 3: Feature image extraction using B-spline multi-scale wavelet transform; in the process of model building, extracting feature images that effectively describe bones of different scales is the basis of model building;

Step 3.1: B-spline multi-scale wavelet transform: The B-spline function quickly converges to a Gaussian function with the increase of the spline order, and its first-order derivative can approach the optimal edge detection operator; use B-spline wavelet for multi-scale Scale edge enhancement; use function:

H h (({e e}^{iω iω})) = = Σ Σ {h h}_{k k} {e e}^{- - ikω ikω} = = \frac{{\overset{^^}{θ θ}}_{m m} ((22 ω ω))}{{\overset{^^}{θ θ}}_{m m} ((ω ω))} = = \frac{{e e}^{- - iϵω iϵω} {((\frac{sin sin ω ω}{ω ω}))}^{m m + + 11}}{{e e}^{- - i i \frac{ϵω ϵω}{22}} {((\frac{sin sin ((ω ω / / 22))}{ω ω / / 22}))}^{m m + + 11}} = = {e e}^{- - i i \frac{ϵω ϵω}{22}} {((cos cos \frac{ω ω}{22}))}^{m m + + 11} - - - - - - ((1818))

G G ((ω ω)) = = \frac{\overset{^^}{{ψ ψ}_{m m}} ((22 ω ω))}{{\overset{^^}{θ θ}}_{m m} ((ω ω))} = = \underset{k k}{Σ Σ} {g g}_{k k} {e e}^{- - ikω ikω} - - - - - - ((2020))

Find B-spline wavelet transform filter coefficients:

In the formula: h _k is the coefficient of the low-pass Gaussian filter, g _k is the coefficient of the high-pass Gaussian filter, and the obtained h _k and g _k are used to perform fast wavelet decomposition on the image, and the rows of the original image are respectively convolved with h _k and g _k , and then downsampling its column, you can get two sub-images, and then filter and down-sample the two sub-images along the direction of the column, you can get four output sub-images of 1/4 size, diagonal detail image, vertical detail image, horizontal detail image and approximate image;

Step 3.2: Gaussian filter local 2-jet feature image extraction, because the feature vector constructed by Gaussian partial derivatives can be used to describe the local features of the image, according to the formula (7), extract the 2-jet of the detail image after wavelet transformation of different scales jet feature (I, I _x , I _y , I _xx , I _xy , I _yy ) graph to build a regression model;

Step 4: Use support vector regression to predict the bone image, and establish the bone model by obtaining the regression function f(x)=ω·x+b; set the sample set as ( _xi , y _i ), i=1,2 ,...,l, _xi ∈ R ⁿ , the sample set S is a linear approximation of the ε-insensitive loss function. When the linear regression function f(x)=ω·x+b fits ( _xi , y _i ), it is assumed that the fitting error accuracy of all training data is ε, namely:

|y _i -f(x _i )|≤εi=1,2,...,l (22)

Let d _i denote the distance of a point ( _xi , y _i ) ∈ S to the hyperplane f(x):

{d d}_{i i} = = \frac{| | < < w w,, x x > > + + b b - - y the y | |}{\sqrt{11 + + {| | | | w w | | | |}^{22}}} - - - - - - ((23 twenty three))

Since the set S is an ε-linear approximation, there is |<w,x>+bf(x _i )|≤ε,i=1,...,m. can get:

\frac{| | < < w w,, x x > > + + b b - - y the y | |}{\sqrt{11 + + {| | | | w w | | | |}^{22}}} ＜ < \frac{ϵ ϵ}{\sqrt{11 + + {| | | | w w | | | |}^{22}}},, i i = = 11,, . . . . . .,, l l - - - - - - ((24 twenty four))

So there is

{d d}_{i i} ＜ < \frac{ϵ ϵ}{\sqrt{11 + + {| | | | w w | | | |}^{22}}},, i i = = 11,, \cdot \cdot \cdot \cdot \cdot \cdot,, l l - - - - - - ((2525))

(23) shows that

The optimal hyperplane is obtained by maximizing

obtained (i.e. minimized

Find the following optimization problem:

min min \frac{11}{22} {| | | | w w | | | |}^{22} - - - - - - ((2626))

The constraints are |<w,x>+by _i |≤ε,i=1,...,m. In addition, considering the fitting error, the relaxation factor is introduced

When there is an error in the division,

Both are greater than 0, and if the error does not exist, take 0;

Therefore, the standard ε-insensitive support vector regression machine can be expressed as:

\{\begin{matrix} \underset{w w,, b b,, ξ ξ}{min min} \frac{11}{22} {| | | | w w | | | |}^{22} + + C C {Σ Σ}_{i i = = 11}^{l l} (({ξ ξ}_{i i} + + {ξ ξ}_{i i}^{* *})) \\ s the s . . t t . . {y the y}_{i i} - - w w \cdot \cdot {x x}_{i i} - - b b ＜ < ϵ ϵ + + {ξ ξ}_{i i} \\ w w \cdot \cdot {x x}_{i i} + + b b - - {y the y}_{i i} ＜ < ϵ ϵ + + {ξ ξ}_{i i}^{* *} \\ {ξ ξ}_{i i} &GreaterEqual; &Greater Equal; 00 \\ {ξ ξ}_{i i}^{* *} &GreaterEqual; &Greater Equal; 00,, i i = = 1,2 1,2,, . . . . . .,, l l \end{matrix} - - - - - - ((2727))

Among them, C>0 is the balance factor. The Lagrangian function is introduced as:

\begin{matrix} L L ((w w,, b b,, ξ ξ,, α α,, {α α}^{* *},, γ γ)) = = \frac{11}{22} {| | | | w w | | | |}^{22} + + C C {Σ Σ}_{i i = = 11}^{n no} (({ξ ξ}_{i i} + + {ξ ξ}_{i i}^{* *})) - - {Σ Σ}_{i i = = 11}^{n no} {α α}_{i i} [[{ξ ξ}_{i i} + + ϵ ϵ - - {y the y}_{i i} + + f f (({x x}_{i i}))]] \\ - - {Σ Σ}_{i i = = 11}^{n no} {α α}_{i i}^{* *} [[{ξ ξ}_{i i} + + ϵ ϵ + + {y the y}_{i i} - - f f (({x x}_{i i}))]] - - {Σ Σ}_{i i = = 11}^{n no} {γ γ}_{i i} (({ξ ξ}_{i i} {γ γ}_{i i} + + {ξ ξ}_{i i}^{* *} {γ γ}_{i i}^{* *})) \end{matrix} - - - - - - ((2828))

Among them, α _i ,

γ _i ≥ 0, i=1,...,n. The extremum of the function L should satisfy the condition Then the following formula is obtained:

w w = = {Σ Σ}_{i i = = 00}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) {x x}_{i i} - - - - - - ((2929))

{Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) = = 00 - - - - - - ((3030))

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

C C - - {α α}_{i i}^{* *} - - {γ γ}_{i i}^{* *} = = 00,, i i = = 11,, . . . . . .,, l l - - - - - - ((3232))

Substituting (29)-(32) into (28), the dual form of the optimization problem is obtained as:

max max - - \frac{11}{22} {Σ Σ}_{i i,, j j = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) (({α α}_{j j} - - {α α}_{j j}^{* *})) < < {x x}_{i i},, {x x}_{j j} > > + + {Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) {y the y}_{i i} - - {Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} + + {α α}_{i i}^{* *})) ϵ ϵ - - - - - - ((3333))

Constraints are:

{Σ Σ}_{i i = = 11}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) = = 0,0 0,0 ＜ < {α α}_{i i},, {α α}_{i i}^{* *} ＜ < C C,, i i = = 11,, . . . . . .,, l l - - - - - - ((3434))

Since the bone grayscale prediction is a nonlinear regression, a nonlinear mapping φ is first used to map the data to a high-dimensional feature space, and then a linear regression is performed in the high-dimensional feature space to build a model. Since only the inner product operation in the high-dimensional feature space is considered in the above optimization process, nonlinear regression can be realized by replacing <φ(x),φ(y)> with a kernel function k(x,y) . Thus, the optimization equation for nonlinear regression is to maximize the following function:

\begin{matrix} W W ((α α,, {α α}^{* *})) = = - - \frac{11}{22} {Σ Σ}_{i i,, j j}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) (({α α}_{j j} - - {α α}_{j j}^{* *})) K K (({x x}_{i i},, {x x}_{j j})) \\ + + {Σ Σ}_{i i,, j j}^{n no} (({α α}_{i i} - - {α α}_{i i}^{* *})) {y the y}_{i i} - - {Σ Σ}_{i i,, j j}^{n no} (({α α}_{i i} + + {α α}_{i i}^{* *})) ϵ ϵ \end{matrix} - - - - - - ((3535))

The constraint condition is (32). After solving the value of α, f(x) can be obtained:

f f ((x x)) = = {Σ Σ}_{i i = = 11}^{N N} (({α α}_{i i} - - {α α}_{i i}^{* *})) K K (({x x}_{i i},, x x)) + + b b - - - - - - ((3636))

Typically, most α _i or

The value of will be zero, non-zero α _i or

The corresponding samples are called support vectors. According to the optimization conditions, there are:

\{\begin{matrix} {α α}_{i i} [[{ξ ξ}_{i i} + + ϵ ϵ - - {y the y}_{i i} + + f f (({x x}_{i i}))]] = = 00 & ,, i i = = 11,, . . . . . . l l \\ {a a}_{i i}^{* *} [[{ξ ξ}_{i i}^{* *} + + ϵ ϵ - - {y the y}_{i i} + + f f (({x x}_{i i}))]] = = 00 & ,, i i = = 11,, . . . . . . l l \\ ((C C - - {α α}_{i i})) {ξ ξ}_{i i} = = 00 & ,, i i = = 11,, . . . . . . l l \\ ((C C - - {α α}_{i i}^{* *})) {ξ ξ}_{i i}^{* *} = = 00 & ,, i i = = 11,, . . . . . . l l \end{matrix} - - - - - - ((3737))

So the calculation formula of b is as follows:

b = {the y}_{i} - ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),

When α _j ∈ (0,C)

b = {the y}_{i} + ϵ - Σ_{i = 1}^{l} (α_{i} - α_{i}^{*}) K (x_{j}, x_{i}),

when

The value of b can be calculated with any support vector, and the regression function f(x) can be obtained; after the prediction model is established, the bone structure image can be predicted according to the gray value distribution of the lung X-ray image;

Step 5: Perform image subtraction on the normal chest radiograph and the predicted bone image to obtain the predicted soft tissue image.

2. the method for bone suppression in a kind of chest X-ray image according to claim 1, is characterized in that: described step 1.2: the specific steps of alignment training shape vector are as follows:

The specific method of mapping the block data to the feature space is as follows:

Step 1.2.1: Rotate, scale and translate each lung region shape to align with the first shape in the training set;

Step 1.2.2: Calculate the average shape according to the aligned shape;

Step 1.2.3: Rotate, scale and translate the average shape to align with the first shape;

Step 1.2.4: Re-align each shape with the current mean shape;

Step 1.2.5: If the process converges or reaches the specified number of cycles, exit; otherwise go to step.

3. the method for bone suppression in a kind of chest X-ray image according to claim 1, is characterized in that: in the described step 3, the feature image of B-spline multi-scale wavelet transform is extracted as 3rd order B-spline multi-scale wavelet Transformation, when m=3, the filter coefficients are k=-2,-1,0,1,2,

4. The method for bone suppression in a chest X-ray image according to claim 1, characterized in that: in said step 3, a three-scale wavelet decomposition is performed on the chest X-ray image.

5. A method for bone suppression in chest X-ray images according to claim 1, characterized in that: in the step 3.2, the scale of the Gaussian filter is σ=2,4.