Brain Image Synthesis with Unsupervised Multivariate Canonical CSC$\ell_{4}$Net

Image synthesis (a.k.a. image-to-image translation) is commonly performed via appearance transformations, such as linear regression or distribution transformation. Previous works with stand-alone image pairs tend to focus on constructing a linear relationship in different contrasts [24]. Limited by few off-the-shelf paired data, Vemulapalli et al. [28] relaxed the invariable supervision by matching similarities across different modalities of image patches, and then jointly maximizing both global mutual information and local spatial consistency. Huang et al. [11] proposed a data-efficient synthesis method by mapping both a few pairs and large amounts of unpaired patches into a high-dimensional space, and then adopting Laplacian eigenmaps for geometric co-regularization. Neural style transfer [4] is another popular strategy for content-fixed image style translation, which computes the distance via the Gram matrix statistics of pre-trained deep features. GAN [6] was introduced to generate images using a random noise vector with discriminator judgment. Subsequently, many improvements and task-oriented generative models have been proposed. CycleGAN [39] uses a simple yet efficient strategy with cycle-consistent adversarial networks for unsupervised image-to-image translation. To synthesize high-resolution images from semantic labels, Wang et al. [30] proposed an adversarial learning objective which leverages GANs in a conditional setting and discards hand-crafted losses or pre-trained networks. FUNIT [16] contains a content and class-leveled encoder and an overall decoder for synthesizing the analogous image in a desirable category from the given class input. Leveraging the U-Net architecture, a data augmentation method [35] was established by learning both spatial deformation fields and intensity transforms to generate samples. TrGAN [29] focus on improving unsupervised image synthesis and representation learning.

The existing studies on medical imaging [12, 37, 13], e.g. synthesis, segmentation and registration, have shown great promise for either research purposes or clinical analysis, mostly toward a macro objective, i.e. computer-aided diagnosis (CAD). In [14], a probabilistic model was proposed for joint registration and synthesis with cross-modality alignment. Uzunova et al. [27] used a multi-scale GAN to generate large amounts of high-quality medical images by learning a growing resolution conditioned on front scales. Shao et al. [25] presented a diagnosis-guided multi-modal feature selection method for prognostic prediction of a specific disease. Ravi et al. [23] employed adversarial learning in their proposed DaniNet, a degenerative adversarial neuroimage network that allows neurodegeneration to be modeled.

One of the fundamental purposes of medical image processing is to implement CAD, which in turn allows clinicians to make accurate decisions or provide treatment. Toward generalization and practicability, our work is complementary to the aforementioned approaches, where we improve the usage of large-scale heterogeneous medical data in an unsupervised manner.

Convolutional sparse coding (CSC) has been successfully used in a wide range of computer vision and medical image processing problems [3, 8, 10]. CSC deals with the suboptimality of conventional local-independent representations by introducing a global shift-invariant filter. Fast CSC model was proposed by Bristow et al. [3], where the quad-decomposition solves the CSC objective in the Fourier domain, resulting in fast training. CSC-SR, introduced by Gu et al. [8], uses CSC with improved consistency to super-resolve natural images. Heide [10] tackled feature learning with fast and flexible CSC. Huang et al. [12] constructed two invertible mappings based on CSC for cross-modality synthesis and super-resolution of brain images. To reconstruct clean images, while avoiding adversarial attacks, Sun et al. [26] presented a stratified CSC algorithm with the benefit of an input transformation-based defense.

