CN103824281A - Bone inhibition method in chest X-ray image - Google Patents

Abstract

The invention belongs to the technical field of X-ray image processing, and specifically relates to a bone inhibition method in a chest X-ray image. A lung contour model is established for a normal chest image in a dual-energy silhouette image, three-order B batten wavelet transformation is then applied to perform three-scale wavelet decomposition on the image, a characteristic image is extracted by use of a 2-jet method from the decomposed image, then a bone model is established by use of a support vector regression machine, the bone model is applied to an X-ray image for extracting a corresponding prediction bone image, and finally, image subtraction is carried out on the X-ray image and the predicted bone image so as to obtain a predicted soft tissue image. Because most of bone inhibition methods employ a nerve network method, the method is liable to local optimum and the training result is not yet stable. Therefore, the method provided by the invention employs the support vector regression machine to carry out bone prediction, and a better result is obtained.

Description

A kind of method that bone suppresses in chest x-ray image

Technical field

The invention belongs to X-ray image processing technology field, and in particular to a kind of method that bone suppresses in chest x-ray image. 

Background technology

Lung cancer is one of current malignant tumour maximum to human health risk.Particularly since nearly half a century, the environment constantly deterioration and the substantial increase of the smoking size of population caused with air pollution, the incidence of disease and case fatality rate of various countries' lung cancer are all steeply rising.Because lung belongs to inside of human body internal organs, most lung cancer are simply growing silently in body at first, and patient does not have any sensation.It is most of to be in middle and advanced stage when patient is medical because of clinical symptoms such as cough, spitting of blood and pectoralgias, miss the best opportunity for the treatment of.So, early diagnosis and the treatment of lung cancer are the keys for improving patients with lung cancer survival rate. 

Medical imaging can not only find the presence of focus by image, and positioning focus is in the position of full lung, moreover it is possible to observe the physical features such as the size, form and density of focus.The diagnostic imaging of thoracopathy is based on chest x-ray photography, but in chest x-ray image, often because illness is located at the rib or clavicle structure occlusion area in image and causes to fail to pinpoint a disease in diagnosis. 

Dual Energy Subtraction (Dual Energy Subtraction, DES) is a kind of newer imaging technique developed on the basis of the photography of digital chest x-ray.It is different to the energy attenuation mode of x linear lights from soft tissue using bone, and the photoelectric absorption effect difference of the material of different atomic weight, the information of two kinds of sink effects is separated using digital photography, the dampening information of selective removal bone or soft tissue, and then obtain pure soft tissue picture and bone image.Therefore three images, a normal rabat, a soft-tissue image and a bone image can be obtained using dual energy outline technology. 

DES soft tissue subtraction image is to chest micronodular lesion, and calcium scoring venereal disease becomes, lung marking lesion, is particularly at the clear display rate of lung rib, clavicle overlapping nodular lesions apparently higher than common chest x light images.This, which is mainly, subtracts the pulmonary parenchyma image in shadow zone and eliminates the influence of the overlapping factor such as bone, and for rib, the equitant focus of clavicle is shown clearly, and this is more beneficial for the more patient of costal cartilage calcification. 

Some research institutions both domestic and external suppress problem to bone in chest x light images and studied.Some foreign scholars estimate skeletal structure model using normal chest x light images, produce bone and suppress the detection performance that image carrys out the high lesion of figure.Ahmed and Rasheed et al. propose to suppress rear side rib using independent component analysis (Independent Component Analysis, ICA) method, and then improve the visibility of tubercle.Jiann et al. proposes rib restrainable algorithms in a kind of non-parametric normal rabat.The algorithm is split by lung, rib segmentation, the estimation of rib gray scale and the part such as suppresses and constitutes, and proposes to estimate the gray scale of rib using Real Coding Genetic Algorithm (Real-coded genetic algorithm, RCGA), and then realize the suppression of rib. 

Also some scholars utilize the normal breast image and bone image in DES images, construct regression model, and then reach that generation bone suppresses the purpose of image.Loog et al. proposes the Image filter arithmetic based on nonlinear regression to obtain skeletal structure and soft-tissue image from x light images.The logical regression model that DES lung images and bone image are set up using k arest neighbors regression algorithms of the algorithm.Kenji Suzuki et al. propose the image processing techniques based on multiresolution large-scale training artificial neural network (Multiresolution massive training artificial neural network, MTANN) to suppress the contrast of bony areas in x light images.The algorithm uses the x light images in dual intensity image, as training data, regression model to be set up using MTANN with bone structure image.In the bone suppression image of generation, the anatomical structure such as tubercle and blood vessel is high-visible. 

In summary, just just risen for the research that bone in chest x light images suppresses both at home and abroad.Original most of researchs are set up and removed based on the bone extraction in common chest x light images, segmentation, model.But, because skeletal structure is complicated in common x light images, and fringe region gray scale contrast it is relatively low the features such as, realize that the segmentation and suppression of bone are relatively difficult merely with common x light images.And bone image can be obtained from the angle of a completely lower grain husk using the bone regression model of dual intensity image, partitioning algorithm should not be studied again.During being set up in model, using a large amount of DES bone images information, therefore the prognostic chart picture precision obtained is higher. 

The content of the invention

In view of the deficienciess of the prior art, the present invention provides the method that bone suppresses in a kind of high chest x-ray image of stability. 

The present invention takes to be carried out as follows： 

Step 1：The determination of model lung initial profile position； 

Step 1.1：Mark the boundary point on every image lung border in training set； 

Step 1.2：Training shapes of aliging vector；Training shapes is alignd by the operation of scaling, rotation and translation, them is alignd as far as possible closely, if x_iIt is the vector of n point in i-th of shape in training set, x_i=(x_i0,y_i0,...x_ik,y_ik,...,x_in-1,y_in-1)^T.Wherein, (x_ik,y_ik) it is k-th point in i-th of shape, to give two similar shape x_iAnd x_j, selection anglec of rotation θ, scaling s, translation (t_x,t_y), then M (s, θ) [x] represents the conversion that the anglec of rotation is θ and scaling is s, makes following weighted sum minimum by x_iIt is mapped as x_j： 

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

Wherein, $M(s,\mathrm{\θ})\left[\begin{array}{c}{x}_{\mathrm{jk}}\\ y_{\mathrm{jk}}\end{array}\right]=\left(\begin{array}{c}\left(s\cos \mathrm{\θ}\right)x_{\mathrm{jk}}-\left(s\sin \mathrm{\θ}\right)y_{\mathrm{jk}}\\ \left(s\sin \mathrm{\θ}\right)x_{\mathrm{jk}}+\left(s\cos \mathrm{\θ}\right)y_{\mathrm{jk}}\end{array}\right),$ t=(t_x,t_y,...,t_x,t_y)^T.In formula：W is the weighting diagonal matrix of a point, it by each boundary point weight w_kConstitute.Make R_klRepresent the distance of and at k-th point at l-th point,

Represent the weights of the change, then of distance in training set at k-th point

If making a_x=scosθ,a_y=ssin θ, then have four linear equalities according to least square method

