Deep Learning with Domain Adaptation for Accelerated Projection-Reconstruction MR

TL;DR

This paper introduces a deep learning method with domain adaptation to reconstruct high-resolution MR images from under-sampled data, significantly reducing acquisition time and outperforming traditional algorithms.

Contribution

The study presents a novel domain adaptation approach that pre-trains on CT or synthesized data and fine-tunes on limited MR data, improving image quality and speed.

Findings

Outperforms compressed sensing algorithms in image quality

Achieves several orders of magnitude faster computation

Pre-training on similar organ data is more effective

Abstract

Purpose: The radial k-space trajectory is a well-established sampling trajectory used in conjunction with magnetic resonance imaging. However, the radial k-space trajectory requires a large number of radial lines for high-resolution reconstruction. Increasing the number of radial lines causes longer acquisition time, making it more difficult for routine clinical use. On the other hand, if we reduce the number of radial lines, streaking artifact patterns are unavoidable. To solve this problem, we propose a novel deep learning approach with domain adaptation to restore high-resolution MR images from under-sampled k-space data. Methods: The proposed deep network removes the streaking artifacts from the artifact corrupted images. To address the situation given the limited available data, we propose a domain adaptation scheme that employs a pre-trained network using a large number of x-ray computed tomography (CT) or synthesized radial MR datasets, which is then fine-tuned with only a few radial MR datasets. Results: The proposed method outperforms existing compressed sensing algorithms, such as the total variation and PR-FOCUSS methods. In addition, the calculation time is several orders of magnitude faster than the total variation and PR-FOCUSS methods.Moreover, we found that pre-training using CT or MR data from similar organ data is more important than pre-training using data from the same modality for different organ. Conclusion: We demonstrate the possibility of a domain-adaptation when only a limited amount of MR data is available. The proposed method surpasses the existing compressed sensing algorithms in terms of the image quality and computation time.

Figures

Equations

$${\hat{P}}_{\theta }(\omega )$$

$$m(x,y)=\frac{1}{2\pi }{\int }_0^{\pi }d\theta {\int }_{-\infty }^{\infty }|\omega |{\hat{P}}_{\theta }(\omega )e^{j\omega }d\omega .$$

$${\hat{P}}_{\theta }(\omega )$$

$${\hat{P}}_{\theta }(\omega )⟶\hat{P}\in {\mathbb{C}}^{N_{\theta }\times N_{\omega }}.$$

$$m(x,y)⟶M\in {\mathbb{C}}^{N_x\times N_y}.$$

$$L_D(y,f)=E_{x\sim D}{\Vert y-f(x;\eta )\Vert }^2,$$

$$L_P(y,f)$$

$${\mathrm{disc}}_{\mathcal{Y}}(P,Q)=\underset{f\in \mathcal{F}}{\sup }\left|L_P(y,f)-L_Q(y,f)\right|.$$

$$L_Q(y,f)\le \underset{\text{empirical risk}}{\underbrace{L_{\hat{Q}}(y,f)}}+\underset{\text{complexity penalty}}{\underbrace{2{\hat{R}}_m(\mathcal{F})}}+3\sqrt{\frac{\ln (2/\delta )}{m}},$$

$$L_{\hat{Q}}(f)=\frac{1}{m}\sum _{i=1}^m\Vert y_i-f(x_i;\eta ){\Vert }^2$$

$${\hat{R}}_m(\mathcal{F})=E_{\sigma }\left[\underset{f\in \mathcal{F}}{\sup }\left(\frac{1}{m}\sum _{i=1}^m{\sigma }_{i}f(x_i)\right)\right],$$

$$L_P(y,f)$$

$$y_i^{\prime \prime }={\hat{y}}_i^{\prime }+ϵ\Delta y_i^{\prime }$$

$${\eta }^{t+1}\leftarrow \arg \underset{\eta }{\min }\sum _{i=1}^b\Vert f(x_i^{\prime };\eta )-y_i^{\prime \prime }{\Vert }^2,$$

$$\Vert f(x_i^{\prime };\eta )-y_i^{\prime \prime }{\Vert }^2\simeq \Vert {\partial }_{\eta }f(x_i^{\prime };{\eta }^t)(\eta -{\eta }^t)-ϵ\Delta y_i^{\prime }{\Vert }^2$$

$${\eta }^{t+1}={\eta }^t+ϵ\Delta {\eta }^t,$$

$${\eta }^{t+1}$$

$$\underset{x}{\min }\frac{1}{2}||y-Ax|{|}_2^2+\lambda TV(x),$$

$$NMSE=\frac{||x-\hat{x}|{|}_2^2}{||x|{|}_2^2},$$

Full text

Methods: The proposed deep network removes the streaking artifacts from the artifact corrupted images. To address the situation given the limited available data, we propose a domain adaptation scheme that employs a pre-trained network using a large number of x-ray computed tomography (CT) or synthesized radial MR datasets, which is then fine-tuned with only a few radial MR datasets.

Results: The proposed method outperforms existing compressed sensing algorithms, such as the total variation and PR-FOCUSS methods. In addition, the calculation time is several orders of magnitude faster than the total variation and PR-FOCUSS methods. Moreover, we found that pre-training using CT or MR data from similar organ data is more important than pre-training using data from the same modality for different organ.

Conclusion: We demonstrate the possibility of a domain-adaptation when only a limited amount of MR data is available. The proposed method surpasses the existing compressed sensing algorithms in terms of the image quality and computation time.

