The deep convolutional GAN (DCGAN, Radford et al. (2016)) architecture was developed to be able to incorporate convolutions into a GAN model, since the original GAN of Goodfellow et al. (2014) used only fully-connected layers, which limited the size of the output images. The DCGAN3D generator used in this work is described in Appendix A, adapted from the publicly available 2D PyTorch implementation²²2https://github.com/pytorch/examples/blob/master/dcgan/main.py, accessed Feb 2020. The total number of parameters was $\sim 24$ million. In summary, the latent code is expanded to a $4\times 4\times 4$ volume via convolution, and then this volume is progressively upsampled with transposed convolution layers until reaching the output size of $32\times 64\times 64$. Each transposed convolution operation is followed by a 3D batchnorm (Ioffe and Szegedy, 2015) and a ReLU activation, except the last layer which has a tanh activation to bound the generated image in $[-1,1]$.

StyleGAN improves on other GAN techniques by using style-based image generation (Karras et al., 2019). The latent code is processed by a set of fully-connected layers and fed into a convolutional generator at multiple resolution scales to control the statistics of the activation maps using adaptive instance normalisation (Huang and Belongie, 2017). This is a more powerful method of altering overall image appearance than learning across multiple convolutional-only layers.

The styleGAN generator also incorporates a mapping network which maps $z$ to a another, warped space, which has distribution $\mathbb{P}_{w}$. The aim is to allow the network to learn a more natural latent space, such that the $w$ vectors are better structured according to semantic image features (Karras et al., 2019).

Figure 1 shows schematically the architecture of the 3D styleGAN (styleGAN3D) patch generator implemented in this work, and Appendix B describes the details of the latent mapping and synthesis networks. The total number of trainable parameters was $\sim 18$ million. Note that in this work the noise injection channels of the original styleGAN architecture were omitted since we would wish to retain control of all relevant details of the generated images for the potential downstream tasks of anomaly detection and image reconstruction.

bigGAN is a competing high-quality GAN model first described by Brock et al. (2019). By combining various architectural components that had previously been shown to aid GAN training, and systematically optimising the resulting architectures, bigGAN was demonstrated to allow stable training for high-resolution natural images. In particular, bigGAN inherits spectral normalisation Miyato et al. (2018) and self-attention layers Zhang et al. (2019) from previous GAN models. We adapted the publicly available bigGAN implementation³³3https://github.com/ajbrock/BigGAN-PyTorch, accessed Mar 2022 to 3D, with the full architecture listed in Appendix C. Note that due to memory constraints the bigGAN3D architecture had only $\sim 4$ million trainable parameters.

We used a simple convolutional feed-forward network as the discriminator backbone for all methods, as is common in both older and more recent GAN models (Radford et al., 2016; Karras et al., 2018, 2019). We also used minibatch discrimination (MD) to increase variability of GAN output. MD refers to a family of techniques which can be employed to avoid mode collapse (where many $x_{f}$ are similar) by allowing the discriminator to make comparisons between samples in a minibatch (Salimans et al., 2016). MD is particularly useful when the discriminator architecture is shallow, since shallow networks are intrinsically less able to detect mode collapse.

In this work we use a minimum-based MD method, denoted MDmin. Full details are given in Appendix D, but in summary the minimum perceptual difference between samples in a minibatch is provided to the discriminator to facilitate inclusion of whole-batch information. When mode-collapse has occurred for a generator model, this measure is small, allowing the discriminator to more easily detect that the minibatch is fake, providing a strong signal to the generator to diversify its samples.

The details of the discriminator architecture, including the location of the MDmin layer are given in Table 1. The number of trainable parameters was $\sim 11.0$ million. Note that the MDmin layer is placed just before the final convolutional layer to allow the discriminator to quickly detect the MD signal in cases of mode collapse.

We also note that the effectiveness of minibatch discrimination is limited by the maximum batch size that is computationally feasible. For this reason we use a heuristic trick to increase the effective batch size by pre-selecting the $k$ most similar (according to the MDmin) samples from a larger batch of $N$ samples, thereby allowing the discriminator to reliably see the most mode-collapsed samples for a particular model. We denote this training heuristic largeEBS (large effective batch size), and a full description is provided in Appendix D.