$\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)$

Wherein, $X_i=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k{x}_{\mathrm{ik}},$ $Y_i=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k{y}_{\mathrm{ik}},$ $Z=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k(x_{2k}^2+y_{2k}^2),$ $C_1=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k(x_{1k}{x}_{2k}+y_{1k}{y}_{2k}),$ ${\mathrm{<figure-callout\; id="2"\; label="C"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">C</figure-callout>}}_2=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k(y_{1k}{x}_{2k}+x_{1k}{y}_{2k}),$ $w=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k.$ Matrix operation method can be used to solve variable a_x,a_y,t_x,t_y。 

Step 1.3：Set up the model of lung's initial profile position；After training shapes vector alignment, the statistical information of change in shape is found out using principal component analytical method, model is set up accordingly； 

Step 2：Split with reference to the pulmonary parenchyma of half-tone information and shape information：The gray scale and shape information of boundary point in multiple characteristic images are utilized simultaneously so that the boundary intensity that searches, shape information are similar to training image； 

Step 2.1：Extract characteristic image；Carry out principal component analysis and dimensionality reduction carried out to data, the gray scale for the point p+ △ p near image I midpoint p can be deployed to obtain by Taylor, $I(p+\mathrm{\&Delta;p})\&ap;\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{n!}(\mathrm{\&Delta;}{p}_{i1},\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Delta;}{p}_{\mathrm{in}}\frac{{\&PartialD;}^{n}I\left(p\right)}{{\&PartialD;i}_{1\&CenterDot;\&CenterDot;\&CenterDot;}{\&PartialD;i}_n})---\left(2\right)$

It can thus be concluded that, image I second order Taylor is expanded into： 

$\begin{array}{c}I(p+\mathrm{\&Delta;p})=\underset{n=0}{\overset{2}{\mathrm{\&Sigma;}}}\frac{1}{n!}(\mathrm{\&Delta;}{p}_{i1},\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Delta;}{p}_{\mathrm{in}}\frac{{\&PartialD;}^{n}I\left(p\right)}{{\&PartialD;i}_1\&CenterDot;\&CenterDot;\&CenterDot;{\&PartialD;i}_n})=I\left(p\right)+\mathrm{\&Delta;}{p}_{i1}\frac{\&PartialD;I\left(p\right)}{{\&PartialD;i}_1}+\frac{1}{2!}\left(\mathrm{\&Delta;}{p}_{i1}\mathrm{\&Delta;}{p}_{i2}\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;i}_1{\&PartialD;i}_2}\right)\\ =I\left(p\right)+\mathrm{\&Delta;}{p}_x\frac{\&PartialD;I\left(p\right)}{\&PartialD;x}+\frac{1}{2!}(\mathrm{\&Delta;}{p}_x^2\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;x}^2}+\mathrm{\&Delta;}{p}_x\mathrm{\&Delta;}{p}_y\frac{{\&PartialD;}^{2}I\left(p\right)}{\&PartialD;x\&PartialD;y})\\ +\mathrm{\&Delta;}{p}_y\frac{\&PartialD;I\left(p\right)}{\&PartialD;y}+\frac{1}{2!}(\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="p\; y"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">p</figure-callout>}}_y^{}\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">y</figure-callout>}}^2}+\mathrm{\&Delta;}{p}_y\mathrm{\&Delta;}{p}_x\frac{{\&PartialD;}^{2}I\left(p\right)}{\&PartialD;y\&PartialD;x})\end{array}---\left(3\right)$

Gaussian smoothing is carried out to image, it is A (p, σ)=(I*G) (p, σ) that the image after gaussian filtering is carried out for image I (p).The n rank partial derivatives of filtering image areThen it can be seen from Convolution Properties： 

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

In formula：Footmark i₁,...i_nRepresent that data seek partial derivative to variable respectively, Gauss partial derivative is regarded as a wave filter, the combination of Gauss partial derivative response constitutes a complete feature description；N-jet features for a given scale parameter σ can be defined as： 

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

It can thus be seen that as N → ∞, in a certain metric space, image I local feature is made up of the combination of all partial derivatives.It is therefore possible to use Gauss partial derivative constructs a characteristic vector to describe local feature.When partial derivative exponent number is higher, the description to image is more accurate.But, if describing image using excessive partial derivative, feature vector dimension can be increased. 

The candidate point of boundary point and the gray scale cost of each candidate point are obtained using local 2-jet characteristic images, local 2-jet is characterized as： 

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

Step 2.2：Choose the candidate point of boundary point：For each point on initial lung border, the similarity degree of the gray scale of all pixels point and respective point gray scale in training characteristics image in the point search region in all characteristic images is calculated, the maximum point of similarity degree is selected, is used as the candidate point of the boundary point.Similarity degree is the mahalanobis distance h of surrounding pixel point gray scale respective point surrounding pixel point gray scale set into training sample characteristic image in all characteristic images_i： 

Step 2.3：Lung region segmentation is carried out using Dynamic Programming；In boundary point region of search, the grey similarity cost of pixel is the similarity degree h of the surrounding pixel point gray scale of surrounding pixel point gray scale boundary point corresponding to training image_iIn boundary point region of search, the grey similarity cost of pixel is the similarity degree h of the surrounding pixel point gray scale of surrounding pixel point gray scale boundary point corresponding to training image_i；The shape similarity cost of i-th of boundary point is defined as in image： 

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

In formula, v_i=p_i+1-p_i, the shape facility of i-th of boundary point is represented,

The average and covariance of i-th of boundary point shape facility in all training images are represented respectively. 

For i-th of boundary point p in test image_i, having m in specified region of search has smaller grey similarity cost candidate point, then n boundary point will produce n × m gray scale cost matrix： 

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

$h_i=\underset{j=1}{\overset{N}{\mathrm{\&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}(k-1)\right)\\ \sin \left(\frac{2\mathrm{\&pi;}}{n_c}(k-1)\right)\end{array}\right]---\left(11\right)$

In formula：

For on characteristic image, with point p_iFor the center of circle, r_cFor the n on the circle of radius_cThe gray scale of individual point, k=1......n_c。 The average and covariance of the surrounding pixel point gray scale of i-th of boundary point respectively in j-th of characteristic image of training image.N is characterized total number of images,

It is to look for an optimal path to search for optimal profile, and often row selects an element in Matrix C, when along Path selection, and the summation of gray scale and shape similarity cost is minimum, i.e.,： 

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

Wherein, γ is gray scale and shape similarity cost coefficient.Adjust γ values so that both costs play roughly the same effect during boundary search； 

Step 3：Extracted using the characteristic image of B-spline multi-scale wavelet transformation；During model foundation, the characteristic image for extracting effectively description different scale bone is the basis that model is set up； 

Step 3.1：B-spline multi-scale wavelet transformation：B-spline function quickly converges on Gaussian function with the increase of batten exponent number, and its first derivative can approach Optimal edge detection operator；Multi-scale edge enhancing is carried out using B-spline small echo；M rank B-spline basic function N (m) may be defined as： 

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