Let $\mathbf{X}=\{\mathbf{x}_{1},...,\mathbf{x}_{S}\}$ be a source domain training set containing $S$ source modality images, and $\mathbf{Y}=\{\mathbf{y}_{1},...,\mathbf{y}_{T}\}$ be a target domain training set containing $T$ target modality examples. We define $\mathbf{F}$ as the filter of the convolutional feature maps $\mathbf{Z}$. When standard training following the independence assumption illustrated in Eq. (1) is applied to two modalities, this leads to two separate filters of shift-invariant atoms $\mathbf{F}^{x}$, $\mathbf{F}^{y}$, with their corresponding feature maps $\mathbf{Z}^{x}$ and $\mathbf{Z}^{y}$. However, when the same modality data are sampled from random measurements (i.e. intra-modal variation), the new features of these variables are no longer a CSC solution. To overcome this problem, we attempt to regularize the intra-modal data to guarantee a unique convolutionally sparse solution. Considering that a unit normalization of data is taken from the random measurements of one modality, one possibility is to normalize the column elements of $\mathbf{Z}^{x}$ and $\mathbf{Z}^{y}$ to unity, i.e., $\left\|\mathbf{Z}^{x}_{i}\right\|^{2}_{2}=1$, $\left\|\mathbf{Z}^{y}_{j}\right\|^{2}_{2}=1$, $\forall i=1,...,S$, $\forall j=1,...,T$. However, using the unit normalization leads to ”erased” modality information, which is especially the case in the context of neuroimaging data. To circumvent this issue, we eliminate the intra-modal scaling ambiguity by first computing the maximum $\ell_{2}$ norm of $\mathbf{Z}^{x}$ and $\mathbf{Z}^{y}$, respectively, and then performing our intra-modal unit normalization (IUN). The IUN (termed as $\Upsilon$) of the convolutional feature maps becomes

The restricted elements of $\mathbf{Z}^{x}$ and $\mathbf{Z}^{y}$ need to be unified by IUN and satisfy the general unit normalization $\left\|\hat{\mathbf{Z}}^{x}_{i}\right\|^{2}_{2}=1$, $\forall i$, and $\left\|\hat{\mathbf{Z}}^{y}_{j}\right\|^{2}_{2}=1$, $\forall j$.

The CSC model of Eq. (1) has the advantage of computational efficiency when compared to the network-relevant methods. However, the pure penalty under $\ell_{1}$-norm ignores the diversity and complexity of data, leading to unsatisfactory results. Recent studies [6, 36] have shown that, by stacking multiple layers on top of each other, the extracted features are deeper and thus the performance of applications (e.g. classification or reconstruction) can be further boosted. Inspired by the success of extensive efforts on networks, we construct a CSC network architecture with a novel multivariate canonical adapted feature mapping layer, which has cross-modal learning power with synthetic potential.

Without loss of generality, we consider the CSC over IUN, which has multiple layers (termed as CSC-Net). We denote a function $f(\cdot,\cdot)$ for the feature mapping of each layer, such that ${\mathbf{Z}^{x}}^{\prime}=f(\mathbf{X},\mathbf{\Psi}^{x})$, ${\mathbf{Z}^{y}}^{\prime}=f(\mathbf{Y},\mathbf{\Psi}^{y})$, where $\mathbf{\Psi}^{x}=\{\mathbf{F}^{x},\Upsilon,\lambda\}$, $\mathbf{\Psi}^{y}=\{\mathbf{F}^{y},\Upsilon,\lambda\}$. Let $l$ be the network layers, where $l\in[1,L]$. In CSC-Net, the feature maps $\hat{\mathbf{Z}}$ can be defined as the representation of the layer $l$ with tensor properties of height $h$ and width $w$: $\hat{\mathbf{Z}}^{x,\left|l\right|}_{i}\in\mathbb{R}^{N^{\left|l\right|}\times h^{\left|l\right|}\times w^{\left|l\right|}}$, $\forall i,l$, and $\hat{\mathbf{Z}}^{y,\left|l\right|}_{j}\in\mathbb{R}^{N^{\left|l\right|}\times h^{\left|l\right|}\times w^{\left|l\right|}}$, $\forall j,l$. To hierarchically approximate the convolutional feature maps, we construct a sparse intermediate representation which imposes the same structure on the upper layer of the representation. Intuitively, the $l$-th layer of $\mathbf{Z}$ for $\mathbf{X}$ and $\mathbf{Y}$ can be estimated as $\hat{\mathbf{Z}}^{x,\left|l\right|}_{i}=f(\hat{\mathbf{Z}}^{x,\left|l-1\right|}_{i},\mathbf{\Psi}^{x,\left|l-1\right|})$ and $\hat{\mathbf{Z}}^{y,\left|l\right|}_{j}=f(\hat{\mathbf{Z}}^{y,\left|l-1\right|}_{j},\mathbf{\Psi}^{y,\left|l-1\right|})$, respectively.