Keywords: Deep learning, convolutional neural network, domain adaptation, projection reconstruction MRI, compressed sensing

Introduction

In MRI, data acquisition is performed by scanning k-space data of objects according to the sampling trajectory. Thus, the quality and the properties of reconstructed images depend on the k-space sampling patterns, such as Cartesian, radial, and spiral sampling patterns [1].

The radial k-space trajectory has many advantages over the Cartesian k-space trajectory, although the Cartesian k-space trajectory is the most widely used sampling pattern. Specifically, each line measured using the radial k-space trajectory contains both center and periphery frequency information. Since the center frequencies are over-sampled for all data lines, the center position is averaged when image reconstruction is conducted such that motion artifacts from the flow or respiration are reduced with the radial k-space trajectory [2, 3, 4, 5]. In addition, the radial k-space trajectory for sub-sampling experiments shows visually better image quality than the under-sampled Cartesian k-space trajectory, as the associated aliasing artifacts appear as less-coherent streaking patterns that are much less disturbing than wrap-around artifact patterns due to the under-sampling in the Cartesian k-space.

Because the increased scan time is a critical weak point, many investigators have employed compressed sensing [6, 7, 8, 9] in an effort to reduce the scan time. In fact, based on the observation that streaking artifacts are less coherent, the radial trajectory was the first sampling pattern used for the ground-breaking compressed sensing (CS) work by Romberg and Candes [7], and many researchers have developed under-sampled radial k-space image reconstruction methods using compressed sensing (CS) algorithms [10, 11, 12, 13]. One of the remarkable aspects of the CS approaches is that although the corrupted image appears to lose fine detail due to the down-sampling, the seemingly lost information can be recovered by iteratively imposing sparsity in the total variation, wavelets, and other transform domains. However, one of the main technical limitations of CS methods is the increased computational complexity due to the iterative reconstruction process, thus necessitating a new method that obviates these limitations.

Recently, deep learning approaches have achieved tremendous success in various fields, such as classification [14], segmentation [15], denoising [16], and super resolution [17, 18]. In the field of medical imaging, most works have focused on image-based diagnostics. Recently, Wang et al.[19] applied deep learning to CS-MRI. They trained a deep neural network from downsampled reconstruction images to learn an instance of fully sampled reconstruction. Subsequently, they used the deep learning result either for initialization or as a regularization term in classical CS approaches. A deep network architecture using an unfolded iterative CS algorithm has also been proposed [20]. Instead of using handcrafted regularizers, the authors attempted to learn a set of optimal regularizers. A multilayer perceptron was also proposed for accelerated parallel MRI [21]. In the x-ray computed tomography (CT) area, Kang et al.[22] provided the first systematic study of a deep CNN for low-dose CT and showed that a deep CNN using directional wavelets can more efficiently remove low-dose-related CT noises. Unlike these low-dose artifacts from reduced tube currents, the streaking artifacts originating from sparse projection views show globalized artifacts that are difficult to remove using conventional denoising CNNs [23, 24, 25]. Jin et al.[26] and Han et al.[27] independently proposed multi-scale residual learning networks using U-Net [15] to remove the global streaking artifacts caused by sparse projection views.

Recall that the projection slice theorem can convert the radial k-space data into a sinogram data format corresponding to parallel beam x-ray CT [28]. Accordingly, the radial k-space sampling data can be reconstructed by a CT reconstruction method, such as filtered back-projection (FBP). This does not require re-gridding of the k-space data to the Cartesian grid, which makes the reconstruction more accurate. This method is often referred to as projection reconstruction (PR) in the MR literature. The degree of similarity between sparse-view CT and accelerated radial acquisition in MR allows us to exploit the synergy from deep learning-based CT reconstruction [27].

In particular, given this similarity between projection-reconstruction MR and CT, this paper proposes a novel means of domain adaptation [29, 30] of a type never been exploited in MR reconstruction approaches which rely on deep learning [19, 20]. Note that a large dataset is usually required when training a deep neural network. In contrast to x-ray CT, the radial trajectory is not widely used in MRI; hence, a large number of PR datasets cannot be easily obtained in most MR sites. Thus, we want to transfer an advanced deep network from a CT dataset and adapt the network parameter to suit MR reconstruction. Although our deep network could be trained with only radial MR data, one of the most important innovations and focuses of this work is to consider this common scenario using domain adaptation, which is an active field of research in deep learning [29, 30]. Given a pre-trained network from CT, the actual MR training only requires very few radial MR datasets, which significantly reduces the training time and expands the applicability of deep learning for MR imaging. By extending this idea furthermore, we also demonstrated that pre-training using synthetic radial MR data from public domain MR images is also an effective approach for domain adaptation when the underlying structure are similar. This finding suggests that the domain adaptation can be widely used to supplement the limited training data set in machine learning approach for MR reconstruction.

Note that our approach can be more accurately classified as a restoration method for artifact removal rather than as a reconstruction method because the trained network is used to restore the image from an image contaminated with streaking artifacts. However, the final goal is to restore high-resolution images from insufficient data, as it is in deep learning-based reconstruction methods [19, 20].

Theory

MR Radial k-space Trajectory versus CT Sinogram

If $m(x,y)$ denotes an MR image, then the radial k-space data obtained at the angle $\theta$ can be expressed as

[TABLE]

Because this is a Fourier transform in the polar coordinates, the inverse transform is given by :

[TABLE]