Wherein, $\left(\begin{array}{c}m\\ t\end{array}\right)=\frac{t!}{m!(t-m)!}.$ It follows that zero order B-spline is： 

As m >=1, m rank B-spline functions N_m(x) it can be represented with convolution iterative relation： 

Obviously, N_m(x) it is non-negative and compact schemes, its support width is suppN_m=[0, m+1], support center is (m+1)/2, and N_mOn the Central Symmetry of (m+1)/2. 

The Fourier transformation of zeroth order B-spline isThen the Fourier formalism of m B-spline is： 

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

If by m B-spline base N_m(x) N is moved to_m[x+ (m+1)/2], then symmetrical centre moves to the origin of coordinates, the linear phase that Fourier transformation will no longer occur in above formula.But, when m is even number, it will appear from half-integer node.To avoid the occurrence of half-integer node, when m is even number, by N_m(x) N is moved to_m[x+ (m+1)/2], then symmetrical centre moves to 1/2.By N_m(x) function after integer translation is designated as θ_m(x), its Fourier transformation is： 

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

Two scaling relation formulas are obtained by multiresolution analysis

It can try to achieve： 

$H\left(e^{\mathrm{i\&omega;}}\right)=\mathrm{\&Sigma;}{h}_k{e}^{-\mathrm{ik\&omega;}}=\frac{{\widehat{\mathrm{\&theta;}}}_m\left(2\mathrm{\&omega;}\right)}{{\widehat{\mathrm{\&theta;}}}_m\left(\mathrm{\&omega;}\right)}=\frac{e^{-\mathrm{i\&epsiv;\&omega;}}{\left(\frac{\sin \mathrm{\&omega;}}{\mathrm{\&omega;}}\right)}^{m+1}}{e^{-i\frac{\mathrm{\&epsiv;\&omega;}}{2}}{\left(\frac{\sin (\mathrm{\&omega;}/2)}{\mathrm{\&omega;}/2}\right)}^{m+1}}=e^{-i\frac{\mathrm{\&epsiv;\&omega;}}{2}}{\left(\cos \frac{\mathrm{\&omega;}}{2}\right)}^{m+1}---\left(18\right)$

In formula, when m is odd number, ε=0；When m is even number, ε=1.Can obtain low pass filter by Euler's formula and binomial expansion is： 

When m is odd number,

When m is even number,

The first derivative of B-spline function is chosen as wavelet function, i.e., m+1 times B-spline function is 2^-1First differential function is m B-spline wavelet function on yardstick： 

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

Its Fourier transformation is

And ψ_m(x) the admissibility condition of small echo is met.By two scaling relations of multiresolution analysis Wavelets and scaling function

It can try to achieve： 

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

When m is odd number, $G\left(\mathrm{\&omega;}\right)={4\mathrm{ie}}^{-\frac{\mathrm{i\&omega;}}{2}}\sin \frac{\mathrm{\&omega;}}{2};$ When m is even number, $G\left(\mathrm{\&omega;}\right)={4\mathrm{ie}}^{\frac{\mathrm{i\&omega;}}{2}}\sin \frac{\mathrm{\&omega;}}{2}.$ It follows that time-domain finite impulse response is： 

As m=3, every filter coefficient is k=- 2, -1,0,1,2,g_k={0,0,-2,2,0}.As can be seen here, the finite length of the wave filter of three rank B-spline small echos is 5. 

Ask for B-spline wavelet transform filter coefficient： 

In formula：h_kFor low-pass filter coefficients, g_kFor high-pass filter coefficient, gained h is utilized_kAnd g_kCarry out Wavelet fast decomposition to image, the row of original image respectively with h_kAnd g_kConvolution is carried out, then arranges it carry out down-sampling, two width subgraphs can be obtained, two width subgraphs are filtered then along the direction of row and down-sampling is carried out, it is possible to the output subgraph of four 1/4 sizes, diagonal detail pictures is obtained, vertical detail image, level detail image and approximate image； 

Step 3.2：Gaussian filtering part 2-jet characteristic images are extracted, and because the characteristic vector that Gauss partial derivative is constructed can be used for describing the local feature of image, extract 2-jet features (I, the I of the detail pictures after different scale wavelet transformation_x,I_y,I_xx,I_xy,I_yy) figure sets up regression model.According in formula (7) formula, the yardstick of gaussian filtering is σ=2,4； 

Step is weighed：The prediction of bone image is carried out using support vector regression, by trying to achieve regression function f (x)=ω x+b, skeleton model is set up：If sample set is (x_i,y_i), i=1,2 ..., l, x_i∈Rⁿ, sample set S is the linear approximation of ε-insensitive loss function.As linear regression function f (x)=ω x+b fittings (x_i,y_i) when, it is assumed that the error of fitting precision of all training datas is ε, i.e.,： 

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

Allow d_iRepresent point (x_i,y_i) ∈ S to hyperplane f (x) distance： 

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

Because S collection is ε-linear approximation, have |<w,x>+b-f(xi)|≤ε,i=1,...,m.It can obtain: 

$\frac{|＜w,x>+b-y|}{\sqrt{1+{\left|\right|w\left|\right|}^2}}＜\frac{\mathrm{\&epsiv;}}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,...,l---\left(24\right)$

Then have

$d_i＜\frac{\mathrm{\&epsiv;}}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,l---\left(25\right)$

(23) formula shows

It is point in S to the upper bound of the distance of hyperplane.The linear approximation collection S of ε-insensitive loss function best fit approximation hyperplane is to obtain hyperplane by maximizing the point in S to the upper bound of hyperplane distance； 

Optimal hyperlane is by maximizing

What is obtained (minimizes

).Therefore, as long as minimizing | | w | |²It can be obtained by best fit approximation hyperplane.Then linear regression problem is just turned to： 

Seek following optimization problem： 

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

It is constrained to |<w,x>+b-y_i|≤ε,i=1,...,m.Additionally, it is contemplated that the situation of error of fitting, introduces relaxation factor

When division has error,

Both greater than 0, error is not present and takes 0； 

Thus, the insensitive support vector regression of criterion epsilon can be expressed as： 

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

Wherein, C>0 is balance factor.Introducing Lagrangian is： 

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

Wherein, α_i,

γ i >=0, i=1 ..., n.Function L extreme value should meet condition

Then following formula is obtained： 

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

$\underset{i=1}{\overset{n}{\mathrm{\&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) are updated in (28), the dual form for obtaining optimization problem is： 

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

It is constrained to： 

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

Because bone gray scale is predicted as nonlinear regression, so data are first mapped to a high-dimensional feature space using a Nonlinear Mapping φ, then carrying out linear regression in high-dimensional feature space sets up model.Replaced due in superincumbent optimization process, only taking into account the inner product operation in high-dimensional feature space, therefore with a kernel function k (x, y)<φ(x),φ(y)>Nonlinear regression can just be realized.Then, the optimization method of nonlinear regression is the following function of maximization： 

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

Its constraints is (32).After the value for solving α, it is possible to obtain f (x)： 

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

Under normal circumstances, most of α_iOr

Value will be zero, the α being not zero_iOrCorresponding sample is referred to as supporting vector.According to optimal condition (KKT conditions), have in saddle point： 

$\left\{\begin{array}{cc}{\mathrm{\&alpha;}}_i[{\mathrm{\&xi;}}_i+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)]=0& ,i=1,...l\\ a_i^{*}[{\mathrm{\&xi;}}_i^{*}+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)]=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)$

