2slbrpshort = 2S-LBRP, long = two-step Laplace-Beltrami regularized projection
arapshort = ARAP, long = as-rigid-as-possible bdtshort = BDT, long = bagged
decision tree bfmshort = BFM, long = Basel face model cnnshort = CNN, long =
convolutional neural network cpdshort = CPD, long = coherent point drift
cishort = CI, long = cephalic index crshort = CR, long = cephalic ratio
ctshort = CT, long = computed tomography cvaishort = CVAI, long = cranial
vault asymmetry index dlshort = DL, long = deep learning dtshort = DT, long =
decision tree ecgshort = ECG, long = electrocardiography ecgishort = ECGI,
long = electrocardiographic imaging fnnshort = FNN, long = feedforward neural
network ganshort = GAN, long = generative adversarial network gpashort = GPA,
long = generalized Procrustes analysis gpmmshort = GPMM, long = Gaussian
process morphable model icpshort = ICP, long = iterative closest points
icpdlbrpshort = ICPD-LBRP, long = iterative coherent point drift with Laplace-
Beltrami regularized projection icpdshort = ICPD, long = iterative coherent
point drift gpushort = GPU, long = graphics processing unit knnshort = kNN,
long = k-nearest-neighbors lbrpshort = LBRP, long = Laplace-Beltrami
regularized projection lbshort = LB, long = Laplace-Beltrami ldashort = LDA,
long = linear discriminant analysis lmshort = LM, long = landmark mapshort =
MAP, long = maximum a posteriori mrishort = MRI, long = magnetic resonance
imaging mlshort = ML, long = machine learning mlpshort = MLP, long = multi
layer perceptron nbshort = NB, long = naïve Bayes nicpashort = N-ICP-A, long =
nonrigid iterative closest points affine nicptshort = N-ICP-T, long = nonrigid
iterative closest point translation nnshort = NN, long = neural network
obbshort = OBB, long = oriented bounding boxes osnicpshort = OS-N-ICP, long =
optimal step nonrigid iterative closest points pcashort = PCA, long =
principal component analysis pdmshort = PDM, long = point distribution model
ppcashort = PPCA, long = probabilistic principal component analysis psmshort =
PSM, long = posterior shape model qdashort = QDA, long = quadratic
discriminant analysis ransacshort = RANSAC, long = random sample consensus
rfshort = RF, long = random forest rmsshort = RMS, long = root mean squared
shapshort = SHAP, long = SHapley Additive exPlanations sidsshort = SIDS, long
= sudden infant death syndrome ssmshort = SSM, long = statistical shape model
ssimshort = SSIM, long = structural similarity index measure ssimccshort =
SSIMcc, long = structural similarity index measure to closest clinical sample
stoshort = STO, long = sellion tragion orientation svdshort = SVD, long =
singular value decomposition svmshort = SVM, long = support vector machine
wpcashort = WPCA, long = weighted principal component analysis

# Impact of Data Synthesis Strategies for the Classification of
Craniosynostosis

Matthias Schaufelberger  Institute of Biomedical Engineering (IBT), Karlsruhe
Institute of Technology (KIT), Kaiserstr. 12, 76137 Karlsruhe, Germany
Reinald Peter Kühle  Department of Oral, Dental and Maxillofacial Diseases,
Heidelberg University Hospital, Im Neuenheimer Feld 400, Heidelberg, Germany
Andreas Wachter  Institute of Biomedical Engineering (IBT), Karlsruhe
Institute of Technology (KIT), Kaiserstr. 12, 76137 Karlsruhe, Germany
Frederic Weichel  Department of Oral, Dental and Maxillofacial Diseases,
Heidelberg University Hospital, Im Neuenheimer Feld 400, Heidelberg, Germany
Niclas Hagen  Institute of Medical Informatics, Heidelberg University
Hospital, Im Neuenheimer Feld 130.3, Heidelberg, Germany  Friedemann Ringwald
Institute of Medical Informatics, Heidelberg University Hospital, Im
Neuenheimer Feld 130.3, Heidelberg, Germany  Urs Eisenmann  Institute of
Medical Informatics, Heidelberg University Hospital, Im Neuenheimer Feld
130.3, Heidelberg, Germany  Jürgen Hoffmann  Department of Oral, Dental and
Maxillofacial Diseases, Heidelberg University Hospital, Im Neuenheimer Feld
400, Heidelberg, Germany  Michael Engel  Department of Oral, Dental and
Maxillofacial Diseases, Heidelberg University Hospital, Im Neuenheimer Feld
400, Heidelberg, Germany  Christian Freudlsperger  Department of Oral, Dental
and Maxillofacial Diseases, Heidelberg University Hospital, Im Neuenheimer
Feld 400, Heidelberg, Germany  Werner Nahm  Institute of Biomedical
Engineering (IBT), Karlsruhe Institute of Technology (KIT), Kaiserstr. 12,
76137 Karlsruhe, Germany

(April 2023)

## Abstract

Introduction: Photogrammetric surface scans provide a radiation-free option to
assess and classify craniosynostosis. Due to the low prevalence of
craniosynostosis and high patient restrictions, clinical data is rare.
Synthetic data could support or even replace clinical data for the
classification of craniosynostosis, but this has never been studied
systematically.

Methods: We test the combinations of three different synthetic data sources: a
ssm, a gan, and image-based pca for a cnn-based classification of
craniosynostosis. The cnn is trained only on synthetic data, but validated and
tested on clinical data.

