Research on BOLD-fMRI Data Denoising Based on Bayesian Estimation and Adaptive Wavelet Threshold

The acquisition of functional magnetic resonance imaging (fMRI) images of blood oxygen level-dependent (BOLD) effect and the signals to be analyzed is based on weak changes in the magnetic field caused by small changes in blood oxygen physiological levels, which are weak signals and complex in noise. In order to model and analyze the pathological and hemodynamic parameters of BOLD-fMRI images effectively, it is urgent to use effective signal analysis techniques to reduce the interference of noise and artifacts. In this paper, the noise characteristics of functional magnetic resonance imaging and the traditional signal denoising methods are analyzed. The Bayesian decision criterion takes into account the probability of the total occurrence of all kinds of references and the loss caused by misjudgment and has strong discriminability. So, an improved adaptive wavelet threshold denoising method based on Bayesian estimation is proposed. By using the correlation characteristics of multiscale wavelet coefficients, the corresponding wavelet components of useful signals and noises are processed differently; while retaining useful frequency information, the noise is weakened to the greatest extent. The new adaptive threshold wavelet denoising method based on Bayesian estimation is applied to the actual experiment, and the results of OEF (oxygen extraction fraction) are optimized. A series of simulation experiments are carried out to verify the effectiveness of the proposed method.


Introduction
OEF (oxygen extraction fraction) refers to the ratio of oxygen uptake from the blood to the total oxygen content of arterial blood when oxygen-rich arterial blood flow passes through capillaries in a region of the brain, which reflects the activity of oxygen metabolism in the brain [1]. It is one of the three physiological indicators of energy metabolism in brain tissue, along with CBF (cerebral blood flow) and CMRO2 (cerebral metabolic rate of oxygen) [2]. OEF is also the difference in oxygen content between arterial and venous blood when tissues take oxygen from the capillary network for utilization. OEF indicates the ability of neurons to utilize oxygen and is an indicator of energy metabolism in brain tissue. Under normal circumstances, the metabolic oxygen supply and oxygen demand of neurons in the brain are in a dynamic equilibrium. When some abnormal conditions such as ischemic stroke and brain injury destroy this equilibrium relationship, physiological metabolic indicators such as cerebral blood flow and OEF will be abnormal. Therefore, the measurement of OEF is of great significance for the prevention and treatment of ischemic stroke. Over the past 25 years, the increase of OEF has been recognized as the gold standard to prove the existence of cerebral hemodynamic damage [3]. In addition, there is no significant difference in OEF values between gray matter and white matter, and the variability in the population is small, which is suitable as a measurement index. At present, the relatively mature OEF measurement technology is to measure the cerebral hemodynamic parameters through PET (positron emission tomography) technology [4]. However, PET technology is very expensive, and it also needs to use 15 O 2 as a radioactive labeling substance. The radioactive substance itself has radiation trauma effect on the tested human body, and the halflife of this radioactive substance is very short, only about 2 minutes, resulting in this technology not being easy to apply in the clinic. With the rapid development of magnetic resonance imaging (MRI) and blood oxygen level-dependent (BOLD) effects, more and more studies have begun to use BOLD for T2 * -weighted imaging to study and process cerebral blood perfusion and metabolism [5][6][7][8]. The difference between this method and other imaging methods is that MRI has the highest spatial resolution for brain lesions and is safer and more comprehensive in tissue structure, physiology, function, metabolism, and other aspects, especially in the diagnosis of neurological diseases.
In 1994, Haacke et al. [9] proposed a theoretical model to measure the signal decline caused by extracellular deoxygenated hemoglobin through a new MRI sequence. In 2000, An and Lin [10] proposed an MRI sequence that combined GRE and spin echo (SE) and solved the brain tissue OEF value through a mathematical model. In 2003, An and Lin continued to improve and reorganize the MRI sequence, trying to measure the OEF value from another perspective [11]. The sequence they proposed does not need to calculate R 2 , which can reduce the measurement error to a certain extent. Based on the research of An et al., researchers have developed a lot of OEF postprocessing software, which can directly output OEF measurement results after inputting the nuclear magnetic resonance images of GESSE sequences [12]. However, there are still few researches on measuring OEF with new nuclear magnetic resonance sequences or algorithms [13].
A mathematical model is established by using the functional magnetic resonance imaging of blood oxygen leveldependent (BOLD-fMRI) effect, which can accurately measure the OEF value of the human brain [14]. However, since the image signal is based on a small change in the physiological level of blood oxygen, the magnetic field is weakly changed, and the noise is complicated. So it is urgent to have effective signal analysis technology to reduce the interference of noise and artifacts [15,16] and improve the accuracy of OEF measurement. In this paper, an adaptive wavelet threshold denoising method based on Bayesian estimation is proposed on the basis of traditional wavelet denoising. After the wavelet decomposition, according to the different wavelet subband characteristics, the improved wavelet scaling function parameter equation is used to determine the optimal threshold for each wavelet decomposition scale level. The simulation results show that compared with the traditional wavelet hard threshold and soft threshold denoising, this method can better improve the SNR and reduce the mean square error of the denoising image and has better denoising effect, making the final measurement result of OEF more accurate.
After reconstruction, the NMR image signal is inherently complex. Two signals are collected through the receiving coil. After orthogonal detection, the two image signals have zero mean value with equal variance and Gaussian noise except the phase difference. Magnetic resonance image signals obey Gaussian distribution in the signal region with high SNR and Rayleigh distribution in the signal region with low SNR [17]. The noise is no longer completely independent but related to the signal.
The noise in NMR image signals obeys Rician distribution, as shown in Figure 1.
As can be seen from the figure, if the signal-to-noise ratio is high, the Rician distribution [18] will be close to the Gaussian distribution; if the signal-to-noise ratio is close to 0 (in this case, only noise exists), the Rician distribution will become a Rayleigh distribution. When the signal-to-noise ratio is low, it is signal dependent, and it is very difficult to remove the random variation and deviation rate of Rician noise generation. How to choose an effective denoising method to effectively separate the signal from the noise is also very meaningful and challenging.
The wavelet threshold denoising method [19,20] is to process the decomposed wavelet coefficients by finding appropriate thresholds in the appropriate method and to remove the wavelet coefficients belonging to the noise, leaving the wavelet coefficients belonging to the signal, and then reconstruct the processed wavelet coefficients to get the denoised signal. Classical threshold processing wavelet denoising has the following types: wavelet hard threshold denoising, soft threshold denoising, and soft and hard threshold denoising [21]. Figure 2, the hard threshold denoising is that after the wavelet coefficients are decomposed, the wavelet coefficients with absolute values greater than T are retained, and the wavelet coefficients with absolute values less than T are set to zero. Set the threshold T to 0.5. The specific expression is as follows [22]:

Wavelet Hard Threshold Denoising. As shown in
It can also be seen from the figure that the wavelet hard threshold denoising has a sudden change when the absolute value of the wavelet coefficient is equal to T, which will lead to excessive denoising and being unnatural [23,24]. Figure 3, the soft threshold denoising is that after the wavelet coefficient is decomposed, shrink the wavelet coefficients whose absolute value is greater than T to the zero point, and set the wavelet coefficient whose absolute value is smaller than T to zero. Set the threshold T to 0.5. The specific expression is as follows [22]:

Soft Threshold Denoising. As shown in
In this way, the coefficient change is more natural and softer than the hard threshold denoising [25,26].

Calculation Method of OEF.
In this paper, the twochamber model [12,27] is used to calculate the relevant parameters, and the OEF value is finally calculated. The two-chamber model treats the magnetized material and the tissue of the magnetized material as two chambers. In MRI 2 Oxidative Medicine and Cellular Longevity analysis of cerebral blood vessels, magnetized substance refers to hemoglobin in the capillaries and blood vessels; the tissue of the magnetized substance is brain substance.
In the BOLD effect fMRI image, the brain parenchyma corresponding to a pixel contains many capillaries. Assuming that the distribution direction of capillaries in the brain parenchyma is random, its magnetic moment can be decomposed into short-term scale S s ðtÞ and long-term scale S l ðtÞ: ρ is a constant, and λ is the relative volume of the magnetized material. δ w in equation (3) can be obtained from Among them, Hct is the hematocrit, B 0 is the intensity of the external magnetic field, Δχ is the magnetic susceptibility shift between deoxyhemoglobin and oxyhemoglobin, usually 0.18 ppm/Hct, γ is the magnetic susceptibility of the substance, and δ w is the frequency deviation of the magnetized substance. Calculate λ by short-term scale fitting of the brightness change curve in the image and then substitute λ into equation (4) to calculate δ w by long-term scale fitting, and finally, substitute δ w into equation (5) to calculate the OEF value [28].

