MRI Scan Synthesis Methods Based on Clustering and Pix2Pix

We consider a missing data problem in the context of automatic segmentation methods for Magnetic Resonance Imaging (MRI) brain scans. Usually, automated MRI scan segmentation is based on multiple scans (e.g., T1-weighted, T2-weighted, T1CE, FLAIR). However, quite often a scan is blurry, missing or otherwise unusable. We investigate the question whether a missing scan can be synthesized. We exemplify that this is in principle possible by synthesizing a T2-weighted scan from a given T1-weighted scan.

Our first aim is to compute a picture that resembles the missing scan closely, measured by average mean squared error (MSE). We develop/use several methods for this, including a random baseline approach, a clustering based method and pixel-to-pixel translation method by Isola et al. [15] (Pix2Pix) which is based on conditional GANs. The lowest MSE is achieved by our clustering-based method.

Our second aim is to compare the methods with respect to the effect that using the synthesized scan has on the segmentation process. For this, we use a DeepMedic model trained with the four input scan modalities named above. We replace the T2-weighted scan by the synthesized picture and evaluate the segmentations with respect to the tumor identification, using Dice scores as numerical evaluation. The evaluation shows that the segmentation works well with synthesized scans (in particular, with Pix2Pix methods) in many cases.

1. 1.

   Note, however, that our work was mostly done before that challenge was posed, so in particular, we are using the BraTS 2019 data set for our experiments.

This work was partly funded by the German Research Foundation (DFG), project numbers 416767905 and 456558332.

Authors and Affiliations

1. Department of Computer Science, Heinrich-Heine-University Düsseldorf, Universitätsstraße 1, 40225, Düsseldorf, Germany

   Giulia Baldini & Melanie Schmidt

2. Department of Diagnostic and Interventional Radiology, University Hospital Aachen, Pauwelsstraße 30, 52074, Aachen, Germany

   Charlotte Zäske

3. Department of Radiology, University Hospital of Cologne, Kerpener Street 62, 50937, Cologne, Germany

   Liliana L. Caldeira

Corresponding authors

Correspondence to Giulia Baldini or Melanie Schmidt .

Editors and Affiliations

1. University of Utah, Salt Lake City, UT, USA

   Joseph Finkelstein

2. Ben-Gurion University of the Negev, Beer Sheva, Israel

   Robert Moskovitch

3. University of Pavia, Pavia, Italy

   Enea Parimbelli

Ethics Approval

Ethical approval was not required for this study, as it exclusively utilizes data made available during the BraTS segmentation challenge.

Appendix

A Detailed Description of BrainClustering

Figure 2 gives an overview of the BrainClustering process, which can be divided into multiple steps.

Training Step 1: Both images are segmented into macro clusters. The segmentation is done by solving a one-dimensional k-means problem optimally by dynamic programming with the implementation \( k \)-means1d  [11, 34, 37]. In the example, the number of macro clusters k is 3. Since the macro clusters are supposed to correspond to the different tissues, they should intuitively be equal to the number of different tissues. However, the process works better if we allow a little more macro clusters to account for varieties of tissues. Empirically, a good range is between three and six clusters.

Training Step 1b: As an intermediate step after finding the macro clusters, we need to identify which cluster represents what tissue. The tissues in T1W can be identified by ordering the clusters according to the average intensity value of their pixels. We could do a similar approach to finding the tissues in T2W, but we instead find them by matching them to the found clusters in T1W as described below in Step 2b (cluster label matching).

Training Step 2: Next, we compute a micro clustering for each tissue which captures the shadings of the tissues and the relation between the shading in T1W and T2W. We again implement this step by optimally solving a one-dimensional k-means problem. The number of micro clusters corresponds to the number of different shades that we allow. In Fig. 2, we use 3 micro clusters per tissue for visualization purposes, while in practice, we use at least 100 micro clusters.

Training Step 2b: Now, there are micro clusters for every tissue, both in T1W and in T2W. They are small patches of the same shading, and we aim to identify how such a patch in T1W is mapped to T2W. For this, we need to match the clusters in T1W and T2W. We call this process cluster label matching. The matching is done by first excluding voxels that are only present in one of the scans (since the scans are registered to the same template, these are only a few voxels), and then finding the label that maximizes the number of voxels that the matched clusters have in common. The problem reduces to maximum weighted matching problem, which we solve with the Hungarian method [16].

After this step, we have a macro clustering and every macro clustering is subdivided into small patches, and every small patch in T1W has its specified counter part in T2W.

Training Step 3: The next step is to capture the relationship between the shadings in T1W and T2W. The idea behind this is that we assume that it is possible to map the shading of a specific tissue in T1W to the shading of the same tissue in T2W. We model this by using a function \( f_{t}\) for every tissue which is supposed to translate intensity values from T1W to intensity values in T2W. We later want to use the function \( f_{{t}}({i}_1) = {i}_2 \) to predict the T2W intensities of a patient whose T2W scan is missing.

For every tissue type t, we have computed a micro clustering and matched the micro clusters between the two scans. Now we compute the average of the points in each micro cluster of T1W, which we call \( a_1 \). Then, we compute the average for each corresponding micro cluster of T2W, which we call \( a_2 \). Finally, we add a new entry \( (a_1, a_2) \) to \( f_{{t}} \), which is the map corresponding to the current macro cluster/tissue type. It may happen that there already is an entry with key \(a_1\) stored in the map. In this case, we compute the average with the moving average formula:

where \(\#_{t}(a_1)\) is the number of times we tried to add an entry with key \(a_1\) to the map \(f_{{t}}\). Note that consequently we have to keep track of \(\#_{t}(a_1)\), i.e., we have to annotate each entry of the map with the corresponding cardinality.

Training Step 1-3: are done for all T1W/T2W scan pairs available in the training data set, and the values found in Step 3 are always inserted into the same tables. So we get as many tables as we have macro clusters, and up to as many rows as we have micro clusters for all training images. Since we process many scans, the full tables get too large, and in reality, we only store a meaningful subset of the rows. When querying for a T2W, we can compute missing values by interpolation.

What is left to describe is how we now synthesize images based on our mapping tables. We implemented two options (the second one is faster).

Synthesizing T2W : For a patient whose T2W is missing, we preprocess the T1W scan and then cluster it with \( k \)-means1d to obtain a macro clustering. We load all tables \(f_{{t}}\) into the memory which were computed with the same number of macro clusters. Let \( {t}\) be one of the macro clusters. We consider each point p belonging to cluster \( {t}\). We find the two rows of \(f_{{t}}\) with the closest intensity values to p and interpolate \(f_{{t}}(p)\) from these two rows. The resulting intensities form a synthesized scan. Once it is computed, it is postprocessed using a \( 3\times 3 \) median filter [9] to remove salt-and-pepper noise which was present in the synthesized images.

Synthesizing T2W in Search Mode: We noticed that having large tables computed by training on many images (e.g. 200) produces noise in the results. Thus, we propose an additional method that does not precompute a large model before querying, but instead computes a small model every time we want to answer a query, which we call Search. Given a T1W query image, we first search in the training data set for the w patients whose T1W has the smallest mean squared error to the query T1W (where w is a small constant, e.g. \( w = 5 \)). Then, we create a small model by performing the training process only on these w patients, and produce the synthesized T2W with this model. This approach has the downside that we have to create a new model for every single query, thus significantly increasing the query time. However, since we choose small values for w, a synthesized T2W image for one input can be computed in around ten minutes, without the need of computing a model beforehand.

BrainClustering (BC), Pix2Pix (P2P) and original T2W scans and corresponding superimposed segmentation generated with DeepMedic (together with the real T1W, T1CE and FLAIR), shown for all cases from the testOUR data set. In the last column, the ground truth segmentation is superimposed on the T1W scan.

Full size image

Continuation of Fig. 3.

Full size image

The segmentations in the second and in the fourth row have been created using the respective T2W and the three original T1W, T1CE, FLAIR. The rightmost image in the “Real T1W” column shows the real T1W image and the ground truth segmentation.

Full size image

[Best viewed in color] Comparison of different approaches with respect to MSE in the brain mask and (if available) in the tumor areas. The scores represent the distribution of the patients. A lower value indicates a better score. The crosses represent the averages, the thick bars are the medians and the dots are the outliers.

Full size image

[Best viewed in color] Comparison with respect to the true tumor segmentation (for testBraTS and testOUR) and the original tumor segmentation (for testNoTruth) using Dice and undirected \( 95^{\text {th}} \) Hausdorff distance (HD95). For testBraTS and testOUR, we also include the segmentation which is produced by using the real T2W as Original to give a positive baseline as to what scores are achievable if T2W was reproduced perfectly. We compute the Dice score and the undirected 95% percent Hausdorff distance on three different regions and average the scores. The box plots represent the distribution of the patients. For the Dice score a higher value indicates a better value whereas for the Hausdorff distance the opposite holds. The crosses represent the averages, whereas the thick bars are the means and the dots are the outliers.

Full size image

We compare the classic Train &Test mode and Search mode with \( w \in \{5, 10\} \) and using 3, 4, 5 and 6 macro clusters, for a total of 12 different models.

B Additional Modifications to Pix2Pix

We implement some modifications to Pix2Pix:

1. 1.

   We modify the loading functionality to accept three-dimensional NIfTI images.

2. 2.

   We extend the data augmentation process with an additional library, TorchIO [28], which implements useful preprocessing and data augmentation routines for medical imaging.

3. 3.

   We add the mixed precision training functionality from NVIDIA APEX (available at https://github.com/NVIDIA/apex). Operations like matrix to matrix multiplication and convolution are then performed in half-precision floating-point format, which results in a speed-up at training time.

4. 4.

   Since a big portion of our scans is background, the L1 loss is only computed on the actual brain voxels.

5. 5.

   Recent studies have shown that the transpose convolution operation during upsampling creates checkerboard artifacts [25, 36]. We implement the solution of Wojna et al. [36]: linear additive upsampling. This method has been successfully used for the generator of Pix2PixHD  [35] and applied on CT scans [12]. We substitute the transpose convolution operation with a two-factor upsampling followed by a 4-factor reduction of the number of channels, and finally by a \( 3\times 3 \) convolution with stride 1. We use bilinear upsampling if the scans are two-dimensional, and trilinear upsampling for the three-dimensional case.

C Additional Figures, Diagrams and Tables

(See Figs. 3, 4, 5, 6, 7 and Table 2).

Reprints and permissions

© 2024 The Author(s), under exclusive license to Springer Nature Switzerland AG

Policies and ethics