Results: The combination of a ssm and a gan achieved an accuracy of more than
0.96 and a F1-score of more than 0.95 on the unseen test set. The difference
to training on clinical data was smaller than 0.01. Including a second image
modality improved classification performance for all data sources.

Conclusion: Without a single clinical training sample, a cnn was able to
classify head deformities as accurate as if it was trained on clinical data.
Using multiple data sources was key for a good classification based on
synthetic data alone. Synthetic data might play an important future role in
the assessment of craniosynostosis.

##  1 Introduction

Craniosynostosis is a group of head deformities affecting infants involving
the irregular closure of one or multiple head sutures and its prevalence is
estimated to be between four and ten cases per 10,000 live births [1]. As
described by Virchow’s law [2], depending on the affected suture distinct
types of head deformities arise. Genetic mutations have been identified as one
of the main causes of craniosynostosis [3, 4], which has been linked to
increased intracranial pressure [5] and decreased brain development [6]. The
most-performed therapy is surgical intervention consisting of resection of the
suture and cranial remodeling of the skull. It has a high success rate [7] and
is usually performed within the first two years of age. Early diagnosis is
crucial and often involves palpation, cephalometric measurements, and medical
imaging. ct imaging is the gold standard for diagnosis, but makes use of
harmful ionizing radiation which should be avoided, especially for very young
infants. Black-bone mri [8] is sometimes performed, but requires sedation of
the infants to impede moving artifacts. 3D photogrammetric scanning enables
the creation of 3D surface models of the child’s head and face and is a
radiation-free, cost-effective, and fast option to quantify head shape. It can
be employed in a pediatrician’s office and has potential to be used with
smartphone-based scanning approaches [9].

Due to the low prevalence, craniosynostosis is included in the list of rare
diseases by the American National Organization for Rare Disorders. Beside the
few data, strict patient data regulations, and difficulties in anonymization
(photogrammetric recordings show head and face), there are no publicly
available clinical datasets of craniosynostosis patients available online.
Synthetic data could potentially be used as a substitute to develop algorithms
and approaches for the assessment of craniosynostosis, but only one synthetic
dataset based on a ssm from our group [10] has been made publicly available so
far. Scarce training data and high class imbalance due to the different
prevalences of the different types of craniosynostosis [4] call for the usage
of synthetic data to support or even replace clinical datasets as the primary
resource for dl-based assessment and classification. The inclusion of
synthetic data could facilitate training due to the reduction of class
imbalance and increase the classifier’s robustness and performance.
Additionally, synthetic data may also be used as a cost-effective way to
acquire the required training material for classification models without
manually labeling and exporting a lot of clinical data. Using synthetic data
for classification studies in a supporting manner or as a full replacement for
clinical data has gained attraction in several fields of biomedical
engineering (e.g. [11, 12]), especially if clinical data is not abundant.
While classification approaches of craniosynostosis on ct data [13], 2D images
[14], and 3D photogrammetric surface scans [15, 16, 17] have been proposed,
the dataset sizes were below 500 samples (e.g. [17], [15], and [13]) and
contained a high class imbalance. The usage of synthetic data is a
straightforward way to increase training size and stratify class distribution.

However, although the need for synthetic data had been acknowledged [15],
synthetic data generation for the classification of head deformities has not
been systematically explored yet. With the scarce availability of clinical
data and multiple options of synthetic data generation available, we aim to
test the effectiveness of multiple data synthesis methods both individually
and as multi-modal approaches for the classification of craniosynostosis.
Using synthetic data as training material facilitates not only the development
of larger and more robust classification approaches, but also makes data
sharing easier and increases data availability. A popular approach for 3D data
synthesis is statistical shape modeling. It describes the approach to model 3D
geometric shape variations by means of statistical analysis. With the
application of head deformities, they have been employed to distinguish
clinical head parameters [18], to evaluate head shape variations [19], to
assess therapy outcome [20], and to classify craniosynostosis [16]. Although
their value in the clinical assessment of craniosynostosis has been shown, the
impact of ssm-based data augmentation for the classification of
craniosynostosis has not been evaluated yet. With the introduction of a
conversion of the 3D head geometry into a 2D image, image-based cnn-based
classification [17] can be applied on low-resolution images. gan [21] have
been suggested as a data augmentation tool [15] and have been able to increase
classification performance for small datasets [22].

The goal of this work is to employ a classifier based on synthetic data, using
three different types of data synthesis strategies: ssm, gan, and image-based
pca. The three modalities are systematically compared regarding their
capability in the classification of craniosynostosis when trained only on
synthetic data. We will demonstrate that the classification of
craniosynostosis is possible with a multi-modal synthetic dataset with a
similar performance to a classifier trained on clinical data. Additionally, we
propose a gan design tailored towards the creation of low-resolution images
for the classification of craniosynostosis. Both the gan, the different ssm,
and pca, were made publicly available along as all the 2D images from the
synthetic training, validation and test sets.

##  2 Methods

###  2.1 Dataset and Preprocessing

