Deep convolutional neural networks for bias field correction of brain magnetic resonance images

As a low-frequency and smooth signal, the bias field has a certain destructive effect on magnetic resonance (MR) images and is the main obstacle for doctors' diagnosis and image processing (such as segmentation, texture analysis, and registration). Before analyzing a damaged MR image, a preprocessing step is required to correct the bias field in the image. Unlike traditional bias field removal algorithms based on signal models and a priori assumptions, deep learning methods do not require precisely modeling signals and bias fields and do not need to adjust parameters. An MR image with the bias field is input and the corrected MR image is output after the deep neural network being trained on a large training set. In this paper, we propose taking the original image and the local feature images of the bias field in multiple frequency bands obtained by a Log-Gabor filter bank as input, correcting the bias field of a brain MR image through a deep separable convolutional neural network. Meanwhile, to speed up the training process and improve bias correction performance, we apply residual learning and batch normalization. We conducted the same test on BrainWeb simulation database and Human Connectome Project real data set, the consistency of qualitative and quantitative evaluation shows that our training model demonstrates better performance than the traditional state-of-the-art N4 and non-iterative multi-scale (NIMS) methods. Especially for the images with high-intensity non-uniformity level, the bias field has been well corrected.


Introduction
Among the many image acquisition technologies, magnetic resonance imaging (MRI), as a noninvasive medical imaging technology, has the characteristics of high resolution and high contrast in soft tissues and plays an important role in diagnosing the pathology of human internal tissues and organs [1]. Especially in brain image research, MRI is superior to other imaging techniques [2]. However, in the MRI imaging process, due to the influence of factors such as the uneven static magnetic field, the eddy current formed by the gradient field, and the uneven sensitivity of the receiving coil, MRI will have uneven intensity. This phenomenon is referred to as MR image artifacts or the bias fields. The bias field is a poor low-frequency smooth signal (as shown in Fig. 1), and it will change the intensity value of the pixel, which makes the same tissue in the image have a different intensity distribution. The existence of the bias field may affect the accuracy of clinical diagnosis and reduce the performance of quantitative image processing (such as segmentation, registration, and classification) algorithms. Therefore, before further processing, an image, the bias field of the MR image must be removed.
In past research, a variety of correction methods for the bias field had been proposed [3]. According to the different sources of the bias field, the correction methods are mainly divided into two types: prospective methods [4,5] and retrospective methods [6,7]. Prospective methods, such as model-based calibration [8], and multicoil imaging [9,10], mainly correct the bias fields caused by specific nuclear magnetic resonance equipment. This method can only handle the bias fields caused by equipment and environmental reasons, but cannot handle the bias fields introduced by the patients. Retrospective methods are image-based post-processing methods. They are correction methods proposed for the bias field generated by the shape, position, direction, and specific magnetic permeability of an imaged object. They do not distinguish the source of the bias field and rely solely on the image intensity and prior knowledge of an imaging object [11]. Retrospective correction methods can be roughly divided into the following five categories [12]: (1) Filtering. The non-uniformity field of the image is obtained by extracting the low-frequency components of the image. For example, Li et al. [13] used the homomorphic filtering method to correct the bias field in the magnetic resonance images. (2) Surface fitting [14,  Deep convolutional neural networks for bias field correction… 15]. Model non-uniformity fields by fitting a parametric smooth surface. For example, Zhao et al. [16] proposed to correct the intensity inhomogeneity of MR images by using a non-local mean pattern-matching algorithm to find voxels with similar structures to fit a smooth polynomial surface. (3) Segmentation based [17][18][19]. Tissue segmentation and non-uniformity correction are performed simultaneously, and the non-uniformity correction in the segmentation process is combined to effectively improve the segmentation performance. For example, Cai et al. [20] proposed an adaptive-scale active contour model based on image entropy and semi-naive Bayesian classifier to segment and estimate the bias field of non-uniform images simultaneously. (4) Statistical modeling. This approach is to incorporate prior knowledge into the energy minimization framework. For example, George and Kalaivani [21] imposed sparse constraints on uncorrupted images in the frame domain and inhomogeneity fields in the Fourier domain, and incorporated L1 norm to estimate the optimal true image and bias field simultaneously. (5) Histogram based [15,22]. The most representative method for correction directly based on image intensity histogram is N3 [23] and its extension N4 [24]. N4 is currently the most popular and widely used bias correction method, which uses a robust B-spline approximation procedure and an improved hierarchical optimization method for bias field correction [25]. In addition, George et al. [26] proposed a new method for correcting intensity non-uniformity using non-iterative multi-scale (NIMS) method in 2017. This method does not need to segment an image, nor does it need to have any prior knowledge of the scanner or the subject, and it extracts the bias field by using the Log-Gabor filter bank to perform multiscale analysis on the input image. These methods above mainly rely on the information of the acquired image and use prior knowledge about the spatial relationships of the imaged anatomical structure and/or intensity probability distribution to extract the hidden bias field in the image [12]. Retrospective methods usually make some assumptions about the acquisition process and can handle the bias field artifacts related to the scanner and anatomy. Compared with prospective methods that can only correct bias fields caused by MR scanning equipment, retrospective methods can also eliminate the bias fields introduced by patients.
The N4 method requires foreground extraction, and the bias removal efficiency may depend on the accuracy of background removal, the initialization of the width of the bias field histogram, and the spline distance [24]. It is time-consuming when it is executed in ITK (The Insight Segmentation and Registration Toolkit). The NIMS method does not require iteration, so the running time is relatively short, but it easily causes filtering artifacts in the high-contrast area of brain images, and it is not effective at removing the biases of brain images with serious pollution. One disadvantage of these traditional methods is that most of them involve complex optimization problems and reasoning can be time-consuming [26]. Recently, deep learning models have been successfully applied to various computer vision tasks and medical image computing tasks, and have achieved the most advanced results in both domains, mainly thanks to task-oriented feature learning and end-to-end efficient inference architectures. Inspired by this, it would be interesting to develop a deep learning-based method for predicting the bias field of each input MR image directly. Han et al. [27] proposed an unsupervised learning method to correct the inhomogeneity of MR images. This method is based on the Deep Image Prior (DIP) algorithm [28,29], which is not only very time-consuming but also unstable in optimization. Therefore, in this paper, we use deep learning to build a bias removal network (BIASNet) to eliminate the intensity non-uniformity in MR images.
In order to obtain more local contrast information on the bias field, before conducting network training, we first use the Log-Gabor filter bank to extract the local feature map of the bias field located in multiple frequency bands from a T1 image and take the local feature map and the original image as the input of the network. Then, in BIASNet, we introduce deep separable convolution to obtain more pixellevel features, residual learning and batch normalization are used to speed up the training process, and BIASNet can automatically conduct deep learning to obtain a more accurate bias field. Our BIASNet has conducted the same tests on the simulated database BrainWeb and the real dataset HCP. The qualitative and quantitative results show that our model significantly improves the accuracy and efficiency compared with the popularly available and state-of-the-art methods N4 and NIMS. Especially for image with high INU level, the bias field has been well corrected.
Overall, the main contributions of this paper are summarized as follows: 1. We propose a bias removal network BIASNet combined with dilated convolution and depthwise separable convolution, which can extract global and local context information and accurately predict the bias field under a significant speed gain; 2. Use the Log-Gabor filter bank to extract the local feature maps of the bias field in 3 frequency bands from the image as prior knowledge, so as to carefully estimate the local bias field of network training and improve the debiasing ability of BIASNet. 3. We design a mixed loss function by integrating mean squared error (MSE) and structural similarity (SSIM) loss functions. It can consider brightness, contrast, and structure indicators to obtain a more intuitive unbiased image; 4. BIASNet uses a residual learning mechanism to prevent gradient vanishing and gradient exploding problems.