In addition to the multi-layer sparse and deeper representation, another core module of the CSC-Net is the multivariate canonical adaptation (MCA). The multivariate formulation is initially normalized by IUN for intra-modal unity, and then we optimize the inter-modal (i.e. cross-modal data) structure using the constructed space, allowing us to utilize the feature maps learned in the previous step. Many domain-related studies [11, 38, 17] are based on the concept of ‘feature shift’, which can be summarized as learning a domain-invariant feature representation between the source and target domains for objective transformation. Motivated by the state-of-the-art domain adaptation works [21, 5], we cross-transfer both intra-modal and inter-modal features to a high-level projective space to handle the multivariate heterogeneous medical data. Specifically, we begin by importing a reproducing kernel Hilbert space (RKHS), where the learned features between the CSC-Net layers and RKHS can be formed to minimize the maximum mean discrepancy (MMD) [7] between the source and target domains optimally, using a kernel $k(x_{i},y_{j})=\left\langle\hat{\mathbf{Z}}^{x}_{i},\hat{\mathbf{Z}}^{y}_{j}\right\rangle_{\mathcal{H}}$ under a Hilbert space $\mathcal{H}$. Note that, $\left\langle\cdot,\cdot\right\rangle_{\mathcal{H}}$ represents the inner product, and $k$ is defined on the vector. Following the virtue of the MMD function in Eq. (3),

we adapt the multiple feature layers using the multilayer MMD penalty (termed as $\mathcal{L}_{\mathcal{H}}(X,Y)$) for adapting the cross-modal CSC-Net, which is defined as

Here, the previously described kernel $k$ is computed by the Gaussian kernel function defined on the vectorization of tensors $\hat{\mathbf{Z}}^{x,\left|l\right|}_{i}\in\mathbb{R}^{N^{\left|l\right|}\times h^{\left|l\right|}\times w^{\left|l\right|}}$ and $\hat{\mathbf{Z}}^{y,\left|l\right|}_{j}\in\mathbb{R}^{N^{\left|l\right|}\times h^{\left|l\right|}\times w^{\left|l\right|}}$ of the layer $l\in L$. The calculation in a RKHS space can be expressed as $e^{{-\left\|\hat{\mathbf{Z}}_{i}^{x,\left|l\right|}-\hat{\mathbf{Z}}_{j}^{y,\left|l\right|}\right\|}^{\frac{2}{p}}}$, where $p$ is a bandwidth parameter reflected in the Gaussian kernel function. $\prod_{l\in L}k^{\left|l\right|}(\hat{\mathbf{Z}}_{i}^{x,\left|l\right|},\hat{\mathbf{Z}}_{j}^{y,\left|l\right|})$ gives the multiplicative results of the inner products.

Generally, CSC consists of minimizing a convolutional model-fitted least-square system and a sparse regularizer by adopting an $\ell_{0}$ or $\ell_{1}$ penalty to promote the expected sparsity for recovering objectives. The $\ell_{0}$-norm is NP-hard when finding a local minimum of a nonconvex function, while the $\ell_{1}$-norm provides a unique solution through a convex program in the polynomial time. Although existing CSC algorithms can be separately optimized by alternating subproblems to update sparse feature maps and filters under the non-smooth $\ell_{1}$ penalty and the $\ell_{2}$ constraint, when applied to neuroimaging (which has a high-dimensional character), this leads to a high computational complexity due to the quasi-polynomial problem.