The original seminal work on GANs introduced by Goodfellow et al. (2014) proposed an intuitive approach where the discriminator aims to correctly classify real and fake images, while the generator aims to force the discriminator to misclassify fake images as real. Concretely, the losses for the two networks are:

and

where $\sigma(\cdot)$ denotes the sigmoid function. Note that expectation values $\mathbb{E}[\cdot]$ are defined over the entirety of the distributions, but in practice these losses are calculatated batch-wise as an approximation.

There are many alternative losses that have been proposed for training GANs, including Wasserstein distance (Arjovsky et al., 2017; Gulrajani et al., 2017), least-squares (Mao et al., 2017), Hinge loss (Lim and Ye, 2017), etc. Evaluating all such loss functions for a given application is a time-consuming exercise that we leave for future work. In this work we focus on one of the relativistic losses of Jolicoeur-Martineau (2019), which has been incorporated successfully into large-scale GAN models of medical images (Quiros et al., 2021). The equations for the relativistic losses are:

and

where:

Intuitively, these relativistic losses reframe the task of the discriminator and generator so that the discriminator now aims to assign higher output to real images than for fake ones. The generator then aims to ‘fool’ the discriminator to do the opposite. Note that in this case the generator learns directly from both the fake and real images, whereas for the standard GAN losses (Equations 1 and 2), the generator only receives direct feedback from the discriminator on the fake images.

We used the publicly available LUNA16 lung CT dataset (Setio et al., 2017), comprising 888 lung CT images, annotated by four radiologists. These annotations include malignancy scores and segmentations for nodules $\geq$3\text{\,}\mathrm{m}\mathrm{m}$$ diameter. We chose to exclude nodules and pathological images from the present study to approximate the training of a GAN for a downstream anomaly detection task.

First, we exclude nodules by selectively rejecting sampled patches that contain any ground truth nodule voxels. Second, we considered all scans that had any single nodule with a median malignancy score (across radiologists) $\geq 4$ as diseased, excluding the whole CT volume from further use. Removing these malignant scans left 636 scans which were split into 509 training scans and 127 held out for future testing.

Using various combinations of the above GAN components results in a number of overall comparison methods for training our GAN. Table 2 summarises the combinations we chose to compare in this work.

The training of GAN model involves sampling from the real image distribution, as well as generating fake images using the generator. As shown above, the GAN output was of size $32\times 64\times 64$ vox³. Therefore patches sampled from the real CT images were also of this size. Patches were sampled from within the lung using the LUNA16 provided lung masks, excluding all voxels annotated as nodules as described above. Figure 2 show a sample of real patches (centre slice only), highlighting the high variability in lung wall position and vessel structure. Note that resampling to a standard voxel size and other data augmentations (translation, rotation, etc) were not performed.

Unless stated otherwise, for a given CT image, 672 patches were sampled at a time, windowed to $[-1000,400]$ HU, and then rescaled to $[-1,1]$. These patches were then used to create minibatches of size 48, with 14 minibatches per image. For each update of the discriminator, this batch of 48 real images is passed through to obtain $D(x_{r})$, followed by a batch of 48 fake images to get $D(x_{f})$. A similar sampling process is followed for the generator updates. A single pass through all 509 training CT scans (i.e. 7126 iterations) is considered an epoch of training. For each epoch, patches are sampled separately, ensuring that with high probability no single patch is repeated during the training.

For all GAN models, the generator and discriminator learning rates were set the same at $0.0001$, using the Adam optimiser (Kingma and Ba, 2015) with $\beta_{1}=0.5$. The discriminator and generator were each updated once per iteration. Training was run for 20 epochs using either an NVIDIA GeForce GTX 1080 Ti or an NVIDIA TITAN Xp GPU.