Then the calculating formula for obtaining b is as follows： 

$b=y_i-\mathrm{\&epsiv;}-\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x_j,x_i),$ Work as α_j∈(0,C) 

$b=y_i+\mathrm{\&epsiv;}-\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x_j,x_i),$ When

∈(0,C) 

B value can be calculated with any one supporting vector, regression function f (x) is obtained；Set up after forecast model, generation bone structure image can be predicted according to the grey value profile of lung's x light images； 

Step 5：The soft-tissue image for doing image subtraction to normal rabat and the bone image predicted to predict. 

The step 1.2：Alignment training shapes are vectorial to be comprised the following steps that： 

The described specific method that the data of piecemeal are mapped into feature space is as follows： 

Step 1.2.1：Rotation, each lung region shape of zooming and panning, make it be alignd with first shape in training set； 

Step 1.2.2：According to alignment shape, average shape is calculated； 

Step 1.2.3：Rotation, zooming and panning average shape make it be alignd with first shape； 

Step 1.2. is weighed：Again each shape is alignd with current average shape； 

Step 1.2.5：If process restrains or to designated cycle number of times, exit；Otherwise step is gone to. 

The decomposition of three multi-scale wavelets is done in the step 3 to chest x light images. 

Advantages of the present invention is：The present invention first sets up lung outlines model to the normal rabat in dual intensity sketch figure picture, reapply three rank B-spline wavelet transformations and the decomposition of three multi-scale wavelets is done to image, image after decomposition is extracted into characteristic image using 2-jet methods, reuse support vector regression and set up skeleton model, skeleton model is applied in x-ray image, so as to extract corresponding prediction bone image, image subtraction finally is done so as to the soft-tissue image predicted with the bone image of x-ray image and prediction.Because suppressing bone method mostly uses neural net method, and the method is easily absorbed in local optimum, training result less stable, so this patent carries out bone prediction using support vector regression, has obtained more excellent result. 

Brief description of the drawings

Fig. 1 is the overview flow chart of the present invention； 

Flow chart is cut in the lung differentiation that Fig. 2 is the present invention； 

Fig. 3 is gaussian filtering part 2-jet extraction characteristic image schematic diagrames in the present invention； 

Fig. 4 is the schematic diagram of best fit approximation hyperplane in the present invention； 

Fig. 5 is support vector regression schematic diagram in the present invention. 

Embodiment

As shown in Figures 1 to 5, the present invention takes is carried out as follows： 

Step 1：The determination of model lung initial profile position； 

Step 1.1：Mark the boundary point on every image lung border in training set； 

Step 1.2：Training shapes of aliging vector；Training shapes is alignd by the operation of scaling, rotation and translation, them is alignd as far as possible closely, if x_iIt is the vector of n point in i-th of shape in training set, x_i=(x_i0,y_i0,...x_ik,y_ik,...,x_in-1,y_in-1)^T.Wherein, (x_ik,y_ik) it is k-th point in i-th of shape, to give two similar shape x_iAnd x_j, selection anglec of rotation θ, scaling s, translation (t_x,t_y), then M (s, θ) [x] represents the conversion that the anglec of rotation is θ and scaling is s, makes following weighted sum minimum by x_iIt is mapped as x_j： 

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

Wherein, $M(s,\mathrm{\&theta;})\left[\begin{array}{c}{x}_{\mathrm{jk}}\\ y_{\mathrm{jk}}\end{array}\right]=\left(\begin{array}{c}\left(s\cos \mathrm{\&theta;}\right)x_{\mathrm{jk}}-\left(s\sin \mathrm{\&theta;}\right)y_{\mathrm{jk}}\\ \left(s\sin \mathrm{\&theta;}\right)x_{\mathrm{jk}}+\left(s\cos \mathrm{\&theta;}\right)y_{\mathrm{jk}}\end{array}\right),$ t=(t_x,t_y,...,t_x,t_y)^TIn formula：W is the weighting diagonal matrix of a point, it by each boundary point weight w_kConstitute.Make R_klRepresent the distance of and at k-th point at l-th point,Represent the weights of the change, then of distance in training set at k-th point

If making a_x=scosθ,a_y=ssin θ, then have four linear equalities according to least square method

$\left(\begin{array}{cccc}{X}_2& -Y_2& W& 0\\ Y_2& {\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">X</figure-callout>}}_2& 0& W\\ Z& 0& X_2& {\mathrm{<figure-callout\; id="2"\; label="Y"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">Y</figure-callout>}}_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)$

Wherein, $X_i=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k{x}_{\mathrm{ik}},$ $Y_i=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k{y}_{\mathrm{ik}},$ $Z=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k(x_{2\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">k</figure-callout>}}^2+y_{2k}^2),$ $C_1=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k(x_{1k}{x}_{2k}+y_{1k}{y}_{2k}),$ ${\mathrm{<figure-callout\; id="2"\; label="C"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">C</figure-callout>}}_2=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k(y_{1k}{x}_{2k}+x_{1k}{y}_{2k}),$ $w=\underset{k=0}{\overset{n-1}{\mathrm{\&Sigma;}}}{w}_k.$ Matrix operation method can be used to solve variable a_x,a_y,t_x,t_y。 

Step 1.3：Set up the model of lung's initial profile position；After training shapes vector alignment, the statistical information of change in shape is found out using principal component analytical method, model is set up accordingly； 

Step 2：Split with reference to the pulmonary parenchyma of half-tone information and shape information：The gray scale and shape information of boundary point in multiple characteristic images are utilized simultaneously so that the boundary intensity that searches, shape information are similar to training image； 

Step 2.1：Extract characteristic image；Carry out principal component analysis and dimensionality reduction carried out to data, the gray scale for the point p+ △ p near image I midpoint p can be deployed to obtain by Taylor, $I(p+\mathrm{\&Delta;p})\&ap;\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{n!}(\mathrm{\&Delta;}{p}_{i1},\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Delta;}{p}_{\mathrm{in}}\frac{{\&PartialD;}^{n}I\left(p\right)}{{\&PartialD;i}_{1\&CenterDot;\&CenterDot;\&CenterDot;}{\&PartialD;i}_n})$

                                                                       (2) 

It can thus be concluded that, image I second order Taylor is expanded into： 