Pioneering research [32] shows that maximizing the element-wise $\ell_{4}$-norm enables an interpretative, stable and robust algorithm for reducing the per-iteration cost. Moreover, the global geometry of the $\ell_{4}$-norm over a unit $\ell_{2}$-sphere guarantees a randomly initialized first-order gradient descent algorithm [33], since each saddle point has negative curvature. Based on the essence of the $\ell_{4}$-penalized heuristic formulation, i.e. replacing $\min\left\|\cdot\right\|_{1}$ as $\max\left\|\cdot\right\|_{4}^{4}$, we attempt to address the shortcomings of CSC by modeling $\ell_{4}$ maximization directly as the sparsity penalty within the objective of CSC$\ell_{4}$Net, which can be formulated as

Here, $\mathcal{S}(X)=\sum_{i=1}^{S}\left\|\hat{\mathbf{Z}}_{i}^{x,\left|l\right|}\right\|_{4}^{4}$ and $\mathcal{S}(Y)=\sum_{j=1}^{T}\left\|\hat{\mathbf{Z}}_{j}^{y,\left|l\right|}\right\|_{4}^{4}$ represent the sparsity cost function for penalizing $\mathbf{Z}^{x}$ and $\mathbf{Z}^{y}$ to be sparse. Mathematically, $\left\|\cdot\right\|_{4}^{4}$ is the element-wise $\ell_{4}$-norm with an expression $\left\|\mathbf{Z}\right\|_{4}^{4}=\sum_{i,j}\mathbf{Z}_{i,j}^{4},\forall\mathbf{Z}\in\mathbb{R}$. $\mathcal{L}_{\mathcal{R}}(X)$ denotes the reconstruction loss of $\mathbf{X}$ with $\frac{1}{2}\left\|\mathbf{X}-\sum_{i=1}^{S}\mathbf{F}_{i}^{x,\left|l\right|}\ast\hat{\mathbf{Z}}_{i}^{x,\left|l\right|}\right\|_{2}^{2}$ and $\mathcal{L}_{\mathcal{R}}(Y)$ is the reconstruction loss of $\mathbf{Y}$ having $\frac{1}{2}\left\|\mathbf{Y}-\sum_{j=1}^{T}\mathbf{F}_{j}^{y,\left|l\right|}\ast\hat{\mathbf{Z}}_{j}^{y,\left|l\right|}\right\|_{2}^{2}$.

Once the training stage is complete, we can get two sets of filters, $\mathbf{F}^{x}$ and $\mathbf{F}^{y}$, and the associator $\mathbf{P}$. Given a test image $\mathbf{T}^{x}$ with the source modality, the desirable target modality version of $\mathbf{T}^{x}$ can be treated as $\mathbf{T}^{y}=\sum(\mathbf{Z}^{tx}\mathbf{P})\ast\mathbf{F}_{y}$, where $\mathbf{Z}^{tx}=f(\mathbf{T}^{x},\mathbf{\Psi}^{x})$.