All data from this study was provided from the Department of Oral and
Maxillofacial Surgery of the Heidelberg University Hospital, in which patients
with craniosynostosis are routinely recorded for therapy planning and
documentation purposes. The recording device is a photogrammetric 3D scanner
(Canfield VECTRA-360-nine-pod system, Canfield Science, Fairfield, NJ, USA).
We used a standardized protocol which had been examined and approved by the
Ethics Committee Medical Faculty of the University of Heidelberg (Ethics
number S-237/2009). The study was carried out according to the Declaration of
Helsinki and written informed consent was obtained from parents.

Figure 1: Landmarks provided in the dataset, used for the alignment for
statistical shape modeling and the coordinate system creation of the distance
maps [17]. The three landmarks on the right exist for both left and right part
of the head. Figure 2: Pie chart of the class ratios in the clinical dataset
(control 56 %, coronal 5 %, metopic 14 %, sagittal 25 %). The legend in the
center shows the absolute number of samples in the dataset (in total 496
samples). Figure 3: The four classes of the dataset with their distinct head
shapes and their resulting distance maps representation. Top row: frontal
view, middle row: top view, bottom row: 2D distance maps.

Each data sample was available as a 3D triangular surface mesh. We selected
the 3D photogrammetric surface scans from all available years (2011–2021). If
multiple scans for the same patient were available, we selected only the last
preoperative scan to avoid duplicate samples of the same patients. All patient
scans had been annotated by medical staff with their diagnosis and 10
cephalometric landmarks. Fig. 1 shows the available landmarks on the dataset.
We retrieved patients with coronal suture fusion (brachycephaly and unilateral
anterior plagiocephaly), sagittal suture fusion (scaphocephaly), and metopic
suture fusion (trigonocephaly), as well as a control group with the dataset
distribution displayed in Fig. 2. Besides healthy subjects, the control group
also contained patients suffering from mild positional plagiocephaly without
suture fusion. Subjects with positional plagiocephaly in the control group
were treated with helmet therapy or laying repositioning. In contrast, all
patients suffering from craniosynostosis required surgical treatment and
underwent remodeling of the neurocranium. The four head shape resulting from
craniosynostosis are visualized in Fig. 3.

We used the open-source Python module pymeshlab [23] (version 2022.2) to
automatically remove some recording artifacts such as duplicated vertices and
isolated parts. We also closed holes resulting from incorrect scanning and
removed irregular edge lengths by using isotropic explicit re-meshing [24]
with a target edge length of 1 mm. In an earlier work [17], we defined a 2D
encoding of the 3D head shape (“distance maps”, displayed in Fig. 3, bottom
row) which was also included in the pre-processing pipeline with the default
parameter of [17].

Figure 4: Data subdivision for the synthetic-data-based classification and the
creation of synthetic data. The test set was separated initially from the
dataset, while the validation set was used to produce the synthetic samples on
which the cnn was trained. Green: data, blue: 3D-2D image conversion, dark
red: generative models.

###  2.2 Data subdivision

We did not use the full clinical dataset (validation and test set according to
Fig. 4) as training data for the data generation models (gan, ssm, and pca)
since the statistical information of the test set would be included in the
synthetic data sources, leading to leakage (an overestimation of the model
performance due to statistical information “leaking” into the test set).
Instead, we chose the schematic displayed in Fig. 4. We used a stratified
50–50 split of the clinical data and used one half of the samples as the
validation set and the other half as the test set.

The test set was separated from the validation set, only to be used for the
final evaluation of the classifier. Following this approach, the test set did
neither have any influence on the synthetic data, nor was it incorporated in
the validation set and should therefore be a true representation of unknown
data to the classifier. The validation set was used to select the best network
during training and for hyperparameter tuning, but not as training material.
Additionally it was used as the original (training) data on which we built the
synthetic image generators. The synthetic training set was then created from
the validation set according to the three data synthesis approaches described
below: ssm, gan, and pca. The three approaches operated on different domains:
While the ssm was applied directly on the 3D surface scans, the gan and the
pca used the 2D distance map images. All images were created as
28$\times$28-sized craniosynostosis distance maps which was sufficient for
good classification in an earlier study [17]. We describe each of the three
individual approaches ssm, gan, and pca below.

###  2.3 Data Synthesis

####  2.3.1 ssm

The pipeline for the ssm creation (similar to [25]) consists of initial
alignment, dense correspondence establishment, and statistical modeling to
extract the mean shape and the principal components from the sample covariance
matrix (see also Fig. 5). For correspondence establishment, we employed
template morphing.

Figure 5: The statistical shape model pipeline employed in this study. The
target scan is colored green with the deforming template in white.

We used the mean shape of our previously published ssm [10] as a template
which would be morphed onto each of the target scans. Procrustes analysis was
employed on the ten cephalometric landmarks to obtain a transformation
including translation, rotation, and isotropic scaling from the template to
each target according to the cephalometric landmarks on the face and ears. For
correspondence establishment, we employed the lbrp approach [26] to morph the
template onto each of the targets. We used two iterations: a high stiffness
fit (providing a now landmark-free transformation from template to the target,
improving the alignment also from the back of the head not covered with the
landmarks) and a low stiffness fit (allowing the template to deform very close
to the targets [27]). The deformed templates were then in dense
correspondence, sharing the same point IDs across all scans and were used for
further processing.