$\begin{array}{c}I(p+\mathrm{\&Delta;p})=\underset{n=0}{\overset{2}{\mathrm{\&Sigma;}}}\frac{1}{n!}(\mathrm{\&Delta;}{p}_{i1},\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Delta;}{p}_{\mathrm{in}}\frac{{\&PartialD;}^{n}I\left(p\right)}{{\&PartialD;i}_1\&CenterDot;\&CenterDot;\&CenterDot;{\&PartialD;i}_n})=I\left(p\right)+\mathrm{\&Delta;}{p}_{i1}\frac{\&PartialD;I\left(p\right)}{{\&PartialD;i}_1}+\frac{1}{2!}\left(\mathrm{\&Delta;}{p}_{i1}\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="p\; i"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">p</figure-callout>}}_i\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;i}_1{\&PartialD;i}_2}\right)\\ =I\left(p\right)+\mathrm{\&Delta;}{p}_x\frac{\&PartialD;I\left(p\right)}{\&PartialD;x}+\frac{1}{2!}(\mathrm{\&Delta;}{p}_x^2\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;x}^2}+\mathrm{\&Delta;}{p}_x\mathrm{\&Delta;}{p}_y\frac{{\&PartialD;}^{2}I\left(p\right)}{\&PartialD;x\&PartialD;y})\\ +\mathrm{\&Delta;}{p}_y\frac{\&PartialD;I\left(p\right)}{\&PartialD;y}+\frac{1}{2!}(\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="p\; y"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">p</figure-callout>}}_y^{}\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HDA0000453943940000022.png"\; state="\{\{state\}\}">y</figure-callout>}}^2}+\mathrm{\&Delta;}{p}_y\mathrm{\&Delta;}{p}_x\frac{{\&PartialD;}^{2}I\left(p\right)}{\&PartialD;y\&PartialD;x})\end{array}---\left(3\right)$

Gaussian smoothing is carried out to image, it is A (p, σ)=(I*G) (p, σ) that the image after gaussian filtering is carried out for image I (p).The n rank partial derivatives of filtering image are

Then it can be seen from Convolution Properties： 

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

In formula：Footmark i₁,...i_nRepresent that data seek partial derivative to variable respectively, Gauss partial derivative is regarded as a wave filter, the combination of Gauss partial derivative response constitutes a complete feature description；N-jet features for a given scale parameter σ can be defined as： 

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

It can thus be seen that as N → ∞, in a certain metric space, image I local feature is made up of the combination of all partial derivatives.It is therefore possible to use Gauss partial derivative constructs a characteristic vector to describe local feature.When partial derivative exponent number is higher, the description to image is more accurate.But, if describing image using excessive partial derivative, feature vector dimension can be increased. 

The candidate point of boundary point and the gray scale cost of each candidate point are obtained using local 2-jet characteristic images, local 2-jet is characterized as： 

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

Step 2.2：Choose the candidate point of boundary point：For each point on initial lung border, the similarity degree of the gray scale of all pixels point and respective point gray scale in training characteristics image in the point search region in all characteristic images is calculated, the maximum point of similarity degree is selected, is used as the candidate point of the boundary point.Similarity degree is the mahalanobis distance h of surrounding pixel point gray scale respective point surrounding pixel point gray scale set into training sample characteristic image in all characteristic images_i： 

Step 2.3：Lung region segmentation is carried out using Dynamic Programming；In boundary point region of search, the grey similarity cost of pixel is the similarity degree h of the surrounding pixel point gray scale of surrounding pixel point gray scale boundary point corresponding to training image_iIn boundary point region of search, the grey similarity cost of pixel is the similarity degree h of the surrounding pixel point gray scale of surrounding pixel point gray scale boundary point corresponding to training image_i；The shape similarity cost of i-th of boundary point is defined as in image： 

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

In formula, v_i=p_i+1-p_i, the shape facility of i-th of boundary point is represented,

The average and covariance of i-th of boundary point shape facility in all training images are represented respectively. 

For i-th of boundary point p in test image_i, having m in specified region of search has smaller grey similarity cost candidate point, then n boundary point will produce n × m gray scale cost matrix： 

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

$h_i=\underset{j=1}{\overset{N}{\mathrm{\&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}(k-1)\right)\\ \sin \left(\frac{2\mathrm{\&pi;}}{n_c}(k-1)\right)\end{array}\right]---\left(11\right)$

In formula：For on characteristic image, with point p_iFor the center of circle, r_cFor the n on the circle of radius_cThe gray scale of individual point, k=1......n_c。 

The average and covariance of the surrounding pixel point gray scale of i-th of boundary point respectively in j-th of characteristic image of training image.N is characterized total number of images,

It is to look for an optimal path to search for optimal profile, and often row selects an element in Matrix C, when along Path selection, and the summation of gray scale and shape similarity cost is minimum, i.e.,： 

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

Wherein, γ is gray scale and shape similarity cost coefficient.Adjust γ values so that both costs play roughly the same effect during boundary search； 

Step 3：Extracted using the characteristic image of B-spline multi-scale wavelet transformation；During model foundation, the characteristic image for extracting effectively description different scale bone is the basis that model is set up； 

Step 3.1：B-spline multi-scale wavelet transformation：B-spline function quickly converges on Gaussian function with the increase of batten exponent number, and its first derivative can approach Optimal edge detection operator；Multi-scale edge enhancing is carried out using B-spline small echo；M rank B-spline basic function N (m) may be defined as： 

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

Wherein, $\left(\begin{array}{c}m\\ t\end{array}\right)=\frac{t!}{m!(t-m)!}.$ It follows that zero order B-spline is： 

As m >=1, m rank B-spline functions N_m(x) it can be represented with convolution iterative relation： 

Obviously, N_m(x) it is non-negative and compact schemes, its support width is suppN_m=[0, m+1], support center is (m+1)/2, and N_mOn the Central Symmetry of (m+1)/2. 

The Fourier transformation of zero order B-spline isThen the Fourier formalism of m B-spline is： 

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

If by m B-spline base N_m(x) N is moved to_m[x+ (m+1)/2], then symmetrical centre moves to the origin of coordinates, the linear phase that Fourier transformation will no longer occur in above formula.But, when m is even number, it will appear from half-integer node.To avoid the occurrence of half-integer node, when m is even number, by N_m(x) N is moved to_m[x+ (m+1)/2], then symmetrical centre moves to 1/2.By N_m(x) function after integer translation is designated as θ_m(x), its Fourier transformation is： 

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

Two scaling relation formulas are obtained by multiresolution analysisIt can try to achieve： 

$H\left(e^{\mathrm{i\&omega;}}\right)=\mathrm{\&Sigma;}{h}_k{e}^{-\mathrm{ik\&omega;}}=\frac{{\widehat{\mathrm{\&theta;}}}_m\left(2\mathrm{\&omega;}\right)}{{\widehat{\mathrm{\&theta;}}}_m\left(\mathrm{\&omega;}\right)}=\frac{e^{-\mathrm{i\&epsiv;\&omega;}}{\left(\frac{\sin \mathrm{\&omega;}}{\mathrm{\&omega;}}\right)}^{m+1}}{e^{-i\frac{\mathrm{\&epsiv;\&omega;}}{2}}{\left(\frac{\sin (\mathrm{\&omega;}/2)}{\mathrm{\&omega;}/2}\right)}^{m+1}}=e^{-i\frac{\mathrm{\&epsiv;\&omega;}}{2}}{\left(\cos \frac{\mathrm{\&omega;}}{2}\right)}^{m+1}---\left(18\right)$

