Magnetic Resonance Imaging (MRI) is a prevalent non invasive imaging technique that produces high contrast anatomical images without ionizing radiation. However, due to long acquisition times, MRI scans are prone to noise and artifacts corruption. In this paper, we address the denoising problem by leveraging the intrinsic nature of the K-space domain where the noise naturally occurs during acquisition. We propose a light, eficient U-Net architecture that is specifically tailored to operate directly in the K-space. The model considers a residual learning-based estimate of the noise component across a range of noise levels and distributions and in particular for additive Gaussian noise. We also propose a SNR degradation based progressive training scheme that greatly improves performance across a wide range of noise levels. The network is computationally cost-efective and can be run on CPU in acceptable inference time, making it suitable for real time or resource constrained applications.
MRI is a powerful medical imaging technique that provides detailed images of the internal structures of the body of a patient. The main benefit of MRI is the possibility to acquire images in multiple planes (sagittal, coronal, and axial) without repositioning the patient. MRI provides a very high contrast images of soft tissues, using their water content and molecular properties. It is also non invasive because it doesn’t require X-rays or other radiations that could be harmful for the patient.
When a patient is placed inside an MRI scanner, the machine generates a magnetic field that aligns the protons in fat and water molecules along the direction of the field. Once the alignment is established, the scanner applies radio frequency pulses at the natural frequency of the protons. The energy stored in the pulses excites the alignment, and the protons precess along the axis of the magnetic field. When the pulses stop, the protons begin to realign themselves to their original configuration. As they realign, they release energy in the form of radio frequency signals, which are picked up by special receiving coils. They do not appear directly as an image; they are in K-space, a domain in which data is recorded in terms of spatial frequencies. To obtain a doctor-readable image, an inverse Fourier transform is then performed on this complex valued data, resulting in high-resolution grayscale images that relate to the anatomical structure of the tissues being examined.
Although magnetic resonance imaging (MRI) is a valuable tool for medical imaging, the resulting scans are inherently noisy and present some challenges. The acquisition process is time-consuming and expensive, so the patient must remain still throughout the scan, otherwise the resulting images will contain motion artifacts.
Additionally, the MRI machine is sensitive to external factors such as temperature and electrical pulses, which can cause additional noise to be introduced. The thermal mobility of protons can potentially introduce noise into the image background, especially in low-signal regions.
Mathematically, certain probability distributions can be used to characterize noise in MRI scans. The equation 1 describes the bell-shaped symmetric probability distribution of Gaussian noise, often known as normal noise. This type of noise is additive, defined by a constant standard deviation across the image, and is typically attributed to thermal motion and electrical noise.
, () =
1 √2 − (− )2 2 2
Rician noise is prominent in the low signal-to-noise ratio (SNR) regions of the scans. It follows a Rician distribution (equation 2) and emerges when the magnitude of the MRI signal is nonnegative and follows a Rayleigh distribution, while the phase is uniformly distributed. This type of noise is particularly dificult to handle compared to other types of noise.
, () = − (− )2 2 2 2 0 ︁( )︁ 2 (1) (2) creating variability in the returned signal. It follows a They reduce noise caused by rapid change or motion but Poisson distribution (equation 3). occasionally perturb content as well.
Non-linear filters are applied when the noise is not () = − (3) uniformly distributed. They rely on adaptive kernels ! or other rules of convolution. They are more advanced
Finally, the T1 and T2 relaxation phases are usually and can function diferently depending on local intensity
subject to exponential noise. This type of noise reflects patterns.
the intrinsic unpredictability of the decay of MRI signals The second is anisotropic difusion filtering, which is a
over time and has an exponential distribution (equation variation of standard spatial filtering in that it adapts its
4). behavior to the local image properties. The algorithm
ad () = − (4) justs the number of filtering iterations depending on the
intensity gradient in the neighborhood: more iterations
A modern and popular approach to deal with long ac- where there are smooth transitions, fewer where there
quisition processes is to take a shorter one and capture are sharp transitions. This technique preserves edges
less details, this involves acquiring only a portion of the more efectively than linear filters, but certain details
K-space data which corresponds to a subset of spatial along boundaries are still lost.
frequencies. With this process, the generated images The second technique is the non-local means filter,
contains more artifacts and more noise, but reliable and which better exploits image redundancy. In this case,
powerful deep learning models can enhance these scans. the algorithm treats larger structures as meaningful
feaThis paper is focused on the implementation of a model tures and smaller patterns as noise that should be
refor denoising that can work with diferent levels of noise; moved. Each pixel is adjusted by averaging values from
in particular the focus is on denoising MRI scans with an a broader neighborhood area weighted by structure
simiadditive Gaussian noise at the K-space level. The model, larity. This technique overcomes some of the deficiencies
which follow a residual approach, should adapt also to of anisotropic difusion but at the cost of more
compuother noise probability distributions. tation, since it involves a search in large neighborhoods
for each pixel.
2. Related Works Finally, methods that combine domain and range
filtering consider both spatial proximity and intensity
similarFrom image pixel analysis to contemporary deep learn- ity. A highly successful algorithm in this category is the
ing models, a number of methods have been developed bilateral filter, initially derived by Tomasi and Manduchi
recently to solve the issue of denoising in magnetic reso- [
Linear filters are ideal if the noise is evenly distributed. denoising by first converting the image from the spaThey are realized with fixed smoothing kernels and are tial domain to a new representation, where noise and not dependent on the content. There is a subgroup that appropriate image structures can be separated more easacts in the spatial domain: the kernel acts on each pixel of ily. Standard filters or operations specifically designed a 2D neighborhood. The goal is to smooth, detect edges, for the operation can then be applied to the new domain. or sharpen. Traditional examples include the Gaussian The main classes within this class difer depending on the iflter, the Sobel operator, and the Laplacian kernel. Linear type of transform used and the character of the resulting iflters remove Gaussian noise, but blur details as well. representation.
The second type is temporal. There, the kernel is used One of the simplest methods is frequency domain dein the time domain. There, the kernel is applied along noising, which applies the Fourier transform to represent time, on the same spatial pixels across two frames. It the image in terms of spatial frequency. Noise, being detects temporal changes, the most widely used ones high frequency by nature, can be minimized in this doare temporal averaging, diferencing, and optical flow. main without losing low-frequency data. This method is usually efective when the noise intensity is moderate, coder presented by Hinton and Salakhutdinov [19]. The but breaks down when the noise increases. main application is to learn a latent and more compact
The wavelet transform is an improvement over this representation of an image by downgrading it in the first method, as it decomposes the image into sets of frequency half of the model and then reconstructing it back in the bands at diferent scales and allows simultaneous local- second half; this process can be adapted to denoising by ization in both space and frequency. The multi resolution passing noisy images as input and expecting as output structure allows for better separation between noise and their denoised version. Similar architectures have also signal at diferent frequency levels. However, the wavelet been applied in other areas of the medical sector, such basis can be challenged when representing curved or as in the processing of speech signals for the automatic highly textured structures, which limits its applicability detection of speech disorders [20]. Some years later, Ronin some medical imaging scenarios. neberger et al. [21] proposed a new architecture, the
This drawback was addressed with the creation of the U-Net, consisting in a fusion between the Autoencoder curvelet transform, a generalization of wavelets specif- and the ResNet, that achieved good results in biomedical ically designed to produce a coarser representation of image segmentation. Also the U-Net can be adapted to curves and edges. It is well suited to handle long and perform general image denoising and so in particular anisotropic features, and is ideal for images with highly MRI denoising. directional textures. Despite this advantage, the trans- Although the deep learning MRI denoising perforform produces artifacts in smooth regions, where the mance was adequate, a new problem arose. MRI scans redundancy of the representation is lower. with noise distributions other than Gaussian could not
The contourlet transform takes this concept a step be accommodated by the networks. This drawback stems further by combining multiscale and directional decom- from the straightforward methodology of the models, position. It starts with a Laplacian pyramid for coarse which learned to produce a clean image from a noisy imimage structures and then uses directional filtering to age. This indicated that the models had only worked on capture edges at diferent orientations. This produces a denoising one distribution, or perhaps a limited number more expressive and flexible representation, especially of related distributions. Zhang et al. [22] proposed a new for contours and geometric shapes. However, the in- approach that they called residual learning in which the creased complexity of this method leads to increased model learns to predict the residual image, i.e. the noise computational costs, especially for the fine-grained edge from the noisy image, with the possibility of focusing on preservation required at multiple scales. more and diferent noise distributions. Residual learning also introduced regularization in the training process and boosted the image denoising performances.
In recent years, machine learning and then deep learning
started to become popular, because of the better
computational performance ofered by modern GPUs and the
brilliant results obtained. The first improvement in the
use of deep learning for images was the introduction of
hybrid neural networks models [
Afterwords, researchers believed that increasing the number of convolutional layers in the networks would have lead to better performances, but it was not true because the models sufered from the vanishing gradient problem. To overcome this issue, He et al. [18] proposed ResNet, a new network with the so called residual layers where the new outputs are computed with an additive update from the previous inputs. Moreover, ResNet allowed to obtain better performances and a more stable training with less convolutional layers.
ResNet was also one of first deep learning models used for general image denoising, together with the Autoen
(IXI) [23], which includes approximately 600 MRI scans in NIFTI format. Data were acquired from healthy subjects using diferent acquisition protocols (T1, T2, PD, MRA, DTI) at three hospitals in London, namely Hammersmith Hospital (Philips 3T), Guy’s Hospital (Philips 1.5T), and the Institute of Psychiatry (GE 1.5T). Each scan is associated with protocol specific parameters, publicly available on the dataset website.
procedure ofers several advantages. 2D image process- niques were not used. This is because the U-Net model ing is faster and more suitable for real-time situations, already tends to adjust the activations thanks to the conas it uses less processing resources than full-volume pro- catenation mechanism between the old and new features cessing. A more focused investigation is also possible, in the upsampling phase. Another reason is related to as regions of interest are often located in a single plane, eficiency, the network remains faster and, in the tests such as axial, coronal or sagittal. This makes it possible carried out, the denoised images with or without normalto focus on the relevant anatomical elements and elimi- ization did not present relevant diferences. PReLU was nates the need to consider the entire volume. 2D images adopted as the activation function, since the presence of are easier to visualise and evaluate and provide a recog- a trainable parameter improves the overall performance nisable and comprehensible representation of the data. of the network according to the metrics used. Another advantage is the reduction in the overall size The network can be described with simple equations of the data, which facilitates its transmission or storage, following its division in layers. especially in resource limited environments. The down part can be seen as the function
Despite these advantages, relying solely on 2D slices DownConv2d () defined as in equation 6, in which the has some disadvantages. It may not accurately represent parameter represents the features computed by the complex structures or dynamic processes and may result previous layer. in the loss of 3D or temporal information. In some situations, combining both modalities and examining selected = PReLU (Conv2d ()) (6) slices in the context of the full volume may be more use- +1 = PReLU (Conv2d ()) ful. For each of the 581 subjects, for a total of 43,575 The up part can be seen as the function images, the 25 central slices in the coronal, axial, and UpConv2d (, − − 1) defined as in equation 7, sagittal planes are selected to avoid completely or mostly in which the parameter represents the features comblack images. A padding technique involving continuous puted by the previous layer and the − − 1 parameter zero-filling is used to uniformly scale all photos to the represents the correspondent features computed by the maximum size, since slices acquired from diferent planes down part to be concatenated. have variable sizes. ingA(f7te5r%p),rveaplriodcaetisosnin(g1,5t%he) adnadtatseestt(i1s0s%pl)i;t into the train- = cat(+P1R=eLUPR(eCLoUnv(TC2odnv(2d)(),))− − 1) (7)
Before loading the images in the dataset they are
normalized in the [
In particular, the levels of added noise is based on the a batch of noisy MRI scans in the K-space and the output Signal-to-Noise Ratio (SNR), with the standard deviation is the estimated batch of noise also in the K-space. computed using equation 5, where |ℱ []| is the average of the magnitude of the Fourier transform of the image. 3.3.1. Loss Function
(5) task, so a well suited loss function is the mean squared error (MSE) defined in equation 8.
of the classic U-Net. The main modifications concern In particular we use the MSE applied to noisy scans several aspects of the network. In both the descending and free noise scans is called pixel loss (PL) (equation 9), and ascending parts, five convolutional layers are present. which is the main part of the total loss function. Although additional layers can be added, the size of the Some works like Zhao et al. [24] or Mustafa et al. preprocessed MRI images is 256×256, so adding more [25] showed that in image regression tasks, the MSE layers would lead to activation maps that are too small. succeeded in performing the task with good results, but it Furthermore, the network is faster and lighter by avoid- blurred the resulting images. A popular way to overcome ing additional layers. this problem is to add other parts in the total loss function,
Convolutions do not include bias, since the introduc- focusing on the features [24], or to design from scratch a tion of bias parameters has shown a drastic worsening specific and more complex loss function [25]. of the network performance, generating images with For the proposed network, we integrates other two anomalous visual artifacts. Similarly, normalization tech- parts in the total loss function, the frequency loss (FL) MSE (, ) =
1 ∑︁ (, − )2 =1
(8) pling and 5 upsampling layers, connected via skip connections. (equation 10) that is the sum of MSE applied to real and imaginary part of the Fourier transform of the noisy scans and the free noise ones, and the edge loss (EL) (equation 11 that is the MSE applied to some edge features of the noisy scans and the free noisy ones; in particular the function used are the Laplacian and/or the Sobel filters.
PL (, ) = MSE (, ) (9)
EL (, ) = MSE ( () , ())
previously introduced parts: the pixel loss (PL), the frequency loss (FL) and the edge loss (EL).
vironment with an NVIDIA T4 GPU (CUDA 11.8), 12 GB
We used ADAM optimizer, the starting learning rate of 0.001 and batch size of 64. A learning rate scheduler was used to reduce the learning rate by a factor of 0.1 at every instance that the validation loss failed to decrease for two successive epochs. Training was also monitored using an early stopping criterion that terminated if validation loss failed to drop for five epochs.
age reconstruction evaluation metrics,namely the Peak
PSNR(, ) = 10 log10 ︂( max()2 )︂ MSE(, ) (13) where and are the reference and reconstructed images, respectively, and MSE is the mean squared error.
SSIM(, ) =
(2 + 1)(2 + 2) ( 2 + 2 + 1)( 2 + 2 + 2) , (14) where , are the local means, 2 , 2 the variances, and the covariance between and , with 1 and 2 being constants to stabilize the division. Two pairs of training experiments can be distinguished. The first uses
type of normalization. The studies vary in how noisy scans are produced within each pair. Specifically, the second uses a progressive technique, while the first uses a random one.
configurations, each evaluated under diferent conditions. The first group of experiments employed the Mean
used the Mean Absolute Error (MAE), along with a difer- reduction. For this reasons, the other two experiments
ent input normalization strategy based on centering the are made with the same loss function, but substituting
data with zero mean and unit standard deviation, instead each instance of MSE with MAE.
of scaling it to the [
In the first experiment, the training set was corrupted 1. This because this type of normalization is widely used by Gaussian noise with a signal to noise ratio (SNR) ran- in lots of works and tends to have better performances domly sampled from a discrete uniform distribution be- in most of the cases. tween -5 and 5. For the validation set, the SNR was fixed at -5 to maintain a challenging evaluation scenario. The 4.3. Results model was trained with early stopping based on validation loss. From Table 1 and Figures 2a–3d, the following trends are
The second experiment introduced noise in a progres- easily observable. The most obvious is that progressive sive manner, the training began with images at SNR = 5 noise injection consistently outperforms random noise and gradually moved down to SNR = -5. At each level, the injection. With both loss functions, MSE and MAE, the model was trained separately using a validation set with progressive setting improves the average by about 4.8 dB matching noise level. This setup was designed to let the for PSNR and from 0.14 to 0.15 for SSIM. In addition to model adapt incrementally to increasing noise intensity. achieving higher scores, this training approach is also
The third and fourth experiments mirrored the struc- more stable, with lower variance at various noise levels. ture of the first and second, respectively, but used the Comparing the two loss functions, the model trained with MAE loss instead of MSE, and applied the aforementioned MSE has slightly better SSIM values. The diference is centered normalization. These variations aimed to evalu- small, about 0.01 on average, but constant. This indicates ate the impact of a loss function less sensitive to outliers that the model trained with MSE is better at keeping and a normalization scheme more common in deep learn- structural details intact. In contrast, the model trained ing workflows. Detailed results for all experiments are with MAE is superior in PSNR by about 0.2 dB on average. reported in Table 1. This holds true for all noise levels and suggests that MAE may be better at removing noise, albeit at the expense of 4.2. The MAE and the centered standard slightly worse structural fidelity These results highlight the behavior of the proposed architecture under diferent normalization training conditions. To further contextualize its performance, we compare it against existing state-of-the-art denoising approaches.
in the training procedure allowing the model to be more robust to diferent levels of noise earlier, the progressive approach has faster convergence time because the model only has to adapt to the next level of noise from the previous without the abrupt shifts.
Moreover, MSE can be sensitive to the presence of noise or outliers in the data because it strongly penalizes large deviations from the true values, due to the squared function. So the resulting denoised images could be overly influenced by noise. Instead, the Mean Absolute Error (MAE) (equation 15) is more robust to noise and outliers.
MAE (, ) =
1 ∑︁ |, − | =1 (15)
It treats errors uniformly, due to the absolute value function making it less susceptible to extreme values in the data. Consequently, the denoised images produced tend to be more noise resistant. Also, MAE tends to produce images that are visually crisper, but with better preservation of fine structures. MSE instead tends to suppress high-frequency details and edges in favor of noise
they are compared with other four denoising methods. In particular, with Optimized NLM by Coupé et al. [26], WSM by Coupé et al. [27], MCDnCNNg and MCDnCNNs by Jiang et al. [28].
It is useful to point out that the four chosen methods work with diferent datasets and have a diferent noise generation procedure. Therefore, the collected metrics to make comparisons are diferent (a kind of average on the level of noise) from the real ones.
Analyzing the table 2, the two fast U-Net trained with the progressive approach have overall very similar metrics with the other methods. In particular they are better on MRI scans with low noise, while they lose some performances on scans with high noise. This is due to the trade-ofs made to keep the network fast.
improve its speed while maintaining good performance.
The noise generation procedure, based on the K-space In this work, a new approach to MRI denoising is pre- and the signal-to-noise ratio (SNR), allows to obtain more sented, introducing several novelties in both the noise realistic noisy scans. This realism is also due to the depengeneration phase and the training strategy, as well as dence of the noise on the mean of the Fourier transform an architectural modification of the classical U-Net to modulus. However, this type of generation makes the (a) MSE, random (b) MSE, SNR=5 (c) MSE, SNR=0 (d) MSE, SNR=-5 comparison with other methods more complex, since in most works a simulated noise based on predefined percentages is preferred.
The progressive training approach, where the SNR varies from 5 to − 5, represents a turning point for the improvement of the results both in terms of PSNR and SSIM, as well as leading to a clear visual quality in the denoised images. Furthermore, the adaptive nature of this strategy allows the network to maintain the denoising performance already acquired along all the considered noise levels.
The proposed U-Net, composed of only five layers and without normalizations, can run in short times even on CPU. This compromise in terms of speed does not lead to a significant degradation of the metrics, which remain comparable with those obtained by other methods based on deep learning. It should also be noted that, although processing time is rarely analyzed in the literature, today in the medical field it is increasingly important to obtain reliable results in short times.
ChatGPT, Grammarly in order to: Grammar and spelling check, Paraphrase and reword. After using this tool/service, the authors reviewed and edited the content as needed and take full responsibility for the publication’s content.