gpa was performed to remove both rotational and translational components on
all the morphed templates so that the mean shape could be determined and
removed. The remaining zero mean data matrix served as a basis for the
principal component analysis. To counterbalance higher point density in the
facial regions, we used weighted pca instead of ordinary pca for the
statistical modeling. The weights were assigned according to the surface area
that each point encapsulated and computed using the area of each triangle of
the surface model. We created one ssm for each class, ensuring that the models
were independent from each other and did not contain influences from the other
classes. We cut off the coefficient vectors after 95 % of the normalized
variance to remove noise and ensured only the most important components were
included in the ssm. The synthesis of the model instances could then be
performed as

${\mathbf{\mathrm{s}}}={\mathbf{\mathrm{\bar{s}}}}+{\mathbf{\mathrm{V}}}{\mathbf{\mathrm{\Lambda}}}^{\frac{1}{2}}{\mathbf{\mathrm{\alpha}}},$
(1)

with ${\mathbf{\mathrm{\bar{s}}}}$ denoting the mean shape,
${\mathbf{\mathrm{V}}}$ the principal components,
${\mathbf{\mathrm{\Lambda}}}$ the sample covariance matrix, and
${\mathbf{\mathrm{\alpha}}}$ the shape coefficient vector. We created 1000
random shapes of each class using a Gaussian distribution of the shape
coefficient vector and created craniosynostosis distance maps for each sample.

####  2.3.2 Image-based pca

We used ordinary pca as the last modality to generate 2D image data. While the
ssm also made use of pca in the 3D domain, image-based pca operated directly
on the 2D images. This was a computationally inexpensive and less
sophisticated alternative to both gan and ssm since neither extensive model
training and hyperparameter tuning, nor 3D morphing and correspondence
establishment was required. We employed ordinary pca for each of the four
classes separately and we again created 1000 samples for each class. Since ssm
is related to pca, the data synthesis could be performed as

${\mathbf{\mathrm{i}}}={\mathbf{\mathrm{\bar{i}}}}+{\mathbf{\mathrm{V}}}{\mathbf{\mathrm{\Lambda}}}^{\frac{1}{2}}{\mathbf{\mathrm{\alpha}}},$
(2)

with ${\mathbf{\mathrm{\bar{i}}}}$ denoting the mean image in vectorized
shape, ${\mathbf{\mathrm{V}}}$ again the principal components,
${\mathbf{\mathrm{\Lambda}}}$ the sample covariance matrix, and
${\mathbf{\mathrm{\alpha}}}$ the coefficient vector of the principal
components. We again drew 1000 random vectors from a Gaussian distribution and
transformed them back into 2D image-shape.

####  2.3.3 gan

The gan combines multiple suggestions from different gan designs and was
designed as a conditional [28] deep convolutional [29] Wasserstein [30] gan
with gradient penalty [31] (cDC-WGAN-GP). The design in terms of the
intermediate image sizes is visualized in Fig. 6. For the full design
including all layers, consult Appendix A.

Figure 6: Visualization of the intermediate image sizes from the used gan
model. Left: generator, right: critic (discriminator). The filter kernel sizes
are described in the Appendix A.

We opted for a design including a mixture between transposed, interpolation,
and normal convolutional filter kernels, which prevented checkerboard
artifacts and large patches. The combination of interpolation layers and
transposed convolutional layers lead to better images than each of the
approaches alone (see also in Appendix B Fig. 12) present in our previous
approach [32]. The conditioning of the gan was implemented as an embedding
vector controlling the image label that we wished to synthesize. We trained
the gan for 1000 epochs using the Wasserstein distance [30] which is
considered to stabilize training [33]. Instead of the originally proposed
weight clipping, we used a gradient penalty [31] of $\lambda=1$. We used 10
critic iterations before updating the generator and a learning rate of
$\alpha=3\cdot 10^{-5}$ for both networks. The loss $L$ can be described as
follows [31]:

$L=\mathbb{E}_{\tilde{x}\sim\mathbb{E}_{D}}D(\tilde{x}|y)-\mathbb{E}_{x\sim\mathbb{E}_{G}}D(x|y)+\lambda(||\nabla_{\hat{x}}D(\hat{x})||_{2}-1)^{2}$
(3)

with $\tilde{x}$ denoting the generator samples $G(z|y)$ and $\hat{x}=\epsilon
x+(1-\epsilon)\tilde{x}$ with $\epsilon$ denoting a uniformly distributed
random variable between 0 and 1 [31].

Figure 7: Image development of the gan generator during different stages of
training visualized as a 2$\times$2 grid.

###  2.4 Image assessment

We used ssimcc as the basis for a metric to assess the similarity of the
synthetic images to the clinical images and defined the
$\mathrm{SSIM}_{\mathrm{cc}}$ for each synthetic sample by using the minimum
ssimcc with respect to each clinical sample of the same class $N$:

$\mathrm{SSIM}_{\mathrm{cc},i}=\min\nolimits_{\forall n\in
N}\mathrm{SSIM}(p_{i,\mathrm{synthetic}},p_{n,\mathrm{clinical}})$ (4)

It has to be noted that the $\mathrm{SSIM}_{\mathrm{cc}}$ itself did not
assess the quality of the synthetic images, but was rather designed to
evaluate the similarity to the clinical images. With this approach, we tried
to quantify a “good” data generator: The data should not be very similar to
the original data (because then we could simply use the original data), but
also not too different (because then they might not be a true representation
of the underlying class anymore). “Good” images should not be “too close” to
1, but also not “too low”.

###  2.5 cnn Training