In formula, when m is odd number, ε=0；When m is even number, ε=1.Can obtain low pass filter by Euler's formula and binomial expansion is： 

When m is odd number,

When m is even number,

The first derivative of B-spline function is chosen as wavelet function, i.e., m+1 times B-spline function is 2^-1First differential function is m B-spline wavelet function on yardstick： 

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

Its Fourier transformation is

And ψ_m(x) the admissibility condition of small echo is met.By two scaling relations of multiresolution analysis Wavelets and scaling function

It can try to achieve： 

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

When m is odd number, $G\left(\mathrm{\&omega;}\right)={4\mathrm{ie}}^{-\frac{\mathrm{i\&omega;}}{2}}\sin \frac{\mathrm{\&omega;}}{2};$ When m is even number, $G\left(\mathrm{\&omega;}\right)={4\mathrm{ie}}^{\frac{\mathrm{i\&omega;}}{2}}\sin \frac{\mathrm{\&omega;}}{2}.$ It follows that time-domain finite impulse response is： 

As m=3, every filter coefficient is k=- 2, -1,0,1,2,

g_k={0,0,-2,2,0}.As can be seen here, the finite length of the wave filter of three rank B-spline small echos is 5. 

Ask for 3 rank B-spline wavelet transform filter coefficients： 

In formula：h_kFor low-pass filter coefficients, g_kFor high-pass filter coefficient, gained h is utilized_kAnd g_kCarry out Wavelet fast decomposition to image, the row of original image respectively with h_kAnd g_kConvolution is carried out, then arranges it carry out down-sampling, two width subgraphs can be obtained, two width subgraphs are filtered then along the direction of row and down-sampling is carried out, it is possible to the output subgraph of four 1/4 sizes, diagonal detail pictures is obtained, vertical detail image, level detail image and approximate image； 

The present embodiment does the decomposition of three multi-scale wavelets to chest x light images, 12 Zhang Xiaobo exploded view pictures is obtained, wherein 3 approximate images, 3 level detail images, 3 vertical detail images, 3 diagonal detail pictures.It is few that the skeletal structures that diagonal is distributed are presented due in chest x light images.Therefore, the feature of diagonal detail pictures is not extracted. 

Step 3.2：Gaussian filtering part 2-jet characteristic images are extracted, and because the characteristic vector that Gauss partial derivative is constructed can be used for describing the local feature of image, extract 2-jet features (I, the I of the detail pictures after different scale wavelet transformation_x,I_y,I_xx,I_xy,I_yy) figure sets up regression model.According in formula (7) formula, the yardstick of gaussian filtering is σ=2,4； 

Step is weighed：The prediction of bone image is carried out using support vector regression, by trying to achieve regression function f (x)=ω x+b, skeleton model is set up；If sample set is (x_i,y_i), i=1,2 ..., l, x_i∈Rⁿ, sample set S is the linear approximation of ε-insensitive loss function.As linear regression function f (x)=ω x+b fittings (x_i,y_i) when, it is assumed that the error of fitting precision of all training datas is ε, i.e.,： 

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

Allow d_iRepresent point (x_i,y_i) ∈ S to hyperplane f (x) distance： 

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

Because S collection is ε-linear approximation, have |<w,x>+b-f(x_i)|≤ε,i=1,...,m.It can obtain: 

$\frac{|＜w,x>+b-y|}{\sqrt{1+{\left|\right|w\left|\right|}^2}}＜\frac{\mathrm{\&epsiv;}}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,...,l---\left(24\right)$

Then have

$d_i＜\frac{\mathrm{\&epsiv;}}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,l---\left(25\right)$

(23) formula shows

It is point in S to the upper bound of the distance of hyperplane.The linear approximation collection S of ε-insensitive loss function best fit approximation hyperplane is to obtain hyperplane by maximizing the point in S to the upper bound of hyperplane distance； 

Optimal hyperlane is by maximizing

What is obtained (minimizes

).Therefore, as long as minimizing | | w | |²It can be obtained by best fit approximation hyperplane.Then linear regression problem is just turned to： 

Seek following optimization problem： 

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

It is constrained to |<w,x>+b-y_i|≤ε,i=1,...,m.Additionally, it is contemplated that the situation of error of fitting, introduces relaxation factorWhen division has error,

Both greater than 0, error is not present and takes 0； 

Thus, the insensitive support vector regression of criterion epsilon can be expressed as： 

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

Wherein, C>0 is balance factor.Introducing Lagrangian is： 

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

Wherein, α_i,

γ i >=0, i=1 ..., n.Function L extreme value should meet condition

Then following formula is obtained： 

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

$\underset{i=1}{\overset{n}{\mathrm{\&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) are updated in (28), the dual form for obtaining optimization problem is： 

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

It is constrained to： 

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

Because bone gray scale is predicted as nonlinear regression, so data are first mapped to a high-dimensional feature space using a Nonlinear Mapping φ, then carrying out linear regression in high-dimensional feature space sets up model.Replaced due in superincumbent optimization process, only taking into account the inner product operation in high-dimensional feature space, therefore with a kernel function k (x, y)<φ(x),φ(y)>Nonlinear regression can just be realized.Then, the optimization method of nonlinear regression is the following function of maximization： 

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

Its constraints is (32).After the value for solving α, it is possible to obtain f (x)： 

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

Under normal circumstances, most of α_iOr

Value will be zero, the α being not zero_iOrCorresponding sample is referred to as supporting vector.According to optimal condition (KKT conditions), have in saddle point： 

$\left\{\begin{array}{cc}{\mathrm{\&alpha;}}_i[{\mathrm{\&xi;}}_i+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)]=0& ,i=1,...l\\ a_i^{*}[{\mathrm{\&xi;}}_i^{*}+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)]=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)$

Then the calculating formula for obtaining b is as follows： 

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

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

B value can be calculated with any one supporting vector, regression function f (x) is obtained；Set up after forecast model, generation bone structure image can be predicted according to the grey value profile of lung's x light images； 

Step 5：The soft-tissue image for doing image subtraction to normal rabat and the bone image predicted to predict. 

The step 1.2：Alignment training shapes are vectorial to be comprised the following steps that： 

The described specific method that the data of piecemeal are mapped into feature space is as follows： 

Step 1.2.1：Rotation, each lung region shape of zooming and panning, make it be alignd with first shape in training set； 

Step 1.2.2：According to alignment shape, average shape is calculated； 

Step 1.2.3：Rotation, zooming and panning average shape make it be alignd with first shape； 

Step 1.2. is weighed：Again each shape is alignd with current average shape； 

Step 1.2.5：If process restrains or to designated cycle number of times, exit；Otherwise step is gone to. 

The decomposition of three multi-scale wavelets is done in the step 3 to chest x light images. 

Claims (5)

1. a kind of method that bone suppresses in chest x-ray image, is comprised the following steps that：

Step 1：The determination of model lung initial profile position；

Step 1.1：Mark the boundary point on every image lung border in training set；

