CN102867294B - Fourier-wavelet regularization-based coaxial phase contrast image restoration method - Google Patents

Abstract

The invention belongs to the field of biomedical engineering and medical imaging, and relates to a fourier-wavelet regularization-based coaxial phase contrast image restoration method. The method comprises the steps: collecting a plurality of images of a tool edge instrument, and acquiring tool edge section curves of each image at different positions to obtain an integrated line spread function h(x) of a coaxial phase contrast imaging system; placing an imaging object, and collecting images to obtain an imaging result y(x); calculating fourier transformation of the y(x) and the h(x) to obtain Y(omega) and H(omega); obtaining fourier regularization evaluation of an ideal phase contrast imaging result through fourier inverse transformation, and effectively evaluating the wavelet coefficient signal wj(x) of the ideal phase contrast f(x); and conducting Db6 orthogonal wavelet transform on the fourier regularization evaluation, and carrying out regularization to obtain wavelet inverse transformation and obtain an evaluation result. According to the fourier-wavelet regularization-based coaxial phase contrast image restoration method, the contrast of the phase contrast image under deterioration effect can be effectively improved, the fidelity of the restoration image can be ensured, and the image restoration contrast and signal to noise ratio can be improved.

Description

Based on the line phase-contrast image recovery method of Fourier and small echo regularization

Art

The invention belongs to biomedical engineering and Medical Imaging field, relate to a kind of image recovery method,

Background technology

Breast cancer is current women first " killer ", and according to World Health Organization's statistics, within 2000, the whole world about has 1,200,000 women to be diagnosed as breast cancer, and 500,000 people die from breast cancer.Breast cancer not only jeopardizes patient vitals, and causes the damage of muliebrity organ, and serious harm women is physically and mentally healthy.Although China is not breast cancer country occurred frequently, average annual growth exceeds national 1-2 percentage point occurred frequently nearly ten years, with the speed increase of annual 3%-4%.Because China is populous, the diagnosis and treatment of breast cancer have become day by day heavy and have needed the social concern of solution badly, and realize the key that early diagnosis is this social concern of solution, raising survival and quality of life.

The Main Means of current mammary gland routine inspection is molybdenum (rhodium) target X ray breast imaging art, but long-term clinical practice shows, this technology in sensitivity, specificity, all there is great shortcoming in the aspect such as security and comfortableness: on the one hand this technology exists the loss up to 10-15%; Be the positive by this diagnosis of technique on the other hand and final biopsy diagnosis rate is 25-29%.Especially seriously, because the breast of young woman is too fine and close, diagnostic accuracy is caused seriously to reduce.

Until last century Mo, X-ray phase contrast imaging theory (X-ray phase contrast imaging, XPCI) proposition, has broken traditional x-ray imaging theory, for the early-stage breast cancer anatomy imaging diagnostic techniques of realizing ideal brings new dawn.Research shows, under identical radiation dose, the contrast resolution of phase contrast imaging improves about 10 times than traditional X-ray line absorption contrast imaging, significantly improves the image visibility of soft-tissue imaging.At present, the researchs of various countries all in Efforts To Develop X-ray phase contrast such as Europe, the United States, Japan, Australia, x-ray phase contrast imaging is extensively thought " can give radiodiagnosis medical science bring dramatic change " micro-imaging technique.

The technology being currently available for X-ray phase contrast can be divided into interference imaging method (inteeferometry) according to image-forming principle difference, diffraction enhanced imaging method (diffraction enhanced imaging, DEI), coaxial imaging method (in-line holography), grating differential phase contrast imaging method (Differential phase contrast, and coded aperture phase contrast imaging method (Coded-aperture based x-ray phase contrast imaging) DPC), these methods are at X source, experimental provision, the aspects such as detector and imaging performance index are very different.Interference imaging method and diffraction enhanced imaging method all need accurate monochromating crystal, cause them at least to there are three problems in actual applications thus: the monochromating crystal 1, in system must carry out high precision alignment, and therefore imaging system is very responsive to environmental perturbation; 2, monochromating crystal means homogeneous X-ray bundle, this just causes most of light source cannot meet the luminous flux requirements of imaging system, solve this problem at present, but synchrotron radiation source is costly mainly through synchrotron radiation source, floor area is large, is unfavorable for clinical practice and popularization; 3, crystal can absorb a part from the X ray after measurand outgoing, is therefore difficult to control optimal imaging dosage.Grating differential phase contrast imaging method is the phase contrast imaging technology that development in recent years is got up, because the method can adopt common X-ray light source, therefore larger clinical practice potentiality are possessed, but still there are some at present and need technical solution problem badly in it, such as, DPC method needs to introduce high-precision parallel grating (light source grating, phase grating and absorption grating), adds the requirement of system accuracy and stability; Due to light receiver limited angle, reduce system luminous flux, inevitably cause system exposure times elongated; In addition, some external factor such as patient temperature also can bring very important impact to grating precision.Coded aperture phase contrast imaging method can adopt relatively large interval, coded aperture (1 ~ 2 order of magnitude larger than the fence interval in DPC method), improve the requirement to system accuracy in DPC method, but it is still faced with the other problems needing solution in DPC method badly.Line phase-contrast imaging does not need to introduce extra optical devices directly can obtain interior of articles fine structure information, its imaging mode is otherwise known as free-space propagation phase contrast imaging (free-space propagation), prove during the people such as Wilkins study in early days that line phase-contrast imaging can adopt heterogeneous X-ray light source, thus avoided the problem of imaging system luminous flux aspect.On the other hand, the investment of line phase-contrast imaging system is little, and floor area is no more than typical X-ray photography, and therefore generalization is strong, is expected the Main Means that the photography of alternative typical X-ray becomes Mass Screening of Breast Cancer.In current all kinds of phase contrast imaging technology, the light source that coaxial imaging method, grating differential phase contrast imaging method and coded aperture phase contrast imaging method adopt is all the feasible X source of engineering, has the relative advantage of clinical practice.Wherein, coaxial imaging method can be considered the basis of grating differential phase contrast imaging method and coded aperture phase contrast imaging method, and latter two method is equivalent in coaxial imaging optical path, introduce extra optical devices (parallel grating or coded aperture).As can be seen here, the technical development of coaxial imaging method effectively must promote the development of grating differential phase contrast imaging method and coded aperture phase contrast imaging method.