Resnet18 was used as a classifier since it showed the best performance on this
type of distance maps [17]. We used pytorch’s [34] publicly available, pre-
trained Resnet18 model and fine-tuned the weights during training. During
training, all images were reshaped to a size of 224$\times$224 to match the
input size of Resnet18. We performed a different run of cnn training on all
seven combinations of the synthetic data. The cnn was trained only on
synthetic data (except for the clinical scenario which was trained on clinical
data for comparison). During training, we evaluated the model on both the
(purely synthetic) training data and the (clinical) validation set (see also
Fig. 8). The best-performing network was chosen according to the maximum
F1-score on the validation set. The test set was never touched during training
and only evaluated in a final run after training.

Figure 8: Classification training using the synthetic data, the validation
data, and the test set. The cnn classifier using clinical data uses the
validation data as a training set. Green: data, blue: violet: classification
models. Figure 9: Number of training samples in each classification scenario.
The clinical scenario has less than 500 samples while all synthetic scenarios
have at 4000, 8000, or 12000 samples.

When multiple data sources were used, the models had a different number of
training samples (see Fig. 9) and all synthetically-trained models were
trained for 50 epochs. Convergence was achieved usually already during the
first ten epochs, indicating that there was sufficient training material for
each model. We used Adam optimizer, cross entropy loss, a batch size of 32
with a learning rate of $1\cdot 10^{-4}$, weight decay of $0.63$ after each 5
epochs. To evaluate the synthetically-trained models against a clinically
trained model, we additionally employed one cnn trained on clinical data,
trained with the same parameters except a higher learning rate of $1\cdot
10^{-3}$.

We used the following types of data augmentation during training: Adding
random pixel noise (with $\sigma=1/255$), adding a random intensity (with
$\sigma=5/255$) across all pixels, horizontal flipping, and shifting images
left or right (with $\sigma=12.44\,\mathrm{pixels}$). All those types of data
augmentation corresponded to real-world patient and scanning modifications:
Pixel noise corresponded to scanning and resolution errors, adding a constant
intensity was equal to a re-scaling of the patient’s head, horizontal flipping
corresponded to the patient as if they were mirrored in real life, and
shifting the image horizontally modeled an alignment error in which the
patient effectively turns their head $20^{\circ}$ left or right during
recording.

All the clinical 2D data, the gan, and the statistical models were made
publicly available111https://github.com/KIT-IBT/craniosource-gan-pca-ssm. We
included a script to create synthetic samples for all three image modalities
to allow users to create a large number of samples. The synthetic and clinical
samples of this study are available on Zenodo [35].

##  3 Results

###  3.1 Image evaluation

Figure 10: Images of all three data modalities and clinical samples. From top
to bottom the image modalities: ssm, gan, pca, clinical. From left to right
the four classes: Control, coronal, metopic, sagittal.

Fig. 10 shows image of each of the different data synthesis types compared
with the clinical images. From a qualitative, visual examination, the
synthetic images had similar color gradients, shapes, and intensities as the
clinical images. gan images appeared slightly noisier than the other images
and did not show the left and right ear visible in the other images.

Figure 11: Boxplots of $\mathrm{SSIM}_{\mathrm{cc}}$ (ssim to closest clinical
sample) of each class for each of the synthetic data generators.

From the quantitative comparison (see Fig. 11), ordinary pca images were
substantially and consistently more similar to the clinical images than the
other two modalities (differences of the medians larger than 0.02), while ssm
and gan images were less similar, with the ssm images being the most
dissimilar for the coronal class.

###  3.2 Classification results

All comparison presented here were carried out on the untouched test set.
According to the classification results for the synthetic training in Tab. 1,
the ssm was the best single source of synthetic data with an F1-score higher
than 0.85. All combinations of synthetic models showed F1-scores higher than
0.8. The classifier on the clinical data scored an accuracy above 0.96, but
was surpassed by the combination of gan and ssm. F1-score was highest for the
clinical classification (0.9533), but the combination of ssm and gan scored a
very close F1-score (0.9518). Including a second data source always improved
the F1-score compared to a model with a single data source (adding pca to gan
by 0.29, adding ssm to pca by 0.16, adding ssm to gan by 0.1).

Table 1: cnn-Classification comparison on the test set trained on different synthetic data sources. Boldface: best results among the data source. Synthetic data source | Accuracy | F1-score
---|---|---
GAN | 0.4274 | 0.4930
PCA | 0.7581 | 0.6997
SSM | 0.9153 | 0.8547
GAN-PCA | 0.8508 | 0.7823
GAN-SSM | 0.9677 | 0.9518
PCA-SSM | 0.9153 | 0.8595
GAN-PCA-SSM | 0.9597 | 0.9445
Clinical | 0.9637 | 0.9533

##  4 Discussion

Without being trained on a single clinical sample, the cnn trained from the
combination of the ssm and the gan was able to correctly classify 95 % of the
data. Classification performance on the synthetic data proved to be equal to
or even slightly better than training on the clinical data, at least for the
data generated using the ssm and the gan (and optionally also pca). This
suggests that certain combinations of synthetic data might be indeed
sufficient for a classification algorithm to distinguish between types of
craniosynostosis. Compared with classification results from other works, the
purely synthetic-data-based classification performs in a similar range and
sometimes even better than other approaches on clinical data [15, 17, 16, 36,
13].