Log-Gabor filter bank
The Gabor filter, as a statistical image local information extraction method, has good multi-channel and multi-resolution characteristics, but it is not a band-pass filter in the strict sense. When the bandwidth is greater than 1-octave, even part of the Gabor filter will produce non-DC components, so the bandwidth of the Gabor filter can only be limited to the 1 octave range [30]. If you want a greater spectrum coverage, you will need to employ many Gabor filters, which will increase the amount of calculations you will have to do. David J. Field introduced the Log-Gabor filter in 1978 to solve the nonzero DC component in the Gabor filter at a higher bandwidth [31]. In addition to the advantages of the Gabor filter, the Log-Gabor filter also compensates for some of its shortcomings [30]. (1) The Log-Gabor filter has no DC component and the bandwidth is not limited. (2) The transfer function of the Log-Gabor filter has an extended tail at the high-frequency end, which is more effective for natural image encoding than the Gabor function that over-represents low-frequency components. (3) The Log-Gabor function's measurement of the logarithmic frequency scale is consistent with that of the human visual system. Due to the singularity of the logarithmic function at the origin, the response function of the Log-Gabor filter cannot be represented in the spatial domain, but the filter can be constructed in the frequency domain. The polar coordinate expression of the two-dimensional Log-Gabor filter function in the frequency domain is [32]: where G is the radial component, G is the direction component, 0 and 0 are the center frequency of the filter and the direction of the filter, respectively. and are the frequency and direction of the image to be filtered, respectively, and 0 takes a value between 0 and 1. When constructing the filter bank, as 0 changes, to obtain a constant shape ratio filter, the value of must be adjusted to keep their ratio constant. The bandwidth B of the filter depends on the value of , and the bandwidth direction B is determined by ,B = 2 √ 2 ln 2 . Figures 2 and 3 are the images of the 2D Gabor function and the 2D Log-Gabor function in the frequency domain and the spatial domain, respectively. Figures 2 and 3 show that the 2D Log-Gabor filter and 2D Gabor filter have relatively similar function images in the spatial domain, which indicates that the 2D Log-Gabor filter not only inherits the multi-resolution and multi-channel advantages of the Gabor filter [33] but also compensates for its lack of bandwidth. By constructing the Log-Gabor filter bank, we can extract the local feature information of the image from filters with different directions and different center frequencies.
The extraction of the bias field has high requirements for frequency and spatial positioning. The Log-Gabor filter can provide any bandwidth to ensure better spatial positioning [34]. The Log-Gabor filter bank can achieve frequency positioning through a (1) set of tuned Log-Gabor filters. Therefore, we choose Log-Gabor function to construct a Log-Gabor filter bank to extract the local feature information of a bias field. The 3-octave bandwidth contains a large range of frequency components, which increases the difficulty of separating the bias field information from the decomposition band, and the 2-octave bandwidth can ensure that a smaller number of filters are used to effectively extract all frequency bands with bias fields [26]. Therefore, each Gaussian filter in the filter bank uses the best bandwidth of 2-octave bands, namely octaves. Therefore, in order to maintain the bandwidth of the 2-octave filter, when 0 changes, the scale factor must be adjusted to keep the ratio 0 constant. We construct a filter bank composed of three tuned discrete Log-Gabor filters, and the transfer function of the filter bank is represented in Formula (2). The filters are arranged in descending order of cutoff frequency. Among the three filters, the one with the smallest cutoff wavelength has the largest spectral radius in the frequency domain. Figure 4 shows the local feature maps obtained after applying the Log-Gabor filter bank to the T1w image from BrainWeb. Deep convolutional neural networks for bias field correction…