The proposed CSC$\ell_{4}$Net is evaluated on three datasets: IXI, ¹¹1http://brain-development.org/ixi-dataset NAMIC Multimodality, ²²2http://insight-journal.org/midas/collection/view/190 and BraTS. ³³3https://www.med.upenn.edu/sbia/brats2018/data.html IXI dataset contains 578 healthy subjects, NAMIC dataset includes 10 normal controls and 10 schizophrenic cases, and BraTS dataset has 220 brain tumor subjects. We conduct extensive experiments on three scenarios: (1) Proton Density (PD) $\rightleftharpoons$ T2 on IXI dataset, (2) T1 $\rightleftharpoons$ T2 on NAMIC dataset, (3) FLAIR $\rightleftharpoons$ T1 on BraTS. In each dataset, the existing well-aligned paired images are removed, with half of them per a side for a strict unsupervised setting. Specifically, 150 unpaired PD-w and T2-w MRI are selected from IXI dataset, 6 unpaired T1-w, T2-w acquisitions are picked from NAMIC dataset, and 70 unpaired T1-w and FLAIR data are chosen from BraTS dataset for training. The split 50 samples from IXI, 2 samples from NAMIC and 30 samples from BraTS are used for validation. The remaining data: 100 (IXI), 4 (NAMIC), and 40 (BraTS) are used for testing. It is worth noting that, in our experimental datasets, IXI only contains healthy images while both NAMIC and BraTS include pathological data. In other words, the first setting in our experiments considers healthy cases, the second setting involves a mix of healthy and pathological examples, and the third setting tests our method on common diseases. We tune the hyper-parameters of our model on the validation set. To verify whether the synthesized results can replace the ground truths with diagnostic acceptability, we feed both real scans and all generations into a classic and commonly used segmentation algorithm, FMRIB software library (FSL) [15]. We classify cerebrospinal fluid (CSF), gray matter (GM), white matter (WM) and optional tumor lesions to obtain the average quantification of the whole brain volume. FSL⁴⁴4https://fsl.fmrib.ox.ac.uk/fsl/fslwiki is a comprehensive library of analysis, specific for functional and structural brain imaging data. We obtain the segmented results under the tissue prior probability templates in a default image space, and therefore there is no guarantee that our segmentation will exactly follow other methods. For evaluation criteria, we use PSNR, SSIM and Dice score (which measures segmentation overlap, with a higher value meaning a better result) to assess the performance of different methods. We demonstrate that the proposed CSC$\ell_{4}$Net exhibits competitive performance with fewer layers (and the corresponding parameters) compared to other deep learning methods.

Inspired by ResNet [9], CSC$\ell_{4}$Net includes seven CSC layers with layered MMD and RM constraints, where the spatial subsampling operation is performed with a stride of 2 in the last two bottleneck layers. Batch normalization is incorporated layer-wise, following each CSC layer, for faster convergence. The last CSC layer lies on top of a global spatial average pooling layer. The bottleneck layer in CSC$\ell_{4}$Net outputs 64, 128, 256 with 3, 2, 2 times repeated stacks, respectively. We use a learning rate of 0.0002, a batch size of 16, a total of 200 epochs, and a sparse regularization parameter of 0.02. We work in 2D slices and employ the Adam optimizer. For the layered MMD parameters, we use the Gaussian kernel with bandwidths equipped as median pairwise squared distances. We use the same network configurations for all datasets.

We compare our CSC$\ell_{4}$Net with the most relevant and state-of-the-art synthesis methods: 1) MIMECS [24]; 2) V-US [28]; 3) cycleGAN [39]; 4) 3D-cGAN [20]; 5) MSGAN [18]. Note that, MIMECS and MSGAN are supervised methods, so we input originally paired data for their training. V-US, cycleGAN and 3D-cGAN are unsupervised approaches, so we input the selected unpaired images for training. For fair comparison, we empirically set all parameters of the compared methods following their recommended data to reach their best performance. Notably, we exactly followed the data processing steps implemented by MIMECS, V-US and 3D-cGAN, where a standard intensity correction and skull-striping were preformed. Therefore, no extra data processes are implemented.