The ssm appeared to be the data source contributing the most to the
improvement of the classifier: Not only did it score highest among the unique
data sources, but it was also present in the highest scoring classification
approaches. One reason for this might be that it was also the least similar
data source for most of the classes. Due to the inherent modeling of the
geometric shape in 3D, the created 2D distance maps are always created from 3D
samples, while pca and the gan could, in theory, create 2D images which do not
correspond to a 3D shape. In contrast, the gan-based classifiers only showed a
good classification performance when combined with a different data modality
and its synthesized images seemed to show less pronounced visual features than
the other two modalities. However, the ssimcc based metric shows no
substantial difference between the gan images and the other two modalities.
However, one possible reason might be that the gan learned features of
multiple classes and the images might still contain features which are derived
from images from other classes. The pca images were neither required, nor
detrimental for a good classification performance. According to the ssimcc,
the pca images were the most similar images to its clinical counterparts.

Overall, a combination of different data modalities seemed to be the key
element for achieving a good classification performance. Both ssm and pca
model the data according to a Gaussian distribution, while the gan uses an
unrestricted distribution model. The different properties of modeling the
underlying statistical distribution of a Gaussian distribution (ssm and pca)
on the one hand, and without an assumed distribution (gan) on the other hand
might have lead to a compensation of their respective disadvantage increasing
overall performance for the combinations. One limitation of this study is the
small dataset. As the clinical classification uses the same dataset for
training and validation, this might make it prone to overfitting. However, the
resulting classification metrics achieved in this study were similar to a
classification study on clinical data alone [17] which suggests that over-
fitting has not been an issue.

##  5 Conclusion

We showed that it is possible to train a classifier for different types of
craniosynostosis based solely on artificial data synthesized by a ssm, pca,
and a gan. Without having seen any clinical samples, a cnn was able to
classify four types of head deformities with an F1-score higher than 0.95 and
performed comparable to a classifier trained on clinical data. The key
component in achieving good classification results was using multiple, but
different data generation models. Overall, the ssm was the data source
contributing most to the classification performance. For the gan, using a
small image size and alternating between transposed convolutions and
interpolations were identified as key elements for suitable image generation.
The datasets and generators were made publicly available along with this work.
We showed that clinical data is not required for the classification of
craniosynostosis paving the way into cost-effective usage of synthetic data
for automated diagnosis systems.