Theory and method
With a deeper network architecture, the deep learning-based method can expand the receptive field of a network to obtain a larger range of contextual information about an image, thereby obtaining a more accurate bias field. Convolutional neural networks (CNNs) have rapidly developed in image classification and segmentation since 2010 and have achieved great success [35]. The denoising convolutional neural networks (DnCNNs) using residual learning strategies proposed by Zhang et al. have achieved good performance [36]. A deep convolutional network for image super-resolution (SRCNN) proposed by Dong et al. realized the conversion from low-resolution images to high-resolution images [37]. At present, few papers use deep learning to remove the bias field in an image. We first introduce a simple and widely used mathematical model of a bias field [26]: is the bias field; and n(x) represents the high-frequency noise introduced due to the defect of the MRI scanner, which is considered to be independent of the bias field. Assuming that n(x) follows a Gaussian probability distribution [38,39], it can be removed by conventional low-frequency filtering. In the current work, no explicit preprocessing of additive noise is considered. Therefore, Formula (3) can be simplified as: u is the desired signal, f is the contaminated signal, and b is the bias field. The basic idea of DL (deep learning) based on a NN (neural network) is to use the following formula to obtain the predicted bias field signal: where Net represents the NN architecture and is a network parameter, including weights and biases. Therefore, the predicted real image can be obtained, and at the same time, the multiplicative noise contained in the image is removed. In this section, we will introduce the network architecture and how to optimize the parameters.