Finally, for the styleGAN3D approach, so-called style-mixing regularisation was used (Karras et al., 2019) which swaps two latent codes $w_{1}$ and $w_{2}$ at a random depth in the convolutional generator. Since the GAN training process still encourages the output to be realistic, the effect is to decouple the effect of the latent code at different scales (see Karras et al. (2019) for further details). Note that in comparison to the standard 2D styleGAN method of Karras et al. (2019), progressive growing of images was not found to be necessary for styleGAN3D due to the small relative size of the output patches compared to the convolutional kernel sizes.

GAN models were evaluated throughout training with the commonly-used Fréchet inception distance (FID) which was proposed and explained in detail by Heusel et al. (2017). In brief, the FID aims to measure the similarity between the distributions of real training images ($\mathbb{P}$) and fake, GAN produced images ($\mathbb{Q}$). To do this, a large number ($\sim 5000$–$10000$) of images from each distribution is sampled. These images are passed through a pre-trained network to obtain the activations at a deep layer, which can be assumed to be Gaussian in distribution. Measuring the Fréchet distance between these two distributions provides a measure of similarity between the original distributions. The FID provides a measure of both average image quality (i.e. ‘realism’) and variability, and so it is sensitive to mode collapse.

As mentioned, FID uses a pre-trained network to provide higher-level representations of the input images. In common practice, this network is an Inception network pre-trained on a natural image task (Heusel et al., 2017), which presents two considerations for use in the current investigation. Firstly, the pre-trained network accepts only 2D input. Therefore, when calculating FID for our 3D GANs, we use only the central slices of generated images. Secondly, since the data domain is different (natural images vs lung CT patches), the transferability of the standard Inception network to lung CT could be questionable. However, we note that calculating the FID between two random sets of real images produces a low FID value of $\leq 4$, suggesting that producing a low FID score is a valid target for GAN models of lung CT patches.

We used FID to measure GAN performance throughout training, with five training runs performed for each method listed in Table 2. For each method we selected the lowest-FID model from across the five training runs as the best. These lowest FIDs were then compared between architectures, and the corresponding models were retained for further investigation and characterisation.

To complement the results of the automated model evaluation, we also performed a human observer study on the best performing 3D GAN models. One human observer (SE) was provided 200 2D image slices per model, split 50:50 into real and fake. The task was to classify real and fake images, and this was repeated three times to assess intra-observer variability. Receiver operating characteristic (ROC) analysis of the results was then performed.

The analyses and model selection described above were performed only in 2D due to the reliance of the standard FID on a 2D Inception network, and the infeasibility of performing large observer studies with 3D patches. To investigate the 3D, domain-specific structure of the generated patches, additional analysis was performed.

Firstly, although the conventional FID measurement uses a trained 2D Inception network, there has been research proposing the use of a 3D network to enable a 3D FID score (Sun et al., 2020). In this case the network from which activations are retrieved is a 3D ResNet pre-trained on different medical imaging segmentation datasets (Chen et al., 2019). Transfer learning experiments showed that the learnt features in such a network were beneficial for a range of tasks, including segmentation of lung tissue in CT images (Chen et al., 2019). We use one model provided by Chen et al. (2019) to calculate a 3D FID (FID3D) for our lung tissue GANs. Specifically, we use a pre-trained ResNet-10 model, and obtain the output activation maps, which have dimensions of $512\times 4\times 8\times 8$ when given $1\times 32\times 64\times 64$ input. To avoid excessive memory usage we spatially average these maps to obtain a $512$ element feature per input patch. Processing $\sim 10000$ real and fake images allows the Fréchet distance to be measured between real and fake datasets as for standard 2D FID.

Secondly, we note that if the distributions of real and fake images are close, then we expect that the distributions of any quantity derived from those images should also be close, i.e. comparing the distributions of any characteristic of real and fake images is valid. Beyond this, using a quantity relevant for the task at hand is preferable. To the best of our knowledge, little research has been performed on the expected appearance of healthy lung tissue in CT image patches. Therefore we opted to analyse the appearance of a prominent feature of healthy lung tissue patches: the vasculature. Specifically we chose to model the complexity of the vasculature by finding the number of branch points in an automated fashion. 3D patches were skeletonised using standard python scikit-image library functions in order to provide a sparse representation of the 3D structure. These 3D skeletons were visually assessed using 3D rotating maximum intensity projections (MIPs). The number of branch points was automatically detected, and the distribution of branch points across $\sim$10k images was recorded. Figure 3 shows an example skeletonisation of a generated patch, with automatically detected branch-points highlighted.