## References

  * [1] L.Ronald French, Ian T. Jackson, and L.Joseph Melton.  A population-based study of craniosynostosis.  Journal of Clinical Epidemiology, 43(1):69–73, January 1990.
  * [2] John A. Persing, John A. Jane, and Mark Shaffrey.  Virchow and the Pathogenesis of Craniosynostosis: A Translation of His Original Work.  Plastic and Reconstructive Surgery, 83(4):738–742, April 1989.
  * [3] Anna K Coussens, Christopher R Wilkinson, Ian P Hughes, C Phillip Morris, Angela van Daal, Peter J Anderson, and Barry C Powell.  Unravelling the molecular control of calvarial suture fusion in children with craniosynostosis.  BMC Genomics, 8(1):458, December 2007.
  * [4] Sheree L. Boulet, Sonja A. Rasmussen, and Margaret A. Honein.  A population-based study of craniosynostosis in metropolitan Atlanta, 1989–2003.  American Journal of Medical Genetics Part A, 146A(8):984–991, April 2008.
  * [5] Dominique Renier, Christian Sainte-Rose, Daniel Marchac, and Jean-François Hirsch.  Intracranial pressure in craniostenosis.  Journal of Neurosurgery, 57(3):370–377, September 1982.
  * [6] Kathleen A. Kapp-Simon, Matthew L. Speltz, Michael L. Cunningham, Pravin K. Patel, and Tadanori Tomita.  Neurodevelopment of children with single suture craniosynostosis: A review.  Child’s Nervous System, 23(3):269–281, January 2007.
  * [7] Jeffrey A. Fearon, Rachel A. Ruotolo, and John C. Kolar.  Single Sutural Craniosynostoses: Surgical Outcomes and Long-Term Growth:.  Plastic and Reconstructive Surgery, 123(2):635–642, February 2009\.
  * [8] Anne Saarikko, Eero Mellanen, Linda Kuusela, Junnu Leikola, Atte Karppinen, Taina Autti, Pekka Virtanen, and Nina Brandstack.  Comparison of Black Bone MRI and 3D-CT in the preoperative evaluation of patients with craniosynostosis.  Journal of Plastic, Reconstructive & Aesthetic Surgery, 73(4):723–731, April 2020.
  * [9] Inés Barbero-García, José Luis Lerma, and Gaspar Mora-Navarro.  Fully automatic smartphone-based photogrammetric 3D modelling of infant’s heads for cranial deformation analysis.  ISPRS Journal of Photogrammetry and Remote Sensing, 166:268–277, August 2020.
  * [10] Schaufelberger, Matthias, Kühle, Reinald Peter, Wachter, Andreas, Weichel, Frederic, Hagen, Niclas, Ringwald, Friedemann, Eisenmann, Urs, Hoffmann, Jürgen, Engel, Michael, Freudlsperger, Christian, and Nahm, Werner.  A statistical shape model of craniosynostosis patients and 100 model instances of each pathology, November 2021.
  * [11] Claudia Nagel, Matthias Schaufelberger, Olaf Dössel, and Axel Loewe.  A Bi-atrial Statistical Shape Model as a Basis to Classify Left Atrial Enlargement from Simulated and Clinical 12-Lead ECGs.  In Statistical Atlases and Computational Models of the Heart. Multi-Disease, Multi-View, and Multi-Center Right Ventricular Segmentation in Cardiac MRI Challenge, volume 13131, pages 38–47. Springer International Publishing, Cham, 2022.
  * [12] Jorge Sánchez, Giorgio Luongo, Mark Nothstein, Laura A. Unger, Javier Saiz, Beatriz Trenor, Armin Luik, Olaf Dössel, and Axel Loewe.  Using Machine Learning to Characterize Atrial Fibrotic Substrate From Intracardiac Signals With a Hybrid in silico and in vivo Dataset.  Frontiers in Physiology, 12:699291, July 2021.
  * [13] Carlos S. Mendoza, Nabile Safdar, Kazunori Okada, Emmarie Myers, Gary F. Rogers, and Marius George Linguraru.  Personalized assessment of craniosynostosis via statistical shape modeling.  Medical Image Analysis, 18(4):635–646, May 2014.
  * [14] Seyed Amir Hossein Tabatabaei, Patrick Fischer, Sonja Wattendorf, Fatemeh Sabouripour, Hans-Peter Howaldt, Martina Wilbrand, Jan-Falco Wilbrand, and Keywan Sohrabi.  Automatic detection and monitoring of abnormal skull shape in children with deformational plagiocephaly using deep learning.  Scientific Reports, 11(1):17970, September 2021.
  * [15] Guido de Jong, Elmar Bijlsma, Jene Meulstee, Myrte Wennen, Erik van Lindert, Thomas Maal, René Aquarius, and Hans Delye.  Combining deep learning with 3D stereophotogrammetry for craniosynostosis diagnosis.  Scientific Reports, 10(1):15346, December 2020.
  * [16] Matthias Schaufelberger, Reinald Kühle, Andreas Wachter, Frederic Weichel, Niclas Hagen, Friedemann Ringwald, Urs Eisenmann, Jürgen Hoffmann, Michael Engel, Christian Freudlsperger, and Werner Nahm.  A Radiation-Free Classification Pipeline for Craniosynostosis Using Statistical Shape Modeling.  Diagnostics, 12(7):1516, June 2022.
  * [17] Matthias Schaufelberger, Christian Kaiser, Reinald Kühle, Andreas Wachter, Frederic Weichel, Niclas Hagen, Friedemann Ringwald, Urs Eisenmann, Jörgen Hoffmann, Michael Engel, Christian Freudlsperger, and Werner Nahm.  3D-2D Distance Maps Conversion Enhances Classification of Craniosynostosis.  IEEE Transactions on Biomedical Engineering, pages 1–10, 2023.
  * [18] J.W. Meulstee, L.M. Verhamme, W.A. Borstlap, F. Van der Heijden, G.A. De Jong, T. Xi, S.J. Bergé, H. Delye, and T.J.J. Maal.  A new method for three-dimensional evaluation of the cranial shape and the automatic identification of craniosynostosis using 3D stereophotogrammetry.  International Journal of Oral and Maxillofacial Surgery, 46(7):819–826, July 2017.
  * [19] Naiara Rodriguez-Florez, Jan L. Bruse, Alessandro Borghi, Herman Vercruysse, Juling Ong, Greg James, Xavier Pennec, David J. Dunaway, N. U. Owase Jeelani, and Silvia Schievano.  Statistical shape modelling to aid surgical planning: Associations between surgical parameters and head shapes following spring-assisted cranioplasty.  International Journal of Computer Assisted Radiology and Surgery, 12(10):1739–1749, October 2017.
  * [20] Pam Heutinck, Paul Knoops, Naiara Rodriguez Florez, Benedetta Biffi, William Breakey, Greg James, Maarten Koudstaal, Silvia Schievano, David Dunaway, Owase Jeelani, and Alessandro Borghi.  Statistical shape modelling for the analysis of head shape variations.  Journal of Cranio-Maxillofacial Surgery, 49(6):449–455, June 2021\.
  * [21] Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative Adversarial Networks, June 2014.
  * [22] Thomas Pinetz, Johannes Ruisz, and Daniel Soukup.  Actual Impact of GAN Augmentation on CNN Classification Performance:.  In Proceedings of the 8th International Conference on Pattern Recognition Applications and Methods, pages 15–23, Prague, Czech Republic, 2019. SCITEPRESS - Science and Technology Publications.
  * [23] Paolo Cignoni, Marco Callieri, Massimiliano Corsini, Matteo Dellepiane, Fabio Ganovelli, and Guido Ranzuglia.  MeshLab: An Open-Source Mesh Processing Tool.  Eurographics Italian Chapter Conference, page 8 pages, 2008.
  * [24] Nico Pietroni, Marco Tarini, and Paolo Cignoni.  Almost Isometric Mesh Parameterization through Abstract Domains.  IEEE Transactions on Visualization and Computer Graphics, 16(4):621–635, July 2010.
  * [25] Hang Dai, Nick Pears, William Smith, and Christian Duncan.  A 3D Morphable Model of Craniofacial Shape and Texture Variation.  In 2017 IEEE International Conference on Computer Vision (ICCV), pages 3104–3112, Venice, October 2017. IEEE.
  * [26] Hang Dai, Nick Pears, and William Smith.  Augmenting a 3D morphable model of the human head with high resolution ears.  Pattern Recognition Letters, 128:378–384, December 2019.
  * [27] Hang Dai, Nick Pears, William Smith, and Christian Duncan.  Statistical Modeling of Craniofacial Shape and Texture.  International Journal of Computer Vision, 128(2):547–571, February 2020.
  * [28] Mehdi Mirza and Simon Osindero.  Conditional Generative Adversarial Nets, November 2014.
  * [29] Alec Radford, Luke Metz, and Soumith Chintala.  Unsupervised Representation Learning with Deep Convolutional Generative Adversarial Networks, January 2016.
  * [30] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein GAN, December 2017.
  * [31] Ishaan Gulrajani, Faruk Ahmed, Martin Arjovsky, Vincent Dumoulin, and Aaron Courville.  Improved Training of Wasserstein GANs, December 2017.
  * [32] Christian Kaiser, Matthias Schaufelberger, Reinald Peter Kühle, Andreas Wachter, Frederic Weichel, Niclas Hagen, Friedemann Ringwald, Urs Eisenmann, Michael Engel, Christian Freudlsperger, and Werner Nahm.  Generative-Adversarial-Network-Based Data Augmentation for the Classification of Craniosynostosis.  Current Directions in Biomedical Engineering, 8(2):17–20, August 2022.
  * [33] Martin Arjovsky and Léon Bottou.  Towards Principled Methods for Training Generative Adversarial Networks, January 2017.
  * [34] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: an imperative style, high-performance deep learning library.  In Advances in neural information processing systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [35] Matthias Schaufelberger, Reinald Kühle, Andreas Wachter, Frederic Weichel, Niclas Hagen, Friedemann Ringwald, Urs Eisenmann, Jürgen Hoffmann, Michael Engel, Christian Freudlsperger, and Werner Nahm.  GAN, PCA, and Statistical Shape Models for the Creaiton of Synthetic Craniosynostosis Distance Maps, July 2023.
  * [36] Saloni Agarwal, Rami R. Hallac, Ovidiu Daescu, and Alex Kane.  Classification of Craniosynostosis Images by Vigilant Feature Extraction.  In Advances in Computer Vision and Computational Biology, pages 293–306. Springer International Publishing, Cham, 2021\.