BIASNet architecture
Contextual information is very important for the reconstruction of pixels destroyed by the bias field. The deep learning-based method mainly captures more global contextual information by expanding the receptive field [40]. Although increasing the depth and width of the deep network can expand the receptive field, when the depth of the network is very large, continuing to increase the depth of the network may lead to the disappearance of the gradient, and the use of widening methods may generate more parameters, resulting in network overfitting. Therefore, we set the dilation rate to 2 in the convolutional layer of the convolution residual (CR) component or separable convolution residual (SCR) component of BIASNet and use the dilated convolution to extract global features, as shown in Fig. 5.
The CR component consists of the dilated convolution, residual learning, batch normalization (BN), and ReLU functions. The SCR component consists of separable expanded convolution, residual learning, batch normalization (BN), and ReLU functions. Residual learning fundamentally breaks the symmetry of the network and improves the ability to characterize the network. Moreover, a residual network is easier to optimize, and better accuracy can be obtained from the greatly increased depth [41]. By setting the BN layer, the generalization ability of the network is improved.
BIASNet does not directly output an unbiased image but is used to predict the bias ⌢ b ij at each pixel f ij of the original image. Therefore, according to Formula (4), the unbiased intensity can be calculated as f ij ⌢ b ij . The process of MR image bias removal under this architecture is shown in Fig. 6.
BIASNet has 4 input layers, which are the original image T1w and the feature images of scale 1, scale 2, and scale 3 from the Log-Gabor filter bank. Each input layer extracts features through a convolutional layer with a 3 × 3 convolution kernel and 64 filters and then uses an addition layer to perform feature fusion. Subsequently, 15 SCR components are cascaded. The convolution kernel of each SCR component is 3 × 3 and there are 64 filters. In order to expand the receptive field, the dilation rate is set to 2. Finally, through a 1 × 1 convolution layer, a bias field image of the same size as the input is generated, and an unbiased image can be obtained by using the formula. The combination of residual learning and batch normalization in the SCR (CR) component speeds up the training speed and improves the bias removal performance. Because the output should have the same size as the input, the network does not have a pooling layer.

Optimization of BIASNet
In order to optimize the parameter of the proposed BIASNet, the total loss function is defined as: where Loss mse is the mean squared error (MSE) function between the desired bias image and the estimated ones from the input contaminated images. It is defined as follows: Loss disssim is the average structural difference function between the estimated bias field image and the desired bias field image and is defined as follows: In Formulas (6) ~ (8), f i is the i-th image contaminated by the bias field, u i is the corresponding real image, represents the network parameter, and and are the hyperparameters that control the ratio of Loss mse and Loss disssim . In experiment, we set = 1 and = 0.01.
SSIM Net(f i ; ), ( f i u i ) is used to measure the structural similarity between the estimated bias field image and the desired bias field image. Its value range is [0, 1], and it is a full reference image quality evaluation index [42]. Its formula is defined as follows: where, C 1 and C 2 are constants. In order to avoid the denominator being 0, they are usually set as C 1 = (K1 * R) 2  (6) L = * Loss mse + * Loss disssim Deep convolutional neural networks for bias field correction…