We evaluate the effectiveness of our method by synthesizing images on different datasets. The results of CSC$\ell_{4}$Net and the five baseline methods are compared in Table 1. As further validation, we also apply the synthesized results to segment major tissues in each dataset and summarize the average performances in Table 1. First, we compare our CSC$\ell_{4}$Net against several state-of-the-art synthesis methods on IXI for transforming PD-w MRI to T2-w MRI and vice versa. For better interpretation, we provide the visualization results in Fig. 2. From both Table 1 (top) and Fig. 2, we observe that CSC$\ell_{4}$Net achieves much better performance (visually and quantitatively) than the other five baseline methods in all cases. The average synthesis performances of CSC$\ell_{4}$Net on IXI dataset are {PSNR: 36.45dB, SSIM: 0.9065} and {PSNR: 38.18dB, SSIM: 0.9596} for T2-w $\rightarrow$ PD-w and PD-w $\rightarrow$ T2-w, respectively. CSC$\ell_{4}$Net gains a significant {PSNR, SSIM} performance improvement of {3.34dB, 0.0508} and {2.89dB, 0.0625}, respectively, compared to the best baseline MSGAN. In addition, we report the segmentation accuracy (Dice scores), which are based on the synthesized results, in Table 1. As we intuitively expected, if the model can provide results with data fidelity, the synthesis task is able to deliver practical usable results for further diagnosis. For medical image synthesis-driven segmentation, our method again outperforms the other methods, i.e., MIMECS, V-US, cycleGAN, 3D-cGAN and MSGAN, by a large margin. Our results are more representative than those of other state-of-the-art image synthesis methods on IXI dataset, showing dice improvements of 3.91% and 7.27% for T2-w $\rightarrow$ PD-w and PD-w $\rightarrow$ T2-w respectively, w.r.t the best baseline.

In addition to healthy cases, we also evaluate our algorithm on the pathological datasets (i.e. NAMIC, BraTS). We adopt the same evaluation criteria as on IXI dataset. We provide the qualitative results of our method in Figs. 3-4. These figures demonstrate how BrainGAN can handle schizophrenic and brain tumor synthesis on NAMIC and BraTS dataset, respectively. The proposed method consistently obtains the best performance, particularly when the lesions in the T1-w data appear with much lower contrast than in the FLAIR brain MRI (shown in Fig. 4). In Table 1 (middle and bottom), similarly for PD $\rightleftharpoons$ T2 synthesis, CSC$\ell_{4}$Net outperforms the compared methods with improvements of {1.41dB, 0.0580} and {1.74dB, 0.0339} for T1 $\rightleftharpoons$ T2 on NAMIC dataset, and {3.67dB, 0.0582} and {2.77dB, 0.0314} for T1 $\rightleftharpoons$ FLAIR on BraTS dataset, respectively. Using the defined settings (refer to Section 5.1) and all synthesized data from both the forward and backward transformation, Table 1 compares the segmentation results. As with the Dice scores reported in Table 1 (middle and bottom), our method again gains significant performance enhancements, i.e., increasing by 13.54% for T1 $\rightleftharpoons$ T2, 9.02% for T2 $\rightleftharpoons$ T1, 9.72% for T1 $\rightleftharpoons$ FLAIR, and 3.76% for FLAIR $\rightleftharpoons$ T1. For clarity, Fig. 5 provides detailed comparison results of all methods distributed across different datasets.

Since CSC$\ell_{4}$Net comprises a combination of several components, we perform an extensive ablation study to better understand why it is able to obtain state-of-the-art results. Considering the number of experiments in the above studies, we focus on the case of PD-w $\rightarrow$ T2-w from IXI dataset and report our ablation results in Table 2. We separately adopt {IUN, $\mathcal{L}_{\mathcal{H}}$, $\mathcal{L}_{\mathcal{M}}$} as the additional module under the baseline CSC, and evaluate the effects in terms of image quality and synthesis-based segmentation performance. When each of the components is separately combined with CSC, the PSNR, SSIM, and Dice scores are improved by {2.42dB, 3.73dB, 3.54dB}, {0.05, 0.07, 0.07}, {0.35%, 2.19%, 7.79%}, respectively. With the assistance of $\mathcal{L}_{\mathcal{H}}$ and $\mathcal{L}_{\mathcal{M}}$, the image quality and synthesis-based segmentation have the greatest impact. CSC with different combinations of two modules achieves {7.42dB. 7.05dB, 9.33dB}, {0.16, 0.16, 0.19}, and {13.44%, 17.7%, 26.43%} improvements over the baseline.