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

Abstract

The invention belongs to the fields of biomedical engineering and medical imaging, and relates to a coaxial phase contrast image restoration method based on Fourier-wavelet regularization, comprising: collecting multiple images of knife-edge instruments and obtaining images of different positions from each image Knife-edge section curve to obtain the overall line spread function h(x) of the coaxial phase contrast imaging system; place the imaging object, collect the image, and obtain the imaging result y(x); calculate the Fourier transform of y(x) and h(x) Transform to obtain Y(ω) and H(ω); obtain Fourier regularized estimation of ideal phase contrast imaging results by inverse Fourier transform Efficient estimation of wavelet coefficient signal w _j (x) of ideal phase contrast f(x); estimation of Fourier regularization Carry out Db6 orthogonal wavelet transform and perform regularization to get Carry out wavelet inverse transform to obtain the estimation result. The invention can effectively improve the contrast of the phase contrast image under the deterioration effect, ensure the fidelity of the restored image, and can improve the image restoration contrast and the signal-to-noise ratio.

Description

On-Axial Phase Contrast Image Restoration Method Based on Fourier and Wavelet Regularization

Technical field

The invention belongs to the fields of biomedical engineering and medical imaging, and relates to an image restoration method.

Background technique

Breast cancer is currently the number one "killer" of women. According to the statistics of the World Health Organization, in 2000, about 1.2 million women in the world were diagnosed with breast cancer, and 500,000 people died of breast cancer. Breast cancer not only endangers the lives of patients, but also causes damage to female sexual organs, seriously endangering women's physical and mental health. Although China is not a country with a high incidence of breast cancer, its average annual growth rate has been 1-2 percentage points higher than that of high-incidence countries in the past ten years, with an annual growth rate of 3%-4%. Due to the large population in our country, the diagnosis and treatment of breast cancer has become an increasingly serious social problem that needs to be solved urgently. Early diagnosis is the key to solve this social problem and improve the survival rate and quality of life of patients.

The current main method of routine breast examination is molybdenum (rhodium) target X-ray mammography, but long-term clinical practice shows that this technology has major shortcomings in terms of sensitivity, specificity, safety and comfort: on the one hand, the technology has The missed detection rate is as high as 10-15%; on the other hand, the final biopsy diagnosis rate of positive diagnosis by this technique is 25-29%. Especially serious, because the breasts of young women are too dense, resulting in a serious loss of diagnostic accuracy.

Until the end of the last century, the theory of X-ray phase contrast imaging (XPCI) was proposed, which broke the traditional concept of X-ray imaging and brought a new dawn for the ideal anatomical imaging diagnosis technology of early breast cancer. . Studies have shown that under the same radiation dose, the contrast resolution of phase contrast imaging is about 10 times higher than that of traditional X-ray absorption contrast imaging, which significantly improves the image visibility of soft tissue imaging. At present, Europe, the United States, Japan, Australia and other countries are vigorously carrying out research on X-ray phase contrast imaging. X-ray phase contrast imaging has been widely regarded as a microscopic imaging technology that can bring about great changes in radiological diagnostic medicine.

Currently available X-ray phase-contrast imaging technologies can be divided into interferometry, diffraction enhanced imaging (DEI), in-line holography, and grating differential phase imaging according to different imaging principles. Contrast imaging (Differential phase contrast, DPC) and coded-aperture based x-ray phase contrast imaging (Coded-aperture based x-ray phase contrast imaging), these methods have great advantages in X light source, experimental equipment, detector and imaging performance indicators. s difference. Interferometric imaging and diffraction-enhanced imaging both require precise monochromatic crystals, which lead to at least three problems in their practical applications: 1. The monochromatic crystals in the system must be aligned with high precision, so the imaging system has no impact on the environment. The perturbation is very sensitive; 2. Monochromatic crystals mean monochromatic X-ray beams, which causes most light sources to fail to meet the luminous flux requirements of the imaging system. At present, synchrotron radiation sources are mainly used to solve this problem, but synchrotron radiation sources are costly. It occupies a large area, which is not conducive to clinical application and promotion; 3. The crystal will absorb part of the X-rays emitted from the measured object, so it is difficult to control the optimal imaging dose. The grating differential phase contrast imaging method is a phase contrast imaging technology developed in recent years. Because this method can use ordinary X-ray light sources, it has great clinical application potential, but there are still some technical problems that need to be solved urgently, such as , the DPC method needs to introduce high-precision parallel gratings (source gratings, phase gratings, and absorption gratings), which increases the requirements for system accuracy and stability; due to the limited light receiving angle, the system luminous flux is reduced, and the exposure time of the system inevitably changes. In addition, some external factors such as patient body temperature will also have a non-negligible impact on the accuracy of the grating. The coded aperture phase-contrast imaging method can use a relatively large coded aperture interval (1 to 2 orders of magnitude larger than the fence interval in the DPC method), which improves the system accuracy requirements in the DPC method, but it still faces the problem of Other issues that need to be addressed urgently. Coaxial phase contrast imaging can directly obtain the fine structure information inside the object without introducing additional optical devices, and its imaging method is also called free-space propagation phase contrast imaging (free-space propagation). Wilkins et al. It is proved that coaxial phase contrast imaging can adopt polychromatic X-ray light source, thus avoiding the problem of luminous flux of imaging system. On the other hand, the coaxial phase-contrast imaging system has a small investment and occupies no more area than conventional X-ray photography. Therefore, it is highly scalable and is expected to replace conventional X-ray photography as the main means of breast cancer screening. Among the current various phase contrast imaging techniques, the light sources used in the coaxial imaging method, grating differential phase contrast imaging method and coded aperture phase contrast imaging method are engineering feasible X light sources, which have relative advantages in clinical application. Among them, the coaxial imaging method can be regarded as the basis of the grating differential phase contrast imaging method and the coded aperture phase contrast imaging method. The latter two methods are equivalent to introducing an additional optical device (parallel grating or coded aperture ). It can be seen that the technical development of the coaxial imaging method can effectively promote the development of the grating differential phase contrast imaging method and the coded aperture phase contrast imaging method.