Quantitative standard for bias removal
In order to evaluate the bias removal performance, in the experiments on the same dataset, our method was quantitatively evaluated and compared with the N4 and NIMS methods in terms of the structural similarity index(SSIM), coefficient of variation(CV), and joint coefficient of variation (CJV). The SSIM is a structural similarity index that measures the similarity of two images, and its definition is shown in Formula (9). It measures the similarity between the real image and the corrected image and judges whether the bias removal is good according to the degree of similarity.
CV is the coefficient of variation of the same tissue. The smaller the CV is, the better the uniformity of an image [12]; and its calculation formula is as follows: where, and respectively represent the variance and mean of the pixels in the area. In each original image and the corresponding corrected image, the CVs of the white matter and gray matter tissue regions were respectively calculated, and the quality of the bias removal in the same tissue region was determined according to the CV.
However, the coefficient of variation alone cannot provide information about the joint variation of the distribution of different tissues. When the CV of one-tissue region changes and the other tissue region remains unchanged, the effect of image non-uniformity correction is still not ideal, and the joint coefficient of variation (CJV) can overcome these limitations. The calculation of the CJV, which is a direct indicator of the degree of overlap between the two tissue regions, requires two different tissue regions; and the smaller the CJV is, the better the image uniformity correction effect is [40]. For two different organizational regions T 1 and T 2 , the calculation formula is as follows: BIASNet was implemented using the publicly available TensorFlow framework and Keras artificial neural network library with Adam optimization and trained using a deep learning acceleration computing server. The server is configured with a 3.20 GHz Core i5-6500H CPU, an NVIDIA GeForce RTX 2080 (11 GB) GPU and 32 GB RAM. During the training, we used Adam for stochastic optimization and the mini batch size was 10. For the training based on 2D slices, we trained 200 epochs for our model; and on the 3D images, we trained our model based on patches for 10 epochs. In the training process, the learning rate was set to 1e-3.
N4, NIMS, and BIASNet were all performed on the same platform. In the N4 method, we set the threshold of the stik.otsu threshold () function to 200.In the NIMS method, we set the minimum cut-off wavelengths corresponding to the BrainWeb and HCP datasets as 24 and 34, respectively, the multiplicative factor = 3 , the indirect measure of filter bandwidth k = 0.55 , and the S-Golay filter kernel sizes for BrainWeb and HCP datasets were chosen to be 52 and 75, respectively.