##  Appendix A Generative adversarial network structure

This is the gan structure of generator and discriminator employed for the
creation of the synthetic data as an output of according to models’ __str__
method called via print(model).

    Generator28(
      (embed): Embedding(4, 100)
      (gen): Sequential(
        (0): Sequential(
          (0): ConvTranspose2d(200, 256, kernel_size=(5, 5), stride=(1, 1),
                    bias=False)
          (1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=True)
          (2): ReLU(inplace=True)
        )
        (1): Sequential(
          (0): Interpolate(size=(8, 8),bilinear,align_corners=True)
          (1): BatchNorm2d(256, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=True)
          (2): ReLU(inplace=True)
        )
        (2): Sequential(
          (0): Conv2d(256, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1),
                    bias=False)
          (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=True)
          (2): ReLU(inplace=True)
        )
        (3): Sequential(
          (0): Interpolate(size=(15, 15),bilinear,align_corners=True)
          (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=True)
          (2): ReLU(inplace=True)
        )
        (4): Sequential(
          (0): ConvTranspose2d(128, 128, kernel_size=(3, 3), stride=(1, 1),
                    bias=False)
          (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=True)
          (2): ReLU(inplace=True)
        )
        (5): Sequential(
          (0): Interpolate(size=(30, 30),bilinear,align_corners=True)
          (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=True)
          (2): ReLU(inplace=True)
        )
        (6): Conv2d(128, 1, kernel_size=(3, 3), stride=(1, 1),
                    bias=False)
        (7): Tanh()
      )
    )

    Discriminator28(
      (net): Sequential(
        (0): Sequential(
          (0): Conv2d(2, 32, kernel_size=(4, 4), stride=(2, 2), padding=(1, 1),
                    bias=False)
          (1): InstanceNorm2d(32, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=False)
          (2): LeakyReLU(negative_slope=0.2)
        )
        (1): Sequential(
          (0): Conv2d(32, 128, kernel_size=(4, 4), stride=(2, 2), padding=(1, 1),
                    bias=False)
          (1): InstanceNorm2d(128, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=False)
          (2): LeakyReLU(negative_slope=0.2)
        )
        (2): Sequential(
          (0): Conv2d(128, 256, kernel_size=(5, 5), stride=(2, 2), padding=(1, 1),
                    bias=False)
          (1): InstanceNorm2d(256, eps=1e-05, momentum=0.1, affine=True,
                    track_running_stats=False)
          (2): LeakyReLU(negative_slope=0.2)
        )
        (3): Conv2d(256, 1, kernel_size=(3, 3), stride=(1, 1))
      )
      (embed): Embedding(4, 784)
    )

##  Appendix B Failed GAN attempts

We show artifacts arising from only using transposed convolutional layers
(ConvTranspose2d), using only up-scaling interpolation layers (Interpolate),
or from large gradient penalties which prohibits training in Fig. 12.

Figure 12: Failed gan images arranged in a 2$\times$2 grid with artifacts
arising from poor network design or bad training conditions. From left to
right: Deconvolution artifacts, interpolation artifacts, and noise artifacts.