which is often referred to as filtered back-projection (FBP). In this equation, $|\omega|$ denotes a ramp filter. In fact, the FBP formula in (2) is the standard reconstruction formula for parallel-beam X-ray CT, where $\hat{P}_{\theta}(\omega)$ is obtained by taking the 1-D Fourier transform along the $t$ direction of the projection data $P_{\theta}(t)$. More specifically, we have the projection slice theorem [28]:

[TABLE]

where $P_{\theta}(t):=\int\int dxdy{~{}m(x,y)\delta(t-x\cos\theta-y\sin\theta)}$ is the projection data at angle $\theta$. This similarity between radial acquisition and CT implies that a deep network trained from CT data can be used synergistically in projection-reconstruction MRI. Since the formulations in (1) and (2) are continuous-domain formulation, in practice the k-space data $\hat{P}_{\theta}(\omega)$ are sampled along angle direction and read-out direction, forming a discrete matrix:

[TABLE]

Similarly, the image $m(x,y)$ is also reconstructed on a discrete grid:

[TABLE]

With any insufficiency with regard to the number of radial k-space scan lines $N_{\theta}$, streaking artifacts occur in the reconstructed image from MR and CT. For example, Figs. 1(a) and (b) show the FBP reconstruction results of CT and MR from 48 view projection and 45 radial k-space lines, respectively. Here, the streaking artifacts were obtained from the magnitude images. In Fig. 1(a), the first and the second columns show very different CT images, but similar streaking artifact patterns are observed from those images. Interestingly, similar phenomena were also observed in the MR images, as shown in Fig. 1(b). Thus, we are interested in transferring a network trained from CT domain for MR image restoration. This can be done using domain adaptation which is the main topic in the following section.

Domain Adaptation

A deep neural network typically requires a large number of datasets. Unfortunately, a Cartesian k-space sampling pattern is much more widely used in MRI than a radial sampling pattern, and for this reason it is often difficult to collect many datasets for radial sampling patterns. On the basis of the observation that streaking artifact patterns are similar regardless of the imaging modalities, we utilize a domain-adaptation technique to overcome the insufficiency of radial MR dataset.

Specifically, $x\in X\subset{\mathbb{R}}^{N_{x}\times N_{y}}$ denotes the real-valued magnitude image of the sparse-view reconstruction, and $y\in Y\subset{\mathbb{R}}^{N_{x}\times N_{y}}$ corresponds to the streaking artifact-free images in the magnitude domain. We define a domain as a pair consisting of a distribution $D$ over $X$ and a regression function $f:X\mapsto Y$ that belongs to the function class ${\mathcal{F}}$ as represented by a neural network. We consider two domains: a source domain and a target domain. In this paper, the source domain is CT data or synthesized radial MR data, while the target domain is in vivo radial MR data.

We denote by $Q$ the source domain and by $P$ the target domain distributions. In the domain-adaption scenario, we receive two samples: a labeled sample of $m$ points from the source domain ${{\mathcal{S}}}=\{(x_{1},y_{1}),\cdots,(x_{m},y_{m})\}$ with $x_{1},\cdots,x_{m}\in{\mathbb{R}}^{N_{x}\times N_{y}}$ drawn according to $Q$, matched with the target sample $y_{1},\cdots,y_{m}\in{\mathbb{R}}^{N_{x}\times N_{y}}$, and an unlabeled sample drawn according to the target distribution $P$. In addition, a small amount of label data from the target domain ${{\mathcal{T}}}^{\prime}=\{(x_{1}^{\prime},y_{1}^{\prime}),\cdots,(x_{s}^{\prime},y_{s}^{\prime})\}$ are available. This is based on the practical acquisition scenario in which there are ample projection datasets from CT or synthesized radial MR, whereas relatively few in vivo radial MR datasets are available.

For any distribution $D$, a loss function is defined as a function that minimizes the risk in the target domain, i.e.,

[TABLE]

where $\eta\in{\mathbb{R}}^{N_{\eta}}$ denotes the network parameters, and $x,y\in{\mathbb{R}}^{N_{x}\times N_{y}}.$ Hence, our goal in the domain adaptation is to minimize the loss in the target (i.e., MR) domain. By means of triangular inequality, we have the equations below

[TABLE]

where what is termed the $\mathcal{Y}$-discrepancy between two domains is defined by [31, 32]

[TABLE]

To minimize the risk in the target domain, we are interested in minimizing the risk in the source domain $L_{Q}(y,f)$ and the $\mathcal{Y}$-discrepancy between the two domains. Note that all the calculation of $L_{Q}(y,f)$ and $\mathcal{Y}$-discrepancy is in the magnitude domain, so we do not need to consider the complex nature of the MR image.

However, there are several technical issues. The first of these is that the associated probability distribution $Q$ is unknown. To address this issue, approaches from statistical learning theory [33] offer the bound of the risk of a learning algorithm with regard to a complexity measure (e.g., VC dimension, shatter coefficients) and empirical risk. Specifically, with a probability of $\geq 1-\delta$ for some small $\delta>0$, for every function $f\in{\mathcal{F}}$, we have [34]

[TABLE]

where the empirical risk is given by

[TABLE]

and the Rademacher complexity $\hat{R}_{n}({\mathcal{F}})$ is defined by

[TABLE]

where $\sigma_{1},\cdots,\sigma_{n}$ are independent random variables uniformly chosen from $\{-1,1\}$. Thus, we can obtain the following generalization error bound for the domain adaptation:

[TABLE]