Finally, the relationship between the latent space and the distribution of branch points was investigated using the UMAP dimensionality reduction technique (McInnes et al., 2020). Using the freely available implementation of UMAP⁴⁴4https://umap-learn.readthedocs.io/en/latest/index.html, accessed Oct 2021, the 512-dimensional latent spaces for each method ($z$ for DCGAN3D and bigGAN3D, $w$ for styleGAN3D) were transformed into a two-dimensional embedded space. The embedded space was calculated using 50000 random samples to mediate the curse of dimensionality. For 1000 of these samples, the original latent vectors were passed through the corresponding GAN model to produce fake images which could then be analysed as above to find the number of branch points. Latent space structure could then be visualised by labelling each point in the embedded latent space with its corresponding number of branch points.

Table 3 shows the minimum observed FIDs for each model, as well as the average and standard deviation (SD) of the minimum over multiple runs. The SD of the minimum FID over the five runs provides a measure of training stability, since it reflects the ability of the method to provide similar results under different initialisations and sequences of training images. The significance of the effects of the MDmin and larger effective batch size techniques was established by performing Welch’s $t$-test on the minimum FID in the five training runs for pairs of methods using the same generator architecture (Figure 4). FID values and network losses throughout the training are provided in Appendix E.

For the DCGAN3D method, including the MDmin layer provides a small, non-significant ($p=0.64$) improvement, reducing the minimum FID from $89.0$ to $80.3$ and the average minimum FID only from $93.1\pm 3.1$ to $91.0\pm 9.0$. The increased SD over five runs suggests that the modest performance improvement comes at the cost of reduced training stability. In addition, we found that the performance of DCGAN3D-base was not significantly affected by batch size (Appendix E).

On the other hand, MDmin was observed to have a significant effect for the styleGAN3D architecture. Without the MDmin layer, the minimum FID observed was $122.3$, with an average minimum of $136.9\pm 10.7$. Including the MDmin layer reduced this to a minimum of $41.0$ with an average of $44.8\pm 4.2$ ($p<<0.05$); a reduction in both mean and SD representing improved training stability. For the bigGAN3D architecture, the use of MDmin reduced FID from $93.8\pm 5.4$ to $58.4\pm 5.8$ which was a significant improvement ($p<<0.05$).

The heuristic method to increase the effective batch size (largeEBS) provided non-significant improvements for the DCGAN3D and bigGAN3D architectures, although for the styleGAN3D architecture, the minimum FID reduced from $41.0$ with just the MDmin layer to $35.3$ when using the increased effective batch size, reflecting a significant reduction in mean minimum FID ($p<0.05$).

Figure 5 shows some samples for the minimum FID models for each of the GAN methods. The DCGAN3D-base method shows a variety of images, but the image quality is poor, with periodic artefacts visible. On the other hand, the styleGAN3D-base and bigGAN3D-base methods produce high-quality images, but with severe mode collapse. The DCGAN3D-MDmin and DCGAN3D-MDmin-largeEBS methods improve on image quality compared to DCGAN3D-base, but some mode collapse has arisen. This seems to contradict the purpose of the MDmin and largeEBS techniques, though it should be noted that the minimum FID for the best DCGAN3D-base training run occurred early in training, whereas the DCGAN3D-MDmin method reached a minimum at a later point. Because of this, a direct comparison of the variability of the two methods is not appropriate, since mode collapse usually occurs later in training.

In contrast to DCGAN3D, the styleGAN3D-MDmin model reduces mode collapse while retaining image quality compared to the styleGAN-base method. For bigGAN3D, the MDmin did not visibly reduce mode collapse to the same extent as for styleGAN, despite the significant drop in FID recorded in Table 3. However, image quality seems to be qualitatively improved, with more structure visible within the lung parenchyma. Similarly to the DCGAN3D methods, this is likely due to the early occurrence of the minimum FID for the bigGAN3D-base method, before image quality was stable. Using the largeEBS trick seems to further reduce mode collapse for the styleGAN architecture, whereas for the bigGAN3D architecture the effect was negligible.