Up to now, coaxial phase-contrast imaging technology has been widely and continuously concerned in the field of radiology at home and abroad. This technology is considered to be one of the most suitable microscopic imaging technologies for clinical medical application transformation under the current conditions. It is worth mentioning that companies engaged in the development of X-ray phase-contrast imaging equipment have emerged abroad. For example, Japan’s Konica Minolta (http://konicaminolta.jp) was the first to use phase-contrast imaging technology for clinical breast diagnosis in the world, and Launched the world's first phase-contrast imaging mammography system based on common X-ray sources. But so far, the advantages of X-ray coaxial phase-contrast imaging in the clinical diagnosis of breast cancer have not yet been shown. Konica Minolta's phase-contrast imaging mammography system has better imaging quality than the existing traditional X-ray mammary digital imaging. The improvement of the system is not large, and the expected effect reported in the relevant literature has not been achieved. The recall rate and cancer detection rate obtained by using this system are not significantly different from those of traditional methods.

By analyzing the real physical process of coaxial phase contrast imaging and the phase contrast imaging model reported in the literature, it can be found that the key issues restricting the clinical application of coaxial phase contrast imaging technology are mainly reflected in the imaging system. Defects, such as the X light source is not an ideal point source, the performance of the detector is limited by factors such as its own resolution and point spread function, and there are various harmful noises in the system.

Aiming at the non-ideality of the detector, Olivo, a scholar from the University of London, UK, made a preliminary attempt by using the Wiener deconvolution method, but the simulation and test results of Olivo are aimed at the synchrotron radiation source system. The synchrotron radiation source system has a very good light source and a very high signal-to-noise ratio, but the system is expensive. For example, the Shanghai Synchrotron Radiation Source was built in 2010, with a total investment of 1.4 billion and an area of about 200,000 square meters. On the other hand, the quality of the restored image obtained by Olivo using the dimensional deconvolution method is still far from the actual demand. The results obtained by Olivo are as follows: ideally, the phase contrast imaging contrast is above 200%; due to the deterioration caused by the non-ideality of the detector, the phase contrast drops to only 3%; through the Wiener deconvolution method, the phase contrast can be Increased to 7.2% and 11.9%, but in the case of 11.9% with a higher phase contrast, the introduction of noise is larger.

Since micro-focus coaxial phase contrast imaging is an imaging technology more suitable for clinical application, it is of more significant research value and significance to improve the quality of phase contrast imaging for this imaging technology. At the same time, exploring a more effective new method for improving the quality of phase contrast has become one of the key issues facing the development and application of coaxial phase contrast imaging technology and systems, which is of great significance.

Contents of the invention

The gist of the present invention is to propose a method that can improve the contrast and signal-to-noise ratio of X-ray coaxial phase contrast imaging results, so as to solve the key problems faced by coaxial phase contrast imaging under the current engineering technology conditions: due to the current engineering conditions , micro-focus X light source, non-ideality of the detector and system noise cause the phase contrast image quality to deteriorate and the phase contrast to decrease. Technical scheme of the present invention is as follows:

A coaxial phase contrast image restoration method based on Fourier-wavelet regularization, comprising the following steps:

① Use the coaxial phase-contrast imaging system to collect multiple images of knife-edge instruments;

② Obtain the knife-edge section curves at different positions from each image,

③ average the knife-edge section curve, and then calculate the derivative of the average curve to obtain the overall line spread function h(x) of the coaxial phase contrast imaging system;

④ Place the imaging object, collect the image, and obtain the imaging result y(x);

⑤ Calculate the Fourier transform of y(x) and h(x) to obtain Y(ω) and H(ω)

⑥ Through the formula Get Fourier Regularized Frequency Domain Results α is the regularization parameter,

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

⑦ Effectively estimate the wavelet coefficient signal w _j (x) of the ideal phase contrast f(x), and obtain the estimated signal

⑧ Estimation of Fourier regularization Perform Db6 orthogonal wavelet transform, the wavelet function is expressed as ψ, the scale function is expressed as φ, and the wavelet coefficient signal is obtained after wavelet transform

⑨For wavelet coefficient signal Use the regularization function Λ _j (x) to carry out wavelet regularization, and obtain the regularization result ${\tilde{w}}_j\left(x\right)={\hat{w}}_j\left(x\right){\mathrm{\Λ}}_j\left(x\right);$

⑩ right Carry out the wavelet inverse transform of the wavelet function ψ, and finally obtain the Fourier-wavelet regularization estimation result of the ideal phase contrast imaging result

As a preferred implementation, in step 7, the Fourier regularization estimation Using Db4 orthogonal wavelet transform, the wavelet function is expressed as ψ', and the scaling function is expressed as φ', to obtain the wavelet coefficient signal right Perform hard localization, that is, ${\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 is the noise equation of the wavelet function ψ′ at the jth scale, and ρ _j is the scale-dependent threshold factor; then the Perform Db4 orthogonal wavelet inverse transform to obtain denoising estimated signal right Perform Db6 orthogonal wavelet transform to obtain wavelet coefficient signal That is, the effective estimation signal of w _j (x).

Adopting the new technology of Fourier-wavelet regularized image restoration for coaxial phase contrast imaging of the present invention can effectively improve the contrast of the phase contrast image under the deterioration effect and ensure the fidelity of the restored image. The invention achieves better results than the Wiener filtering method and the Tikhonov method in terms of improving image restoration contrast and signal-to-noise ratio. On the basis of coaxial phase contrast imaging, the present invention realizes the effective quality improvement of unsatisfactory imaging results, so the research results of the present invention can effectively make up for the coaxial phase contrast imaging system actually clinically used for early diagnosis of breast cancer. Insufficient imaging performance. The application of this method will provide technical support for effectively realizing the diagnosis of small lesion tissues of early breast cancer, and provide strong support for the in-depth clinical practice and research of coaxial phase contrast imaging of breast cancer.

Description of drawings

Figure 1. The performance transfer function of the detector.

Figure 2. Imaging results of 300 micron fiber in the case of distance from light source to object/distance from object to detector = 80cm/80cm, (a) two-dimensional image, (b) cross-sectional curve.

Figure 3. The cross-sectional view of the ideal imaging result of the 300 micron fiber under the condition that the distance from the light source to the object/distance from the object to the detector = 80cm/80cm.

Figure 4. Cross-sectional view of fiber image obtained after Wiener deconvolution (a) α = 0.1 (b) α = 0.07.

Figure 5. (a) The corresponding L curve when the regularization matrix is selected as the Laplacian operator, the optimal regularization parameter λ=2.6088 (b) corresponds to the regularization matrix and regularization parameters in Figure 5 (a), and the obtained The phase contrast image restoration results, the contrast is 15.9%.

Figure 6. Cross-sectional view of fiber image obtained with Fourier-wavelet regularization.

Figure 7 (a) Phase contrast images of six 0.5mm pencil leads obtained in the experiment (b) Wiener filtering results (α=0.2) (c) Wiener filtering results (α=0.1) (d) TikHonov regularization results ( e) Fourier-wavelet regularization results.

Figure 8 Average cross-sectional view of the first 2 pencil lead phase-contrast images.

Detailed ways

Aiming at the deterioration of phase contrast imaging results caused by detector, light source and system noise in the process of X-ray coaxial phase contrast imaging, the present invention proposes a Fourier-wavelet regularization method for coaxial phase contrast imaging, which can effectively suppress the above-mentioned The influence of the deterioration effect on the imaging result can improve the image contrast of the details of the phase contrast result, and at the same time ensure the fidelity of the restored image. The present invention will be described in detail below from several aspects in conjunction with the accompanying drawings and embodiments.

1 Digital X-ray Imaging System

The experimental imaging system is Pixarray 100 small animal digital radiography system, manufactured by BIOPTICS Company of the United States. The detector of the system is a 1024×1024 CCD array with a pixel size of 50 μm×50 μm and 14 gray levels. The horizontal and vertical spatial resolutions are 20 pixels per millimeter. The focal spot size of the X-ray tube is 50 μm. The full width at half maximum of the point spread function of the detector is 110 μm. In the experiment, the working voltage of the X-ray source is 33kVp, and the working current is 0.5mA. During phase contrast imaging, the distance from the X-ray source to the object and the distance from the object to the detector are both set to 80cm. Under this setting, the focal spot image formed by the light source on the detector has a full width at half maximum of 50 μm. The point spread function of the overall system is the convolution result of the point spread function of the detector and the point spread function of the light source. We adopted the method suggested by Olivo, a scholar at the University of London, to place the knife-edge object on the plane of the imaging object, so that we can directly obtain the line spread function curve of the entire system, and then rotate through the center of the line spread function to obtain the overall system point spread function.

For the convenience of explanation, the present invention first gives the phase contrast image deterioration under the one-dimensional model, and based on this, extends to the phase contrast image deterioration under the two-dimensional model, and gives the corresponding image restoration results.

2 Wiener Deconvolution Technology for Coaxial Phase Contrast Imaging

Taking the one-dimensional model as an example, for the micro-focus X-ray coaxial phase contrast imaging system, due to the non-ideality of the detector and the deterioration effect caused by system noise, the coaxial phase contrast imaging result y(x) we actually obtained can be Expressed by the following formula:

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

Among them, y(x) is the phase contrast imaging result obtained in the actual experiment, * is the convolution operator, x is the spatial position coordinate, f(x) is the coaxial phase contrast imaging result under ideal conditions, and LSF(x) is Detector line spread function, s(x) is the line spread function imaged by the light source on the detector, h(x)=LSF(x)*s(x) is the overall line spread function of the system, it is the detector line The convolution result of the spread function and the light source line spread function, n(x) is the system noise.

The deconvolution method can remove the deterioration effect of the imaging result caused by the system point spread function, so as to approach the ideal coaxial phase contrast imaging result. But the deconvolution method requires the system point spread function to be non-singular. When the point spread function of the system is singular, the deconvolution method will amplify the high-frequency part of the system noise, resulting in blurred image details and failing to obtain the required image restoration results.

Therefore, in the presence of system noise, our goal is to find an efficient estimate of the ideal imaging result f(x), namely:

$\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; *</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

here is the effective estimation result of f(x) under the minimum mean square error criterion.

Wiener deconvolution ends by providing the filter φ(x) required by equation (2). Mathematically, the Wiener deconvolution method is the application of Wiener filtering to the deconvolution problem in the presence of noise terms. The expression formula of Wiener deconvolution in the frequency domain is as follows:

$\mathrm{<span\; class="notranslate">\; \&Phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span>}{{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; /</span>{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

Where Φ(ω), F(ω), H(ω) and N(ω) are the Fourier transforms of φ(x), f(x), h(x) and n(x) respectively, H ^* ( ω) is the complex conjugate of H(ω).

Thus, an effective estimate of the ideal phase-contrast imaging result in the frequency domain can be obtained, namely

$\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; Y</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span>\mathrm{<span\; class="notranslate">\; \&Phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; Y</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span>\frac{{\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span>}{{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; /</span>{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$

here and Y(ω) is and the Fourier transform result of y(x).

Since F(ω) and N(ω) are generally unknown and cannot be measured in the actual system, the formula of Wiener deconvolution is often simplified, and the system characteristic related item |N(ω) ² /F(ω )| ² is replaced by a constant α, which is the regularization parameter:

$\mathrm{<span\; class="notranslate">\; \&Phi;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; h</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span>}{{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$

The choice of the optimal α value depends on the noise of the acquired image: a large α value better suppresses noise but leads to signal distortion. Although a small value of α can obtain a more accurate signal, the price is that more noise is introduced. Olivo used this Wiener deconvolution-based image restoration technique in its research.

3 Tikhonov regularization technique for coaxial phase contrast imaging

In order to illustrate the Tikhonov regularization technique for coaxial phase contrast imaging, we first transform the expression of formula (1) from the convolution form to the matrix form:

y=Hf+n (6)

Here, the lowercase bold letters y, f and n represent the vector form of y(x), f(x) and n(x) in formula (1), respectively, and the bold capital letter H represents h from the convolutional form (formula 1). (x) Convolution kernel matrix converted to matrix form (Equation 6). The specific form of H is:

$\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; =</span>\left[\begin{array}{cccc}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& <span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>\\ \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; )</span>& <span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>\\ <span\; class="notranslate">\; \&CenterDot;</span>& & & \\ <span\; class="notranslate">\; \&CenterDot;</span>& & & \\ <span\; class="notranslate">\; \&CenterDot;</span>& & & \\ \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>& \mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; N</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>& <span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>\end{array}\right]<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

The image y obtained by actually imaging the object is the result of the ideal result being deteriorated by the system transfer function and including the noise term n.

Both Equation (6) and Equation (1) can be referred to as deterioration models. Once the deterioration has been modeled, the next step is to develop a solution to it. As mentioned earlier, due to the existence of noise items, direct deconvolution cannot obtain the required image restoration results. Usually, this type of problem can be called an ill-conditioned problem. Using regularization methods can convert an ill-conditioned problem into a well-conditioned problem, so as to obtain the desired image restoration results.

The purpose of the regularization method is to introduce relevant information about the desired recovery result, so as to stabilize the problem and obtain an effective stable solution. For formula (6), various regularization methods can be adopted to realize the image restoration problem. Among them, the Tikhonov regularization method has good characteristics in terms of solution stability and effectiveness.

Based on the basic principle of the Tikhonov regularization method, considering that H in formula (6) is an ill-conditioned or singular matrix, in order to obtain the solution of the required characteristics, we adopt the regularization of the following minimization criterion:

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

Among them, the Tikhonov regularization matrix L and the regularization parameter λ need to be selected according to actual needs. In many cases, when the original signal is considered to be a substantially continuous signal, the Tikhonov regularization matrix generally uses a high-pass filter operator to improve the smoothness of the restored signal. Because the actual phase contrast image is smooth and continuous, we use the Laplacian high-pass filter operator as the Tikhonov regularization matrix, which has great advantages in image boundary restoration and restoration signal fidelity, so it is compatible with the present invention The method of comparison has better persuasiveness.

Corresponding to the selected regularization matrix L, we first need to determine the regularization parameter λ. There are many ways to choose regularization parameters, including cross-validation method, constrained maximum likelihood method and L-curve method. Among them, the L-curve method can most vividly express the trade-off selection process between the error residual ||Hf-y|| and the high-frequency noise ||Lf|| in the signal. In the L-curve method, the abscissa is ||Hf-y||, and the ordinate is ||Lf||, from which many points (||Hf-y||,||Lf||), different The points correspond to different regularization parameters λ, and a curve is obtained after fitting. Choosing a large regularization parameter λ can obtain a smaller half-norm ||Lf|| of the solution, but at the cost of introducing a lot of noise. While selecting a small regularization parameter λ can obtain a higher signal-to-noise ratio, but the degree of signal distortion increases. Therefore, for the trade-off between signal recovery accuracy and signal noise, choose the best balance point between the two, that is, the point with the largest curvature on the positioning curve, which is the position of the regularization parameter λ we need.

Using the Tikhonov regularization method for the deterioration results of coaxial phase contrast imaging, we can obtain numerical solutions, and the results are as follows:

$\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; f</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; h</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; L</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; h</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; )</span>$

By comprehensively restoring the image contrast and signal-to-noise ratio, the Tikhonov regularization method can obtain better results than the Wiener filtering method.

4 Fourier-wavelet regularization technique

The present invention aims to study the Fourier-wavelet regularization technique suitable for coaxial phase contrast image restoration. The technique includes two steps of Fourier regularization and wavelet regularization denoising post-processing. In the first step, the deterioration effect in the phase contrast imaging results is alleviated by inverting the operator operation and Fourier regularization operation, but this step also brings noise amplification effect. In the second step, the noise introduced in the first step is effectively reduced by the wavelet regularization method.

In the first step of the Fourier-wavelet regularization technique, a regularized sharpening filter, such as a Wiener filter (formula (5)), is usually used in the inversion operator operation and the Fourier regularization operation. In order to recover more phase contrast fringe information of the object boundary, a small regularization parameter α is generally adopted in the Fourier-wavelet regularization algorithm, which also leads to the cost of noise amplification.

In the second step of the Fourier-wavelet regularization technique, the wavelet regularization operation is performed on the boundary sharpening signal obtained in the first step, and the non-linear filtering based on wavelet transform is used to remove the sharpening signal caused by the regularization in the first step. noise amplified by the filter. This nonlinear filtering operation is based on adaptive thresholding in the wavelet domain. What is used here is redundant wavelet operation. Our estimated signal obtained by Wiener filtering Carry out wavelet transformation, the redundant wavelet scale function is represented by the symbol φ, and the wavelet function is represented by the symbol ψ. right After wavelet transform, a series of scale coefficient signals can be obtained and wavelet coefficient signal The subscript j here represents the scale of the wavelet transform. The purpose of wavelet regularization is to signal from wavelet coefficients In order to remove noise artifacts, the present invention adopts the wavelet domain Wiener filtering method to realize. Using wavelet domain Wiener filtering, the wavelet coefficient signal The ideal wavelet regularization implementation method is as follows:

${\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; \&Lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 10</span>}<span\; class="notranslate">\; )</span>$

${\mathrm{<span\; class="notranslate">\; \&Lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 11</span>}<span\; class="notranslate">\; )</span>$

Here Λ _j (x) represents the ideal regularization of wavelet coefficients, is the noise variance at position x at the j-th wavelet scale, and w _j (x) is the wavelet coefficient signal of the coaxial phase-contrast imaging result f(x) under ideal conditions.

Due to the spatial invariance of Wiener filtering, the amplified residual noise in the first step of the Fourier-wavelet regularization algorithm is stationary. At the same time, because the redundant wavelet transform also has space invariance, the residual noise in a certain wavelet space should also be stable in the whole signal. Therefore, wavelet regularization can be rewritten as

${\mathrm{<span\; class="notranslate">\; \&Lambda;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; w</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 12</span>}<span\; class="notranslate">\; )</span>$

here Represents the noise variance at the jth wavelet scale, which can be calculated by the following formula

${\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}\frac{{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; \&Psi;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{<span\; class="notranslate">\; (</span>{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; /</span>{<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; W</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; \&omega;</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 13</span>}<span\; class="notranslate">\; )</span>$

Here σ ² is the overall noise estimate, Ψ _j (ω) is the Fourier transform result of the wavelet function ψ at the jth wavelet scale, W _j (ω) is the Fourier transform of w _j (x). In the present invention, the results of Fourier regularization Take the smallest scale of wavelet transform, and then estimate the median value of the scale coefficient to obtain σ ² .

Λ _j (x) and cannot be obtained directly, because the ideal wavelet coefficient w _j (x) corresponding to the theoretically predicted phase-contrast imaging result f(x) cannot be obtained (formula (12), formula (13)). Therefore, we effectively estimate to replace w _j (x). The results can be regularized by Fourier The wavelet coefficients are obtained by hard domain value. In this step, another wavelet function needs to be used. In order to distinguish it from the previously used scaling function φ and wavelet function ψ, the scaling function used here is represented by the symbol φ′, while the wavelet function is represented by the symbol ψ′. Considering the wavelet function ψ′ _j at the jth scale, the wavelet coefficient signal Indicates the result of Fourier regularization The result of wavelet transform using wavelet function ψ′. Then, in the ψ′ wavelet function space, the wavelet coefficient signal Perform domain shrinkage, that is,

${\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span>& <span\; class="notranslate">\; |</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; ></span>{\mathrm{<span\; class="notranslate">\; \&rho;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}\\ \mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; ,</span>& <span\; class="notranslate">\; |</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; w</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&le;</span>{\mathrm{<span\; class="notranslate">\; \&rho;</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}}_{\mathrm{<span\; class="notranslate">\; j</span>}}\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 14</span>}<span\; class="notranslate">\; )</span>$

here is the noise equation of the wavelet function ψ′ at the jth scale, and ρ _j is the scale-dependent threshold value factor, which can be realized according to the Mallat wavelet selection method. The wavelet coefficient signal after shrinking the threshold value Perform the inverse wavelet operation to obtain the denoised effective estimation signal of the phase contrast imaging result f(x)

have We can use the scale function φ and the wavelet function ψ to perform wavelet transformation on it to obtain an effective estimate of the ideal wavelet coefficient w _j (x) Instead of w _j (x), Λ _j (x) and Thus, using formulas (10) and (11) to obtain the coefficient of wavelet regularization Finally passed to Perform wavelet inverse transform to obtain the Fourier-wavelet regularization estimation result of the ideal phase contrast imaging result f(x)

In order to compare with the Wiener deconvolution results adopted by Olivo and illustrate the advantages of the present invention, we adopt the criteria for judging the quality of phase contrast restored images in Olivo's research results. The phase contrast calculation formula used by Olivo is as follows

$\mathrm{<span\; class="notranslate">\; Comtrast</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; max</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; min</span>}}}{{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; background</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 15</span>}<span\; class="notranslate">\; )</span>$

Here I _max is the maximum value of the positive overshoot peak in the phase contrast result, I _min is the minimum value of the negative overshoot valley in the phase contrast result, and I _background is the background intensity outside the phase contrast overshoot details.

5 Application process of Fourier-wavelet regularization technology for coaxial phase contrast imaging

The flow process of the new method for phase contrast image recovery based on Fourier-wavelet regularization of the present invention is described as follows:

1. X-ray coaxial phase contrast imaging parameter setting: In the present invention, the distance from the light source to the object is set to 80 cm, and the corresponding distance from the object to the detector is 80 cm.

2. Set the exposure parameters of the digital radiography system, place the knife-edge instrument at a distance of 80cm from the detector surface (that is, it is located on the plane of the imaging object), continuously collect 15 images, and obtain 50 knife-edge section curves at different positions from each image, and then 15*50 knife-edge cross-section curves are averaged, and then the derivative of the average curve is calculated to obtain the overall linear diffusion function h(x) of the system.

3, place imaging object, adopt the pencil lead of 300 micron polyethylene fiber and 0.5mm diameter among the present invention. Under the setting of imaging parameters (distance from light source to object/distance from object to detector=80cm/80cm), image the object and obtain the imaging result y(x).

4. Calculate the Fourier transform of y(x) and h(x) to obtain Y(ω) and H(ω). by formula Get Fourier Regularized Frequency Domain Results Then the Fourier regularized estimation of the ideal phase contrast imaging result f(x) is obtained by inverse Fourier transform

5. Effectively estimate the wavelet coefficient signal w _j (x) of the ideal phase contrast f(x), that is, take the estimated signal to replace w _j (x): for Fourier regularization results Using Db4 orthogonal wavelet transform, the wavelet function is expressed as ψ', and the scaling function is expressed as φ', to obtain the wavelet coefficient signal right Perform hard localization, that is, ${\tilde{w}}_j^{\&prime;}\left(x\right)=\left\{\begin{array}{cc}{\hat{w}}_j^{\&prime;}\left(x\right),& \left|{\hat{w}}_j^{\&prime;}\left(x\right)\right|>{\mathrm{\&rho;}}_j{\widehat{\mathrm{\&sigma;}}}_j\\ 0,& \left|{\hat{w}}_j^{\&prime;}\left(x\right)\right|\&le;{\mathrm{\&rho;}}_j{\widehat{\mathrm{\&sigma;}}}_j\end{array}\right.,$ then to Perform Db4 orthogonal wavelet inverse transform to obtain denoising estimated signal right Perform Db6 orthogonal wavelet transform to obtain wavelet coefficient signal That is, the effective estimation signal of w _j (x).

6. Estimation of Fourier regularization Carry out Db6 orthogonal wavelet transform, the wavelet function is expressed as ψ, and the scaling function is expressed as φ. Obtain wavelet coefficient signal after wavelet transform

7. For wavelet coefficient signals Carry out wavelet regularization, the regularization function adopts Λ _j (x), and obtains the regularization result

right Carry out the wavelet inverse transform of the wavelet function ψ, and finally obtain the Fourier-wavelet regularization estimation result of the ideal phase contrast imaging result

The present invention adopts the Pixarray 100 digital radiography device for small animals produced by BIOPTICS Company of the United States to construct a coaxial phase contrast imaging system. Firstly, the transfer function of the system detector is obtained through the knife-edge device. Figure 1 shows the detector transfer function measured by the knife-edge method, and the full width at half maximum of the system transfer function obtained from the measurement is 120 microns.

Fig. 2(a) shows the imaging result of the 300 micron fiber obtained by the imaging system under the condition that the distance from the light source to the object/the distance from the object to the detector=80cm/80cm. Due to the poor observability of the details in the two-dimensional image, we give the cross-sectional curve corresponding to Fig. 2(a) (see Fig. 2(b)). It can be seen from Figure 2(b) that due to the deterioration effect of the system transfer function, the actual contrast is only about 2.8%, and there is relatively serious system noise (the signal-to-noise ratio of the calculated signal is 16.2dB).

For the convenience of analysis, we also simulated the one-dimensional cross-sectional view of the ideal result of coaxial phase-contrast imaging of 300-micron polyethylene fibers under the above system settings, as shown in Figure 3. It can be seen from Figure 3 that for 300 micron polyethylene fiber imaging, the contrast of the ideal phase contrast imaging result can reach more than 200% under the condition that the distance from the light source to the object/the distance from the object to the detector = 80cm/80cm. Compared with Fig. 2, due to the deterioration effect of the detector, the contrast of the fiber image drops to only about 2.8%.

We first adopted the Wiener deconvolution method in Olivo's research, and set the Wiener filter parameters α=0.1 and α=0.07. Through Wiener deconvolution, we obtain a cross-sectional view of the image restoration results shown in Fig. 4. After Wiener filtering, when using α=0.1, the image contrast increases from the original 2.8% to 8.7% (as shown in Figure 4(a)), when using α=0.07, the image contrast increases from the original 2.8% to 16.6% . But when α=0.07 is used, the background noise of the restored image is much higher than when α=0.1. This is consistent with the results given by the Olivo study.

Then we use the Tikhonov regularization method for coaxial phase contrast imaging for phase contrast image restoration. For the above-mentioned deteriorated image (as shown in Figure 2), the Laplacian operator is selected as the regularization matrix, and then the corresponding optimal regularization parameter λ is selected by using the L-curve method. Figure 5(a) shows the L curve obtained when the regularized matrix identity matrix is selected, the error residual norm ||Hf-y|| is the abscissa, and the solution norm ||f|| is the ordinate , the curves are given in log-logarithmic form for easy observation and analysis. It can be seen from the figure that the optimal regularization parameter λ corresponds to the inflection point of the L curve, that is, λ=2.6088. Fig. 5(b) shows a cross-sectional view of the corresponding image restoration result, and it can be obtained by calculation that the contrast after image restoration is increased to 15.9%.

Figure 6 shows the image restoration results obtained by using the Fourier-wavelet regularization method, and the image contrast is increased from the original 2.8% to 18.5%. In the first step of the Fourier-wavelet regularization method, during Fourier regularization, we adopted a Fourier regularization parameter α=0.03, which is smaller than α=0.07 used in Wiener filtering. Compared with Wiener filtering and Tikhonov regularization methods, it can be seen that the Fourier-wavelet regularization method can obtain the highest restoration contrast and the best signal-to-noise ratio.

In order to further compare the feasibility and effectiveness of the above three methods in phase contrast image restoration, we performed phase contrast imaging on another object, and then used the above three methods to obtain image restoration results. Figure 7(a) is the original phase-contrast image of six 0.5mm-diameter pencil leads placed in sequence, and the system settings are consistent with Experiment 1. By calculation, the signal-to-noise ratio of the image is 18.7dB. Because the pencil lead itself has a certain absorption, the dark fringes of the phase contrast fringes are covered by the absorption contrast, and only the bright fringes can be seen. Figure 7(b) and Figure 7(c) are the results obtained when Wiener filtering is used, and the filtering parameters are set to α=0.2 and α=0.1 respectively, and the signal-to-noise ratios of the corresponding images are 25.4dB and 21.3dB respectively. Figure 7(d) is the phase contrast image deconvolution result obtained by Tikhonov regularization method. The optimal regularization parameter λ=6.2357 is obtained by the L-curve method, and the signal-to-noise ratio of the finally obtained image is 24.2dB. The image restoration result obtained by Fourier-wavelet regularization is shown in Fig. 7(e), and the signal-to-noise ratio of this figure is 30.1dB.

In order to effectively compare the improvement of image contrast by the above three methods, we give the average cross-sectional view of the first two pencil leads in Figure 7. The average contrast of the original phase contrast image of the pencil lead before image restoration was 10.7%. It should be pointed out here that the total contrast of 10.7% includes two parts: 8.9% of absorption contrast and 1.8% of phase contrast. Through Wiener filtering, corresponding to α=0.2 and α=0.1, the average contrast is increased to 11.6% and 12.0%. The absorption contrast did not change before and after image restoration, and the phase contrast increased from the original 1.8% to 2.7% (α=0.2) and 3.1% (α=0.1). The Tikhonov regularization method can obtain 11.8% contrast, including 8.9% absorption contrast and 2.9% phase contrast. While for the Fourier-Wavelet regularization method, the average contrast increases to 12.7%, including 3.8% phase contrast. It can be seen that the phase contrast (3.8%) obtained by the Fourier wavelet regularization method is greater than that obtained by the Tikhonov regularization method (2.9%), and is also greater than the phase contrast obtained by the Wiener filter method (2.7% and 3.1%).

The final results show that, aiming at the deterioration benefit of coaxial phase contrast imaging system imaging results under the current engineering conditions, a new technology of Fourier-wavelet regularized image restoration oriented to coaxial phase contrast imaging of the present invention can effectively improve The contrast of the phase contrast image under the deterioration effect, and the fidelity of the restored image is guaranteed. This technology has achieved better results than the Wiener filtering method and the Tikhonov method in terms of improving image restoration contrast and signal-to-noise ratio. On the basis of coaxial phase contrast imaging, the present invention realizes the effective quality improvement of unsatisfactory imaging results, so the research results of the present invention can effectively make up for the coaxial phase contrast imaging system actually clinically used for early diagnosis of breast cancer. Insufficient imaging performance. The application of this method will provide technical support for effectively realizing the diagnosis of small lesion tissues of early breast cancer, and provide strong support for the in-depth clinical practice and research of coaxial phase contrast imaging of breast cancer.

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{\&Lambda;}}_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^{\&prime;}\left(x\right)=\left\{\begin{array}{cc}{\hat{w}}_j^{\&prime;}\left(x\right),& \left|{\hat{w}}_j^{\&prime;}\left(x\right)\right|>{\mathrm{\&rho;}}_j{\widehat{\mathrm{\&sigma;}}}_j\\ 0,& \left|{\hat{w}}_j^{\&prime;}\left(x\right)\right|\&le;{\mathrm{\&rho;}}_j{\widehat{\mathrm{\&sigma;}}}_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).