Step 1.2：Training shapes of aliging vector；Training shapes is alignd by the operation of scaling, rotation and translation, them is alignd as far as possible closely, if x_iIt is the vector of n point in i-th of shape in training set, x_i=(x_i0,y_i0,...x_ik,y_ik,...,x_in-1,y_in-1)^T, wherein, (x_ik,y_ik) it is k-th point in i-th of shape, to give two similar shape x_iAnd x_j, selection anglec of rotation θ, scaling s, translation (t_x,t_y), then M (s, θ) [x] represents the conversion that the anglec of rotation is θ and scaling is s, makes following weighted sum minimum by x_iIt is mapped as x_j：

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

Wherein,

t=(t_x,t_y,...,t_x,t_y)^T, in formula:W is the weighting diagonal matrix of a point, it by each boundary point weight w_kConstitute.Make R_klRepresent the distance of and at k-th point at l-th point,Represent the weights of the change, then of distance in training set at k-th point

Step 1.3：Set up the model of lung's initial profile position；After training shapes vector alignment, the statistical information of change in shape is found out using principal component analytical method, model is set up accordingly,

It is characterized in that：Also comprise the following steps：

Step 2：Split with reference to the pulmonary parenchyma of half-tone information and shape information：The gray scale and shape information of boundary point in multiple characteristic images are utilized simultaneously so that the boundary intensity that searches, shape information are similar to training image；

Step 2.1：Extract characteristic image；Carry out principal component analysis and dimensionality reduction carried out to data, the gray scale for the point p+ △ p near image I midpoint p can be deployed to obtain by Taylor, $I(p+\mathrm{\&Delta;p})\&ap;\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}\frac{1}{n!}(\mathrm{\&Delta;}{p}_{i1},\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Delta;}{p}_{\mathrm{in}}\frac{{\&PartialD;}^{n}I\left(p\right)}{{\&PartialD;i}_{1\&CenterDot;\&CenterDot;\&CenterDot;}{\&PartialD;i}_n})---\left(2\right)$

It can thus be concluded that, image I second order Taylor is expanded into：

$\begin{array}{c}I(p+\mathrm{\&Delta;p})=\underset{n=0}{\overset{2}{\mathrm{\&Sigma;}}}\frac{1}{n!}(\mathrm{\&Delta;}{p}_{i1},\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{\&Delta;}{p}_{\mathrm{in}}\frac{{\&PartialD;}^{n}I\left(p\right)}{{\&PartialD;i}_1\&CenterDot;\&CenterDot;\&CenterDot;{\&PartialD;i}_n})=I\left(p\right)+\mathrm{\&Delta;}{p}_{i1}\frac{\&PartialD;I\left(p\right)}{{\&PartialD;i}_1}+\frac{1}{2!}\left(\mathrm{\&Delta;}{p}_{i1}\mathrm{\&Delta;}{p}_{i2}\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;i}_1{\&PartialD;i}_2}\right)\\ =I\left(p\right)+\mathrm{\&Delta;}{p}_x\frac{\&PartialD;I\left(p\right)}{\&PartialD;x}+\frac{1}{2!}(\mathrm{\&Delta;}{p}_x^2\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;x}^2}+\mathrm{\&Delta;}{p}_x\mathrm{\&Delta;}{p}_y\frac{{\&PartialD;}^{2}I\left(p\right)}{\&PartialD;x\&PartialD;y})\\ +\mathrm{\&Delta;}{p}_y\frac{\&PartialD;I\left(p\right)}{\&PartialD;y}+\frac{1}{2!}(\mathrm{\&Delta;}{p}_y^2\frac{{\&PartialD;}^{2}I\left(p\right)}{{\&PartialD;y}^2}+\mathrm{\&Delta;}{p}_y\mathrm{\&Delta;}{p}_x\frac{{\&PartialD;}^{2}I\left(p\right)}{\&PartialD;y\&PartialD;x})\end{array}---\left(3\right)$

Gaussian smoothing is carried out to image, it is A (p, σ)=(I*G) (p, σ) that the image after gaussian filtering is carried out for image I (p), and the n rank partial derivatives of filtering image are

Then it can be seen from Convolution Properties：

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

In formula：Footmark i₁,...i_nRepresent that data seek partial derivative to variable respectively, Gauss partial derivative is regarded as a wave filter, the combination of Gauss partial derivative response constitutes a complete feature description；

The candidate point of boundary point and the gray scale cost of each candidate point are obtained using local 2-jet characteristic images, local 2-jet is characterized as：

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

Step 2.2：Choose the candidate point of boundary point：For each point on initial lung border, the similarity degree of the gray scale of all pixels point and respective point gray scale in training characteristics image in the point search region in all characteristic images is calculated, the maximum point of similarity degree is selected, is used as the candidate point of the boundary point.Similarity degree is the mahalanobis distance h of surrounding pixel point gray scale respective point surrounding pixel point gray scale set into training sample characteristic image in all characteristic images_i：

Step 2.3：Lung region segmentation is carried out using Dynamic Programming；In boundary point region of search, the grey similarity cost of pixel is the similarity degree h of the surrounding pixel point gray scale of surrounding pixel point gray scale boundary point corresponding to training image_iIn boundary point region of search, the grey similarity cost of pixel is the similarity degree h of the surrounding pixel point gray scale of surrounding pixel point gray scale boundary point corresponding to training image_i；The shape similarity cost of i-th of boundary point is defined as in image：

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

In formula, v_i=p_i+1-p_i, the shape facility of i-th of boundary point is represented,

The average and covariance of i-th of boundary point shape facility in all training images are represented respectively,

For i-th of boundary point p in test image_i, having m in specified region of search has smaller grey similarity cost candidate point, then n boundary point will produce n × m gray scale cost matrix：

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

$h_i=\underset{j=1}{\overset{N}{\mathrm{\&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}(k-1)\right)\\ \sin \left(\frac{2\mathrm{\&pi;}}{n_c}(k-1)\right)\end{array}\right]---\left(11\right)$

In formula：

For on characteristic image, with point p_iFor the center of circle, r_cFor the n on the circle of radius_cThe gray scale of individual point, k=1......n_c。

The average and covariance of the surrounding pixel point gray scale of i-th of boundary point respectively in j-th of characteristic image of training image, N are characterized total number of images,

It is to look for an optimal path to search for optimal profile, and often row selects an element in Matrix C, when along Path selection, and the summation of gray scale and shape similarity cost is minimum, i.e.,：

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

Wherein, γ is gray scale and shape similarity cost coefficient.Adjust γ values so that both costs play roughly the same effect during boundary search；

Step 3：Extracted using the characteristic image of B-spline multi-scale wavelet transformation；During model foundation, the characteristic image for extracting effectively description different scale bone is the basis that model is set up；

Step 3.1：B-spline multi-scale wavelet transformation：B-spline function quickly converges on Gaussian function with the increase of batten exponent number, and its first derivative can approach Optimal edge detection operator；Multi-scale edge enhancing is carried out using B-spline small echo；Utilize function：