Finally, Figure 6 shows interpolations between two randomly sampled latent codes for each of the minimum FID models. Note that DCGAN3D and bigGAN3D methods were spherically interpolated in $z$-space, whereas styleGAN methods were linearly interpolated in $w$-space, as described in the original description of the method (Karras et al., 2019). A smooth latent space is a desirable quality of GAN models, which should be reflected by smooth transitions between sampled images. The DCGAN3D methods provide smooth transitions despite their poor image quality. The bigGAN3D methods also exhibit smoothness, with a superior image quality. In general the styleGAN3D interpolations in $w$-space are smooth, although some discontinuities/non-smooth regions are present, possibly due to interpolations traversing low-density regions of $w$-space.

Figure 7 shows the results of the observer study. The DCGAN3D methods, suffering from periodic artefacts and poor image quality, were easily distinguishable from the real images, with $\textrm{FPR}<0.025$ and $\textrm{TPR}>0.88$ for all models. Slightly improved performance was observed using the bigGAN3D models, with $\textrm{FPR}=0.04$ and $\textrm{TPR}=0.93$ for the bigGAN3D-MDmin method, but some remaining convolution artefacts made the fake images conspicious compared to real samples. The styleGAN3D-base method performed better than all the DCGAN3D and bigGAN methods, but the mean observer performance was still high at $\textrm{FPR}=0.14$ and $\textrm{TPR}=0.89$. When including the MDmin layer, observer performance dropped to $\textrm{FPR}=0.29$ and $\textrm{TPR}=0.66$, representing improved model performance. The large effective batch size heuristic worsened styleGAN3D performance slightly, yielding $\textrm{FPR}=0.22$ $\textrm{TPR}=0.69$.

The FID3D score for each of the models is shown in Table 4. Similarly to the conventional 2D FID scores in Table 3, the best performing method overall is the styleGAN3D-MDmin-largeEBS model which produced a value of 0.57 for the FID3D metric. Interestingly, the effect of the largeEBS method is different in 2D FID compared to 3D FID. For example, the largeEBS reduces 2D FID for the bigGAN3D, compared to using MDmin alone, whereas in terms of FID3D, it rises considerably. We hypothesise that this is because the bigGAN3D architecture has fewer trainable parameters, which means that the MDmin method alone may be fully utilising network capacity. Pushing further with the largeEBS method may cause the training to become more unstable, resulting in worse 3D structure. We leave a full comparison of 2D and 3D FID metrics for future work.

Figure 8 shows some example skeletons calculated from random outputs of each GAN model. The DCGAN3D models provide noisy skeleton images, due to the poor quality of the generated samples. In contrast the styleGAN3D and bigGAN3D outputs produce cleaner skeletons with clearer structure, more closely resembling real images.

To analyse the distributions of 3D branch points in real and fake images, ROC curves were produced for each method by varying a threshold number of branch points and recording the true and false positive rates (TPR and FPR) of the resulting binary classification. Methods that generate fake images with a similar distribution of branch points to the real data will produce areas under the ROC curve (AUCs) closer to 0.5, and vice versa. Figure 9 shows the ROC curves each of the methods. The DCGAN methods provide AUCs of $0.750\pm 0.003$ (DCGAN-base), $0.693\pm 0.004$ (DCGAN3D-MDmin), and $0.784\pm 0.003$ (DCGAN3D-MDmin-largeEBS), demonstrating that the DCGAN architecture is unable to produce a realistic distribution of branch points. In contrast, the styleGAN3D architecture produces distributions of branch points closer to the real one, reflected by AUC values of $0.535\pm 0.004$ for styleGAN3D-base, $0.671\pm 0.004$ for styleGAN3D-MDmin, and $0.608\pm 0.004$ for styleGAN3D-MDmin-largeEBS. Overall, the bigGAN3D methods provided the best performance, with a lowest AUC of $0.524\pm 0.004$.