So far, line phase-contrast imaging technique is subject to extensively giving more sustained attention of domestic and international irradiation image field always, and this technology is the most applicable under being considered to conditions present realizes one of micro-imaging technique of clinical medicine application conversion.It is worth mentioning that, there is the company being engaged in X-ray phase contrast equipment development abroad, as Japanese Konica Minolta company (http://konicaminolta.jp) is used for clinical breast disease diagnosis phase contrast imaging technology in the world the earliest, and be proposed the phase contrast imaging mammography system of global first item based on common X source.But up to now, the X ray line phase-contrast advantage be imaged in breast cancer clinical diagnosis is also far from showing, the phase contrast imaging mammography system of Konica Minolta company, image quality improves little relative to existing Traditional x-ray mammary gland digital imaging system, do not reach the Expected Results of pertinent literature report, the recall rate adopting this system to obtain and cancer recall rate and classic method are not significantly distinguished.

By the actual physical process of analysis line phase-contrast imaging and the phase contrast imaging model of bibliographical information, can find, restricting current line phase-contrast imaging technique is mainly reflected in picture system aspects in the key issue of clinical application, picture system self also exists defect, such as X source not desirable point source, detector performance is subject to the restriction of self factor such as resolution and point spread function, and system exists all kinds of hazardous noises etc.

For detector imperfection problem, the scholar Olivo of London university adopts Wiener deconvolution method to carry out preliminary trial, but the simulation of Olivo and test result for be synchrotron radiation origin system.Synchrotron radiation origin system has extraordinary light source and very high signal to noise ratio (S/N ratio), but this systematic cost is expensive, the Chinese Shanghai synchrotron radiation light source built up for such as 2010, gross investment 1,400,000,000, floor area about 200,000 square metres.On the other hand, the quality of the Recovery image that Olivo adopts dimension deconvolution method to obtain, also has suitable gap with actual demand.It is as follows that Olivo obtains result: ideally, and phase contrast imaging contrast is more than 200%; Because detector imperfection causes deterioration, phase contrast contrast drops to only has 3%; By Wiener deconvolution method, phase contrast can bring up to 7.2% and 11.9%, but in 11.9% situation that phase contrast is higher, introduces noise larger.

Because the imaging of Microfocus X-ray line phase-contrast is the imaging technique being more suitable for realizing clinical application, therefore, the raising realizing phase contrast imaging quality for this imaging technique has more significant researching value and meaning.Meanwhile, explore a kind of new method of more effective phase contrast Quality advance, become one of key issue applied towards line phase-contrast imaging technique and System Development, tool is of great significance.

Summary of the invention

Purport of the present invention proposes a kind of method that can improve X ray line phase-contrast imaging results contrast and signal to noise ratio (S/N ratio), under solving current engineering condition with this, the key issue that line phase-contrast imaging faces: due under current engineering specifications, Microfocus X-ray X source, the imperfection of detector and system noise etc. cause phase contrast image deterioration, and phase contrast reduces.Technical scheme of the present invention is as follows:

Based on a line phase-contrast image recovery method for Fourier-small echo regularization, comprise the following steps:

1. the multiple image of line phase-contrast imaging system acquires edge of a knife utensil is utilized;

2. from the edge of a knife cross section curve of every width Image Acquisition diverse location,

3. edge of a knife cross section curve is averaged, then averaged curve is differentiated, obtain overall wire spread function h (x) of line phase-contrast imaging system;

4. place imaging object, gather image, obtain imaging results y (x);

5. calculate the Fourier transform of y (x) and h (x), obtain Y (ω) and H (ω)

6. formula is passed through obtain Fourier's regularization frequency-domain result α is regularization parameter,

H ^*(ω) be H (ω) conjugate complex number; Fourier's regularization estimation of desirable phase contrast imaging result is then obtained by inverse Fourier transform

7. to the wavelet coefficient signal w of desirable phase contrast f (x) _jx () is effectively estimated, obtain estimated signal

8. Fourier's regularization is estimated carry out Db6 orthogonal wavelet transformation, wavelet function is expressed as ψ, and scaling function is expressed as φ, obtains wavelet coefficient signal after wavelet transformation

9. to wavelet coefficient signal utilize Regularization function Λ _jx () carries out small echo regularization, obtain regularization result ${\tilde{w}}_j\left(x\right)={\hat{w}}_j\left(x\right){\mathrm{\Λ}}_j\left(x\right);$

10. right carry out the wavelet inverse transformation of wavelet function ψ, the final Fourier-small echo regularization estimated result obtaining desirable phase contrast imaging result

Preferably, in the 7th step, Fourier's regularization is estimated adopt Db4 orthogonal wavelet transformation, wavelet function is expressed as ψ ', and scaling function is expressed as φ ', obtains wavelet coefficient signal right carry out hard domains value, namely ${\tilde{w}}_j^{\′}\left(x\right)=\left\{\begin{array}{cc}{\hat{w}}_j^{\′}\left(x\right),& \left|{\hat{w}}_j^{\′}\left(x\right)\right|>{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\\ 0,& \left|{\hat{w}}_j^{\′}\left(x\right)\right|\≤{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\end{array}\right.,$ Here the noise equation of wavelet function ψ ' under jth yardstick, ρ _jthat yardstick is correlated with the thresholding factor; Then right carry out the inverse transformation of Db4 orthogonal wavelet, obtain denoising estimated signal right carry out Db6 orthogonal wavelet transformation, obtain wavelet coefficient signal be w _jeffective estimated signal of (x).

Adopt a kind of Fourier towards line phase-contrast imaging of the present invention-small echo regularized image to recover new technology, effectively can improve the contrast of phase contrast image under degradation effects, and ensure that the fidelity of Recovery image.The present invention, in raising Postprocessing technique contrast and signal to noise ratio (S/N ratio), all obtains than Wiener Filtering and the better result of Tikhonov method.The present invention is on the basis of line phase-contrast imaging, achieve and improve the effective mass of undesirable imaging results, therefore achievement in research of the present invention effectively can make up the deficiency for imaging performance in the line phase-contrast imaging system of early diagnosing mammary cancer in actual clinical.The application of the method, provides technical support by the diagnosis for effectively realizing early-stage breast cancer minute lesion tissue, provides powerful support for for the clinical practice and research of carrying out the line phase-contrast imaging of breast cancer in a deep going way provide.

Accompanying drawing explanation

The performance transport function of Fig. 1 detector.

Fig. 2 .300 micrometer fibers at light source to object distance/object to detector distance=80cm/80cm situation under imaging results, (a) two dimensional image, (b) cross-section curve.

Fig. 3 .300 micrometer fibers at light source to object distance/object to detector distance=80cm/80cm situation under the cross-sectional view of ideal image result.

Fig. 4. after Wiener deconvolution, fibre image cross-sectional view (a) α=0.1 (b) α=0.07 of acquisition.

Fig. 5. (a) selects regularization matrix to be the corresponding constantly L curve of Laplace operator, the regularization matrix of optimum regularization parameter λ=2.6088 (b) corresponding diagram 5 (a) and regularization parameter, the phase contrast image restoration result obtained, contrast is 15.9%.

Fig. 6. adopt the fibre image cross-sectional view that Fourier-small echo regularization obtains.

Phase contrast figure (b) Wiener filtering result (α=0.2) (c) Wiener filtering result (α=0.1) (d) TikHonov regularization result (e) Fourier-small echo regularization result of 6 0.5 millimeter of pencil-leads that Fig. 7 (a) experiment obtains.

The average cross section figure of front 2 the pencil-lead phase contrast image of Fig. 8.

Embodiment

For in X ray line phase-contrast imaging process, due to detector, light source and system noise cause the deterioration problem of phase contrast imaging result, the present invention proposes the Fourier-Wavelet Regularization towards line phase-contrast imaging, the above-mentioned degradation effects influence on RT of effective suppression, improve the image contrast of phase contrast result details, ensure the fidelity of Recovery image simultaneously.The present invention is elaborated from several aspect below in conjunction with drawings and Examples.

1 digital X-ray imaging system

Experiment imaging system is Pixarray 100 toy digital radiation imaging system, is manufactured by BIOPTICS company of the U.S..The detector of this system is the ccd array of 1024 × 1024, and pixel size is 50 μm × 50 μm, 14 grades of gray scales.Horizontal and longitudinal spatial resolution is every millimeter of 20 pixels.The focal spot size of X-ray tube is 50 μm.The halfwidth of detector point spread function is 110 μm.In experiment, the operating voltage of x-ray source is 33kVp, and working current is 0.5mA.During phase contrast imaging, x-ray source is set and is 80cm to the distance of object and object to the distance of detector.Under this is arranged, the focal spot picture that light source becomes on the detector is 50 μm of halfwidths.The point spread function of total system is the convolution results of detector point spread function and light source point spread function.The method that the scholar Olivo that we have employed London university advises, edge of a knife object is placed in imaging object plane, we just directly can obtain the line spread function curve of whole system like this, then by line spread function central rotation, can obtain the point spread function of entire system.

For convenience of explanation, the phase contrast image of the present invention first under one-dimensional model worsens, and based on this, the phase contrast image expanding to two dimensional model worsens, and provides corresponding Postprocessing technique result.

2 towards the Wiener deconvolution technology of line phase-contrast imaging

For one-dimensional model, for microfocus X-ray line phase-contrast imaging system, due to the degradation effects impact that imperfection and the system noise of detector cause, line phase-contrast imaging results y (x) of our actual acquisition can represent with formula below:

y(x)＝f(x)*LSF(x)*s(x)+n(x)＝f(x)*h(x)+n(x) (1)

Wherein, y (x) is the phase contrast imaging result that actual tests obtains, * be convolution algorithm symbol, x is locus coordinate, f (x) is line phase-contrast imaging results ideally, LSF (x) is detector lines spread function, s (x) is the line spread function of light source imaging on the detector, h (x)=LSF (x) * s (x) is the line spread function of entire system, it is the convolution results of detector lines spread function and light source line spread function, and n (x) is system noise.

Inverse convolution method can remove the imaging results degradation effects that system point spread function causes, thus close to desirable line phase-contrast imaging results.But inverse convolution method requires that system point spread function is nonsingular.When system point spread function is unusual, the HFS of system noise will amplify by inverse convolution method, causes image detail fuzzy, cannot obtain required Postprocessing technique result.

Therefore, when there is system noise, our object finds one of ideal image result f (x) effectively to estimate, that is:

$\hat{f}\left(x\right)=\mathrm{\φ}\left(x\right)*y\left(x\right)---\left(2\right)$

Here it is the effective estimated result to f (x) under minimum mean square error criterion.

Wiener deconvolution terminates to provide wave filter φ (x) required for formula (2).Mathematically, Wiener deconvolution method applies for the Wiener filtering of the deconvolution problem that there is noise item.Wiener deconvolution is as follows in the expression formula of frequency field:

$\mathrm{\Φ}\left(\mathrm{\ω}\right)=\frac{H^{*}\left(\mathrm{\ω}\right)}{{\left|H\left(\mathrm{\ω}\right)\right|}^2+{\left|N\left(\mathrm{\ω}\right)\right|}^2/{\left|F\left(\mathrm{\ω}\right)\right|}^2}---\left(3\right)$

Wherein Φ (ω), F (ω), H (ω) and N (ω) are φ (x), f (x) respectively, the Fourier transform of h (x) and n (x), H ^*(ω) be H (ω) conjugate complex number.

Thus, the effective estimation of desirable phase contrast imaging result in frequency field can be obtained, namely

$\hat{F}\left(\mathrm{\ω}\right)=Y\left(\mathrm{\ω}\right)\mathrm{\Φ}\left(\mathrm{\ω}\right)=Y\left(\mathrm{\ω}\right)\frac{H^{*}\left(\mathrm{\ω}\right)}{{\left|H\left(\mathrm{\ω}\right)\right|}^2+{\left|N\left(\mathrm{\ω}\right)\right|}^2/{\left|F\left(\mathrm{\ω}\right)\right|}^2}---\left(4\right)$

Here with Y (ω) be with the Fourier transform result of y (x).

Be generally the unknown in systems in practice due to F (ω) and N (ω) and cannot measure, so often the formula of Wiener deconvolution is simplified, by system performance continuous item | N (ω) ²/ F (ω) | ²substitute with a constant α, α is regularization parameter:

$\mathrm{\Φ}\left(\mathrm{\ω}\right)=\frac{H^{*}\left(\mathrm{\ω}\right)}{{\left|H\left(\mathrm{\ω}\right)\right|}^2+\mathrm{\α}}---\left(5\right)$

The noise gathering image is depended in the selection of optimum a-value: large α value still can cause distorted signals by restraint speckle better.And although little α value can obtain signal more accurately, its cost is the introduction of more noises.What Olivo adopted in its research is exactly this image-recovery technique based on Wiener deconvolution.

3 towards the Tikhonov Regularization Technique of line phase-contrast imaging

For convenience of explanation towards the Tikhonov Regularization Technique of line phase-contrast imaging, the expression of formula (1) is first converted into matrix form from convolution form by us:

y＝Hf+n (6)

Here, small letter bold-faced letter y, f and n represents y (x) in formula (1) respectively, the vector form of f (x) and n (x), bold race capital H representative is transformed into the convolution kernel matrix after matrix form (formula 6) from the h (x) of convolution form (formula 1).The concrete form of H is:

$H=\left[\begin{array}{cccc}h\left(1\right)& h\left(N\right)& h(N-1)& \·\·\·h\left(2\right)\\ h\left(2\right)& h\left(1\right)& h\left(N\right)& \·\·\·h\left(3\right)\\ \·& & & \\ \·& & & \\ \·& & & \\ h\left(N\right)& h(N-1)& h=(N-2)& \·\·\·h\left(1\right)\end{array}\right]---\left(7\right)$

The actual image y obtained image objects is that desired result is worsened by ssystem transfer function, and contains the result of noise item n.

Formula (6) and formula (1) all can be called as deterioration model.Set up after worsening model, next step sets up solution for it exactly.Mention above, due to the existence of noise item, directly inverse convolution cannot obtain required Postprocessing technique result, under normal circumstances, this kind of problem can be called ill-conditioning problem.Adopt regularization method an ill-conditioning problem can be converted to good state problem, thus the Postprocessing technique result required for obtaining.

The object of regularization method introduces the relevant information of required restoration result, with this stable problem, obtains effective stable solution.For formula (6), multiple regularization method can be taked to realize Postprocessing technique problem.Wherein Tikhonov regularization method has good characteristic in stability of solution and validity.

Based on the ultimate principle of Tikhonov regularization method, consider that in formula (6), H is morbid state or singular matrix, in order to obtain the solution of required characteristic, we take the regularization minimizing criterion as follows:

||Hf-y|| ²+λ||Lf|| ²(8)

Wherein Tikhonov regularization matrix L and regularization parameter λ needs to select according to actual needs.Under many circumstances, when thinking that original signal is the signal of a basic continous, Tikhonov regularization matrix generally adopts high-pass filtering operator to improve the slickness recovering rear signal.Because actual phase contrast image is smooth continuous print, we adopt Laplce's high-pass filtering operator as Tikhonov regularization matrix, this operator recovers in image boundary and has very large advantage in restoring signal fidelity, therefore does one with method of the present invention and compares and have better cogency.

Regularization matrix L selected by correspondence, first we need to determine regularization parameter λ.Select the method for regularization parameter can have multiple, comprise cross validation method, constraint maximum likelihood method and L curve method.Wherein L curve method image can state error residue amount || Hf-y|| and signal high-frequency noises || weigh selection course between Lf||.In L curve method, horizontal ordinate is || Hf-y||, and ordinate is || Lf||, can obtain a lot of point (|| Hf-y||, || Lf||) thus, the regularization parameter λ that different points is corresponding different, obtains a curve through over-fitting.Select large regularization parameter λ can obtain the semi-norm of less solution || Lf||, but cost introduces very large noise.Although and select little regularization parameter λ can obtain higher signal to noise ratio (S/N ratio), distorted signals degree improves.Therefore, the balance of signal recuperation precision and signal noise, the optimal balance point of both selections, that point that namely on auditory localization cues, curvature is maximum is exactly the position of the regularization parameter λ needed for us.

Worsen result to line phase-contrast imaging and adopt Tikhonov regularization method, we can obtain numerical solution, and result is as follows:

$\hat{f}=(H^{T}H+{\mathrm{\λ}}^2{L}^{T}L)H^{T}y---\left(9\right)$

The contrast of comprehensive Recovery image and signal to noise ratio (S/N ratio), Tikhonov regularization method can obtain the result more excellent than Wiener Filtering.

4 Fourier-small echo Regularization Technique

The present invention is intended to study the Fourier-small echo Regularization Technique being applicable to line phase-contrast Postprocessing technique.This technology comprises Fourier's regularization and these two steps of small echo regularization de-noising aftertreatment.The first step, by reversion operator operation and Fourier's regularization computing to alleviate the degradation effects in phase contrast imaging result, but this step brings noise enlarge-effect simultaneously.Second step, effectively reduces the noise introduced in the first step by Wavelet Regularization.

In the first step of Fourier-small echo Regularization Technique, in reversion operator operation and Fourier's regularization computing, usually adopt regularization sharpening filter, such as S filter (formula (5)).In order to recover more object boundary phase contrast stripe information, in Fourier-small echo regularization algorithm, generally taking less regularization parameter, result also in the cost that noise amplifies thus.

In the second step of Fourier-small echo Regularization Technique, small echo regularization computing is taked, by removing the noise amplified by regularization sharpening filter in a first step based on the nonlinear filtering of wavelet transformation to the boundary sharpening signal that the first step obtains.This nonlinear filtering computing is the adaptation field value process based on wavelet field.What adopt here is redundant wavelet computing.The estimated signal that we obtain Wiener filtering carry out wavelet transformation, redundant wavelet scaling function symbol φ represents, wavelet function symbol ψ represents.Right after wavelet transformation, a series of scale coefficient signal can be obtained and wavelet coefficient signal here subscript j represents the yardstick of wavelet transformation.Small echo regularization object is from wavelet coefficient signal middle removal noise artefact, have employed Wavelet Domain Wiener Filtering method to realize in the present invention.Adopt Wavelet Domain Wiener Filtering, to wavelet coefficient signal take desirable small echo regularization implementation method as follows:

${\tilde{w}}_j\left(x\right)={\hat{w}}_j\left(x\right){\mathrm{\Λ}}_j\left(x\right)---\left(10\right)$

${\mathrm{\Λ}}_j\left(x\right)=\frac{1}{1+{\mathrm{\σ}}_j^2\left(x\right)/w_j^2\left(x\right)}---\left(11\right)$

Here Λ _jx () represents the desirable regularization to wavelet coefficient, the noise variance being positioned at position x under jth wavelet scale, w _jx () is the wavelet coefficient signal of line phase-contrast imaging results f (x) ideally.

Due to the space invariance of Wiener filtering, the residual noise amplified in Fourier-small echo regularization algorithm first step is stable.Simultaneously because redundant wavelet transformation also possesses space invariance, also should be so stable in whole signal at the residual noise of certain wavelet space.Therefore, small echo regularization can be rewritten as

${\mathrm{\Λ}}_j\left(x\right)=\frac{1}{1+{\mathrm{\σ}}_j^2/w_j^2\left(x\right)}---\left(12\right)$

Here represent the noise variance under jth wavelet scale, can be calculated by formula below

${\mathrm{\σ}}_j^2=\underset{\mathrm{\ω}}{\mathrm{\Σ}}\frac{{\mathrm{\σ}}^2{\left|H\left(\mathrm{\ω}\right){\mathrm{\Ψ}}_j\left(\mathrm{\ω}\right)\right|}^2}{{({\left|H\left(\mathrm{\ω}\right)\right|}^2+{\mathrm{\σ}}^2/{\left|W_j\left(\mathrm{\ω}\right)\right|}^2)}^2}---\left(13\right)$

Here σ ²overall noise estimated value, Ψ _j(ω) be for the Fourier transform result of wavelet function ψ under jth wavelet scale, W _j(ω) be w _jthe Fourier transform of (x).To Fourier's regularization result in the present invention get the thinnest yardstick of wavelet transformation, then mediant estimation is got to this scale coefficient and can obtain σ ².

Λ _j(x) and can not directly obtain, because the desirable wavelet coefficient w that phase contrast imaging result f (x) of theoretical prediction is corresponding _jx () is (formula (12), the formula (13)) that cannot obtain.Therefore, we are effectively to estimate carry out alternative w _j(x). can pass through Fourier's regularization result wavelet coefficient carry out hard domains value acquisition.Need in this step to adopt another wavelet function, in order to distinguish with the scaling function φ adopted before and wavelet function ψ, the scaling function adopted here symbol φ ' represents, and wavelet function symbol ψ ' represents.Consider the wavelet function ψ ' under jth yardstick _j, wavelet coefficient signal represent Fourier's regularization result wavelet function ψ ' is adopted to carry out the result of wavelet transformation.Then, we are in ψ ' wavelet function space, to wavelet coefficient signal carry out thresholding contraction, namely

${\tilde{w}}_j^{\′}\left(x\right)=\left\{\begin{array}{cc}{\hat{w}}_j^{\′}\left(x\right),& \left|{\hat{w}}_j^{\′}\left(x\right)\right|>{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\\ 0,& \left|{\hat{w}}_j^{\′}\left(x\right)\right|\≤{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\end{array}\right.---\left(14\right)$

Here the noise equation of wavelet function ψ ' under jth yardstick, ρ _jbe that yardstick is correlated with the thresholding factor, can realize according to Mallat small echo choosing method.Wavelet coefficient signal after thresholding is shunk carry out inverse small echo computing, the effective estimated signal to phase contrast imaging result f (x) through denoising can be obtained

Have we just can adopt scaling function φ and wavelet function ψ to carry out wavelet transformation to it, obtain desirable wavelet coefficient w _jeffective estimation of (x) substitute w _jx (), namely obtains Λ by formula (12) and (13) _j(x) and thus, the coefficient that formula (10) and (11) obtains small echo regularization is recycled finally by right carry out wavelet inverse transformation, the Fourier-small echo regularization estimated result of desirable phase contrast imaging result f (x) can be obtained

In order to the Wiener deconvolution result taked with Olivo compares, advantage of the present invention is described, we adopt the criterion to phase contrast Recovery image quality in Olivo achievement in research.The phase contrast computing formula that Olivo adopts is as follows

$\mathrm{Comtrast}=\frac{I_{\max }-I_{\min }}{I_{\mathrm{background}}}---\left(15\right)$

Here I _maxthe maximal value of positive overshoot crest in phase contrast result, I _minthe minimum value of negative overshoot trough in phase contrast result, I _backgroundit is the background intensity outside phase contrast overshoot details.

5 towards the Fourier-small echo Regularization Technique application flow of line phase-contrast imaging

The flow process that phase contrast image based on Fourier-small echo regularization of the present invention recovers new method is described below:

1, X ray line phase-contrast imaging parameters is arranged: in the present invention, and arranging light source to the distance of object is 80cm, and object is correspondingly 80cm to the distance of detector.

2, the exposure parameter of digital radiation imaging system is set, place edge of a knife utensil range finder surface 80cm (being namely positioned in imaging object plane), continuous acquisition 15 width image, from the edge of a knife cross section curve 50 of every width Image Acquisition diverse location, then 15*50 bar edge of a knife cross section curve is averaged, again averaged curve is differentiated, obtain line spread function h (x) of entire system.

3, place imaging object, in the present invention, adopt the pencil-lead of 300 micron polyethylene fibers and 0.5mm diameter.Under imaging parameters is arranged (light source is to object distance/object to detector distance=80cm/80cm), to image objects, obtain imaging results y (x).

4, calculate the Fourier transform of y (x) and h (x), obtain Y (ω) and H (ω).Pass through formula obtain Fourier's regularization frequency-domain result fourier's regularization estimation of desirable phase contrast imaging result f (x) is then obtained by inverse Fourier transform

5, to the wavelet coefficient signal w of desirable phase contrast f (x) _jx () is effectively estimated, namely get estimated signal carry out alternative w _j(x): to Fourier's regularization result adopt Db4 orthogonal wavelet transformation, wavelet function is expressed as ψ ', and scaling function is expressed as φ ', obtains wavelet coefficient signal right carry out hard domains value, namely ${\tilde{w}}_j^{\′}\left(x\right)=\left\{\begin{array}{cc}{\hat{w}}_j^{\′}\left(x\right),& \left|{\hat{w}}_j^{\′}\left(x\right)\right|>{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\\ 0,& \left|{\hat{w}}_j^{\′}\left(x\right)\right|\≤{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\end{array}\right.,$ Then right carry out the inverse transformation of Db4 orthogonal wavelet, obtain denoising estimated signal right carry out Db6 orthogonal wavelet transformation, obtain wavelet coefficient signal be w _jeffective estimated signal of (x).

6, Fourier's regularization is estimated carry out Db6 orthogonal wavelet transformation, wavelet function is expressed as ψ, and scaling function is expressed as φ.Wavelet coefficient signal is obtained after wavelet transformation

7, to wavelet coefficient signal carry out small echo regularization, Regularization function adopts Λ _jx (), obtains regularization result

right carry out the wavelet inverse transformation of wavelet function ψ, the final Fourier-small echo regularization estimated result obtaining desirable phase contrast imaging result

The Pixarray 100 toy digital radiation imaging device that the present invention adopts BIOPTICS company of the U.S. to produce builds line phase-contrast imaging system.First the transport function of system detector is obtained by knife edge device.Fig. 1 gives and measures by knife-edge method the detector transport function obtained, and measuring the halfwidth obtaining ssystem transfer function is thus 120 microns.

Fig. 2 (a) give imaging system obtain 300 micrometer fibers at light source to object distance/object to detector distance=80cm/80cm situation under imaging results.Because the observability in details of two dimensional image is poor, we illustrate the cross-section curve (see Fig. 2 (b)) of corresponding diagram 2 (a).From Fig. 2 (b), due to the degradation effects of ssystem transfer function, the actual contrast obtained only has about 2.8%, and there is more serious system noise (signal to noise ratio (S/N ratio) that can be calculated signal is 16.2dB).

For the ease of analyzing, we simulate under said system facilities simultaneously, the one dimension cross-sectional view of the line phase-contrast imaging desired result of 300 micron polyethylene fibers, as shown in Figure 3.As can be seen from Figure 3, to 300 micron polyethylene fiber imagings, when light source is to object distance/object to detector distance=80cm/80cm, the contrast of obtainable desirable phase contrast imaging result can reach more than 200%.And comparison diagram 2, due to deterioration detector effect, the contrast of fiber picture drops to only has about 2.8%.

Wiener deconvolution method during first we adopt Olivo to study, and Wiener filtering parameter alpha=0.1 and α=0.07 are set.By Wiener deconvolution, we obtain the cross-sectional view of Postprocessing technique result as shown in Figure 4.After Wiener filtering, when adopting α=0.1, image contrast brings up to 8.7% (as Fig. 4 (a)) from original 2.8%, and when adopting α=0.07, image contrast brings up to 16.6% from original 2.8%.But when adopting α=0.07, after recovering, the ground unrest of image is than high a lot of during α=0.1.This studies with Olivo the result provided is consistent.

The Tikhonov regularization method that then we have employed towards line phase-contrast imaging carries out phase contrast image recovery.For above-mentioned deterioration image (as Fig. 2), select Laplace operator as regularization matrix, then adopt L curve method to choose the regularization parameter λ of corresponding optimum.Fig. 5 (a) gives when selecting regularization matrix formula unit matrix, error residue amount norm || and Hf-y|| is the norm of horizontal ordinate, solution || and f|| is the L curve of ordinate gained, gives curve by log-log mode, so that observation and analysis.As can be seen from the figure, the corresponding L point of inflexion on a curve of optimum regularization parameter λ, i.e. λ=2.6088.Fig. 5 (b) gives the cross-sectional view of corresponding Postprocessing technique result, and by can be calculated, after Postprocessing technique, contrast brings up to 15.9%.

Fig. 6 gives the Postprocessing technique result adopting Fourier-Wavelet Regularization to obtain, and image contrast brings up to 18.5% from original 2.8%.In the first step of Fourier-Wavelet Regularization, during Fourier's regularization, we have employed Fourier regularization parameter=0.03, and this is less than α=0.07 adopted during Wiener filtering.Compare with Wiener filtering and Tikhonov regularization method, can find out, Fourier-Wavelet Regularization can obtain the highest recovery contrast and best signal to noise ratio (S/N ratio).

In order to feasibility in phase contrast image recovery of more above-mentioned 3 kinds of methods further and validity, we carry out phase contrast imaging to another object, then adopt method in above-mentioned 3 to obtain Postprocessing technique result respectively.Fig. 7 (a) is the original phase contrast image of the pencil-lead to the 0.5mm diameter that 6 orders are put, and Operation system setting is consistent with experiment 1.By calculating, the signal to noise ratio (S/N ratio) of this image is 18.7dB.Because pencil-lead itself has certain absorption, therefore the dark fringe of phase contrast striped is absorbed contrast covering, only to see bright fringe.Fig. 7 (b) and Fig. 7 (c) adopts Wiener filtering, and the result obtained when filtering parameter is set to α=0.2 and α=0.1 respectively, the signal to noise ratio (S/N ratio) of correspondence image is 25.4dB and 21.3dB respectively.Fig. 7 (d) is the phase contrast image deconvolution result adopting Tikhonov regularization method to obtain, and obtains optimum regularization parameter λ=6.2357 by L curve method, and the final signal to noise ratio (S/N ratio) obtaining image is 24.2dB.Adopt the Postprocessing technique result of Fourier-small echo regularization acquisition as shown in Fig. 7 (e), the signal to noise ratio (S/N ratio) of this figure is 30.1dB.

In order to effectively contrast above-mentioned 3 kinds of method raisings to image contrast, we illustrate the average cross section figure of front 2 pencil-leads in Fig. 7.Before carrying out Postprocessing technique, the average contrast of the original phase contrast figure of pencil-lead is 10.7%.Here it may be noted that total contrast 10.7% comprises two parts: the absorption contrast of 8.9% and the phase contrast of 1.8%.By Wiener filtering, corresponding α=0.2 and α=0.1 respectively, average contrast brings up to 11.6% and 12.0%.Absorb contrast not change before and after Postprocessing technique, phase contrast brings up to 2.7% (α=0.2) and 3.1% (α=0.1) from original 1.8%.Tikhonov regularization method can obtain the contrast of 11.8%, comprising the absorption contrast of 8.9% and the phase contrast of 2.9%.And for Fourier-Wavelet Regularization, average contrast brings up to 12.7%, comprising the phase contrast of 3.8%.As can be seen here, the phase contrast (3.8%) that Fourier's Wavelet Regularization obtains, be greater than the phase contrast (2.9%) that Tikhonov regularization method obtains, be also greater than the phase contrast (2.7% and 3.1%) that Wiener Filtering obtains.

Net result shows, for under current engineering specifications, the deterioration benefit of line phase-contrast imaging system images result, a kind of Fourier towards line phase-contrast imaging of the present invention-small echo regularized image is adopted to recover new technology, effectively can improve the contrast of phase contrast image under degradation effects, and ensure that the fidelity of Recovery image.This technology, in raising Postprocessing technique contrast and signal to noise ratio (S/N ratio), all obtains than Wiener Filtering and the better result of Tikhonov method.The present invention is on the basis of line phase-contrast imaging, achieve and improve the effective mass of undesirable imaging results, therefore achievement in research of the present invention effectively can make up the deficiency for imaging performance in the line phase-contrast imaging system of early diagnosing mammary cancer in actual clinical.The application of the method, provides technical support by the diagnosis for effectively realizing early-stage breast cancer minute lesion tissue, provides powerful support for for the clinical practice and research of carrying out the line phase-contrast imaging of breast cancer in a deep going way provide.

Claims (2)

1., based on a line phase-contrast image recovery method for Fourier and small echo regularization, comprise the following steps:

1. the multiple image of line phase-contrast imaging system acquires edge of a knife utensil is utilized;

2. from the edge of a knife cross section curve of every width Image Acquisition diverse location;

3. edge of a knife cross section curve is averaged, then averaged curve is differentiated, obtain overall wire spread function h (x) of line phase-contrast imaging system;

4. place imaging object, gather image, obtain imaging results y (x);

5. calculate the Fourier transform of y (x) and h (x), obtain Y (ω) and H (ω);

6. formula is passed through obtain Fourier's regularization frequency-domain result α is regularization parameter, H ^*(ω) be H (ω) conjugate complex number; Fourier's regularization estimation of desirable phase contrast imaging result is then obtained by inverse Fourier transform $\hat{f}\left(x\right);$

7. to the wavelet coefficient signal w of desirable phase contrast f (x) _jx () is effectively estimated, obtain estimated signal

8. Fourier's regularization is estimated carry out Db6 orthogonal wavelet transformation, wavelet function is expressed as ψ, and scaling function is expressed as φ, obtains wavelet coefficient signal after wavelet transformation

9. to wavelet coefficient signal utilize Regularization function Λ _jx () carries out small echo regularization, obtain regularization result ${\tilde{w}}_j\left(x\right)={\hat{w}}_j\left(x\right){\mathrm{\Λ}}_j\left(x\right);$

10. right carry out the wavelet inverse transformation of wavelet function ψ, the final Fourier-small echo regularization estimated result obtaining desirable phase contrast imaging result

2. the line phase-contrast image recovery method based on Fourier and small echo regularization according to claim 1, is characterized in that, in the 7th step, estimate Fourier's regularization adopt Db4 orthogonal wavelet transformation, wavelet function is expressed as ψ ', and scaling function is expressed as φ ', obtains wavelet coefficient signal right carry out hard domains value, namely ${\tilde{w}}_j^{\′}\left(x\right)=\left\{\begin{array}{cc}{\hat{w}}_j^{\′}\left(x\right),& \left|{\hat{w}}_j^{\′}\left(x\right)\right|>{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\\ 0,& \left|{\hat{w}}_j^{\′}\left(x\right)\right|\≤{\mathrm{\ρ}}_j{\widehat{\mathrm{\σ}}}_j\end{array}\right.,$ Here the noise equation of wavelet function ψ ' under jth yardstick, ρ _jthat yardstick is correlated with the thresholding factor; Then right carry out the inverse transformation of Db4 orthogonal wavelet, obtain denoising estimated signal right carry out Db6 orthogonal wavelet transformation, obtain wavelet coefficient signal be w _jeffective estimated signal of (x).