For magnetic resonance (MR) images, we can view image synthesis as learning an intensity transformation between two differing contrast images, e.g., from T1-weighted (T1-w) to T2-weighted (T2-w) or FLuid Attenuated Inversion Recovery (FLAIR). Synthesis can generate contrasts not present in the data set—i.e., image imputation—which are useful for image processing applications such as registration and segmentation [1, 2]. The transformation need not be limited to MR images; an example application is MR to computed tomography (CT) registration where it has been shown to improve accuracy when the moving image is synthesized to match the target image’s contrast . Other examples include multi-contrast skull-stripping for MR brain images , which performs better with synthesized T2-w images when the original T2-w images are unavailable.
Methods to carry out image synthesis include sparse recovery-based methods 
, random forest regression[6, 7], registration [8, 9], and deep learning [10, 11]. Evidence suggests that accurate synthesis is heavily dependent on a standard intensity scale across the sample of images used in the training procedure. That is to successfully train a synthesis algorithm the training and testing data must have similar intensity properties (e.g., the mean intensity of white matter should be the same for all input images). This is a problem in MR synthesis since MR images do not have a standard intensity scale.
2 New Work to be Presented
In this paper, we explore seven methods to normalize the intensity distribution of a sample of MR brain images within each of three contrasts (T1-w, T2-w, and FLAIR). We then quantitatively compare their performance in the task of synthesizing T2-w and FLAIR images from T1-w contrasts using three synthesis algorithms. We show results that suggest intensity normalization as a preprocessing step is crucial for consistent MR image synthesis. In this abstract, we only present T1-w to FLAIR synthesis results over six subjects. In the final paper we will include both T1-w to FLAIR and T1-w to T2-w synthesis over a cohort of 20 subjects. We also provide a terse summary of the normalization techniques here and will provide detailed descriptions in the final paper.
We first outline the seven intensity normalization algorithms considered in this paper, then describe the three different synthesis methods. We denote an MR brain image with for . Let and , denote masks for the white matter (WM) and brain mask, respectively, with . ( is segmentation dependent and described below.) and
are the mean and standard deviation, respectively, of the intensities withinwhere is one of , , or .
In this section, we first describe the seven intensity normalization algorithms considered in this paper, namely: 1) Z-score, 2) Fuzzy C-Means (FCM)-based, 3) Gaussian mixture model (GMM) based, 4) Kernel Density Estimate (KDE) based, 5) Piecewise linear histogram matching (HM)[12, 13], 6) WhiteStripe , and 7) RAVEL 
. We then describe three different synthesis routines: 1) polynomial regression, 2) random forest regression, and 3) deep neural network based synthesis. For the following subsections, letbe the MR brain image under consideration where for , the dimensions of , and let be the corresponding brain mask (i.e., the set of indices corresponding to the location of the brain in ).
In Table 1, we provide an overview of the normalization techniques included in this paper. For segmentation based normalization we use three different segmentation algorithms: 1) three class fuzzy c-means (), 2) three class Gaussian mixture model (), and 3) a kernel density estimate () based approach. Each of which provides a definition of WM which is used to define and compute and in the normalization formulation. WhiteStripe 
uses WM values within 5% of the WM peak with respect to the cumulative distribution function (CDF). RAVEL normalization builds upon WhiteStripe by removing unwanted technical variation from a sample of images with a voxel specific correction. RAVEL requires registration to a normalized space, for which we used SyN .
Piecewise linear histogram
matching [12, 13]
based on deciles.
|WhiteStripe||where of the WM CDF|
|RAVEL||where is a voxel dependent correction factor for unwanted variation.|
In the following sections, we will briefly overview the intensity normalization algorithms used in this experiment. Code for the following intensity normalization algorithms is at: https://github.com/jcreinhold/intensity-normalization.
Z-score normalization uses the brain mask for the image to determine the mean and standard deviation of the intensities inside the brain mask. Then the Z-score normalized image is
FCM-based normalization uses fuzzy c-means to calculate a white matter (WM) mask of the image . This WM mask is then used to normalize the entire image to the mean of the WM. The procedure is as follows. Let be the WM mask for the image , i.e., is the set of indices corresponding to the location of the WM in the image . Then the WM mean is . and the FCM-based normalized image is
where is a constant that determines the WM mean after normalization. In this experiment, we use three-class fuzzy c-means to get a segmentation of the WM over the brain mask for the T1-w image and we arbitrarily set .
GMM-based normalization fits a mixture of three normal distributions to the histogram of intensities inside the brain mask. The meanof the mixture component associated with the WM is then used in the same way as the FCM-based method, so the GMM-based normalized image is
where is a constant that determines the WM mean after normalization. The WM mean is determined by picking the mixture component with the maximum intensity mean for T1-w images, the middle intensity mean for FLAIR images, and the minimum intensity mean for T2-w images.
3.1.4 Kernel Density Estimate-based
KDE-based normalization estimates the empirical probability density function (pdf) of the intensities ofover the brain mask using the method of kernel density estimation. In our experiment, we use a Gaussian kernel. The kernel density estimate provides a smooth version of the histogram which allows us to more robustly pick the maxima associated with the WM via a peak finding algorithm. The found WM peak is then used to normalize the entire image, in much the same way as FCM-based normalization. Namely,
where is a constant that determines the WM peak after normaalization. The WM peak is determined in T1-w and FLAIR by picking the peak associated with the greatest intensity (for FLAIR, this is due to the inability to distinguish between the WM and GM peaks) and for T2-w images the WM peak is determined by the highest peak.
3.1.5 Piecewise Linear Histogram Matching
Piecewise linear histogram matching (which we denote as HM for brevity) addresses the normalization problem by learning a standard histogram for a set of contrast images and linearly mapping the intensities of each image to this standard histogram. The standard histogram is learned through averaging pre-defined landmarks of interest on the histogram of a set of images. In Shah et al., the authors demonstrate good results with this method by defining landmarks as intensity percentiles at percent (where the intensity values below 1% and above 99% are extrapolated from the [1,10] and [90,99] percent intervals). We use these landmarks in our method and arbitrarily set the range of the standard scale to . The intensity values of the set of images are then mapped piecewise linearly to the learned standard histogram along the landmarks. For further detail into the method see Nyúl et al. and Shah et al..
WhiteStripe intensity normalization  performs a Z-score normalization based on the intensity values of normal appearing white matter (NAWM). The NAWM is found by smoothing the histogram of the image and selecting the highest intensity peak for T1-w images (the peaks for the other contrasts are determined in the same way as described in the KDE section). Let be the intensity associated with this peak. The “white stripe” is then defined as the 10% segment of intensity values around . That is, let be the cdf of the specific MR image inside its brain mask , and define . Then, the white stripe is defined as the set
Let be the sample standard deviation associated with . Then the WhiteStripe normalized image is
RAVEL normalization  adds an additional normalization step to WhiteStripe by removing unwanted technical variation (defined below) from a sample of images. Following the notation in the original paper, the method assumes that cerebrospinal fluid (CSF) is associated with technical variation, and—after WhiteStripe normalization—the CSF intensities can be written as
where is an matrix of CSF intensities, represents the unknown technical variation, and is a matrix of the residuals. The CSF intensity values in are determined by deformably co-registering the images, finding a CSF mask for each deformably registered image, and taking the intersection across all the masks.
We then use singular value decomposition to write. Then is an
matrix of right singular vectors and we can useright singular vectors to form a linear basis for the unwanted factors , where is the unknown true rank of . That is, we use as a surrogate for , where is the subset of columns of
collected into a matrix. We then do voxel-wise linear regression to estimate the coefficients. The RAVEL normalized image is then defined as
where are the coefficients of unwanted variation associated with the voxel found through linear regression. In our experiments, we follow the original paper and fix 111The first right singular vector is highly correlated (¿95%) with the mean intensity of the CSF.. For deformable registration, we use SyN to register all images to one image in the data set.
Image synthesis can be thought of as a regression on image intensities, i.e., learning a parametric or non-parametric mapping from one contrasts intensity distribution to another contrasts intensity distribution. We explore three image synthesis methods: 1) polynomial regression (PR), 2) random forest (RF) regression, and 3) deep neural network (DNN)-based synthesis.
PR samples 100,000 voxels in and train a regressor from the intensities of and its six connected neighbors in the source image to the intensity of in the target image. We use a third-order polynomial regressor to learn the mapping from the source intensity to the target intensity. Our RF regression is based on a simplified model of the work of Jog et al. . Training is done from 100,000 sampled voxels with the image patch using voxels in the six primary directions at 3, 5, and 7 indices away from the center voxel and no multi-scale framework. This simplified RF model also only synthesizes the corresponding center voxel of the patch under consideration.
Our DNN approach uses a 4-level U-net [ronneberger2015miccai], extracting patches from axial, sagittal, and coronal orientations to learn the synthesis, in a similar architecture to Zhao et al. 
. We use instance normalization and leaky ReLUs with parameter 0.2 as the activation function since Z-score, WhiteStripe, and RAVEL allow for negative values in the images. It is a good approximation for the state-of-the-art in image synthesis. Complete detail of all three synthesis approaches will be included in the final paper.
Image synthesis can be described as a regression on the intensities of the images, i.e., learning a parametric or non-parametric mapping from one contrasts intensity distribution to another contrasts intensity distribution. In this section we describe three methods of image synthesis: 1) polynomial regression (PR), 2) random forest regression (RF), and 3) deep neural network (DNN)-based synthesis.
3.2.1 Polynomial Regression
For polynomial regression, we randomly select 100,000 voxels inside the brain mask. For the source images, we extracted patches around each of these voxels where the patches include the center voxel and its six neighbors. For the target images, we extract only the corresponding center voxel. We extract the patches in this way across all images, so for images we have an feature matrix for the source images and an
feature matrix for the target images. We use a third-order polynomial as the regressor to learn the mapping from the source feature matrix to the target feature matrix. We use this naïve model to provide a low-variance baseline for image synthesis methods.
3.2.2 Random Forest Regression
Similar to polynomial regression, in random forest regression—inspired by Jog, et al.—we randomly select 100,000 voxels inside the brain mask. For the source images, we extracted patches that comprise the center voxel, its six neighbors, and the voxels in the six primary directions at 3, 5, and 7 voxels away from the center. For the target images, we extract only the corresponding center voxel. We extract the patches in this way across all images, so for images we have an feature matrix for the source images and an feature matrix for the target images. For the random forest regressor that learns the mapping between the source feature matrix and the target feature matrix, we set the number of trees to 60 and the number of samples in a leaf node to 5.
We use a 4-level U-net and extract patches from axial, sagittal, and coronal orientations to learn the synthesis. Patches are extracted in this fashion for data augmentation. We use instance normalization and leaky ReLUs with parameter 0.2 as the activation function since Z-score, WhiteStripe, and RAVEL allow for negative values in the images. The architecture follows Zhao, et al.
who used a similar structure for a synthesis task. We trained the network for 100 epochs for all sets of normalized images excluding the unnormalized images with which we trained the DNN for 400 epochs. This discrepancy in the number of epochs used is due to a failure of convergence observed in the first 100 epochs for the unnormalized images.
3.3 Quality Assessment
We use three different metrics to quantitatively determine the performance of the synthesis result. Note that all three metrics compare the result to the ground truth images which were not used in training any synthesis methods. The metrics are: 1) normalized cross-correlation (NCC), 2) mean structural similarity (MSSIM) , and 3) mutual information (MI). We use these metrics as opposed to MSE or PSNR as the data have been scaled to different ranges, making MSE and PSNR not easily comparable across normalization routines.
For evaluation, we use 18 data sets from the Kirby-21 data set . All of the subjects for the data sets are verified to be healthy subjects. From these 18 data sets, we use the T1-w, T2-w, and FLAIR images. All the images are resampled to 1mm, bias field corrected with N4 , and each T2-w and FLAIR image is affinely registered to the corresponding T1-w image with the ANTs package . The brain mask for the images are found with ROBEX  and the mask is used during normalization and applied to the images before synthesis such that the background is zero in all the images.
show the mean and the bootstrapped 95% confidence interval of the T1-to-FLAIR and T1-to-T2 synthesis, respectively, for the quality metrics averaged over all testing data sets, for every normalization scheme and synthesis algorithm. We use the Wilcoxon signed-rank test to compare the distributions of each normalized method, for all metrics, against the corresponding unnormalized results per synthesis algorithm. We use a statistical significance level ofand show that this threshold is met in Figs. 1 and 2 with an asterisk above the corresponding bar. Figures 3 and 4 show results for the various synthesis algorithms with unnormalized training data (denoted raw) and normalized training data use the FCM approach.
The experiments show that synthesis results are robust to the choice of normalization algorithm, which are stable around the same levels across all metrics. This qualitative result, observed in Figs. 1 and 2, is reinforced with statistical tests. We use the Wilcoxon signed-rank test (with Bonferroni correction) to show statistically significant difference between any of the presented normalization algorithms for each metric (); however, no normalization algorithm consistently met this threshold for any metric with any synthesis algorithm in either T1-to-FLAIR or T1-to-T2 synthesis. An interesting finding is that—in T1-to-T2 synthesis—the random forest regressor qualitatively performs more robustly on unnormalized data, but both the DNN and polynomial regression methods fail; in terms of NCC, the DNN synthesis has zero mean because of negative correlation in some of the testing results.
Failure cases of synthesis in T1-to-FLAIR and T1-to-T2 unnormalized images are shown in Figs. 3 and 4, respectively, which can be compared to the successfully synthesized FCM-normalized images in the same figures. We discuss in the following section.
5 Discussion and Conclusion
We have shown that: 1) synthesis methods are substantially improved with the addition of an intensity normalization pre-processing step, especially DNN synthesis; 2) synthesis is robust to the choice of normalization method as we see no statistically significant difference in the presented normalization methods.
results from the histogram of a particular input T1-w image being different than the majority of T1-w images the model was trained on. In this case, the problem histogram is compressed such that the grey matter peak was nearly aligned with the average location of the WM peaks for all but one of the training set (where the outlier on the training set also has the grey matter peak in the vicinity of the WM peak average for the training set).
The fact that we fail to synthesize unnormalized images correctly in the best case scenario—all of our training and testing images came from the same cohort acquired on the same scanner with the same pulse sequence and all of the images are of healthy patients—points to the importance of intensity normalization as a preprocessing step in any synthesis pipeline. While the highlighted failure case is remarkable, the synthesized versions of the remaining images also exhibit more subtle failure. Specifically, we see poor correspondence in intensities between slices. That is, if you scan through the images on the plane through which the image was synthesized (in this case axial), the result appears like a reasonable synthesis; however, when the image is viewed in the saggital plane we see significant variation in the intensities of neighboring slices and this variation is not observed in the synthesis results of normalized images (see Fig. 5 for an example). While this slice-to-slice variation is partly due to using a 2D synthesis method, 2D synthesis is commonly used in state-of-the-art synthesis methods[2, 10, 11]. Since the DNN performs better across all metrics when the images are normalized, normalization is suggested as a pre-processing step before training or testing any sort of patch-based DNN.
Acknowledgements.This work was supported in part by the NIH/NINDS grant R01-NS070906 and by the National MS Society grant RG-1507-05243.
-  Iglesias, J. E., Konukoglu, E., Zikic, D., Glocker, B., Leemput, K. V., and Fischl, B., “Is synthesizing MRI contrast useful for inter-modality analysis?,” in [MICCAI ], (8149), 631–638 (2013).
-  Huo, Y., Xu, Z., Bao, S., Assad, A., Abramson, R. G., and Landman, B. A., “Adversarial synthesis learning enables segmentation without target modality ground truth,” in [2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018) ], 1217–1220 (April 2018).
-  Roy, S., Carass, A., Jog, A., Prince, J. L., and Lee, J., “MR to CT registration of brains using image synthesis,” SPIE Medical Imaging 9034 (2014).
-  Roy, S., Butman, J. A., and Pham, D. L., “Robust skull stripping using multiple MR image contrasts insensitive to pathology,” NeuroImage 146, 132–147 (2017).
-  Roy, S., Carass, A., and Prince, J. L., “Magnetic resonance image example-based contrast synthesis,” IEEE Transactions on Medical Imaging 32(12), 2348–2363 (2013).
-  Jog, A., Carass, A., Roy, S., Pham, D. L., and Prince, J. L., “Random forest regression for magnetic resonance image synthesis,” Medical Image Analysis 35, 475–488 (2017).
-  Zhao, C., Carass, A., Lee, J., Jog, A., and Prince, J. L., “A supervoxel based random forest synthesis framework for bidirectional MR/CT synthesis,” Simulation and Synthesis in Medical Imaging (SASHIMI) 10557 LNCS(1), 33–40 (2017).
-  Lee, J., Carass, A., Jog, A., Zhao, C., and Prince, J. L., “Multi-atlas-based CT synthesis from conventional MRI with patch-based refinement for MRI-based radiotherapy planning,” in [SPIE Medical Imaging (SPIE-MI 2017) ], 10133 (2017).
-  Cardoso, M. J., Sudre, C. H., Modat, M., and Ourselin, S., “Template-based multimodal joint generative model of brain data,” in [Information Processing in Medical Imaging ], 9123, 17–29 (2015).
-  Wolterink, J. M., Dinkla, A. M., Savenije, M. H., Seevinck, P. R., van den Berg, C. A., and Išgum, I., “Deep MR to CT synthesis using unpaired data,” in [Lecture Notes in Computer Science ], 10557 LNCS, 14–23 (2017).
-  Chartsias, A., Joyce, T., Giuffrida, M. V., and Tsaftaris, S. A., “Multimodal MR Synthesis via Modality-Invariant Latent Representation,” IEEE Transactions on Medical Imaging 37(3), 803–814 (2018).
-  Nyúl, L. G., Udupa, J. K., and Zhang, X., “New Variants of a Method of MRI Scale Standardization,” IEEE Transactions on Medical Imaging 19(2), 143–150 (2000).
-  Shah, M., Xiao, Y., Subbanna, N., Francis, S., Arnold, D. L., Collins, D. L., and Arbel, T., “Evaluating intensity normalization on MRIs of human brain with multiple sclerosis,” Medical Image Analysis 15(2), 267–282 (2011).
-  Shinohara, R. T., Sweeney, E. M., Goldsmith, J., Shiee, N., Mateen, F. J., Calabresi, P. A., Jarso, S., Pham, D. L., Reich, D. S., and Crainiceanu, C. M., “Statistical normalization techniques for magnetic resonance imaging,” NeuroImage: Clinical 6, 9–19 (2014).
-  Fortin, J. P., Sweeney, E. M., Muschelli, J., Crainiceanu, C. M., and Shinohara, R. T., “Removing inter-subject technical variability in magnetic resonance imaging studies,” NeuroImage 132, 198–212 (2016).
-  Leek, J. T. and Storey, J. D., “Capturing heterogeneity in gene expression studies by surrogate variable analysis,” PLoS Genetics 3(9), 1724–1735 (2007).
-  Avants, B. B., Epstein, C. L., Grossman, M., and Gee, J. C., “Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain,” Medical Image Analysis 12(1), 26–41 (2008).
-  Ronneberger, O., Fischer, P., and Brox, T., “U-Net: Convolutional Networks for Biomedical Image Segmentation,” in [MICCAI 2015 ], Lecture Notes in Computer Science 9351, 234–241, Springer Berlin Heidelberg (2015).
-  Zhao, C., Carass, A., Lee, J., He, Y., and Prince, J. L., “Whole brain segmentation and labeling from CT using synthetic MR images,” Machine Learning in Medical Imaging (MLMI) 10541, 291–298 (2017).
-  Wang, Z., Bovik, A. C., Sheikh, H. R., and Simoncelli, E. P., “Image quality assessment: From error visibility to structural similarity,” IEEE Transactions on Image Processing 13(4), 600–612 (2004).
-  Landman, B. A., Huang, A. J., Gifford, A., Vikram, D. S., Lim, I. A. L., Farrell, J. A., Bogovic, J. A., Hua, J., Chen, M., Jarso, S., Smith, S. A., Joel, S., Mori, S., Pekar, J. J., Barker, P. B., Prince, J. L., and van Zijl, P. C., “Multi-parametric neuroimaging reproducibility: A 3-T resource study,” NeuroImage 54(4), 2854–2866 (2011).
-  Tustison, N. J., Avants, B. B., Cook, P. A., Zheng, Y., Egan, A., Yushkevich, P. A., and Gee, J. C., “N4ITK: Improved N3 bias correction,” IEEE Transactions on Medical Imaging 29(6), 1310–1320 (2010).
-  Avants, B. B., Tustison, N., and Song, G., “Advanced normalization tools (ANTS),” Insight j 2, 1–35 (2009).
-  Iglesias, J. E., Liu, C. Y., Thompson, P. M., and Tu, Z., “Robust brain extraction across datasets and comparison with publicly available methods,” IEEE Transactions on Medical Imaging 30(9), 1617–1634 (2011).