The UMAP visualisation of the latent spaces, and their relationship to the number of branch points are shown in Figure 10. For the DCGAN3D and bigGAN3D methods, the UMAP-embedded space shows no structure since the $z$-space is by design distributed normally. Furthermore, no correspondence is apparent between the position of a point in the embedded space and the number of branch points in the corresponding generated image. Conversely, the styleGAN3D methods show both structure in the shape of the $w$-space, and a correspondence between position and number of branch points. This suggests that higher level concepts are contributing to the structure of the full latent space, agreeing with previous literature on the use of styleGAN for medical image generation (Quiros et al., 2021). This is an interesting result, opening the possibility for further investigation into the structure of the latent space with respect to other semantically meaningful features in future work.

It is also important to emphasise that the lack of correspondence between UMAP embeddings and the number of branch points for the bigGAN3D and DCGAN3D methods does not suggest that the full $512$-dimensional latent spaces for these methods are not ordered with respect to the number of branch points, just that any such correspondence is not preserved under the UMAP transformation. Additional analysis to explore whether such correspondence exists for the DCGAN3D and bigGAN3D methods remains for future work.

The following approach to MD was adopted, calculating the L₁ distance, $d_{ij}$, between two images $x_{i}$ and $x_{j}$ within a minibatch as:

where $D_{l}\left(x\right)$ denotes the features of the discriminator network after layer $l$ with $x$ as the input image. In the traditional minibatch discrimination approach (Salimans et al., 2016), these distances (or some transformation of them) are then either summed or averaged over the minibatch to get a per-image score, i.e. $\textrm{MD}_{i}\propto\sum_{j}^{N}d_{ij}$, which is used as input for the subsequent discriminator layer.

We note that taking the sum of $d_{ij}$ allows partial mode collapse to occur if one mode is sufficiently far from the rest of the samples. To overcome this, we propose to take instead the minimum within-batch distance for each image, denoted MDmin, as follows: $\textrm{MDmin}_{i}=\min_{j}{(d_{ij})}$. This represents the maximum similarity between two samples in the minibatch, providing a consistently clearer signal to the discriminator if these images are mode-collapsed. The scalar values per image, $\textrm{MDmin}_{i}$, are then expanded to the correct size and appended to the current feature maps, to provide an extra feature to the subsequent layers of the discriminator.

The MDmin layer described above is able to reduce mode collapse by allowing the discriminator to detect when a minibatch contains similar samples, and feed this information back to the generator. However, mode collapse can clearly still occur, as long as similar samples are unlikely to occur within a minibatch. Size the maximum batch size is bounded by computational limitations (GPU memory), this issue cannot be solved completely by simply increasing the batch size.

We propose to use a heuristic trick to increase the effective batch size, mitigating the effect described above. First, a large batch of $N$ latent codes is sampled and passed through the generator without tracking gradients. This allows a larger batch size than would be normally used in training, since much of the computational load of network training is in the storing of gradients. These generated images are then passed through the discriminator without tracking gradients, and the activation maps at depth $l$ are calculated, similarly as for the MDmin layer. The similarity between all pairs of samples is then calculated, as for the MDmin layer, and the minimum for each sample is taken as a measure of how mode-collapsed that particular sample is within the larger minibatch. The $k$ most similar images are noted, and the corresponding samples are then used for the subsequent training iteration, ensuring that the networks see the most mode-collapsed samples for training purposes. Figure A1 shows this process schematically. Note that with this method, the upper limit to the effective batch size is no longer dependent on memory, since each pair of images can in theory be loaded to memory as required and do not need to be stored otherwise. Real samples are still selected entirely at random; choosing the $k$ most similar samples from a larger batch of real images was seen to degrade performance. Since the learned discriminator representations will not be relevant early in training, we train the networks for $5$ epochs as normal before increasing the effective batch size thereafter. Similar heuristic tricks to select the most relevant/useful samples for GAN training are established in the literature, for example in Sinha et al. (2020), although not for the specific case of reducing mode collapse.