Dataset
To fully validate the bias field correction performance of the proposed BIASNet, we have applied our training and testing on two public datasets: BrainWeb and HCP.
The first dataset is the Simulated Brain Database (SBD, http:// brain web. bic. mni. mcgill. ca/ brain web/) [43]. SBD is a collection of realistic MRI data volumes generated by MRI simulators that have been widely used to evaluate the performance of debiasing methods [26,28,42].
The second dataset is The Human Connectome Project (HCP, https:// www. human conne ctome. org/) [44]. The HCP is an open-access dataset that contains multimodal brain imaging data for healthy young-adults subjects, and specific details can be found in reference [45]. We introduced this dataset to verify the performance of the proposed BIASNet on real human brain MR images.

Results from the simulation database BrainWeb
In this paper, we first studied the performance of the proposed method on the Simulated Brain Database (SBD) [43]. The size of each MRI brain image is 181 × 217 × 181 , and the resolution is 1 × 1 × 1 mm 3 . Here, 2D MRI slices were extracted from the volumes with different INU levels (20%, 40%, and 60%) on the axial plane and were trained and tested respectively. We use rotation, mirroring, translation, and other methods to expand the samples. Finally, for the volume of each INU level, 1359 slices were obtained for training, 151 slices were used for verification and 30 slices were used for testing. In order to further speed up the training process and obtain less redundant regions, we cropped the edge of each slice, and the slice size was 210 × 180.
In order to prove the effectiveness of the method proposed in this paper, we compared our method with the N4 and NIMS methods and performed experiments on replacing the SCR component in BIASNet with the CR component.
In Fig. 7, we drew a comparison of the intensity distribution of the simulated T1w image with 1% noise and 40% INU before and after correction in the 90th column (along the turquoise line). Figure 7a  Deep convolutional neural networks for bias field correction… Fig. 7c, d. This shows that the proposed correction strategy is effective for intensity non-uniformity correction while removing the multiplicative noise contained in the image and retaining the edges and finer details in the input image.
Next, we will quantitatively evaluate the performance of BIASNet (CR), BIAS-Net, and the traditional state-of-the-art N4 and NIMS methods on the BrainWeb database using SSIM, CV, and CJV. Figure 8 shows the SSIMs with different methods on datasets with different INU and noise levels from the BrainWeb database.
In Fig. 8, ni − j(i = 0, 1, 3;j = 20, 40, 60) refers to the dataset with noise level i% and INU j%. Figure 8 shows that the SSIM of the NIMS method is almost the lowest on all the tested datasets. This may be because the NIMS algorithm uses the Log-Gabor filter bank to perform multiscale analysis on the input image to extract the bias field, and filtering may affect the detailed characteristics of the image; therefore the SSIM of the NIMS algorithm is low. The SSIMs of BIASNet (CR) and BIASNet are not much different, and both are better than those of N4 and NIMS. As the INU level increases, this advantage becomes more obvious.
In the experiments, the performance of the proposed scheme was compared quantitatively with those of N4 and NIMS in terms of the CV and CJV. T1w images with three different noise levels (0%, 1%, and 3%) and three different INU amplitudes (20%, 40%, and 60%) from the BrainWeb database are considered for nonuniformity correction. Tables 1 and 2 give the comparison of the CVs of the gray matter region and the white matter region before and after the T1w image correction, respectively. Table 3 provides a comparison of the CJVs of gray matter and white matter regions before and after image correction. The lower the CV is, the smaller the scattering intensity of the segmented region, which in turn indicates that the uniformity of the same tissue region is higher. Tables 1 and 2 show that the performance of the proposed scheme is significantly better than that of the N4 and NIMS methods. Especially in the case of a high noise level and high INU amplitude, the performance advantage of our method is more obvious. Of course, effective bias field correction should increase the contrast of the image, which can be reflected by the appropriate separation of the peaks of the white matter and gray matter distribution or characterized by a sharp drop in the CJV. Table 3 shows that under high noise levels, the CJVs of the BIASNet (CR) and BIASNet are significantly reduced, indicating that the separation of the intensity level of the white matter and gray matter tissue regions is improved. The quantitative analysis of the SSIM, CV, and CJV shows that the scheme proposed in this paper is effective at correcting the bias field under the condition of a higher noise level and INU amplitude. Figure 9 shows the correction results obtained after applying the BIASNet method to BrainWeb simulation data with different noise levels and INU amplitudes. Figure 9a shows the T1w axial image degraded by 20% INU (1% noise), Fig. 9d shows the T1w axial image degraded by 40% INU (no noise), and Fig. 9b, c and Fig. 9e, f correspond to the predicted bias fields and corrected images, respectively. Through visual inspection, it can be clearly seen that the BIASNet method can effectively remove the bias field, and the intensity change of white matter tissue area is significantly reduced in the corrected image.

Results from the real database HCP
The model proposed in this paper is verified on real human brain MR images of the HCP [44] dataset. The dataset contains T1w and T2w structure MR images of healthy adults from a Siemens Skyra 3 T scanner and contains the bias fields of the corresponding images. In the experiment, we randomly selected the T1w structure images of 10 patients in the dataset. The size of each MRI brain image was 260 × 311 × 260 voxels and the resolution was 0.7 × 0.7 × 0.7 mm 3 . We randomly selected 1000 slices from these 10 brain image volumes. In order to train the model better, we expanded the data using rotation and mirroring methods, and finally obtained 2880 training samples, 320 verification samples, and 200 test samples. Similarly, in order to speed up the training process and obtain fewer redundant regions, we cropped the edges of the images, and the size of each sample was 300 × 240. Figure 10 shows the SSIMs obtained after correction with different methods. It can be seen from Fig. 10 that the proposed scheme in this paper also shows good performance on the real human brain dataset. Figure 11 clearly shows that in the original image, the intensity of the white matter region along the selected line highlighted in turquoise has changed very obviously due to the bias field. After correction using the BIASNet method, the white matter intensity variation of approximately 0.25 in the normalized original image is reduced to less than 0.15.
In the test samples, the CVs before and after correction were calculated with the proposed model to evaluate the intra-class changes in white matter and gray matter tissue regions. Table 4 shows the average CV of 200 test samples. Although both the N4 and NIMS algorithms reduce the CVs of gray and white matter in the images, the model proposed in this paper is more effective than these two algorithms. The average CJVs between white matter and gray matter are also listed in Table 5. Table 5 shows that after correction with the BIASNet model, the CJV between white matter and gray matter in the image drops to 64% of the initial value. In the same situation, the N4 and NIMS methods can only be reduced to approximately 77% of the initial value. Obviously, on the real human brain image dataset, the performance of the scheme proposed in this paper is also better than the N4 and NIMS methods in the three aspects of the SSIM, CV, and CJV. Figure 12 shows the corrected results of the BIASNet method on two randomly selected slices of different patients. It can be observed from Fig. 12 that the predicted bias field is smooth, and the intensity of each tissue region of the corrected image looks more consistent, which is obvious in the white matter regions in Fig. 12c, f Therefore, the effectiveness of the proposed scheme for the correction of intensity inhomogeneity in clinical data is visually verified.

Experiment results on 3D images
On the HCP [44] dataset, we further verified the model proposed in this paper on 3D data. From the dataset, we randomly selected T1w structure images of 179 patients, of which 139 samples were used for training, 20 samples were used for verification, and 20 samples were used for testing. In order to speed up the training process and obtain fewer redundant regions, we cropped each MRI brain image volume to 256 × 304 × 256 , and divided the image into patches. In order to better reduce the memory burden, we used 64 × 64 × 64 voxel patches. By using a sliding window strategy with a stride of 48 × 48 × 48 , a total of 20,850 patches were obtained for training, and 3000 patches were used for verification. In the test stage, in order to eliminate the patch boundary artifacts, we changed the stride of the sliding window to 24 × 24 × 24 , so that 17,820 patches could be obtained, and the trained network was applied to the patches of the test set. We performed an averaging operation on the predicted values of the overlapping regions of the patches. We trained our two different models for 10 epochs, with a learning rate of 1e-3, and optimized the models with Adam. Table 6 shows the experimental results of the quantitative evaluation of BIAS-Net based on two-dimensional slices and three-dimensional patches. The corresponding boxplots of the CV and CJV are shown in Fig. 13. It can be seen from Table 6 that in terms of CVs, BIASNet (2D) and BIASNet (3D) both show superior correction performance, especially for the visible high-intensity changes in the white matter tissue regions in the image that were effectively quantitatively corrected by reducing the CVs in the white matter. The CJVs between white matter and gray matter are also listed in Table 6. The CJV can better reflect the gray uniformity of an entire image. BIASNet (3D) considers more spatial information and shows better performance in the quantitative evaluation of the CJV. The boxplot in Fig. 13 depicts the same situation. Since BIASNet (3D) is trained based on  . 13 The boxplot of the CV for gray matter and white matter, and the boxplot of the CJV between gray matter and white matter 1 3 Deep convolutional neural networks for bias field correction… patches, there may be boundary effects after the patches are stitched, resulting in the SSIM of BIASNet (3D) being lower than that of BIASNet (2D). Figure 14 shows a comparison of the results of a patient after correction using BIASNet (2D) and BIASNet (3D). From a visual point of view, there is little difference between Fig. 14c, d, but from the intensity distribution of the 90th column, the intensity distribution of BIASNet (3D) fits better with that of the real image.

The influence of Log-Gabor features on bias removal
On the n1-20, n1-40, and n1-60 datasets of the BrainWeb database, we removed the log -Gabor features of three scales and used the BIASNet model for training. The experimental results of the test set are shown in Table 7. Table 7 shows that after adding Log-Gabor features of different scales in training, the CVs, CJVs, and SSIMs have been improved to varying degrees. The greater the INU amplitude is, the more obvious this advantage. When the INU amplitude is 20%, after the introduction of the Log-Gabor feature, the CJV decreases by 0.2; and when the INU amplitude is 60%, after the introduction of Log-Gabor feature, the value of CJV decreases from 51.6 ± 1.45 to 46.3 ± 2.67. Therefore, in the process of extracting the bias field, Log-Gabor features of different scales can provide additional local information on the bias field.

Discussion
In this work, we propose a bias removal network (BIASNet) combined with dilated convolution and depthwise separable convolution, which can extract global and local context information and accurately predict the bias field. Local feature maps of bias fields located in multiple frequency bands are extracted from contaminated images by a Log-Gabor filter bank as a prior knowledge to aid BIASNet in estimating reliable local bias fields. We design a mixed loss function based on mean squared error and structural similarity, which takes into account luminance, contrast, and structure metrics to obtain a more intuitive unbiased image. In addition, BIASNet employs residual learning mechanisms to avoid gradient vanishing and gradient exploding problems. Moreover, we designed SCR and CR components in BIASNet. Through comparison, we found that the SCR component can extract more local features and obtain more accurate bias fields.
The performance of BIASNet was tested on the BrainWeb database and HCP dataset, and the results showed that the performance of BIASNet was better than the performances of the traditional state-of-the-art N4 and NIMS methods. In the white matter and gray matter regions, the performance of the method proposed in this paper was quantitatively evaluated on the simulated BrainWeb database using the CV. The intraclass variation of the decrease in intensity level observed in the qualitative analysis was verified by experiments. Regardless of the noise level and INU amplitude, the experimental results show that the CJV is also significantly reduced. The quantitative evaluations of the SSIM, CV, and CJV given in Fig. 8, Tables 1, 2, and 3 reveal that the performance of BIASNet is better than those of the N4 and NIMS methods, especially in the case of higher levels of noise and INU in the simulated data. On the HCP real dataset, after image correction, the significant decrease in the ratio of the CV and CJV indicates that the interclass separation of white matter and gray matter tissue regions has been improved and the intra-class intensity changes have been reduced. The small improvement of the SSIM shows that the images corrected by our method are more structurally similar to the real images. Compared with N4 and NIMS, BIASNet produces higher SSIMs and lower CVs and CJVs. In addition, we also performed 3D experiments on the HCP dataset. Table 6 shows that compared with the slice-based 2D experimental results, the CJV of BIASNet (3D) decreases significantly, which indicates that the gray matter and white matter separation is better on the 3D level. However, the SSIM has decreased due to the boundary effect that occurs when the patches are stitched.
To evaluate the efficiency of the proposed BIASNet in testing, we randomly selected one-fold of the SBD database and the HCP dataset and performed N4, NIMS, and BIASNet to compare their execution times. The average execution times (± standard deviation) of N4, NIMS and the proposed BIASNet on these two datasets are listed in Table 8. N4 needs to be executed in ITK, which is time-consuming. NIMS does not require iteration and takes relatively short time. With a speedup factor of approximately ten times, BIASNet is clearly more efficient than N4 and NIMS approaches. As a result, BIASNet is more appropriate and efficient for large datasets than N4 and NIMS approach.

Conclusion
In this paper, we propose a new method for the bias removal of MR brain images based on deep convolutional neural networks. This bias field correction method based on deep learning does not require segmentation of the scanner or the subject of any prior knowledge, and it is a multi-input network combined with residual learning. In this work, one of our limitations is that the experiment on 3D images is performed based on patches, which are influenced by zero-padding in edge part during the convolution process, resulting in grid artifacts after patch stitching. Next, we will study setting different weights for each pixel in the patch to reduce grid artifacts during patch stitching, so as to improve SSIM values on 3D images. In addition, another key limitation of our method is the need for high-quality unbiased ground-truth images from real datasets, which are difficult to obtain in practical applications. How to improve the image analysis performance by combining prior knowledge of organ shapes and locations is not obvious in recent deep learning techniques [46]. Oktay et al. used a new regularization model to integrate anatomical prior knowledge into a CNN for end-to-end learning in 2018. They showed that the method can not only be easily adapted to different analysis tasks (such as image segmentation and enhancement) but can also improve the accuracy of the prediction model [46] Therefore, to overcome the limitations of our method, in future work, we will combine prior knowledge such as anatomy, residual parts, etc., and use special learning techniques such as transfer learning, to reduce future unsupervised bias removal tasks.
In conclusion, the multi-input bias removal convolutional neural network method based on deep learning proposed in this paper can not only accurately estimate the bias fields in MR images, but also effectively retain the detailed MR brain image structure, which can provide high-quality MR images with uniform intensity distribution for subsequent image analysis tasks (such as segmentation, registration, etc.).