In general, the Rademacher complexity term is determined by the structure of the network. In order to reduce the generalization error of domain adaptation for MR reconstruction, we should therefore minimize the empirical risk in the source domain and the discrepancy between the two domains. Empirical risk minimization in the source domain is a standard neural network training process. In our case, we use CT data or synthesized radial MR data so that the network can learn the mapping from artifact-corrupted image to artifact-free image. To minimize discrepancies, we use a small set of label data from radial MR data ${{\mathcal{T}}}^{\prime}=\{(x_{1}^{\prime},y_{1}^{\prime}),\cdots,(x_{s}^{\prime},y_{s}^{\prime})\}$ with $x_{i}^{\prime},y_{i}^{\prime}\in{\mathbb{R}}^{N_{x}\times N_{y}}$ so that the refined network output follows the target distribution. In particular, we use an incremental approach, in which the network parameter is iteratively adjusted to change the network outputs $f(x_{s}^{\prime})$ to reduce the discrepancy.

Specifically, we denote the estimated network parameter at the $t$-th iteration of our algorithm $\eta^{t}\in{\mathbb{R}}^{N_{\eta}}$: each iteration of our method takes a batch of random input $\{x_{i}^{\prime}\}$ and calculates its corresponding output $\hat{y}_{i}^{\prime}=f(x_{i}^{\prime};\eta^{t})$ based on $\eta^{t}$. To realize a better approximation of the target distribution $P$, we restrict the path of evolution such that

[TABLE]

Therefore, we should update the network parameter such that

[TABLE]

where $b$ denotes the batch size. Using the Taylor series expansion, we have

[TABLE]

Accordingly, we have

[TABLE]

The one-step stochastic gradient descent of (6) is then given by

[TABLE]

where $y_{i}^{\prime}\in{\mathbb{R}}^{N_{x}\times N_{y}}$ denotes the MR artifact-free image in the magnitude domain and $f(x_{i}^{\prime};\eta^{t})$ is the $t$-th iteration network output for a given MR input data $x_{i}^{\prime}\in{\mathbb{R}}^{N_{x}\times N_{y}}$. This is then back-propagated across all network levels to update the network parameters.

Note that there are different approaches for domain adaptation, such as fine-tuning [35] and adversary domain adaptation [36]. The difference in these approaches is essentially the definition of $\Delta y_{i}^{\prime}$ in (5) and (7). The aforementioned approach using (5) corresponds to network fine-tuning using the target domain data, but the algorithm can be improved using more advanced domain-adaptation techniques.

Network Architecture

As shown in Fig. 1, the streaking artifacts are displayed as globally distributed artifacts. In order to capture globally spread streaking artifacts, the receptive field of the network should be sufficiently large and the network depth should be increased. Thus, multi-scale learning network using U-net [15] with residual path is useful due to its large receptive field.

Specifically, as illustrated in Fig 2, in the first half of the proposed network, each stage is followed by a max pooling layer, whereas in the latter half of the network, an average unpooling layer is used. In our previous work [27], we showed that this results in an enlarged effective receptive field, which is more advantageous when removing streaking artifacts. In addition, scale-by-scale contracting paths are used to concatenate the results from the front part of the network to the latter part of the network. The numbers of channels for each of the convolution layers are shown in Fig 2. Note that this number is doubled after each pooling layer. In our network, zero-padding was used for the convolution so that the image size does not decrease while passing through the convolution layer.

In addition, the proposed residual network consists of a convolution layer, batch normalization [37], a rectified linear unit (ReLU) [14], and the contracting path connection with concatenation [15]. Specifically, each stage contains four sequential layers composed of convolution with $3\times 3$ kernels, batch normalization, and ReLU layers. The last stage has two sequential layers, and the last layer contains only a convolution layer with a $1\times 1$ kernel. Finally, the proposed U-net has residual path inspired by the Residual Net [38]. This additional residual path facilitates the training by providing additional path for the back-propagation [39].

Data Type

In neural network training of MR images, care should be taken owing to the mismatch in the data type. Note that the input and output data in our pre-trained network are the real data type, since projection data from CT is real-valued. The situation is the same when the projection data is synthesized from real-valued MR public data set. On the other hand, the data obtained from a real MR system is a complex data type. Therefore, we need to convert the MR data to the real-valued image.

Specifically, for the projection-reconstruction MRI, there are two types of back projection methods: magnitude backprojection and complex backprojection [40]. In both methods, to perform back-projection, the complex radial k-space data is first zero-padded along the radial direction, after which 1-D FFT is conducted to generate complex-valued sinogram data. For magnitude back-projection, we then compute the magnitude of the sinogram data and perform FBP. In complex back-projection, the real and imaginary channels are filtered and back-projected separately. For parallel imaging using multiple coils, the magnitude or complex back-projection is performed separately for each coil data. Then, we calculate the square root of the square-sum-of-square (SSOS) image to combine the complex-valued and/or multi-coil images. Thus, the final reconstructed image from both magnitude and complex backprojection becomes a real-valued image, and our network fine-tuning is performed to learn the real-valued artifact-free magnitude images from a real-valued artifact-corrupted image in the magnitude domain. Thus, we do not need to be concerned about the data-type mismatch between the pre-training and fine-tuning.

Materials and Methods

Data Acquisition

Recall that the proposed method consists of two parts: empirical risk minimization using pre-training in the source domain and the discrepancy minimization by means of fine-turning in the target domain. Thus, we need two data sets: one for the source domain and the other for target domain.