Adaptive Wavelet Threshold Denoising Method Based on
Bayesian Estimation. Research on the denoising method of wavelet threshold has been active recently. Many new wavelet threshold denoising methods have been derived for the improvement of soft and hard thresholds and the different calculation methods of thresholds. Based on the Bayesian estimation, LkahwinderKaur proposed the NormalShrink wavelet threshold denoising method for the problem that the Donoho threshold method [29] cannot maximize the separation of image and noise at each level of wavelet. However, this method has a disadvantage; the image after denoising loses some details, and the edges appear blurry. In this paper, the wavelet threshold denoising method of NormalShrink is further improved. A wavelet adaptive threshold method based on Bayesian estimation [30] is used to denoise the image. After wavelet decomposition, according to different wavelet subband characteristics, the improved wavelet scale function parameter equations are used to determine the optimal threshold T that is appropriate for each wavelet decomposition scale.
In the process of denoising the image, if the prior information of the image can be considered and then the optimal threshold is obtained for the risk function, the result error will be reduced. The coefficients of each subband of image wavelet decomposition are basically symmetrically distributed near the zero point, forming a spike distribution, which can be described by a zero-mean generalized Gauss distribution GGD.
The general Gauss distribution is [31,32] where σ X is the standard deviation and β is the morphological parameter.
Suppose that X obeys a mean of zero and the variance is a Gaussian distribution of σ 2 X , that is, X ∼ Nð0, σ 2 X Þ, β = 2, then the Bayesian risk estimation function is For a given set of parameters, the threshold T is searched so that the result of γðTÞ in (4) is the minimum. T * ðσ X , βÞ = arg min rðTÞ T is used to represent the optimized threshold function.
Next, the optimization threshold T * is obtained: The standard density function is So  Oxidative Medicine and Cellular Longevity T B ðσ X Þ is an approximation of T * ðσ X , 2Þ, and its maximum deviation is not more than 0.01. The noise variance σ 2 in (10) can be estimated by using the median absolute variance for the highest frequency subband [28]: Finally, a data-driven, subband-based adaptive threshold is obtained as shown in [33] T

Experimental Method.
There were 12 young normal volunteers in this experiment, all aged 22-30 years old, with no gender selectivity. The 12 normal volunteers were scanned with GESSE sequence of brain parenchyma, and the scanning level was located above the lateral ventricle. Scanning parameters are as follows: tr = 1:5 s, TE = 56 ms, bandwidth 62.5 kHz, image matrix 256 * 256, gradient echo number 32, echo gap 1.5 ms, and layer thickness 7.5 mm. 32 DICOM images obtained from the GESSE sequence are used as the original data input.

Wavelet Threshold Denoising Method.
In this paper, three evaluation indexes of image display, SNR (signal-tonoise ratio), and MSE (mean square error) are considered in the comparison of various denoising methods. Each indicator verifies one of the advantages and disadvantages of the denoising method. Among them, the image display shows that the denoising image can be visually observed; the signal-to-noise ratio is aimed at retaining the signal while suppressing how much noise; the mean square error is suitable for expressing the sharpness of a picture.

SNR (Signal-to-Noise Ratio)
. SNR (signal-to-noise ratio) is one of the more traditional methods for measuring the amount of noise in a signal [34,35]. It is often used as an indicator for denoising effect evaluation. The unit is decibel (dB), which is defined as We can know that the larger SNR is after denoising, the better the denoising effect will be.

MSE (Mean Square Error).
The MSE (mean square error) shows the sharpness of the image [36] and the root mean square error between the original signal A and the denoised estimated signal B is defined as We can know that the smaller MSE is after denoising, the better the denoising effect will be.