$H\left(e^{\mathrm{i\&omega;}}\right)=\mathrm{\&Sigma;}{h}_k{e}^{-\mathrm{ik\&omega;}}=\frac{{\widehat{\mathrm{\&theta;}}}_m\left(2\mathrm{\&omega;}\right)}{{\widehat{\mathrm{\&theta;}}}_m\left(\mathrm{\&omega;}\right)}=\frac{e^{-\mathrm{i\&epsiv;\&omega;}}{\left(\frac{\sin \mathrm{\&omega;}}{\mathrm{\&omega;}}\right)}^{m+1}}{e^{-i\frac{\mathrm{\&epsiv;\&omega;}}{2}}{\left(\frac{\sin (\mathrm{\&omega;}/2)}{\mathrm{\&omega;}/2}\right)}^{m+1}}=e^{-i\frac{\mathrm{\&epsiv;\&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}{\mathrm{\&Sigma;}}{g}_k{e}^{-\mathrm{ik\&omega;}}---\left(20\right)$

Ask for B-spline wavelet transform filter coefficient：

In formula：h_kFor low-pass Gaussian filter coefficient, g_kFor high pass Gaussian filter coefficient, gained h is utilized_kAnd g_kCarry out Wavelet fast decomposition to image, the row of original image respectively with h_kAnd g_kConvolution is carried out, then arranges it carry out down-sampling, two width subgraphs can be obtained, two width subgraphs are filtered then along the direction of row and down-sampling is carried out, it is possible to the output subgraph of four 1/4 sizes, diagonal detail pictures is obtained, vertical detail image, level detail image and approximate image；

Step 3.2：Gaussian filtering part 2-jet characteristic images are extracted, because the characteristic vector that Gauss partial derivative is constructed can be used for describing the local feature of image, according in formula (7) formula, 2-jet features (I, the I of the detail pictures after different scale wavelet transformation are extracted_x,I_y,I_xx,I_xy,I_yy) figure sets up regression model；

Step 4：The prediction of bone image is carried out using support vector regression, by trying to achieve regression function f (x)=ω x+b, skeleton model is set up；If sample set is (x_i,y_i), i=1,2 ..., l, x_i∈Rⁿ, sample set S is the linear approximation of ε-insensitive loss function.As linear regression function f (x)=ω x+b fittings (x_i,y_i) when, it is assumed that the error of fitting precision of all training datas is ε, i.e.,：

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

Allow d_iRepresent point (x_i,y_i) ∈ S to hyperplane f (x) distance：

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

Because S collection is ε-linear approximation, have |<w,x>+b-f(x_i)|≤ε,i=1,...,m.It can obtain:

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

Then have

$d_i＜\frac{\mathrm{\&epsiv;}}{\sqrt{1+{\left|\right|w\left|\right|}^2}},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,l---\left(25\right)$

(23) formula shows

It is point in S to the upper bound of the distance of hyperplane.The linear approximation collection S of ε-insensitive loss function best fit approximation hyperplane is to obtain hyperplane by maximizing the point in S to the upper bound of hyperplane distance；

Optimal hyperlane is by maximizing

What is obtained (minimizes

).Therefore, as long as minimizing | | w | |²It can be obtained by best fit approximation hyperplane.Then linear regression problem is just turned to：

Seek following optimization problem：

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

It is constrained to |<w,x>+b-y_i|≤ε,i=1,...,m.Additionally, it is contemplated that the situation of error of fitting, introduces relaxation factor

When division has error,

Both greater than 0, error is not present and takes 0；

Thus, the insensitive support vector regression of criterion epsilon can be expressed as：

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

Wherein, C>0 is balance factor.Introducing Lagrangian is：

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

Wherein, α_i,

γ_i>=0, i=1 ..., n.Function L extreme value should meet conditionThen following formula is obtained：

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

$\underset{i=1}{\overset{n}{\mathrm{\&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) are updated in (28), the dual form for obtaining optimization problem is：

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

It is constrained to：

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

Because bone gray scale is predicted as nonlinear regression, so data are first mapped to a high-dimensional feature space using a Nonlinear Mapping φ, then carrying out linear regression in high-dimensional feature space sets up model.Replaced due in superincumbent optimization process, only taking into account the inner product operation in high-dimensional feature space, therefore with a kernel function k (x, y)<φ(x),φ(y)>Nonlinear regression can just be realized.Then, the optimization method of nonlinear regression is the following function of maximization：

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

Its constraints is (32), after the value for solving α, it is possible to obtain f (x)：

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

Under normal circumstances, most of α_iOr

Value will be zero, the α being not zero_iOr

Corresponding sample is referred to as supporting vector, according to optimal condition, has in saddle point：

$\left\{\begin{array}{cc}{\mathrm{\&alpha;}}_i[{\mathrm{\&xi;}}_i+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)]=0& ,i=1,...l\\ a_i^{*}[{\mathrm{\&xi;}}_i^{*}+\mathrm{\&epsiv;}-y_i+f\left(x_i\right)]=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)$

Then the calculating formula for obtaining b is as follows：

$b=y_i-\mathrm{\&epsiv;}-\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x_j,x_i),$ Work as α_j∈(0,C)

$b=y_i+\mathrm{\&epsiv;}-\underset{i=1}{\overset{l}{\mathrm{\&Sigma;}}}({\mathrm{\&alpha;}}_i-{\mathrm{\&alpha;}}_i^{*})K(x_j,x_i),$ When

B value can be calculated with any one supporting vector, regression function f (x) is obtained；Set up after forecast model, generation bone structure image can be predicted according to the grey value profile of lung's x light images；

Step 5：The soft-tissue image for doing image subtraction to normal rabat and the bone image predicted to predict.

2. the method that bone suppresses in a kind of chest x-ray image according to claim 1, it is characterised in that：The step 1.2：Alignment training shapes are vectorial to be comprised the following steps that：

The described specific method that the data of piecemeal are mapped into feature space is as follows：

Step 1.2.1：Rotation, each lung region shape of zooming and panning, make it be alignd with first shape in training set；

Step 1.2.2：According to alignment shape, average shape is calculated；

Step 1.2.3：Rotation, zooming and panning average shape make it be alignd with first shape；

Step 1.2.4：Again each shape is alignd with current average shape；

Step 1.2.5：If process restrains or to designated cycle number of times, exit；Otherwise step is gone to.

3. the method that bone suppresses in a kind of chest x-ray image according to claim 1, it is characterised in that：The characteristic image of B-spline multi-scale wavelet transformation is extracted as 3 rank B-spline multi-scale wavelet transformations in the step 3, as m=3, and every filter coefficient is k=- 2, and -1,0,1,2,

g_k={0,0,-2,2,0}.As can be seen here, the finite length of the wave filter of three rank B-spline small echos is 5.

4. the method that bone suppresses in a kind of chest x-ray image according to claim 1, it is characterised in that：The decomposition of three multi-scale wavelets is done in the step 3 to chest x light images.

5. the method that bone suppresses in a kind of chest x-ray image according to claim 1, it is characterised in that：In the step 3.2, the yardstick of gaussian filtering is σ=2,4.

Images