For pre-training with CT data, rebinned fan-beam CT data in our previous study [27] were used to generate the artifact-corrupted and artifact-free images. The CT data was from 3602 slices of actual patient CT data $512\times 512$ in size provided by the 2016 AAPM low-dose CT Grand Challenge (http://www.aapm.org/GrandChallenge/LowDoseCT/). Specifically, the data were acquired from Somatom Definition AS+ and Somatom Definition Flash (manufactured by Siemens Healthcare, Germany) operated in single-source mode. The X-ray tube potential and mAs value were 120 kVp and 200 mAs, respectively. Here, the flying focal spot (FFS) technique [41, 42] served as the means of data acquisition. Because a rebinning process is required to reconstruct images with FFS data, the FFS data was rebinned into fan-beam CT data, and fan-beam FBP is used to obtain the artifact-free images and artifact-corrupted images from full-views and sparse-view data.

In addition, another pre-training were performed using a synthetic radial data from public MR data set. To generate the synthetic radial data, T2 weighted images from Human Connectome Project (HCP) MR dataset (https://db.humanconnectome.org) were used. The T2 weighted brain images contained within the HCP were acquired by Siemens 3T MR system using a 3D spine-echo sequence. The repetition time (TR) and echo time (TE) were 3200 ms and 565 ms, respectively. The number of coils was 32. The field of view (FOV) defined $224\times 224~{}\rm{mm}^{2}$ and the voxel size was 0.7 mm. From the HCP data, the parallel beam projection data are obtained using radon() function in MATLAB.

Once a pre-trained network is given, fine-tuning can be performed using only a few in vivo radial MR datasets to minimize the discrepancy between the source and target domains. In our work, fine-tuning was performed for brain and abdomen MR datasets separately using both pre-trained networks from CT and HCP data set. Here, we aim to investigate the effect of imaging modality (CT vs. MR) and organ (brain vs. abdominal area) in the domain adaptation.

Specifically, the real radial MR data of brain were acquired from a radial spin-echo sequence by a 3.0T MR system manufactured by ISOL technology in Korea. The TR and TE were 2000 ms and 90 ms, respectively. There were 20 slices in total, and the thickness of each slice was 4 mm. The FOV was $220\times 220~{}\rm{mm}^{2}$, and one coil was used. In total, 256 points k-space data instances were obtained for each projection. Fully sampled projection MR data with 180 projection views were acquired from subject 0 for training and validation. Fig. 3 shows the brain MR dataset collected from subject 0 for both training and validation. As shown in Fig. 3(a), 1, 3, 6, 9, and 15 slices from 15 slice images were used in the training dataset. Validation was performed using five distinct slices, as shown in Fig. 3(b). In order to observe the effect of the number of MR datasets, the network parameter was fine-tuned with various numbers of training datasets composed of 1, 3, 6, 9, and 15 MR slices. If only one MR slice was used, the yellow box image in Fig. 3(a) was used. If 3, 6 or 9 slices were used for domain adaptation, then the datasets in rows $(3)$, $(2,4)$, and $(1,3,5)$ in Fig. 3(a) were used. When using 15 slices in the training datasets, all of the datasets in Fig. 3(a) were used. In addition, fully sampled and in vivo accelerated radial MR data consisting of 180, 90, and 45 projection views were obtained from subjects 1 and 2. This dataset was used during the test phase to evaluate the generalization performance of the trained network.

In addition, liver MRI scans of seven healthy volunteers were acquired on the axial plane using a 3D stack-of-radial sequence with golden angle ordering. The experiments were performed on a 3T MRI system (PrismaFit, Siemens Healthcare GmbH, Erlangen, Germany) after the subjects provided informed consent according to a protocol approved by the local Institutional Review Board (IRB) of UCLA. Each subject was instructed to take a 19 second breath-hold, and the acquisition was performed during the breath-hold. The acquisition parameters were as follows: flip angle = 10*∘*, TR = 2.35 ms, TE = 1.12 ms, and FOV = $340\times 340~{}\rm{mm}^{2}$. There were 20 coils. Due to the golden angle ordering, we retrospectively under-sampled k-space data at various acceleration factors by simply taking a portion of the fully sampled radial spokes (302 projections views and 192 points per projection), and we removed the first and last slices out of 20 slices due to artifacts in the fully sampled data. Similar to the brain MRI data, subject 0 was used for training and validation with 15 and 3 slices as the training and validation datasets, respectively. The data for the remaining patients were used during the test phase.

The magnitude and complex back-projection were used to obtain both under-sampled and fully-sampled reconstruction images in the brain and abdominal dataset, respectively. Subsequently, we conducted experiments with several downsampled factors to learn the streaking artifact from various downsampling factors, where the artifact-contaminated images from all downsampling factors were used for training. To perform the back-projections, the complex radial k-space data was first zero-padded with 729 samples along the radial direction, after which 1-D FFT was performed to generate complex-valued sinogram data. The resulting size of the reconstructed real-valued images was $512\times 512$ for both the brain and abdomen images. As described before, the images are finally obtained as the magnitude image from complex-valued reconstruction, and for parallel imaging, the images are obtained as the square root of the sum of square (SSOS) image by combining the complex-valued and/or multi-coil images.

Here, it is important to note that one of the most important advantages of domain adaptation is that we do not need additional pre-training for our MR reconstruction problem if there exists an available pre-trained network [27]. Hence, the computational cost of pre-training for the domain-adaptation scheme can be often considered as zero.

Network Training

The pre-trained network from CT data was trained by means of the stochastic gradient descent (SGD) method. The regularization parameter was $\lambda=10^{−4}$. The learning rate ranged from $10^{-3}$ to $10^{-5}$, and was gradually reduced at each epoch.

Another pre-training was also performed with the synthetic radial MR dataset. T2 weighted magnitude images of HCP dataset were used as artifact-free images, from which synthetic projection MR data were generated using the $radon$ operator in MATLAB. The input images with global artifact caused by insufficient radial acquisition were re-generated with $iradon$ operator using 36, 45, 60, and 90 views projection MR data. As a training dataset, we used 100 subject data, including 3600 slices images. The network parameters were the same as those for the CT pre-trained network.

For given pre-trained networks, fine-tuning was then performed using a small number of MR slices to adjust for the difference between the two domains. The fine-tuning procedure was identical to that in pre-training with a stochastic gradient algorithm using the same parameters except the learning rate. Specifically, when CT pre-trained network was used, the learning rate is applied to $10^{-4}$ to $10^{-6}$. On the other hand, the learning rate is $10^{-5}$ to $10^{-7}$ if synthetic radial MR pre-trained network was used. Other slices with very different structures were then used as the test set for validation. We used 1000 and 3000 epochs for the brain data and the abdominal data, respectively, including the pre-training. Mini-batch data was used using an image patch, and the size of image patch was $256\times 256$. Since CNN learns the convolutional filters which are spatially invariant, the same filter can be used for reconstructing $512\times 512$ images. Using the $256\times 256$ patch size significantly reduces network training time.

The proposed network is designed for reconstruction from 2, 3, 4, and 5 accelerations. To learn various acceleration factors simultaneously, the network was trained using several down-sampling images together. Specifically, for the brain dataset, reconstructions from 180 views were considered as artifact-free data; consequently, the input images used to train the brain dataset consist of batches of reconstructed images using 90, 60, 45, and 36 views. For the abdominal data, 302 views were used to reconstruct the artifact-free image. Thus, the training dataset is composed of downsampled images with 151, 100, 75, and 60 views. The simultaneous training using multiple acceleration factors was helpful in that it reduced over-fitting during the network training process.

The network was implemented using the MatConvNet toolbox (ver.20) [43] in the MATLAB 2015a environment (MathWorks, Natick). We used an NVidia GTX 1080 Ti graphic processor and an Intel i7-7770 CPU (3.60GHz).

Performance Evaluation

In order to demonstrate the comparative advantage of our method, we have compared with CS algorithms including total variation (TV) [44]-based CS [7] and PR-FOCUSS [10]. The TV result was calculated by solving the following TV minimization equation,

[TABLE]

where $A$ is the subsampled projection operator, $y$ is the sub-sampled measurement, $TV(x)$ denotes the TV penalty, and $\lambda$ is the regularization parameter. The TV minimization problem was solved using the alternating direction method of multiplier (ADMM) algorithm with variable splitting with a regularization parameter of $\lambda=10^{-2}$ and an iteration number of 100. The PR-FOCUSS reconstruction scheme utilized the parameter settings from our earlier work [10].

For a quantitative evaluation, we used the normalized mean square error (NMSE). The NMSE is defined by

[TABLE]

where $x$ denotes the reconstruction results from full views and $\hat{x}$ is the reconstruction results from the sub-sampled projection data. The average NMES values were calculated by averaging all slice reconstruction results.

Results

To evaluate the effective domain adaptation schemes, we have evaluated two pre-trained networks: one from CT data and the other from synthetic radial MR data from HCP data set. Note that the CT data from AAPM low-dose CT challenges are mainly from abdominal imaging, whereas the HCP data are for brain imaging. Thus, the source and target domains are different in their imaging modalities and target organs, so we are interested in investigating important factors for domain adaptation.

Fig. 4(c) shows the brain result from the CT pre-trained network, when an under-sampled MR image as shown in Fig. 4(b) is used as the network input. Because the CT pre-trained network is focused on restoring CT global deficient components, the network recognized the detailed MR structures as CT global artifacts and erroneously removed them. Moreover, the pre-trained network was trained using fan-beam CT geometry instead of parallel-beam geometry, as in projection-reconstruction MR, such that the resulting global artifacts were statistically different. Therefore, the detailed structures of MR were not fully recovered, and the restored images appear as CT images. However, domain adaptation using a fine-tuning with radial MR data significantly improves the reconstruction performance. Specifically, Fig. 4(d) shows a restored image from the fine-tuned network with only one MR dataset. The detailed structure was better conserved in the fine-tuned network than in the pre-trained network, as compared in Fig. 4(c) and Fig. 4(d). On the other hand, Fig. 4(e) illustrates the baseline result from the HCP pre-trained network. The detailed MR structure was conserved, whereas the result in 4(c) from CT pre-trained network did not restore the fine detail. Moreover, with additional fine tuning with one MR slice, Fig. 4(f) shows MR reconstruction images more realistic than Fig. 4(d). The HCP pre-trained network can successfully removes global artifact, but the network preserves the detailed MR structures in contrast with the CT pre-trained network.

For a quantitative evaluation of the network performance with respect to the different pre-training strategy, 90, 60, 45, and 36 projection views were retrospectively sub-sampled from 180 views of projection data of two other subjects (subjects 1 and 2). Assuming that 180 views correspond to fully sampled data, Fig. 5(a) shows the corresponding normalized mean square error (NMSE) values. The average NMSE values were calculated by averaging twenty slice restoration results. The results clearly showed that the two proposed networks exceeded CS methods such as the TV and PR-FOCUSS[10]. Among the two domain adaptation schemes using CT and HCP pre-trained network, the HCP pre-trained network has better performance than the CT pre-trained network.

For a qualitative comparison, Fig. 6(a) shows the reconstruction results from retrospective sub-sampling data with 45 views using TV and PR-FOCUSS as well as the proposed network fine-tuned with only one MR dataset from subject 0. The subjective image quality as well as the NMSE value of the slice in Fig 6(a) clearly demonstrated that the proposed method restores the global artifact by eliminating streaking artifacts and preserving the detailed brain structure for all of the views. Next, in vivo downsampled data with 90 and 45 view projections were acquired from subjects 1 and 2. Fully sampled projection data with 180 views were also acquired separately. The quantitative results in Fig 5(a) and the subjective image quality in Fig 6(b) from 90 views also confirmed that the proposed network can be generalized well to in vivo accelerations with significantly improved image quality. In these in vivo experiments, the NMSE value for the reconstruction was calculated from 180 views acquired in a separate scan.

To evaluate the diagnostic quality of the reconstruction, two radiologists (HHK and HJS) performed a blind evaluation of the various reconstruction methods from 90 view data from subject 1 and subject 2. As shown in Fig. 5(b), in both the in vivo and retrospective downsampling cases, the proposed method outperformed the TV and PR-FOCUSS reconstruction. Moreover, radiologists found that the results from the proposed method with 15-slice fine-tuning did not differ from those of the fully sampled case and that the proposed method with 1-slice fine-tuning still offered very good diagnostic quality. On the other hand, the results from PR-FOCUSS and TV showed many artifacts, which affected the diagnosis.

In addition to the quality improvement, the computation time was significantly reduced using the proposed method. The proposed method required only 0.05 seconds for restoration, while the computational times for the TV and PR-FOCUSS methods were approximately 24 to 38 seconds and nearly 29 to 60 seconds, respectively.

Fig. 7(a) shows the abdomen results for the representative slices from subject 1 when the CT pre-trained network was fine-tuned network with the 15 radial MR dataset. The ground truth with 302 views was obtained by the golden-angle method and the downsampled data were generated by collecting 75 views from the first part of the golden-angle projection data. Although the detailed structure of abdominal images is more complex than that of the brain image, the proposed method exhibited stable restoration performance. Fig. 7(b) shows a quantitative evaluation of the two different domain adaptation scheme when 302 views projection data were retrospectively sub-sampled with 75, and 60 views projection data. In contrast with Fig. 5(c), CT pre-trained network is quantitatively better than HCP pre-trained network. This implies that the domain adaptation from similar organ is more important than using the same imaging modality.

Discussions

Advantages of domain adaptation

In order to verify the importance of the domain adaptation, we conducted various comparative studies. First, using the brain dataset, we compared the performance of the proposed network with that of another baseline network trained only with a few radial MR data. For a fair comparison, the network parameters for the MR-only network were equated with those of the proposed network architecture, except for the learning rate. The number of epochs used was 500, equal to that of the proposed method. The convergence plot is shown in Fig 8. The objective and PSNR values were calculated for the labels and restored images with a patch size of $256\times 256$ when the network is trained using retrospective down-sampled images of 36, 45, 60, and 90 views. The dash-line and solid-line define objective functions for train phase and validation phase, respectively. Although the networks were trained with up to 1,000 epochs, they converged reliably after the 500th epoch without over-fitting or under-fitting. Hence, the networks for the 500th epoch were used to restore the MR images. Compared with the MR-only residual network, the domain adaptation scheme showed superior performance.

The restoration results from 60 view-projection data instances are shown in Fig. 9. Domain adaptation resulted in improved restoration, as confirmed in the difference images. Moreover, for the case of brain data, the domain adaptation using HCP pre-training has less residual signals in the edges, indicating that the reconstruction performance improved compared to the CT pre-trained network. In addition, as revealed from the quantitative results in Fig. 5(c) for the retrospective and in vivo sub-sampling cases, the proposed networks significantly outperform the MR-only residual network, as using only a few MR datasets was not enough to allow the learning of the strong global artifact pattern caused by severely under-sampled projection views.

While the radial trajectory is not widely used, one area where they are increasingly used is for phase contrast flow imaging. In this case, the recovering phase information is useful. In this case, we could train the phase network separately similar to [47] . In this case, the magnitude network we have described so far can guide the training and inference of phase network. However, this is beyond the scope of this paper.

Readers may wonder whether we need additional data for source domain ${{\mathcal{S}}}$, such as CT dataset and synthetic MR dataset to use the proposed method. However, this is not the case, as the pre-trained network is feasible as a starting point. The situation is analogous to existing classification studies in the computer vision community using domain adaptation. For example, AlexNet [14] is a classical pre-trained network for the ImageNet dataset [45]. If one is interested in designing an image classification network for a microscopy dataset, for instance, it would be feasible to start with AlexNet and fine-tune it with a microscopy dataset. The designer does not need to run additional pre-training for the ImageNet dataset; therefore, the only additional complexity and datasets come from the fine-tuning stage. Similarly, in our MR network, as a pre-trained network is used as a starting point, the proposed network is quickly trained and converges to a good solution, even with a small amount of MR data. This is an important advantage of the proposed network compared to the existing deep learning approach for MR, which requires a considerably large dataset. In fact, our pre-trained network will be publicly available on the authors’ homepage upon publication so that readers can use it for fine-tuning.

One may be concerned about whether the fine-tune network can deteriorate the pre-trained network. However, the good news is that users should not be concerned about the effect on the pre-trained network from the fine-tuning process because if the user is interested in the CT reconstruction problem, the original pre-trained network can be used for their purpose. Note that our goal is not to design a network that is optimal for both CT and MR data. Rather, our goal of domain adaptation is to design a network that is optimal for MR data.

Dependency on the fine-tuning data size

In order to investigate the dependency on the fine-tuning dataset, we fine-tuned the network with a various numbers of MR datasets. Figs. 10(a) and (b) show in vivo experimental results for brain and abdominal data set with various size of fine tuning data set. Using more radial MR data improved the quality of the restored images. This subjective image quality coincides with the quantitative results shown in Fig. 5(c), showing reduced average NMSE values with more radial MR datasets for fine-tuning. Note that the abdominal image has more complex structures when compared with the brain image. Owing to the image details, the fine-tuning steps for the abdominal image require more MR datasets and more epochs, and we found that the 15 MR dataset and 3,000 epochs are reasonable. Again, with this small amount of data, it is remarkable to observe the superior image reconstruction quality in Fig. 7 relative to the outcomes of all existing methods.

In addition, as shown in Fig 5(c), an increase in the number of MR datasets improves the performance of the MR-only network, and it is impressive to note that the MR-only residual network trained with more than three MR slices outperforms the CS-based approaches (see Figs. 5(a) and 5(b)). This again demonstrates the advantages of the proposed learning architecture. While the MR-only residual network shows a marked improvement with the number of MR datasets, the proposed domain-adaptation method is capable of more stable reconstruction than the MR-only network.

Acceleration Ratio

The current method for the brain data and abdomen dataset requires at least 36 and 75 projection views, respectively. Some readers may wonder whether this is a significant acceleration factor considering the latest radial CS works [12, 13], which require far fewer radial spokes. However, it is important to note that those works [12, 13] are dynamic CS works that exploit temporal redundancies, whereas our work is for static MR reconstruction, which only exploits the redundancies within the frame rather than adjacent temporal frames. In the aforementioned dynamic contrast-enhancement MRI studies [12, 13], the structures remain the same except for breathing-related motions. Therefore, there are many temporal redundancies that can be exploited to achieve this higher acceleration factor. However, to the best of our knowledge, we are not aware of any static CS MR works that can achieve such a high acceleration factor in actual datasets such as those from the brain and liver.

In fact, an extension of the proposed method for dynamic MR images may be feasible and a significantly higher acceleration factor may be achievable using additional temporal redundancy. This is a very important topic in its own right, and it will be reported elsewhere.

Generalization performance

As discussed earlier, the higher the acceleration factor is, the fewer the radial k-space sampling trajectories are in the Fourier domain, which leads to more severe streaking artifacts. However, one of the most important advantages of the proposed deep learning approach is its generalization power. More specifically, streaking artifacts during the training phase do not need to be identical in the test phase. For example, our CT network is trained using a fan-beam trajectory with different acceleration factors, whereas projection-reconstruction MR is a type of parallel-beam geometry. However, even without applying domain adaptation, most of the streaking artifacts can be removed (see Fig. 4). Such remarkable generalization power of a deep network has been an intensive research interest in the statistical learning theory community, and many studies attribute this to the exponential representation power of a deep network [46]. Therefore, as long as we can avoid over-fitting, the same network can be generalized well for various acceleration factors. Interestingly, to avoid over-fitting, we found that training with various acceleration factors was helpful, as described earlier. More specifically, when the proposed network was trained using datasets from various acceleration factors, it was better trained with little over-fitting and was able to reconstruct corrupted images at various acceleration factors.

Conclusion

We developed a novel deep learning approach with domain-adaptation to reconstruct high quality images from sub-sampled k-space data in MRI. The proposed network employed a pre-trained network using CT datasets or synthetic radial MR data, with fine-tuning using a small number of radial MR dataset. If there is a sufficient amount of radial MR data, the deep network with MR data alone showed good restoration results. However, for an insufficient training dataset, the proposed network showed the best performance when combined with pre-training and fine-tuning steps. Among the various factors such as organ and imaging modality, our experiments showed that the similar organ structure is more important, so it is better to use a pre-trained network from similar organ structures. Moreover, the computation speed was much faster than the speeds of conventional compressed sensing methods because the proposed method does not use projection or back-projection, both of which require heavy computational complexity.

The proposed domain-adaptation approach demonstrated the potential for mixing different medical systems when the image artifacts are similar and topologically simple. The principles of medical imaging systems such as computed tomography (CT), magnetic resonance image (MRI), and optical diffraction tomography (ODT) are closely related in the Fourier space according to the Fourier slice theorem and the Fourier diffraction theorem. Based on this relationship, the domain-adaptation approach can be applied to various modalities to remove various artifacts.

Acknowledgment

The authors would like to thanks Dr. Cynthia MaCollough, the Mayo Clinic, the American Association of Physicists in Medicine (AAPM), and grant EB01705 and EB01785 from the National Institute of Biomedical Imaging and Bioengineering for providing the Low-Dose CT Grand Challenge data set. This work is supported by National Research Foundation of Korea, Grant number NRF-2016R1A2B3008104 and NRF-2013M3A9B2076548. The data for a pre-trained network were provided by the Human Connectome Project, MGH-USC Consortium. Thus, this work acknowledges funding from the National Institutes of Health, NIH Blueprint Initiative for Neuroscience Research grant U01MH093765, National Institutes of Health grant P41EB015896, and NIH NIBIB grant K99/R00EB012107