{'fraction_non_alphanumeric': 0.03441063888560041, 'fraction_numerical': 0.007813424884574404, 'mean_word_length': 4.747788614198231, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Information-based Disentangled Representation Learning for Unsupervised MR Harmonization

Accuracy and consistency are two key factors in computer-assisted magnetic resonance (MR) image analysis. However, contrast variation from site to site caused by lack of standardization in MR acquisition impedes consistent measurements. In recent years, image harmonization approaches have been proposed to compensate for contrast variation in MR images. Current harmonization approaches either require cross-site traveling subjects for supervised training or heavily rely on site-specific harmonization models to encourage harmonization accuracy. These requirements potentially limit the application of current harmonization methods in large-scale multi-site studies. In this work, we propose an unsupervised MR harmonization framework, CALAMITI (Contrast Anatomy Learning and Analysis for MR Intensity Translation and Integration), based on information bottleneck theory. CALAMITI learns a disentangled latent space using a unified structure for multi-site harmonization without the need for traveling subjects. Our model is also able to adapt itself to harmonize MR images from a new site with fine tuning solely on images from the new site. Both qualitative and quantitative results show that the proposed method achieves superior performance compared with other unsupervised harmonization approaches.


Introduction
Magnetic resonance (MR) imaging is a commonly used non-invasive imaging modality due to its flexibility and good tissue contrast. For the purposes of describing MR imaging analytically, we can think of an MR image as a function (i.e., arXiv:2103.13283v1 [eess.IV] 24 Mar 2021 imaging equation) of the anatomy being imaged and the associated acquisition parameters [3]. By changing the acquisition parameters or underlying imaging equations, MR images with different contrasts can be generated. To take advantage of this flexibility, MR images of the same anatomy with different contrasts are often acquired in a single session. For example, T 1 -weighted (T 1 -w) images are typically used to achieve balanced contrast between T 2 -weighted (T 2 -w) images [3]. However, a consequence of this flexibility is that there is no standardization when it comes to MR contrasts. For example, both magnetization-prepared rapid gradient echo (MPRAGE) and spoiled gradient echo (SPGR) are commonly used T 1 -w images with very different visual appearances. This lack of standardization makes machine learning (ML) models trained on MPRAGE images often fail on SPGR images and underperform on MPRAGE images acquired by different scanners or with slightly different parameters [19]. Scanner software and calibration differences can also contribute to this effect.
The issue of contrast variation is commonly seen in multi-site studies, where a trained model degrades in performance when tested on data from another site (i.e., the domain shift problem). This is because ML based methods assume the training and testing are conducted on data drawn from the same distribution (domain). This is not the case for MR images acquired from different sites, scanners, or with differing imaging parameters. For example, T 1 -w images acquired from two scanners with different configurations should obviously be treated as two domains. However, T 1 -w and T 2 -w images acquired from the same scanner should also be considered as coming from two domains. MR image harmonization [6] alleviates domain shift by treating the problem as an image-to-image translation (IIT) (or synthesis) task, where the goal is to translate image contrasts between domains (e.g., T 1 -w images acquired from different scanners). MR harmonization can be separated into two categories: supervised and unsupervised. In the supervised setting, MR images of the same anatomy across multiple sites are available; these are known as traveling subjects or inter -site paired data. These images are used to train intensity transformations between sites. However, traveling subjects are impractical in large-scale multi-site harmonization tasks. Unsupervised harmonization methods do not require intersite paired data. Instead, these methods often rely on domain-specific models (e.g., intensity transformations and discriminators). We outline recent related work in IIT and unsupervised domain adaptation (UDA), below.
IIT learns a transformation of images between different domains, e.g., MR to CT [24,26] or T 1 -w to T 2 -w [28]. In both supervised and unsupervised IIT, the goal is to approximate the joint distribution drawn from the (two) domains. Supervised IIT methods use pixel-to-pixel reconstruction error during model training. Recent unsupervised IIT work has explored learning disentangled representations [10,14,25], the idea being to tease apart the domain-invariant and domain-specific information in the representation. As an unsupervised IIT method, unsupervised harmonization faces four challenges. First, the lack of inter-site paired data along with the coupling theory [15] tells us there are infinitely many possible joint distributions given the two marginal distributions. The proposed method consists of a single θ-encoder, a β-encoder, a decoder, and a β-discriminator that work on all domains.
x and x share the same contrast but have different anatomy. The same networks work on all sites.
Therefore, to learn a meaningful harmonization model (joint distribution), further constraints are required. Cycle-consistency is commonly assumed in unsupervised IIT [10,16,27]. However, there is no theory that supports the validity of this assumption. Second, the lack of inter-site paired data means that pixel-to-pixel regularization cannot be easily achieved. Domain-specific discriminators are commonly used in many unsupervised IIT methods [10,16,27]. For harmonization, performance will heavily rely on the discriminators' "judgement" during training, and geometry shift is a common drawback in unsupervised harmonization. Third, the use of site specific-models means that the size of the harmonization model grows with the number of sites. Lastly, most existing harmonization approaches are not able to work on domains not seen in the training data. When testing and training domains differ, most methods require retraining with images from all domains, which is not practical. In general, the goal of UDA is to learn a model from a source domain with labeled data and apply it to a target domain with unlabeled data during testing [9,12,20,23]. Without special design, domain shifts between training and testing can cause a performance drop. Different from IIT, which aims at mapping image data across domains, a UDA tries to adjust the model during testing. A UDA allows the model to detect a domain and then provide self-supervision for model adjustment during testing. UDAs are especially helpful in medical imaging, where training and testing data are likely to come from different sources.
We propose an unsupervised harmonization approach, CALAMITI (contrast anatomy learning and analysis for MR intensity translation and integration), which integrates the merits of both IIT and UDA. Building upon the recent work in [7], we use the routinely acquired multi-contrast MR images within each site (called intra-site paired data) during the same imaging session. However, as we discuss in Section 2.2, this technique alone does not provide a globally disentangled latent space and cannot be easily generalized to data from a new site. CALAMITI is an improved, theoretically grounded, unsupervised harmonization approach based on an information bottleneck (IB) [22] that learns a global, disentangled latent space of anatomical and contrast information and can be easily adapted to a new testing site using only the new data. To our knowledge, this is the first work that overcomes the four challenges in unsupervised harmonization. First, by taking advantage of the intra-site paired data, the proposed method solves an unsupervised IIT problem in a supervised way, avoiding introducing any extra constraint (e.g., cycle-consistency) on the model and achieving better pixel-to-pixel regularization. Second, it has a unified structure for multi-site harmonization, which means that model size does not grow with the number of sites. Third, it provides a global latent space for all training data by encouraging a consistent description of the anatomy. Finally, it is able to adapt to a new site without any retraining on the original data. For all of this work, we also provide a theoretical explanation of the disentangled latent space using IB theory.

The disentangling framework
Our method uses multi-contrast MR images of the same subject within each site (intra-site paired data) to train a cross-site harmonization model. Here, we emphasize the relationship between "site", "domain", and "MR contrast". As shown in Fig. 1(a), given T 1 -w and T 2 -w images from Sites A and B, there are four domains X A1 , X A2 , X B1 , and X B2 , where an alphabetical index indicates site and a numerical index represents contrast. Our goal is to learn a disentangled representation that captures anatomical and contrast information from the input images. Following the notation in [7], the anatomical representation (β) is domain-invariant and the contrast representation (θ) has some domain-specific information. Thus, combining the β from one site with the θ from another allows harmonization across sites. To learn the disentangled representation, we solve the inter-site unsupervised IIT problem based on training from intra-site supervised data. Figure 1(b) outlines our framework, which consists of a θ-encoder, a β-encoder, a decoder, and a β-discriminator that work on all domains. Here, we outline the high-level training strategy using the proposed framework, and we highlight the role of our β-discriminator in Section 2.2. Each site has paired T 1 -w and T 2 -w images-with different imaging parameters at each site-which train a disentangled network in a supervised IIT way. For example, intra-site paired images x A1 and x A2 of the same subject imaged at Site A (in our case T 1 -w and T 2 -w images from Site A) are sent to a β-encoder to extract anatomical information. These images have the same anatomy, so the extracted anatomical representations β A1 and β A2 should be the same. To encourage similarity of β, we randomly shuffle between β A1 and β A2 before decoding as well as introduce a small l 1 loss between the two β's. To prevent contrast representation θ from capturing anatomical information, we provide the θ-encoder with an image of different anatomy (but the same contrast), x A1 . This is achieved by selecting a different slice from the same volume as x A1 . The decoder takes the randomly selected anatomical representation (β A ), concatenated with a θ to generate a synthetic image. The contrast of the synthetic image depends on which θ has been chosen. The same β-encoder, θ-encoder, decoder, and β-discriminator are used for all training sites to achieve a unified structure.
Our β-encoder and decoder both have a U-Net like architecture with four downsampling layers, while the θ-encoder is four convolutional layers followed by three fully connected layers. β is one-hot encoded with multiple channels and the same spatial extents as the input image. For gradients to backpropagate through the one-hot encoded β, we adopt and implement the trick introduced in [7,18], wherein β's are calculated using a Gumbel-softmax layer. One-hot encoding β restricts its capacity, encouraging β to capture only anatomical information.

Creating a consistent anatomical space
To learn a consistent anatomical space for all sites, we introduce a β-discriminator to our framework. Because our training strategy only uses supervised IIT within each site-with no supervision between sites-the β-encoder could possibly learn a distinct β space for each site. In this case, the β's and θ's are disentangled within each site, and we refer to it as a locally disentangled latent space. This is not desirable in harmonization, as combining these β's and θ's across sites would not be ideal. To avoid this, we must encourage the learned β's of all sites to be from the same distribution (i.e., β's and θ's are globally disentangled). This leads us to use a one-class discriminator on β space to encourage distribution similarity. No matter which site an input β comes from, the β-discriminator learns to distinguish whether the β is from Site A or not, further pushing θ to describe the difference between sites as well as different MR contrasts.
The proposed framework solves a number of outstanding problems. First, it performs unsupervised IIT by using supervised IIT during training. This avoids geometry issues inherent in unsupervised IIT by penalizing pixel-to-pixel error during training our framework. This is more effective than introducing a cycle consistency constraint, as cycle consistency still allows a model to learn "circlesquare-circle". Second, our unified harmonization structure means we have one β-encoder, one θ-encoder, one decoder, and one β-discriminator that work on all domains; Section 3 includes an 8-site experiment using this structure. This saves a significant number of parameters when there are many sites-e.g., the CycleGAN [27] requires O(N 2 ) image translation models and N discriminators, for N sites. Third, our discriminator functions on the latent variables β instead of the harmonized images. This encourages our decoder to act like a universal imaging equation, generating synthetic MR images based on any β and θ within the distribution of the training data. This strategy combined with the unified structure, makes our model more robust when there are more sites involved - ----Dewey et al. [7] --CALAMITI during training. Lastly, our discriminator makes a one-class decision: whether an input β is from Site A or not. We show in Section 2.4, that this allows our model to adapt to a new testing site after fine tuning. Table 1 provides a summary comparison of the proposed method with other unsupervised IIT approaches.

Learning from an information bottleneck
By providing the θ-encoder with an image of a different anatomy (but the same contrast) as provided to the β-encoder, we create a conditional variational autoencoder (CVAE) [21] even though the condition variable β is not connected to the θ-encoder. All of these strategies help us to limit the information that can be passed through each of the β and θ channels, which we now show theoretically forms an IB given the model design.
IB theory [22] describes a constrained optimization problem with the goal of learning a compressed latent representation Z such that the mutual information (MI) between Z and the task variable Y is maximized while Z captures minimum information about the input variable X. Mathematically, this can be formulated as Z * = arg min Z I(Z; X) − λI(Z; Y ), where I(·; ·) is the MI and λ is a hyper-parameter. IB theory is closely related to the variational autoencoders (VAEs) and disentangled representation learning (cf. [4,5]). [2] showed that IB is a more general case of the VAE objective.
Our network structure forms a CVAE. To better illustrate the IB in a general setting, we slightly modify the notation (i.e., remove site index), and highlight the CVAE structure. Theorem 1. It can be shown that optimizing our network structure is equivalent to solving a conditional IB problem, i.e., θ * = arg min θ I(X ; θ) − λI(X; θ|β). (1) The proof of Theorem 1 is similar to that in [2], despite the fact that we are solving a CVAE problem with conditionβ. The reason why the first term of Equation 1 is free from conditionβ is becauseβ is disentangled from θ and will be ignored in calculating θ. An intuitive understanding of Equation 1 is that the proposed method learns a contrast representation θ that captures minimum information about the input variable X , while the (conditional) MI between θ and the target variable X is maximized. Since the shared information between variables X and X is the contrast, we would expect θ to capture only contrast information about X after training. Equation 1 can be re-organized as a KL divergence term and a reconstruction term (similar to the CVAE loss [21]), and directly optimized as network loss functions, i.e., where p(θ) is a zero mean unit variance Gaussian distribution. p(θ|x ) and p(x|θ,β) can be modeled by a probabilistic θ-encoder and decoder, respectively. This KL divergence term encourages a bounded value for θ, which even if lightly weighted restricts its possible expression. Accordingly, our network loss functions include a reconstruction loss for supervised IIT (l 1 and perceptual loss [11]), a KL divergence term on θ, an adversarial loss between β-discriminator and β-encoder, and a similarity loss on β between contrasts 4 .

Domain adaptation
Suppose the proposed model was pretrained on Sites A and B, and the goal is to harmonize a new site, Site C, to Site A or B without a retraining that includes data from all sites. As in regular CALAMITI training, the supervised IIT is conducted on images from Site C. However, the decoder and β-discriminator weights are frozen, and only the last few layers of the β-and θ-encoders are updated. We rely on the assumption that our decoder is well-generalized in previous training to produce a variety of contrast images. Our β-discriminator guides the β-encoder to generate β's that follow the previously learned distribution of β; avoiding a Site C specific β space. Thus our β-discriminator acts as a domain shift detector on β space like other UDA methods [9,23]. Until the β-encoder generates β for Site C that is less distinguishable from previously learned β, the β-discriminator will produce a loss for mismatching β's. As the decoder weights are frozen in fine tuning, once β's from Site C match previously learned β distribution, we are able to harmonize images between new and previous sites. Thus by combining β C1 with θ A1 , we can harmonize anatomy imaged at Site C with the corresponding contrast from Site A. Thus Site C can be harmonized to any site included in the original training.

Datasets and preprocessing
MR images, both T 1 -w and T 2 -w, acquired from eight different sites (scanners) were used to evaluate the proposed method. Data source, scanner details, and  imaging parameters are provided in Table 2. Sites A thru F are healthy controls, Sites G and H imaged multiple sclerosis (MS) subjects. Sites D and E have similar scanners and sequences and thus similar contrast. Images underwent preprocessing including N4 inhomogeneity correction, super-resolution for 2D acquired scans, registration to 1 mm 3 MNI space, and white matter (WM) peak normalization. The center 60 axial slices with spatial dimension of 224 × 192 were extracted for each subject.

Qualitative and quantitative evaluation
For Sites C, D, E, and F , there are 10 subjects (600 axial slices) used from each site for training. For the remaining sites, 20 subjects from each site are used in training. There are longitudinal scans in the OASIS3 [13] dataset acquired by different scanners with a short period between visits. These scans are held-out and used as traveling subjects for quantitative evaluation in testing. Specifically, there are seven traveling subjects between Site C and D, and ten traveling subjects between Site E and Site F . The average days between two visits for Sites C/D and Sites E/F are 162 days and 13 days, respectively. In our experiments, β is a four-channel one-hot encoded map, with spatial dimension the same as the image, while θ is a two-dimensional vector. Figure 2 shows harmonized MR images from the eight sites. The mean θ value of all testing images at each site was used to harmonize images between sites. Our analysis has focused on the T 1 -w images, as these represent the images with the greatest disparity across the imaging sites and the primary contrast for neuroimaging analysis. With regard to the T 2 -w images, we achieve similar image quality as the T 1 -w images. Figure 3 shows θ values of the held-out T 1 -w testing images. We observe that Sites D and E overlap, which is good as the sites have identical scanner configurations. Second, images acquired using different methods (MPRAGE, and MEMPRAGE) are separated. Third, the distance in θ space corresponds to human perception; sites with visually similar contrast have closer clusters in θ space. For example, although images in Sites D, E, F , and H are all MPRAGE images, θ points are more closely clustered in Sites D, E and F , than Site H. In Table 3 and Fig. 4, we show qualitative and quantitative comparison of different unsupervised IIT methods. The traveling subjects are used in the quantitative comparison. Specifically, our baseline is the MR images without harmonization (No har), and we compare the structural similarity index measurement (SSIM) and peak signal-to-noise ratio (PSNR) of histogram matching (Hist), Fig. 4. Visual comparison of different harmonization approaches. An MR image from Site G is harmonized to Site A (see Fig. 2 for reference images). Yellow arrows indicate geometry change. CycleGAN [27], Dewey et al. [7], and CALAMITI. Histogram matching is a nontraining method, while the other approaches are ML-based. To select a reference image for histogram matching, we first randomly chose a volume, then selected the same slice number as our source image. For a fair comparison, we consider two training scenarios for CALAMITI: only include two sites (the source and target site in harmonization) or include all eight sites during training. Paired Wilcoxon signed rank tests were conducted between CALAMITI (two sites) and each comparison method under each performance measurement. Results show that CALAMITI has significantly (p < 0.001, N = 420 for Sites C and D, N = 600 for Sites E and F ) better performance over all comparison methods, except for the PSNR of Site D → C. The null hypothesis is that the difference of SSIM or PSNR between the two sites is from a distribution with zero median. Interestingly, CALAMITI has slightly better performance when more sites are used in training. We suggest two possible reasons for this. First, CALAMITI has a unified structure, which makes the whole model less likely to overfit. Second, our β-discriminator is a single class discriminator, so more sites should improve the robustness of the discriminator.

Domain adaptation
We provide both qualitative and quantitative results to test domain adaptation (DA). For qualitative comparison, we trained our harmonization model using Sites G and H and applied the trained model on Site B in testing. The task is to harmonize Site B to Site G or Site H without using data from Sites G or H. As shown in Fig. 5, after fine tuning only on a subset of Site B images, the model is able to adjust itself to produce reasonable harmonization results. Table 4 shows quantitative DA results. In each experiment, the proposed method was trained on the corresponding target site and Site G, while testing and fine tuning were conducted solely on the source site. For example, when evaluating DA in Site C → D, the model was trained on Site D and G, and Site C was used as a previously unseen site for testing and fine tuning. Results show that the proposed method achieves significant (p < 0.001, N = 420) improvements after DA.

Discussion and Conclusion
Both qualitative and quantitative results from our eight-site experiment show the potential of the proposed method in large-scale multi-site studies. There are some limitations. First, the requirement of intra-site paired images in training could potentially restrict some applications-pediatric data for example-where acquiring multi-contrast images is not practical. Second, in our experiments, we only used paired T 1 -w and T 2 -w images. However, the proposed method can be extended to include more contrast MR images such as fluid-attenuated inversion recovery (FLAIR) images to achieve a better disentanglement. Third, our experiments on MS patients show that the proposed method does not produce a satisfactory harmonization result on WM lesion areas. We hypothesize that inclusion of FLAIR images would improve this. Fourth, although satisfactory results have been observed in Section 3.3, the way we used our β-discriminator to update the β-encoder during domain adaptation is theoretically flawed. According to Goodfellow et al. [8], the generator and discriminator must be updated jointly to achieve optimal performance. We view all these limitations as opportunities for future improvements.
In conclusion, we propose an unsupervised MR harmonization approach, CALAMITI, which integrates merits from both unsupervised IIT and UDA, and is grounded in information bottleneck theory. Our model learns a universal imaging equation and a disentangled latent space without inter-site paired data. In contrast to many unsupervised harmonization methods, our model takes advantages of the intra-site paired data to prevent the geometry shift problem. Experiments show that the proposed approach achieves state-of-the-art harmonization performance both visually and in terms of SSIM and PSNR.