Generation of Rician Noise.
Rician noise is not additive, but data dependent. First, a noiseless MR image A is defined on a discrete grid I, A = fa i | i ∈ Ig, and a set of random numbers is used as the brightness of image A. Taking σ as the standard deviation of Gaussian noise, two sets of Gaussian random numbers X = fx i | i ∈ Ig and Y = fy i | i ∈ Ig are formed, and the average of the two sets of numbers is 0 and has the same standard deviation σ. Then, the following M = fm i | i ∈ Ig is the Rician distribution:

Oxidative Medicine and Cellular Longevity
Bayesian estimation adaptive threshold wavelet denoising in MRI image signal denoising.

Generation of Rician Noise
(1) Load the Original MRI Signal. A signal containing Rician noise is generated using MATLAB simulation software, as shown in Figure 4. The left picture shows the original fMRI signal image, and the right picture shows the noisy signal image with Rician noise. And the standard deviation of noise is 0.009. It can be observed that the noisy image is slightly blurred. Since the standard deviation of noise is small, no obvious blur phenomenon is observed.
(2) Denoising Using Gaussian Filtering. The noise image is denoised by Gaussian filtering, which uses a twodimensional operator of size 5 * 5. The image after denoising by Gaussian filtering is shown in Figure 5.
(3) Denoising Using Traditional Wavelet Hard and Soft Thresholds. In the MATLAB simulation program, the original signal is decomposed by the "haar" wavelet function, and the signal is denoised by the soft threshold and the hard threshold, respectively. The result is shown in Figure 6.  (14) and (15), and the SNR and MSE of the denoised image are calculated according to the formula.

Analysis of Results.
Firstly, by observing the noise reduction comparison diagram of Gaussian filtering and three kinds of wavelet transform ( Figures 5-7), it can be seen that the signals obtained after noise reduction by wavelet transform have better similarity and smoothness. Compared with the traditional Gaussian filter denoising method, the signal denoising method based on wavelet transform can remove the noise in the signal more effectively.
The comparison of the advantages and disadvantages of each denoising method is shown in Tables 1 and 2:  Tables 1 and 2 show the SNR and root mean square error of the original signal after denoising the image with different standard deviation noises by Gauss filtering, traditional wavelet threshold transform, and new Bayesian estimation       7 Oxidative Medicine and Cellular Longevity adaptive threshold wavelet transform. According to the evaluation criteria of denoising performance, the signal-to-noise ratio of the denoised signal by wavelet transform is higher than that of the denoised signal by Gaussian filter, and the root mean square error is lower. The data show that the denoised signal by wavelet transform has better similarity with the original signal, and the denoised signal by wavelet transform retains more energy of the original signal.
In order to make a more intuitive comparison, this paper makes a broken line chart based on the table data, as shown in Figures 8 and 9.
It can be concluded that Bayesian estimation of adaptive wavelet denoising is superior to other denoising methods, and its application to fMRI image denoising can improve the degree of signal noise separation and retain useful signals to the maximum extent.
3.5. Simulation Experiment. The new adaptive threshold wavelet denoising method based on Bayesian estimation was applied to the actual experiment to try to optimize the results of mathematical modeling and measurement of OEF [15]. Firstly, the original postprocessing program was used to analyze the experimental data of 12 groups of volunteers, and the OEF results were obtained. Then, the optimized postprocessing program was used to analyze the experimental data of the volunteers to obtain new results. According to the gold standard measured by PET [37,38], the two groups of results were compared and analyzed, as shown in Table 3 and Figure 10.

Conclusions
In the postprocessing process of OEF mathematical modeling, the adaptive threshold wavelet denoising method based on Bayesian estimation retained more useful signals than the original Gaussian filter. The 12 groups of optimized OEF results in the experiment all have a certain degree of increase compared to the OEF values before optimization, making the results closer to the gold standard value of 35% in PET measurement. The experiment preliminarily verified the feasibility and superiority of applying the adaptive threshold wavelet denoising method based on Bayesian estimation to fMRI signal processing and analysis.

Data Availability
The dataset that supports the findings and conclusion of this study are available from the corresponding author on reasonable request. The data are not publicly available due to privacy.