# Beyond optimization - supervised learning applications in relativistic
laser-plasma experiments

Jinpu Lin<EMAIL_ADDRESS>Qian Qian  Jon Murphy  Center for Ultrafast
Optical Science, University of Michigan, Ann Arbor, Michigan 48109, USA
Abigail Hsu  Department of Applied Mathematics and Statistics, Stoney Brook
University, Stony Brook, NY 11794, USA  Yong Ma  Center for Ultrafast Optical
Science, University of Michigan, Ann Arbor, Michigan 48109, USA  Alfred Hero
Department of Electrical Engineering and Computer Science, University of
Michigan, Ann Arbor, Michigan 48109, USA  Alexander G.R. Thomas  Karl
Krushelnick  Center for Ultrafast Optical Science, University of Michigan, Ann
Arbor, Michigan 48109, USA

###### Abstract

We explore the applications of machine learning techniques in relativistic
laser-plasma experiments beyond optimization purposes. We predict the beam
charge of electrons produced in a laser wakefield accelerator given the laser
wavefront change caused by a deformable mirror. Machine learning enables
feature analysis beyond merely searching for an optimal beam charge, showing
that specific aberrations in the laser wavefront are favored in generating
higher beam charges. Supervised learning models allow characterizing the
measured data quality as well as recognizing irreproducible data and potential
outliers. We also include virtual measurement errors in the experimental data
to examine the model robustness under these conditions. This work demonstrates
how machine learning methods can benefit data analysis and physics
interpretation in a highly nonlinear problem of relativistic laser-plasma
interaction.

††preprint: AIP/123-QED

##  I Introduction

High-repetition-rate laser systems have been widely used with evolutionary
algorithms to solve optimization problems in the field of relativistic laser-
plasma interactions, including laser wakefield acceleration He  _et al._
(2015); Dann  _et al._ (2019); Lin  _et al._ (2019), ion acceleration Smith
_et al._ (2020); Noaman-ul Haq  _et al._ (2018), x-ray production Streeter
_et al._ (2018), terahertz generation Hah  _et al._ (2017), laser
filamentation Englesbe  _et al._ (2016); Lefebvre  _et al._ (2018), and laser
focus optimization Nayuki  _et al._ (2005); Lin  _et al._ (2018). However,
evolutionary algorithms usually provide little information other than a local
optimum, which can be hard to interpret. Instead, machine learning methods can
generate predictive models that reveal more information in the dataset and
help understand the physical processes that occur in an experiment.

The broader discipline of plasma physics has adopted various machine learning
methods in recent years. For instance, supervised learning regression
algorithms have been applied to inertial confinement fusion (ICF) experiments
with growing interests, such as Deep Jointly-Informed Neural Networks Gaffney
_et al._ (2019); Hsu, Cheng, and Bradley (2020); Humbird, Peterson, and
McClarren (2019); Maris and Khan (2019) and Gaussian Process regressor
Hatfield  _et al._ (2019). Another popular machine learning technique called
Random Forest has found success in magnetic confinement fusion experiments for
both classification and regression problems Rea  _et al._ (2018); Piccione
_et al._ (2020). In space physics, Gaussian Processes are used to classify
solar wind plasmas into categories Camporeale, Carè, and Borovsky (2017), and
deep neural networks are used to predict solar flares from sunspot data Jiao
_et al._ (2020); Chen  _et al._ (2019). Beyond supervised learning, plasma
physicists find interest in other powerful and increasingly popular machine
learning methods, such as transfer learning Humbird  _et al._ (2019) and
reinforcement learning Witman  _et al._ (2019); Kain  _et al._ (2020). The
laser-plasma community is starting to embrace machine learning techniques as
well. Artificial neural networks are employed to analyze features in high-
order-harmonic spectra Gonoskov  _et al._ (2019) and laser-induced-breakdown
spectra Li  _et al._ (2020). Our work explores the capability of machine
learning techniques in the field of laser-wakefield acceleration using all the
supervised learning methods mentioned above.

Laser wakefield accelerators (LWFAs), first proposed by Tajima and Dawson
Tajima and Dawson (1979), provide a possible alternative to conventional
particle accelerators at a substantially lower size and cost. Taking advantage
of the electric fields in laser-produced plasmas, LWFA can reach acceleration
gradients of 100 GeV/m, which is three orders of magnitude greater than those
produced in conventional accelerators. Extensive studies have been performed
to understand LWFA mechanisms and to experimentally demonstrate energetic
electron beams Malka  _et al._ (2002); Mangles  _et al._ (2004); Geddes  _et
al._ (2004); Faure  _et al._ (2004); Esarey, Schroeder, and Leemans (2009),
and the highest electron energy achieved so far is 7.8 GeV Gonsalves  _et al._
(2019). While the highest energy facilities usually fire a few shots a day,
there has been growing interests in having higher repetition rate lasers with
slightly lower peak power and corresponding plasma targets Roso (2018);
Prencipe  _et al._ (2017); Salehi (2019); Salehi  _et al._ (2017); Feister,
Poole, and Heuer (2019). One of the rationales to increase the repetition rate
is to compare LWFAs to conventional accelerators that operate at tens of Hz or
above. In addition, it allows meeting the requirements for statistical methods
to have better control over experiments He  _et al._ (2015); Dann  _et al._
(2019); Lin  _et al._ (2019); Kim  _et al._ (2017); Liu  _et al._ (2014); Tsai
_et al._ (2018); Shalloo  _et al._ (2020).

In this work, we have predicted the electron beam charges in LWFA given the
laser wavefront modification caused by a deformable mirror using four
supervised learning regression methods: Random Forest, Neural Network, Deep
Jointly-Informed Neural Network, and Gaussian Process. We benchmark these
models and demonstrate applications beyond optimization. We show that
generating higher beam charges favors specific wavefront aberrations, which is
revealed by ranking the feature importance in form of the Zernike
decomposition. We analyze the experimental data quality by evaluating the
model performance on every measured data point. We also characterize the model
robustness against a range of virtual measurement error bars assigned to the
experimental data.

##  II Data and Methods

###  II.1 Data

####  II.1.1 Experimental

The dataset is taken from the recent experiments performed in the Gerard
Mourou Center for Ultrafast Optical Science at the University of Michigan. The
experiments were performed using the Lambda Cubed laser system which provides
35 fs, 20 mJ laser pulses at 480 Hz. The laser wakefield accelerated electron
beam data used here comes from four separate experiment days when running
optimization algorithms. The laser wavefront change was induced by a
deformable mirror and recorded as the 37 actuators’ voltages on the mirror
surface. The electron beam was captured via a scintillator screen that was
imaged by an electron-multiplying CCD camera (Andor Luca-R, 14-bit). The
optics which deliver the laser to the experimental chamber were unchanged
throughout the duration of the four days, and the alignment procedures are a
routine that is carried out the same way each day.

####  II.1.2 Pre-processing

The dataset is pre-processed before the regression modeling. The dataset
contains the information of the optical wavefront change of the driving laser
caused by a deformable mirror as well as the electron beam charge of the
accelerated electrons. The wavefront changed by the deformable mirror has 37
dimensions in space and can be described mathematically by a polynomial, known
as the Zernike polynomial Born and Wolf (2013), to reduce the dimension. In
this case, the coefficients of the first 5 layers (15 terms) in the Zernike
polynomial can accurately reproduce the wavefront. The 15-dimensional vectors
consisting of the Zernike coefficients are used as the input to our supervised
learning models, while the electron beam charges are the output, normalized to
the range (0, 1]. In the context of machine learning, the input is called a
feature and the output is called a label. We have 208 data samples in total.
The dataset is split into two subsets: $80\%$ of the data points are used to
train the models while $20\%$ are for testing. The feature matrix in the
training set is sphered so that its rows have zero sample mean and unity
sample variance. The feature matrix in the test set is updated accordingly
with the same transformation.

####  II.1.3 Correlation

In statistics, the correlation coefficient measures the dependence between two
variables and its value falls in the range of [-1,1]. The absolute value
represents the strength of the dependence while the sign represents the
direction. We calculate the correlation between each feature (Zernike
coefficients) and the electron beam charge in the test dataset, as is
illustrated in Table. 1. Among all the Zernike coefficients, while $z_{10}$
has the largest magnitude of correlation. The correlation matrix will be
compared to the machine learning model predictions in Table. 3, 4 in Section.
III.

Table 1: Correlation between Zernike terms and beam charge.

$z_{0}$ | $z_{1}$ | $z_{2}$ | $z_{3}$ | $z_{4}$ | $z_{5}$ | $z_{6}$ | $z_{7}$
---|---|---|---|---|---|---|---
0.51 | -0.33 | -0.28 | 0.29 | 0.14 | 0.21 | -0.15 | 0.37
$z_{8}$ | $z_{9}$ | $z_{10}$ | $z_{11}$ | $z_{12}$ | $z_{13}$ | $z_{14}$ | charge
-0.42 | 0.20 | -0.55 | 0.35 | 0.47 | -0.42 | -0.019 | 1

###  II.2 Machine learning methods

Four supervised learning regression methods are used to predict the electron
beam charges based on laser wavefront changes. Supervised learning is a branch
of machine learning, which learns a function that maps an input (feature) to
an output (label) based on example input-output pairs in a training
sampleRussell and Norvig (2002). In each of the following supervised learning
methods, the model is trained on the training dataset recursively until it can
accurately predict the labels using the features. The model performance is
then characterized by the test dataset.

####  II.2.1 Random Forest

The Random forest (RF) regressorLiaw, Wiener  _et al._ (2002) is a popular
bagged algorithm for high-dimensional and nonlinear regression. It is based on
the concept of a decision tree, which splits the dataset along some dimensions
and recursively divides the space into regions with similar labels. Being the
most popular bagging (Bootstrap Aggregation) ensemble algorithm, random forest
samples, with replacement at uniform probability, the original dataset $D$
into $m$ datasets ($D_{1},D_{2},...,D_{m}$) with the same size as $D$. For
instance, if the original dataset contains three samples $D=[a,b,c]$, then
$D_{1}$ could be $[a,c,c]$. For each dataset $D_{j}$ in the forest, we train a
full decision tree by splitting the data to $k<n$ dimensions. Only k features
are to be considered when looking for the best split. Since the trees become
much more different as they select different features, we have to increase the
number of trees and average over individual regressors. This bagging process
helps reduce variance effectively. Denote $h_{i}$ as the regressor of the
$i_{th}$ tree, then the bagged regressor is Weinberger (2018):

$h_{i}=\frac{1}{m}\sum_{i=1}^{m}h_{i}(x)$ (1)

In this study, we implement the algorithm using the
Sklearn.ensemble.RandomForestRegressor library in Scikit-learn Pedregosa  _et
al._ (2011). The hyper-parameters to tune are the number of trees, the maximum
depth in a tree, and the maximum number of features when splitting.

####  II.2.2 Deep Neural Network

A Deep Neural Network (DNN) is a feed-forward artificial neural network with
multiple hidden layers. The goal of a feed-forward network is to approximate
some function $f$ Goodfellow, Bengio, and Courville (2016). A typical one
hidden layer can be mathematically described with weight $w$ and bias $w_{0}$
by Eq. (2):

$y=f(w_{0}+w_{1}^{T}x)$ (2)

According to the Universal Approximation Theorem, a feed-forward network with
even one hidden layer can approximate any continuous function from one finite
dimensional space to another under some conditions Hornik (1991). Although
practically it may lead to an infeasibly large layer and fail to generalize,
DNNs are powerful function approximators when learned properly. Compared to
shallow models, DNNs usually can extract better features and learn more
effectively.

In this work, we build a fully-connected five-layer DNN using the
Tensorflow.Keras library Chollet  _et al._ (2015) based on Google’s deep
learning software TensorFlow Abadi  _et al._ (2016). When constructing the
network, we use the rectified linear unit (ReLU) function and the Sigmoid
function as the activation functions for different layers. The cost function
is the mean squared error loss governed by the Adam optimizer Kingma and Ba
(2014) to update the network weights. A $L_{2}$ norm regularization is added
to the loss function to reduce overfitting. The main tuning parameters are the
number of layers, the number of neurons in each layer, the epoch size, and the
initialization of the weight matrix.

####  II.2.3 Deep Jointly-Informed Neural Networks

Deep Jointly-Informed Neural Networks (DJINN) is a machine learning algorithm
that constructs and initializes the deep feedforward neural networks based on
decision trees. It was developed by Humbird et al.Humbird, Peterson, and
McClarren (2018) and it has shown success in training ICF datasets Humbird,
Peterson, and McClarren (2019); Gaffney  _et al._ (2019); Hsu, Cheng, and
Bradley (2020) as well as standard regression datasets such as Boston housing
prices, California housing prices, and diabetes disease progression Humbird,
Peterson, and McClarren (2018). The algorithm starts by constructing a
decision tree or an ensemble of trees, where the number of trees will be the
number of networks in the later stage. It then maps the decision trees to deep
neural networks by taking the decision paths as guidance for the network
architecture and weight initialization. The networks are trained with
backpropagation using TensorFlow Abadi  _et al._ (2016). In the network
architecture, the activation function is the rectified linear unit and the
cost function is governed by Adam optimizer Kingma and Ba (2014). Without
optimizing the architecture of the neural networks, DJINN displays comparable
performance to optimized architectures at a significantly lower computational
cost using their datasets.

The DJINN regression source code is accessible at the LLNL/DJINN github
directory. The main tuning parameters for this study are the maximum depth of
trees, and the number of trees (nets) in ensemble.

####  II.2.4 Gaussian Process

The Gaussian process (GP) is a non-parametric Bayesian algorithm for
supervised learning problems Rasmussen (2003). While most machine learning
algorithms fit the dataset into a model function with weight parameters and
use that function to make predictions, Bayesian methods avoid the intermediate
step and make predictions directly from the dataset. This is achieved by
integrating all possible weight functions in the universe Weinberger (2018):

$P(y|x,D)=\int_{w}P(y|x,w)\cdot P(w|D)\;dw$ (3)

where P is the prediction, w is the weight matrix, and D is the dataset. In
Gaussian process regression, we assume that the data can be fit by some model
with weight function w and a Gaussian distributed noise $\epsilon$:
$y=f(x)=w^{T}x+\epsilon$. Thus the term $P(y|x,w)$ in Eq. 3 is a Gaussian
distribution. The second term on the right-hand-side can also be proved to be
Gaussian using the Bayes’ rule:

$P(w|D)=\frac{P(D|w)\;P(w)}{P(D)}$ (4)

where $P(w)$ is the prior distribution, $P(w|D)$ is the posterior distribution
after D is considered, and $P(D)$ is a normalization. By choosing a Gaussian
prior, Bayes’ rule leads to a Gaussian posterior $P(w|D)$. Marginalizing out w
in Eq. 3, we know $P(y|x,D)$ is also in a Gaussian distribution. The Gaussian
process has a closed form posterior distribution that can be used to quantify
uncertainy of the estimate (posterior mode) through a confidence interval.

Another advantage of Gaussian process regression is that we can specify prior
information about the shape of the model by selecting certain kernel
functions. As is introduced previously, the probability function of the
prediction can be expressed by a Gaussian distribution:
$P(y|x,D)=N(\mu,\Sigma)$, where $\mu$ is the mean and $\Sigma$ is the
covariance matrix. The covariance matrix can be kernelized using selected
kernel functions. In this project, we implement the algorithm using the
$Sklearn.gaussian\\_process$ library in Scikit-learn Pedregosa  _et al._
(2011) with a combination of Matern kernel and Rational Quadratic kernel. The
hyper-parameters to tune are the smoothness, the length-scale, and the scale
mixture parameter in the kernels.

##  III Results

All codes are written in Python. After training the models, we predict the
electron beam charge using the laser wavefront change in the test dataset.
Predicted electron beam charges using the above models are shown in Fig. 1
against measured electron beam charges. A reference line at $45^{\circ}$ is
included and data points closer to the reference line are considered better
predictions. The bottom left corner of the plot is zoomed in and shown to the
right. Detailed statistical evaluations are summarized in Table. 2, in which
we report the mean-square-error (MSE), mean-absolute-error (MAE), R-squared
($R^{2}$) and explained variance score (ExVar) based on the predicted charge
and the measured charge. MSE measures the average squared difference between
the predictions and the real values, which contains information of both
variance and bias. It is the most popular metric when evaluating machine
learning models and we use MSE as the target for the hyperparameter tuning
process. The downside of MSE is that it can be sensitive to outliers, which
MAE handles better by measuring the absolute error instead of the squared
error. $R^{2}$ is the proportion of variance of the measured value from the
prediction. It tells how likely a new sample (out of the dataset) can be
predicted by the model. Explained variance considers bias on top of $R^{2}$.
It is the same as $R^{2}$ if the mean of error is 0. In general, we’d like to
have MSE and MAE close to 0 while $R^{2}$ and ExVar close to unity.

Figure 1: Predicted vs. measured electron beam charge in test dataset. Model | MSE | MAE | $R^{2}$ | ExVar
---|---|---|---|---
Random Forest | 0.00132 | 0.0268 | 0.986 | 0.987
Neural network | 0.00162 | 0.0292 | 0.983 | 0.984
DJINN | 0.00154 | 0.02741 | 0.98403 | 0.98404
Gaussian Process | 0.00185 | 0.0305 | 0.981 | 0.981
Table 2: Evaluation matrix.

All four models demonstrate similar statistics in Table. 2, though RF performs
slightly better and GP gives the largest MSE and MAE scores while the smallest
$R^{2}$ and ExVar scores. It also predicts negative values when the electron
beam charges are small. However, it does not necessarily mean that Random
Forest is the best model and gaussian Process is the worst. The results in the
evaluation matrix are sensitive to the way we split the training set and test
set. We will show in the next section that training and evaluating the model
on different data points yields different results. We will also present more
analyses such as the model consistency against measurement errors and
overfitting-related issues.

####  III.0.1 Data quality

Figure 2: Prediction error of every data point.

Experimental measurements in LWFA can suffer from a lack of reproducibility
and may have outliers in the dataset. Possible sources include the pointing
instabilities and shot to shot fluctuations in high power laser systems, as
well as the irreproducibility in the plasma density profile from gas jets. To
justify the quality of our measured data, we make a prediction using the
wavefront change of each data point and compare that to the corresponding
measured electron beam charge. Instead of splitting the dataset into $80\%$
for training and $20\%$ for testing, here we test on only one data point while
all the other data points are used to train the model. We then compare the
predicted beam charge to the measured electron beam charge of this particular
data point and calculate their difference. This process is looped over the
entire dataset. Fig. 2a plots the prediction error
($\sigma=|y_{predicted}-y_{measured}|$) at each data point in random order
using GP, RF, and DJINN. The prediction errors from DNN are shown in Fig. 2b.
It is observed that the three models in Fig. 2a have similar performance while
the results from DNN in Fig. 2b is very different. We interpret three messages
from these two plots. 1. Prediction errors vary across seven orders of
magnitude, i.e., some data points can be accurately predicted ($\sigma\sim
1e-8$) while some data points can hardly be predicted ($\sigma\sim 0.1$).
Therefore selecting different data points into the training or test set can
lead to different evaluation matrices from the one in Table. 2. 2. This huge
variation can be caused by either the inconsistency of the model across data
points, or this specific data being very different from all the other data
points in the dataset that are used to train the model. We are satisfied with
the reliability of the models, as is characterized in Table. 2. In addition,
the similarity among the three models’ performance in Fig. 2a suggests that
the models are less likely to be inconsistent but in a similar manner. Thus it
provides a potential characterization of the quality of each data point, i,e.,
a data point that has a very large prediction error in all three models can be
considered as having poor data quality and may be dropped as an outlier. 3.
DNN performs differently from the other models, suggesting that it overfits
the data set and it is less reliable in this scenario. It matches the learning
curves in Fig. 4, and a detailed discussion on the overfitting issue can be
found later.

This could have some practical use in analyzing experimental results. For
example, if we drive the same laser pulse into the wakefield accelerator
twice, we might observe different electron beam charges due to a lack of
reproducibility in LWFA. A typical solution would be to calculate the mean
value and the error bar. However, this could include misinformation if there
is a "bad" measurement that should be dropped. By performing the above
analysis we would be able to tell which one of the measured beam charges is
more reliable. Moreover, it can help identify outliers not only in repeated
measurements but anywhere in the parameter space that the experiment scans
across. The labeled data point in Fig. 2a is an example of possibly poor data
quality, where all three models have large prediction errors at this point.

####  III.0.2 Robustness against measurement errors

Figure 3: Model performance against virtual measurement errors.

Studying science with error bars has been a rising trend as experimental
measurements often come with such. In this section, we investigate the
performance of these models against measurement errors of the electron beam
charges. Since measurement errors were not recorded during our experiments, we
include some virtual error bars to every measured electron beam charge. At
each measurement, the true value is assumed to lie in the range of measured
value $\cdot(1\pm X\%)$, where $X=1,3,5,10)$, and 1000 points are drawn
randomly from a normal distribution within this range. Therefore we get 1000
copies of the original dataset with the same wavefront but different electron
beam charges. The reason to have 1000 datasets is to generate enough
statistics to justify the model performance against unsure measurements.
Results are presented in Fig . 3.

Fig. 3a shows the distribution of the test MSE using RF. Each colored line is
generated from 1000 MSEs. During the training process, the model configuration
was kept the same among the 1000 datasets but the weight learning was updated
in each dataset. Measurement errors that define the range of the dataset
fluctuation are $1\%,\;3\%,\;5\%$ and $10\%$, while the corresponding test MSE
distributions are plotted in red, green, blue and black, respectively.
Statistical analysis is summarized in the adjacent table to the right,
illustrating the mean value and standard deviation ($\sigma$) of the test MSE
distribution as well as the percentage of points that fall within one, two, or
three standard deviations around the mean value. Fig. 3b-d present the results
using GP, DJINN, and DNN. The four models share some common performances. The
mean test MSE value increases with the measurement error, which means it is
more likely to make a less accurate prediction when the measurement itself is
less accurate, as expected. The standard deviation also increases with the
measurement error, suggesting a less consistent or less precise model
prediction at larger measurement errors. There are noticeable differences in
the last three table columns of the four models. Remember that the virtual
measurement errors are drawn from a perfect normal distribution, where the
percentage of values that lie within one, two, or three standard deviations
around the mean value are $68.3\%$, $95.5\%$, and $99.7\%$, respectively. As
is shown in the tables, RF and GP retain almost normal distributions in the
sample prediction, DNN gives normal-like distributions and the results from
DJINN are far from normal distributions.

####  III.0.3 Learning curve

Figure 4: Learning curves in Neural Network

In machine learning, overfitting often occurs when a model has learned the
training data so well that it also learns the statistical noise or random
fluctuations in the data. The learning curve is an intuitive tool to visualize
the degree of overfitting. Fig. 4 shows the learning curves in our Neural
Network model, which plots MSE at each epoch for both training data (black)
and test data (red). The training curve tells how well the model learns while
the test curve tells how well the model generalizes. Since the red curve does
not decrease as much as the black curve does, the model overfits. In other
words, the overfitted model performs worse outside the training dataset. Fig.
4 also plots learning curves considering measurement errors in blue and orange
to better reproduce experimental conditions. Measurement errors are included
in a similar way to the ones in the previous section, where the true value is
assumed to lie in the range of measured value$\cdot(1\pm 3\%)$. Instead of
generating 1000 copies of datasets at a time, here we generate only one
dataset with measurement errors and update the measurement errors at every
epoch in the learning process. Obtained learning curves are noisier but
demonstrate less overfitting, as is shown in Fig. 4. After applying
measurement errors, the training curve moves higher but the test curve does
not shift much. Namely, the model finds it harder to learn but it is still
able to make equally accurate predictions. Therefore, including virtual
measurement errors is beneficial as it not only represents practical
experimental conditions better but also decreases overfitting. Note that we
have also tried this using the DJINN model but the learning curves with and
without measurement errors look almost the same. It is not surprising since
the DJINN model does not overfit as much as the neural network model does.

####  III.0.4 Feature importance

We have trained models to predict the electron beam charge upon laser
wavefront modification represented by the first 15 Zernike coefficients. It is
natural to ask how sensitive the beam charge is to these features. We evaluate
the feature importance using our four models and compare them to the
correlation ranking, summarized in Table. 3 and 4. In each row, we list the
four most important features decided by that model, while the numbers are the
orders of the Zernike coefficients. The importance of a feature is measured by
calculating the increase in the model’s prediction error after setting the
feature values to a constant. A feature is said to be important if the
prediction error increases significantly, and less important when the
prediction error does not change much. When evaluating each feature, the
values of this feature of all test data are set to their mean value. Note that
model training is performed prior to this process and the training data are
not modified. It is debatable whether the importance should be computed on the
test data or the training data. The former tells how much the model relies on
each feature for making predictions while the latter tells how much the
feature contributes to the performance of the model on unseen data Molnar
(2020). We show the feature importance computed from both ways. In Table. 3 we
split the dataset to a training set ($80\%$) and a test set ($20\%$) and
evaluate the feature importance on the test set. In Table. 4 we train the
model using the entire dataset and measure the feature importance also on the
entire dataset.

Depending on the model and the evaluation data, the feature importance
rankings in Table. 3 and 4 are slightly different. However, the $0^{th}$,
$1^{st}$, $6^{th}$ and $10^{th}$ Zernike terms are generally believed to be
the more important ones, physically representing the piston, tilt, vertical
trefoil and oblique quadrafoil in wavefront aberration, respectively. It
suggests that controlling these features would be more effective at producing
high electron beam charges in this scenario.

Model | $1^{st}$ | $2^{nd}$ | $3^{rd}$ | $4^{th}$
---|---|---|---|---
Random Forest | 0 | 6 | 1 | 10
Gaussian Process | 1 | 10 | 0 | 6
Neural Network | 0 | 10 | 1 | 3
DJINN | 1 | 10 | 0 | 6
Correlation | 0 | 10 | 1 | 13
Table 3: Evaluate feature (Zernike coefficient) importance on the test dataset. Model | $1^{st}$ | $2^{nd}$ | $3^{rd}$ | $4^{th}$
---|---|---|---|---
Random Forest | 0 | 6 | 1 | 12
Gaussian Process | 10 | 1 | 0 | 6
Neural Network | 0 | 10 | 1 | 8
DJINN | 1 | 10 | 0 | 6
Correlation | 10 | 0 | 12 | 13
Table 4: Evaluate feature (Zernike coefficient) importance on the entire
dataset.

##  IV Discussion

There has been growing interest in the laser-plasma community to discuss what
machine learning can bring us and what it can not. Machine learning methods
are not expected to offer some generalized predictive models that save
experimentalists from carrying out every experiment. This is due to the lack
of reproducibility in high-power laser-plasma experiments, where laser systems
usually suffer from shot-to-shot fluctuations and a given plasma density
profile is hard to duplicate. However, using machine learning techniques can
help us better understand the experiment we did and improve the design of
next-step experiments. It enables deeper physics interpretation of the data
since the predictive accuracy of the regression models is determined by the
data quality. For example, by ranking the feature importance we are able to
identify the Zernike terms that are most sensitive to noise. Some of them
turned out to be the high order terms (vertical trefoil and oblique
quadrafoil) that locate at the edge of the wavefront. This can be explained
from an experimental point of view as follows. Deformable mirrors are
manufactured to have actuators forming into a hexagon matrix, while a
wavefront is usually defined in a circular or rectangular shape, leading to
lack of information on the edge. If the wavefront is measured directly with
wavefront sensors, it is usually required to manually draw a circle that
covers most lighted pixels on the detector as the region of interest. As a
result, uncertainty arises at the pixels on the edges of the region of
interest. If the wavefront is reconstructed using the actuator displacement on
the deformable mirror surface, the phase in these unknown edges also needs to
be defined and is set to be nan in our case. Another possible source of noise
is the imperfect overlap of the laser and the deformable mirror surface:
either the laser beam clipped off the mirror edge or it did not fulfill the
whole mirror surface.

Including virtual measurement uncertainties can be useful even if the
experimental data come with some uncertainties. It may be able to narrow down
the range of uncertainty when the measurement uncertainty is large. For
instance, if a data point in the dataset has some large uncertainty and we
would like to know or at least narrow down the true measurement value, machine
learning regression methods can provide us a possible solution. We can start
by randomly sampling N points within the range of uncertainty of this
measurement value. The next step is to make N copies of the measurement
dataset and replace the point of large uncertainty with one of our sampled
point in each dataset. Therefore, we obtain N similar datasets with difference
only at one point. We then train and test the machine learning models on each
dataset. Those datasets who lead to less accurate predictions (large test MSE)
are less likely to contain true measurement at the uncertain point if the
models make sense.

It is also worth discussing the kernel functions in the Gaussian Process
method. The common way of kernel selecting is either to have expert knowledge
about the dataset or to compare candidate kernels for the best performance. In
this project, we have tried different kernels and decide that a combination of
Matern kernel and Rational Quadratic kernel works best. We can thus infer some
knowledge about our dataset based on the rationale of these kernels. Matern is
a generalization of the popular Gaussian radial basis function (RBF) with
tunable smoothness Rasmussen (2003). The smoothness of the model can be
controlled by a parameter $\mu$ while the $\mu$ value in our case is as low as
0.3, suggesting that the resulting function is far from smooth. The other
kernel that fits our model, the Rational Quadratic kernel, is a sum of RBF
kernels with different length scales. Note that the length scale decides the
safe distance to extrapolate when modeling, and discontinuity can be handled
with a short length scale Duvenaud (2014). According to the optimized
parameters in these two kernels, our model function is neither smooth nor
continuous, which is not surprising as the high-dimensional dataset was taken
from a highly-nonlinear physical process.

In summary, we demonstrate several applications of machine learning in
relativistic laser-plasma experiments beyond optimization purposes. We have
built four supervised learning regression models to predict the electron beam
charge using the laser wavefront change in an LWFA experiment. All four models
present similar statistics in the evaluation matrix, although Random Forest
performs slightly better and Gaussian Process performs slightly worse. To
justify the data quality affected by the irreproducibility in experiments, we
characterize the model prediction on every single data point. Three of the
models show similar performance, providing a potential way of recognizing
outliers without repeated measurements. The Deep Neural Network is easy to
overfit our dataset and thus not the best candidate for analyzing data
quality. We include virtual measurement errors to the measured electron beam
charges and DJINN is found to be the most sensitive to measurement
fluctuations. Having virtual measurement errors is beneficial as it not only
represents experimental conditions better but also decreases overfitting. All
four models give similar feature importance rankings, which are also in line
with the statistical correlation of the dataset. The Deep Neural Network
requires the most computational cost, followed by DJINN, while Gaussian
Process and Random Forest consume the least. Therefore, Random Forest is
recommended when working with datasets from similar laser-plasma experiments.

###### Acknowledgements.

J. Lin is supported by Air Force Office of Scientific Research (AFSOR)
FA9550-16-1-0121. A. Hero is supported by the Consortium for Enabling
Technologies and Innovation under Dept of Energy grant DE-NA0003921. The
authors would like to thank Tony Zhang, Gang Yang, and Yirui Liu for their
contribution in the preparation stage. The authors would like to thank Andrew
Maris for the insightful discussion.

## References

  * He  _et al._ (2015) Z.-H. He, B. Hou, V. Lebailly, J. Nees, K. Krushelnick, and A. Thomas, “Coherent control of plasma dynamics,” Nature communications 6, 1–7 (2015).
  * Dann  _et al._ (2019) S. Dann, C. Baird, N. Bourgeois, O. Chekhlov, S. Eardley, C. Gregory, J.-N. Gruse, J. Hah, D. Hazra, S. Hawkes, _et al._ , “Laser wakefield acceleration with active feedback at 5 hz,” Physical Review Accelerators and Beams 22, 041303 (2019).
  * Lin  _et al._ (2019) J. Lin, Y. Ma, R. Schwartz, D. Woodbury, J. Nees, M. Mathis, A. Thomas, K. Krushelnick, and H. Milchberg, “Adaptive control of laser-wakefield accelerators driven by mid-ir laser pulses,” Optics express 27, 10912–10923 (2019).
  * Smith  _et al._ (2020) J. R. Smith, C. Orban, J. T. Morrison, K. M. George, G. K. Ngirmang, E. A. Chowdhury, and W. M. Roquemore, “Optimizing laser-plasma interactions for ion acceleration using particle-in-cell simulations and evolutionary algorithms,” arXiv preprint arXiv:2005.07746 (2020).
  * Noaman-ul Haq  _et al._ (2018) M. Noaman-ul Haq, T. Sokollik, H. Ahmed, J. Braenzel, L. Ehrentraut, M. Mirzaie, L.-L. Yu, Z. Sheng, L. Chen, M. Schnürer, _et al._ , “Controlling laser driven protons acceleration using a deformable mirror at a high repetition rate,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 883, 191–195 (2018).
  * Streeter  _et al._ (2018) M. Streeter, S. Dann, J. Scott, C. Baird, C. Murphy, S. Eardley, R. Smith, S. Rozario, J.-N. Gruse, S. Mangles, _et al._ , “Temporal feedback control of high-intensity laser pulses to optimize ultrafast heating of atomic clusters,” Applied Physics Letters 112, 244101 (2018).
  * Hah  _et al._ (2017) J. Hah, W. Jiang, Z. He, J. Nees, B. Hou, A. Thomas, and K. Krushelnick, “Enhancement of thz generation by feedback-optimized wavefront manipulation,” Optics Express 25, 17271–17279 (2017).
  * Englesbe  _et al._ (2016) A. C. Englesbe, Z. He, J. A. Nees, A. G. Thomas, A. Schmitt-Sody, and K. Krushelnick, “Control of the configuration of multiple femtosecond filaments in air by adaptive wavefront manipulation,” Optics Express 24, 6071–6082 (2016).
  * Lefebvre  _et al._ (2018) S. Lefebvre, H. Belmouaddine, D. Morris, and D. Houde, “Phase control algorithms and filamentation of ultrashort laser pulses in a scattering medium,” Applied Physics B 124, 210 (2018).
  * Nayuki  _et al._ (2005) T. Nayuki, T. Fujii, Y. Oishi, K. Takano, X. Wang, A. A. Andreev, K. Nemoto, and K.-i. Ueda, “Production of a mev proton with 30 mj laser energy by optimizing the focusing spot using a deformable mirror,” Review of scientific instruments 76, 073305 (2005).
  * Lin  _et al._ (2018) J. Lin, J. H. Easter, K. Krushelnick, M. Mathis, J. Dong, A. Thomas, and J. Nees, “Focus optimization at relativistic intensity with high numerical aperture and adaptive optics,” Optics Communications 421, 79–82 (2018).
  * Gaffney  _et al._ (2019) J. A. Gaffney, S. T. Brandon, K. D. Humbird, M. K. Kruse, R. C. Nora, J. L. Peterson, and B. K. Spears, “Making inertial confinement fusion models more predictive,” Physics of Plasmas 26, 082704 (2019).
  * Hsu, Cheng, and Bradley (2020) A. Hsu, B. Cheng, and P. A. Bradley, “Analysis of nif scaling using physics informed machine learning,” Physics of Plasmas 27, 012703 (2020).
  * Humbird, Peterson, and McClarren (2019) K. D. Humbird, J. L. Peterson, and R. G. McClarren, “Parameter inference with deep jointly informed neural networks,” Statistical Analysis and Data Mining: The ASA Data Science Journal 12, 496–504 (2019).
  * Maris and Khan (2019) A. Maris and S. Khan, “Finding correlations in inertial confinement fusion experimental data using machine learning,” Tech. Rep. (Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States), 2019).
  * Hatfield  _et al._ (2019) P. Hatfield, S. Rose, R. Scott, I. Almosallam, S. Roberts, and M. Jarvis, “Using sparse gaussian processes for predicting robust inertial confinement fusion implosion yields,” IEEE Transactions on Plasma Science 48, 14–21 (2019).
  * Rea  _et al._ (2018) C. Rea, R. Granetz, K. Montes, R. Tinguely, N. Eidietis, J. Hanson, and B. Sammuli, “Disruption prediction investigations using machine learning tools on diii-d and alcator c-mod,” Plasma Physics and Controlled Fusion 60, 084004 (2018).
  * Piccione  _et al._ (2020) A. Piccione, J. Berkery, S. Sabbagh, and Y. Andreopoulos, “Physics-guided machine learning approaches to predict the ideal stability properties of fusion plasmas,” Nuclear Fusion 60, 046033 (2020).
  * Camporeale, Carè, and Borovsky (2017) E. Camporeale, A. Carè, and J. E. Borovsky, “Classification of solar wind with machine learning,” Journal of Geophysical Research: Space Physics 122, 10–910 (2017).
  * Jiao  _et al._ (2020) Z. Jiao, H. Sun, X. Wang, W. Manchester, T. Gombosi, A. Hero, and Y. Chen, “Solar flare intensity prediction with machine learning models,” Space Weather 18, e2020SW002440 (2020).
  * Chen  _et al._ (2019) Y. Chen, W. B. Manchester, A. O. Hero, G. Toth, B. DuFumier, T. Zhou, X. Wang, H. Zhu, Z. Sun, and T. I. Gombosi, “Identifying solar flare precursors using time series of sdo/hmi images and sharp parameters,” Space Weather 17, 1404–1426 (2019).
  * Humbird  _et al._ (2019) K. D. Humbird, J. L. Peterson, B. Spears, and R. McClarren, “Transfer learning to model inertial confinement fusion experiments,” IEEE Transactions on Plasma Science 48, 61–70 (2019).
  * Witman  _et al._ (2019) M. Witman, D. Gidon, D. B. Graves, B. Smit, and A. Mesbah, “Sim-to-real transfer reinforcement learning for control of thermal effects of an atmospheric pressure plasma jet,” Plasma Sources Science and Technology 28, 095019 (2019).
  * Kain  _et al._ (2020) V. Kain, S. Hirlander, B. Goddard, F. M. Velotti, G. Z. Della Porta, N. Bruchon, and G. Valentino, “Sample-efficient reinforcement learning for cern accelerator control,” Physical Review Accelerators and Beams 23, 124801 (2020).
  * Gonoskov  _et al._ (2019) A. Gonoskov, E. Wallin, A. Polovinkin, and I. Meyerov, “Employing machine learning for theory validation and identification of experimental conditions in laser-plasma physics,” Scientific reports 9, 1–15 (2019).
  * Li  _et al._ (2020) L.-N. Li, X.-F. Liu, W.-M. Xu, J.-Y. Wang, and R. Shu, “A laser-induced breakdown spectroscopy multi-component quantitative analytical method based on a deep convolutional neural network,” Spectrochimica Acta Part B: Atomic Spectroscopy , 105850 (2020).
  * Tajima and Dawson (1979) T. Tajima and J. M. Dawson, “Laser electron accelerator,” Physical Review Letters 43, 267 (1979).
  * Malka  _et al._ (2002) V. Malka, S. Fritzler, E. Lefebvre, M.-M. Aleonard, F. Burgy, J.-P. Chambaret, J.-F. Chemin, K. Krushelnick, G. Malka, S. Mangles, _et al._ , “Electron acceleration by a wake field forced by an intense ultrashort laser pulse,” Science 298, 1596–1600 (2002).
  * Mangles  _et al._ (2004) S. P. Mangles, C. Murphy, Z. Najmudin, A. G. R. Thomas, J. Collier, A. E. Dangor, E. Divall, P. Foster, J. Gallacher, C. Hooker, _et al._ , “Monoenergetic beams of relativistic electrons from intense laser–plasma interactions,” Nature 431, 535–538 (2004).
  * Geddes  _et al._ (2004) C. Geddes, C. Toth, J. Van Tilborg, E. Esarey, C. Schroeder, D. Bruhwiler, C. Nieter, J. Cary, and W. Leemans, “High-quality electron beams from a laser wakefield accelerator using plasma-channel guiding,” Nature 431, 538–541 (2004).
  * Faure  _et al._ (2004) J. Faure, Y. Glinec, A. Pukhov, S. Kiselev, S. Gordienko, E. Lefebvre, J.-P. Rousseau, F. Burgy, and V. Malka, “A laser–plasma accelerator producing monoenergetic electron beams,” Nature 431, 541–544 (2004).
  * Esarey, Schroeder, and Leemans (2009) E. Esarey, C. Schroeder, and W. Leemans, “Physics of laser-driven plasma-based electron accelerators,” Reviews of modern physics 81, 1229 (2009).
  * Gonsalves  _et al._ (2019) A. Gonsalves, K. Nakamura, J. Daniels, C. Benedetti, C. Pieronek, T. De Raadt, S. Steinke, J. Bin, S. Bulanov, J. Van Tilborg, _et al._ , “Petawatt laser guiding and electron beam acceleration to 8 gev in a laser-heated capillary discharge waveguide,” Physical review letters 122, 084801 (2019).
  * Roso (2018) L. Roso, “High repetition rate petawatt lasers,” in  _EPJ Web of Conferences_ , Vol. 167 (EDP Sciences, 2018) p. 01001.
  * Prencipe  _et al._ (2017) I. Prencipe, J. Fuchs, S. Pascarelli, D. Schumacher, R. Stephens, N. Alexander, R. Briggs, M. Büscher, M. Cernaianu, A. Choukourov, _et al._ , “Targets for high repetition rate laser facilities: needs, challenges and perspectives,” High Power Laser Science and Engineering 5 (2017).
  * Salehi (2019) F. Salehi,  _High Repetition Rate Laser-Driven Electron Acceleration to Mega-Electron-Volt Energies_ , Ph.D. thesis (2019).
  * Salehi  _et al._ (2017) F. Salehi, A. Goers, G. Hine, L. Feder, D. Kuk, B. Miao, D. Woodbury, K.-Y. Kim, and H. Milchberg, “Mev electron acceleration at 1 khz with< 10 mj laser pulses,” Optics letters 42, 215–218 (2017).
  * Feister, Poole, and Heuer (2019) S. Feister, P. L. Poole, and P. V. Heuer, “Development considerations for high-repetition-rate hedp experiments,” arXiv preprint arXiv:1906.11777 (2019).
  * Kim  _et al._ (2017) H. T. Kim, V. Pathak, K. H. Pae, A. Lifschitz, F. Sylla, J. H. Shin, C. Hojbota, S. K. Lee, J. H. Sung, H. W. Lee, _et al._ , “Stable multi-gev electron accelerator driven by waveform-controlled pw laser pulses,” Scientific reports 7, 1–8 (2017).
  * Liu  _et al._ (2014) C. Liu, J. Zhang, S. Chen, G. Golovin, S. Banerjee, B. Zhao, N. Powers, I. Ghebregziabher, and D. Umstadter, “Adaptive-feedback spectral-phase control for interactions with transform-limited ultrashort high-power laser pulses,” Optics letters 39, 80–83 (2014).
  * Tsai  _et al._ (2018) H.-E. Tsai, K. K. Swanson, S. K. Barber, R. Lehe, H.-S. Mao, D. E. Mittelberger, S. Steinke, K. Nakamura, J. van Tilborg, C. Schroeder, _et al._ , “Control of quasi-monoenergetic electron beams from laser-plasma accelerators with adjustable shock density profile,” Physics of Plasmas 25, 043107 (2018).
  * Shalloo  _et al._ (2020) R. Shalloo, S. Dann, J.-N. Gruse, C. Underwood, A. Antoine, C. Arran, M. Backhouse, C. Baird, M. Balcazar, N. Bourgeois, _et al._ , “Automation and control of laser wakefield accelerators using bayesian optimisation,” arXiv preprint arXiv:2007.14340 (2020).
  * Born and Wolf (2013) M. Born and E. Wolf,  _Principles of optics: electromagnetic theory of propagation, interference and diffraction of light_ (Elsevier, 2013).
  * Russell and Norvig (2002) S. Russell and P. Norvig, “Artificial intelligence: a modern approach,” (2002).
  * Liaw, Wiener  _et al._ (2002) A. Liaw, M. Wiener, _et al._ , “Classification and regression by randomforest,” R news 2, 18–22 (2002).
  * Weinberger (2018) K. Weinberger, “Machine learning for intelligent systems,” lecture notes, Cornell University, www.cs.cornell.edu/courses/cs4780/2018fa/ (2018).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, _et al._ , “Scikit-learn: Machine learning in python,” the Journal of machine Learning research 12, 2825–2830 (2011).
  * Goodfellow, Bengio, and Courville (2016) I. Goodfellow, Y. Bengio, and A. Courville,  _Deep learning_ (MIT press, 2016).
  * Hornik (1991) K. Hornik, “Approximation capabilities of multilayer feedforward networks,” Neural networks 4, 251–257 (1991).
  * Chollet  _et al._ (2015) F. Chollet _et al._ , “Keras,” https://keras.io (2015).
  * Abadi  _et al._ (2016) M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, _et al._ , “Tensorflow: Large-scale machine learning on heterogeneous distributed systems,” arXiv preprint arXiv:1603.04467 (2016).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980 (2014).
  * Humbird, Peterson, and McClarren (2018) K. D. Humbird, J. L. Peterson, and R. G. McClarren, “Deep neural network initialization with decision trees,” IEEE transactions on neural networks and learning systems 30, 1286–1295 (2018).
  * Rasmussen (2003) C. E. Rasmussen, “Gaussian processes in machine learning,” in  _Summer School on Machine Learning_ (Springer, 2003) pp. 63–71.
  * Molnar (2020) C. Molnar,  _Interpretable Machine Learning_ (Lulu. com, 2020).
  * Duvenaud (2014) D. Duvenaud,  _Automatic model construction with Gaussian processes_ , Ph.D. thesis, University of Cambridge (2014).


# A benchmark for 2D foetal brain ultrasound analysis

Mariano Cabezas  Brain and Mind Centre, University of Sydney, Sydney,
Australia  corresponding author(s): Mariano Cabezas
<EMAIL_ADDRESS> Yago Diez  Faculty Of Science, Yamagata
University, Yamagata, Japan  Clara Martinez-Diago  Hospital Universitari de
Girona Doctor Josep Trueta, Girona, Spain  Anna Maroto  Hospital Universitari
de Girona Doctor Josep Trueta, Girona, Spain

###### Abstract

Brain development involves a sequence of structural changes from early stages
of the embryo until several months after birth. Currently, ultrasound is the
established technique for screening due to its ability to acquire dynamic
images in real-time without radiation and to its cost-efficiency. However,
identifying abnormalities remains challenging due to the difficulty in
interpreting foetal brain images. In this work we present a set of 104 2D
foetal brain ultrasound images acquired during the 20th week of gestation that
have been co-registered to a common space from a rough skull segmentation. The
images are provided both on the original space and template space centred on
the ellipses of all the subjects. Furthermore, the images have been annotated
to highlight landmark points from structures of interest to analyse brain
development. Both the final atlas template with probabilistic maps and the
original images can be used to develop new segmentation techniques, test
registration approaches for foetal brain ultrasound, extend our work to
longitudinal datasets and to detect anomalies in new images.

## Background & Summary

Foetal brain development is a complex sequence of events ocurrying throughout
gestation. From early stages of embryonic development, the brain undergoes
structural changes until several months after birth [1, 2]. Therefore,
understanding normal brain development is essential to identify potential
deviations which may lead to neurological disability. Specifically, various
studies have described normal milestones within a specific chronology [3, 4,
5]. These changes can be observed by experts in prenatal diagnosis with skills
in foetal neurology using ultrasound (US) and magnetic resonance imaging
(MRI). However, identifying neurodevelopmental deviations is challenging due
to the difficulty in image interpretation [5].

Currently, US is the established technique for screening due to its ability to
acquire dynamic images in real-time without radiation and cost-efficiently.
The international guidelines recommend the acquisition of a routine first
trimester US within the 11th - 14th weeks and a mid-trimester scan within the
19th - 22th weeks of gestation for anatomical evaluation of the foetus [6].
The standard foetal assessment includes the evaluation of planes acquired
using 2D-US. Even though 3D-US is considered useful for prenatal diagnosis of
some disorders (mainly involving the face, the skeleton, the cardiovascular
system or the brain), the major obstacles for 3D-US implementation worldwide
as the main routine acquisition type are related to foetal motion artefacts
and acoustic shadowing during volume acquisition. Another widely used prenatal
imaging technique is MRI. However, MRI is not used as a primary screening tool
and is instead used as a complementary acquisition when foetal abnormalities
are suspected (even in selected high-risk cases) [7, 8]. Moreover, foetal MRI
performed before the 18th – 22nd weeks does not usually provide additional
information to that obtained by US. Generally, MRI provides a detailed
visualization of structures between the 26th and 32nd weeks, being superior to
US and less susceptible to limitations from maternal body conditions and
foetal presentation (bones do not produce occlusion artefacts in MRI) [8].
Within this context, a foetal brain US atlas based on 3D scans of healthy
foetuses was recently published [9]. This atlas demonstrates the feasibility
of assessing structural changes in the cortex and in the subcortical grey
matter by US. Furthermore, this chronological US-based atlas of the foetal
brain at each gestational age has the potential to become a useful tool to
detect development abnormalities when combined with the standard planes of
mid-trimester routine scans. However, considering that the primary diagnostic
tool in pregnancy remains 2D-US, studies aimed at improving the identification
of fetal structures should focus on 2D-US.

As there is a degree of variability in the shapes of anatomical regions
between individuals, atlases are typically built by taking into account a
number of images considered to fall within normal parameters. A crucial step
in this process is finding corresponding regions between different images and
warping them to a common space, a process known as image registration. In the
case of US atlases, the lower signal-to-noise ratio (SNR), the differences in
location for the foetuses due to the lack of a robust localization technique
and the presence of image artefacts that may blur or shadow salient features
and structure boundaries makes this process particularly challenging.

In this paper we present a set of foetal brain 2D-US images with manual
landmark annotations of structures of interest and soft probabilistic maps for
those structures based on these landmarks in a common space. The aim of this
dataset is to provide a set of US images as a starting point to study
registration and segmentation of foetal brain scans and provide tools for
researchers focusing on foetal brain development. Furthermore, in combination
with the recently published 4D atlas, this dataset can be also used to study
2D to 3D US registration and techniques to determine the “real age” of
gestation and potential abnormalities when comparing to the atlas. To define a
common space, a registration algorithm based on fitting an ellipse to the
skull segmentation of a convolutional neural network (CNN) was used to align
all the images to a common space. In order to avoid misalignements caused by
image quality, a subject ellipse was used to estimate an affine transformation
to a common (atlas) space. We also took advantage of the automatic
differentiation capabilities of the pytorch package in three fundamentally
different optimization settings (training a CNN for segmentation, fitting an
ellipse to a segmentation boundary, and estimating an affine transformation).

## Methods

In order to create the dataset presented, the data collected was processed as
follows: First, the skull was automatically segmented using a UNet network,
then an ellipse was automatically fitted to the shape of the skull. These
ellipses where used to register each image to a reference image. Soft
probabilistic maps where built using the set of registered images. The rest of
this section includes details on each of these steps and a summary of other
publicly available datasets. A visual summary can be found in Figure 1.

### Ultrasound Data

A prospective cohort of low-risk pregnant women was recruited at routine mid-
trimester foetal ultrasound scan. All participants initiated antenatal care
before the 12th weeks of gestation, underwent the first trimester ultrasound
scan between 12th and 14th weeks and had a low risk for aneuploidies in the
first trimester combined screening. Written informed consents were obtained
from participants. A private dataset of 70 pregnant women with a routine mid-
trimester foetal ultrasound scan at $\left[20\sim 20.6\right]$ weeks without
detected abnormalities was acquired, totalling 104 scans (8 women were scanned
three times, 18 women were scanned twice and the remaining 44 women were only
scanned once). The median of the maternal age was 31 (range 18-42). Images
were acquired using high-frequency transabdominal probe (C2-9) of Voluson E10
ultrasound system. For each subject, a transverse view of the foetal head
demonstrating a standard normative transcerebellar scanning plane was manually
annotated by 2 trained clinical experts (Figure 1 top right). Both experts
were part of the Prenatal Diagnosis Unit, one with over 10 years of experience
and the other with 5 years of experience. The experts jointly participated in
the annotations of each image to highlight common structures related to brain
development as follows:

  * •

Skull: 4 landmarks from the inner line of the skull were annotated: 2 at the
level of the middle line and the other 2 in a perpendicular imagined line at
the level of the posterior corners of the cavum septi pellucidi (CSP).

  * •

Thalami: 1 landmark marking the edges of both thalami at the middle line and
the outer edges of the concavity shape were annotated (3 total points).

  * •

Cerebellum: 8 landmarks on the perimeter the cerebellum were annotated.
Specifically, 2 points from the midline, 2 points from the cerebellum external
edges and 4 points in the middle of each cerebellum hemisphere.

  * •

Cavum septi pellucidi (Cavum): 4 points, each marking one corner of its
rectangular shape, were annotated.

  * •

Sylvian Fissure (Sylvius): 2 landmarks, one for each sylvian fissure edge, and
1 landmark in the inflection point of the fissure were annotated. For all the
images, only the inferior fissure was visible.

  * •

Midline: 1 landmark in the upper edge of the midline and 1 landmark in the
upper edge of the CSP were annotated.

An important aspect of acquisition is that the sylvian fissure was always
scanned on the lower part of the image, irrespective of the head orientation
(left to right or right to left). This phenomena has important implications
for registration. When aligning all images to a common 3D space (for example a
3D atlas template) in order to have all images facing the same direction, the
transformation can be modeled with 180$\degree$ rotations over the y axis. If
the transformation is limited to the 2D space (image coordinates), a mirroring
operation has to be applied to align all images. In the code we provide to
process the images, the second option is used for simplicity. Furthermore, our
common space is oriented from left to right (anterior to posterior).

Data acquisition was approved by the ethics committee “Comitè d’Ètica
d’Investigació amb Medicaments CEIM GIRONA” with the code 2023.067.
Furthermore, the subjects were informed and consented to the open publication
of the data.

### Skull segmentation database

The segmentation dataset used to train a skull segmentation CNN was downloaded
from the HC18 grand challenge on “automated measurement of foetal head
circumference using 2D ultrasound images” [10, 11]. The challenge comprises a
set of $800\times 540$-pixel 2D US images with a pixel sizes
$\in(0.052,0.326)$ mm. The data set was split into 999 images for training and
335 for testing. For each image in the training set, an ellipse was manually
fitted to the HC by a trained sonographer but precise segmentations of the
true skull boundary were not provided. No pre-processing techniques were
applied to the images apart from computing the z-score of the intensity values
with respect to the mean and standard deviation of the whole image (non-0
intensities) before feeding them to the segmentation network. The network
trained with all the training set images was then used to roughly segment the
skull of our own 2D US images and provide a starting point to fit an ellipse.

### Other public databases

To the best of our knowledge, only 3 other large datasets focusing on 2-D
ultrasound have been made publicly available, including HC18. As mentioned on
our brief description of HC18, the annotations are limited (only an ellipse
roughly representing the skull is given), the subjects present a large
distribution of gestational ages and no further information on the acquired
plan e is provided. For studies focusing on the segmentation of structures of
interest or registration, new annotations would be needed. Furthermore, the
gestational age and acquisition plane for each subject can have a large impact
on the appearance of the structures. Another public dataset presented by
Burgos-Artizzu et al [12] includes a large set of images (12400) ranging from
the 18th to the 40th week of pregnancy. Images were acquired using six
different machines and labelled at the image level. Similarly to HC18, the
dataset presents a large range of gestational ages and limited annotations
that could only be used to develop classification algorithms. Finally,
Alzubaidi et al [13] released a public dataset of 3832 high resolution images.
In contrast with the other two datasets, no mention of gestational age is
provided and once again rough annotations in the form of bounding boxes are
provided. Moreover, the authors also highlight image resampling as an
additional shortcoming of their dataset.

In our dataset, we focus on a specific acquisition plane and gestational age
as defined by international guidelines and raw images are provided.
Furthermore, landmarks for the most salient points of each structure of
interest are provided with software tools to estimate finer-grained masks and
bounding boxes around the landmarks. In that sense, our dataset provides a
useful tool to address multiple image analysis problems including registration
that could not be easily approached with other available datasets and no
additional annotations.

### Registration method

Medical image registration [14, 15] is a necessary initial step for many
medical image processing applications that rely on group-wise analysis.
Typically, images are registered in pairs: one of them is defined as the
“fixed” image (or reference) and the other as the “moving” image. The moving
image is then warped using a transformation function to generate the final
“moved” image. The transformation is commonly optimised using a predefined
metric that computes image similarity between the reference and the moved
image after transformation. Here, we present a coarse registration method to
roughly align different foetal ultrasound images as described in the following
sections and illustrated in Figure 1. We chose one of the images as the common
“reference” image (image 10) and registered the remaining images to it.

#### Automatic Skull segmentation using a Unet

The Unet architecture [16] is one of the most common CNN architectures for
image segmentation. Due to its encoder-decoder structure and the use of skip
connections, precise segmentations based on multi-scale features can be
obtained for a variety of segmentation applications [17]. Consequently, we
used a 2D Unet trained on the HC18 challenge dataset to provide a rough
segmentation of the skull for all the images in our registration dataset (see
Fig. 1 for an example). Specifically, the encoder and decoder blocks are
comprised of 6 convolutional layers with a residual connection [18] of 32, 32,
128, 256, 256 and 1024 features each (inverse order for the decoder) and a
bottleneck of 1024 features. To optimise the weights of the Unet, the Adam
algorithm with default initial learning rate was used to minimise the binary
cross-entropy loss.

#### Ellipse Registration

An ellipse is a planar curve representing the locus of the points with
constant added distances to two "focal points", as expressed by the quadratic
equation:

$f^{E}(X,Y)=AX^{2}+BXY+CY^{2}+DX+EY+F=0,(X,Y)\in\mathcal{R}^{2}$ (1)

With points that satisfy $f^{E}(X,Y)\neq 0$ being inside the ellipse perimeter
($f^{E}(X,Y)<0$) or outside of it ($f^{E}(X,Y)>0$). The general equation’s
coefficients can be obtained from known semi-major axis $a$ (represented by
the magnitude of the turquoise vector in Fig. 1), semi-minor axis $b$
(represented by the magnitude of the cyan vector in Fig. 1), centre
coordinates $(x_{0},y_{0})$ (represented by the point where the two vectors
meet in Fig. 1) and rotation angle $\theta$ (the angle from the positive
horizontal axis to the ellipse’s major axis as observed in Fig. 1) using the
formulae:

$A=a^{2}\sin^{2}\theta+b^{2}\cos^{2}\theta$ (2)
$B=2(b^{2}-a^{2})\sin\theta\cos\theta$ (3)
$C=a^{2}\cos^{2}\theta+b\sin^{2}\theta$ (4) $D=-2Ax_{0}-By_{0}$ (5)
$E=-Bx_{0}-2Cy_{0}$ (6) $F=Ax^{2}_{0}+Bx_{0}y_{0}+Cy^{2}-a^{2}b^{2}$ (7)

These expressions can be derived from the canonical equation
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1$ by an affine transformation of the
coordinates $(x,y)$ (with a translation $(-x_{0},-y_{0})$ and an angle
$\theta$).

$x=(X-x_{0})\cos\theta+(Y-y_{0})\sin\theta$ (8)
$y=-(X-x_{0})\sin\theta+(Y-y_{0})\cos\theta$ (9)

Given the set of pixel coordinates of the skull segmentation
($(X_{skull},Y_{skull})=\left[(x_{1},y_{1}),\cdots,(x_{N},y_{N})\right]$) we
can fit an ellipse using its parameters ($a$, $b$, $x_{0}$, $y_{0}$ and
$theta$) by minimising the following objective function:

$\hat{a},\hat{b},\hat{x_{0}},\hat{y_{0}},\hat{\theta}=\operatorname*{arg\,min}_{a,b,x_{0},y_{0},\theta}\sum_{i=1}^{N}f^{E}(x_{i},y_{i})^{2},$
(10)

where the coefficients A to F are substituted in $f^{E}$ by their definitions
in equations 2 to 7 and the estimated ellipse parameters.

This process is repeated 5 times, removing erroneous points of the skull
segmentation mask that are one standard deviation away from the mean ellipse
error according to Eq. 10. With the parameters of the skull ellipse estimated
through optimisation, we can now define an affine transformation (referred to
as Ellipse from now on) to move the brain to the centre of the image as:

$\alpha_{e}=\begin{bmatrix}\frac{w}{2a}\cos\theta&\frac{w}{2a}\sin\theta&-w/2\\\
-\frac{h}{2b}\sin\theta&\frac{h}{2b}\cos\theta&-h/2\\\
0&0&1\end{bmatrix}\times\begin{bmatrix}1&0&-x_{0}\\\ 0&1&-y_{0}\\\
0&0&1\end{bmatrix}$ (11)

#### Affine image registration

For comparison, a regular rigid registration of 6 unrestricted parameters was
performed with different initialisations. From the most basic identity
initialisation (referred to as Affine), to an initialisation using the ellipse
parameters of the reference image (referred to as Affine (Reference ellipse))
and a refinement of the ellipse-based affine transformation from Eq. 11
(referred to as Ellipse + Affine).

#### Probabilistic maps

Once the images are co-registered to a common space based on their ellipse, a
two-step process is performed to generate probabilistic maps for all the
structures with more than 2 landmarks.

First, the concave hull of each structure is computed using the alphashapes
package. This concave hull represents a rough segmentation of the location of
the structure. Second, the segmentations for all the subjects are averaged per
structure to determine the probability of each pixel to belong to that
structure.

This approach has some limitations. Namely, some structures have a polygonal
shape, even though they are actually curves (e.g. sylvius) and the final masks
are only a rough representation of the real boundaries (e.g. cerebellum).
However, these polygonal maps can be still used to determine growth milestones
and to provide a rough location of the structure of interest.

## Data Records

The original images with manual landmark annotations (Gimp image editing tool
format) and the co-registered images and probabilistic maps used within this
paper can be found on figshare [19], and are organised with subject id number
(1-52) and an additional number for multiple scans (e.g. 36 and 36.1). Co-
registered images are saved in JPEG format with the “_registered” suffix. The
final probability maps estimated using a combination of co-registered
landmarks and the alphashape package are provided as JPEG images with the name
of the structure (e.g. sylvius.jpeg), while comma separated value (CSV) files
with the point landmarks of all registered subjects are compressed into a
single zip file.

In total, the released dataset consists of 104 annotated 2D US images of
foetal brains on the 20th week of pregnancy. The manual annotations are
described in Ultrasound dataset section.

## Technical Validation

To validate the techniques used for co-registration to a common space, we
focused on common medical imaging metrics for registration using landmarks. In
this section we describe these metrics and provide some qualitative and
quantitative results of the alignment (including a visualisation of the
probabilistic maps of one of the structures of interest).

Regarding the quality of the US images, all images were acquired with a high
frequency ultrasound probe (2-9 MHz) that provides high resolution images
following the ISUOG recommendations. The most constraining factor for the
quality of the ultrasound images is maternal obesity. However, for this study
we excluded women with maternal morbid obesity (body mass index > 40) as it is
one of the factors for high-risk pregnancy. To further illustrate that point,
we provide a comparative example between an image from the dataset and a lower
quality one in figure 2(b).

### Evaluation metrics

For this study we chose to use the anatomical knowledge defined by expert
annotations as the main directive to evaluate the quality of the registration
and avoid focusing exclusively on common pixel-metrics that might be
unreliable and disconnected from physical properties [20]. Specifically, we
used the point annotations described in the Ultrasound Data section with two
point-based metrics and one area-based metric for every anatomical structure
with more than 2 points (all the structures except the midline). For
completeness sake, we also included the structural similarity index metric
(SSIM) pixel-based metric.

Point-based metrics In order to penalise partial matches between anatomical
structures, we considered the Hausdorff distance ($d_{H}$), that computes the
worst possible Euclidean distance ($d(p_{i},p_{j})$) between two sets of
points $P_{f}$ ($|P_{f}|=N$) and $P_{m}$ ($|P_{m}|=M$). We also considered the
average of the minimum Euclidean distances ($d_{E}$) to express global
similarity between structures defined by landmarks.

Area-based metric To evaluate the superposition between two anatomical
structures defined as 2D landmarks (points), we first computed their concave
hulls and then considered their Dice similarity coefficient (polygon DSC).

Image similarity metric For completeness, we also considered the structural
similarity index measure (SSIM) as a pixel-intensity-based metric.

### Comparison of coarse registration approaches

Figures 3, 4 and 5 summarise the results with boxplots and Wilcoxon signed-
rank tests for all the registration methods and metrics considered. Wherever
possible, the results for different anatomical areas are presented separately.
For the Euclidean and Hausdorff metrics, lower values indicate better
registration, while for the SSIM and Dice metric higher results indicate
better registration. Moreover, a qualitative example to illustrate
misalignments between the reference points and the registered landmarks is
provided in Figure 6.

Regarding point metrics, the ellipse (E) and ellipse with affine methods (E+A)
obtain overall better results than the other methods. In general, the
differences observed were found to be statistically significant for both point
metrics and most anatomical structures. Exceptions to this are the thalami and
the cerebellum where the pixel based affine registration method initialised
using the reference ellipse (AFF+I), obtains results that appear worse but are
not significantly different. Comparing the E and E+A methods, small (and not
statistically significant) differences can be observed. Using a refinement
registration after ellipse-based method slightly worsens the metrics when
applied to the skull but improves them in all other anatomical structures.
This is an expected result as the main focus of the first method is to co-
register the skulls (ellipses). The uninitialised affine transformation (AFF)
obtains significantly worse results than these two methods and has a higher
variance of values as illustrated by its wide boxplots. Both methods using
exclusively pixel-wise affine registration (AFF, AFF+I) achieve results that
are worse than the metrics of the original moving image in some cases. This
illustrates the disconnect between pixel-based metrics guiding these methods
and point-based metrics targeting distances between real anatomical structures
(especially for noisy sequences).

Regarding the differences between the Euclidean and Hausdorff metrics
(especially the median values shown by the central line in boxplots), slightly
higher Hausdorff values indicate that some point pairs are further than the
average euclidean distance mean value. This is especially noticeable on the
skull for E and A+E and the thalami for AFF.

The other two metrics show the same general tendencies, even though they focus
on different aspects of the registration. The low SSIM values observed for all
methods (even though E and A+E outperform all methods) illustrats how
challenging this registration scenario is. The variations between individuals
and acquisitions and low SNR make the intensity values particularly
unreliable. On the other hand, the high anatomical DSC results observed for
the ellipse-based methods (includding AFF+I) validate our geometric approach
that relies on a rough initial skull estimate to determine the general shape
and orientation of the brain. Similarly high results obtained for the
cerebellum indicate that a correct skull placement is a crucial important step
towards the registration of all brain structures.

Finally, Figure 7 shows the final heatmap generated from the registered
concave hulls (E+A) of the cerebellum for all output images. Even though faint
outlines of some incorrect registration results outside of the cerebellum
region (delimited by orange points) can be observed, the higher probability
regions in the heatmap clearly correspond to the cerebellum region of the
reference image delimited by the manually annotated orange landmarks.

## Usage Notes

The code provided to analyse the images and perform a rough initial alignment
to a common space has been developed using python. A set of Jupyter notebooks
detailing the use of each step is also provided in the repository with
visualisation examples of each step of the processing pipeline. Furthermore,
we provide the trained weights for the skull segmentation network as part of
the data records (file unet.pt [19]).

Regarding the use of python packages, the code heavily relies on pytorch
(version 1.12.0 with CUDA 10.2) to do the heavy lifting. However, numpy
(version 1.21.6), scipy (version 1.7.3), scikit-learn (version 1.0.2), scikit-
image (version 0.18.1), pandas (version 1.3.4), persim (version 0.3.1) and
alphashape (version 1.3.1) are also used for some of the processing steps or
to analyse different metrics related to the registration. Furthermore, opencv-
python (version 4.5.2.52) and gimpformats (version 2021.1.4) where used to
open the images with python. In particular, manual landmark annotation where
done using the Gimp software and were stored (together with the ultrasound
image) using the xcf format. Finally, matplotlib (version 3.7.1), seaborn
(version 0.11.0) and statannot (version 0.2.3) were used for results
visualisation inside the Jupyter notebooks.

## Code availability

All the code used in the study to generate the final atlas and to co-register
the images is publicly available at https://github.com/marianocabezas/fetal-
brain.

## References

  * [1] Monteagudo, A. & Timor-Tritsch, I. E.  Normal sonographic development of the central nervous system from the second trimester onwards using 2d, 3d and transvaginal sonography.  _Prenatal Diagnosis_ 29, 326–339 (2009).
  * [2] Namburete, A. I. _et al._ Learning-based prediction of gestational age from ultrasound images of the fetal brain.  _Medical Image Analysis_ 21, 72–86 (2015).
  * [3] Cohen-Sacher, B., Lerman-Sagie, T., Lev, D. & Malinger, G.  Sonographic developmental milestones of the fetal cerebral cortex: a longitudinal study.  _Ultrasound in Obstetrics & Gynecology_ 27, 494–502 (2006).
  * [4] Kline-Fath, B. M.  Ultrasound and mr imaging of the normal fetal brain.  _Neuroimaging Clinics of North America_ 29, 339–356 (2019).
  * [5] Guimaraes, C. V. & Dahmoush, H. M.  Fetal brain anatomy.  _Neuroimaging Clinics of North America_ 32, 663–681 (2022).
  * [6] Salomon, L. J. _et al._ Isuog practice guidelines (updated): performance of the routine mid-trimester fetal ultrasound scan.  _Ultrasound in Obstetrics & Gynecology_ 59, 840–856 (2022).
  * [7] Griffiths, P. D. _et al._ Mri in the diagnosis of fetal developmental brain abnormalities: the meridian diagnostic accuracy study.  _Health Technology Assessment (Winchester, England)_ 23, 1 (2019).
  * [8] Prayer, D. _et al._ ISUOG practice guidelines (updated): performance of fetal magnetic resonance imaging.  _Ultrasound in Obstetrics & Gynecology_ 61, 278–287 (2023).
  * [9] Namburete, A. I. _et al._ Normative spatiotemporal fetal brain maturation with satisfactory development at 2 years.  _Nature_ 1–9 (2023).
  * [10] van den Heuvel, T. L. A., de Bruijn, D., de Korte, C. L. & van Ginneken, B.  Automated measurement of fetal head circumference using 2d ultrasound images.  _PloS one_ 13, e0200412 (2018).
  * [11] van den Heuvel, T. L. A., de Bruijn, D., de Korte, C. L. & van Ginneken, B.  Automated measurement of fetal head circumference using 2d ultrasound images [data set].  Zenodo, http://doi.org/10.5281/zenodo.1322001 (2024).
  * [12] Burgos-Artizzu, X. P. _et al._ Evaluation of deep convolutional neural networks for automatic classification of common maternal fetal ultrasound planes.  _Scientific Reports_ 10, 10200 (2020).
  * [13] Alzubaidi, M. _et al._ Large-scale annotation dataset for fetal head biometry in ultrasound images.  _Data in Brief_ 51, 109708 (2023).
  * [14] Fu, Y. _et al._ Deep learning in medical image registration: a review.  _Physics in Medicine & Biology_ 65, 20TR01 (2020).
  * [15] Gholipour, A., Kehtarnavaz, N., Briggs, R., Devous, M. & Gopinath, K.  Brain functional localization: A survey of image registration techniques.  _IEEE Transactions on Medical Imaging_ 26, 427–451, 10.1109/TMI.2007.892508 (2007).
  * [16] Ronneberger, O., Fischer, P. & Brox, T.  U-net: Convolutional networks for biomedical image segmentation.  In _Proceedings of the International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI)_ , 234–241 (2015).
  * [17] Isensee, F., Jaeger, P. F., Kohl, S. A. A. & Maier-Hein, K. H.  nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation.  _Nature: Methods_ (2021).
  * [18] He, K., Zhang, X., Ren, S. & Sun, J.  Deep residual learning for image recognition.  In _IEEE Conf. Comput. Vision Pattern Recog._ , 770 – 778 (2016).
  * [19] Cabezas, M., Diez, Y., Martinez-Diago, C. & Maroto, A.  A benchmark for 2d foetal brain ultrasound analysis [data set].  Figshare, https://doi.org/XXXX (2024).
  * [20] Rohlfing, T.  Image similarity and tissue overlaps as surrogates for image registration accuracy: Widely used but unreliable.  _IEEE Transactions on Medical Imaging_ 31, 153–163 (2012).

## Author contributions statement

All authors conceived and designed the experiments, C.M. and A.M collected the
images and performed the landmark annotations, M.C. and Y.D. conducted the
experiments and analysed the results. All authors reviewed the manuscript.

## Competing interests

The authors declare no competing interests.

## Figures & Tables

---

(a)

(b)

Figure 1: Methodology: 1) The skull is automatically segmented using a UNet
network (top-left). 2) An ellipse is fitted to the skull segmentation (top-
right) 3) to estimate an affine transformation to a reference image (bottom-
right). The axes of the fitted ellipse are warped to the image coordinate axes
and are re-scaled to fit the ellipse in the reference image (bottom-left).

(a) Lower quality image

(b) Image from the dataset

Figure 2: Qualitative comparison between a low quality image where structures
are not clearly visible and and image from the dataset.

|
---|---
|
|

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Quantitative results for all the methods on each structure
(Hausdorff distance, $d_{H}$, lower values indicate better registration). The
upper part of each boxplot figure indicates the results of pairwise
statistical Wilcoxon tests: (ns: 5.00e-02 $<$ p $\leq$ 1.00e+00, *: 1.00e-02
$<$ p $\leq$ 5.00e-02, **: 1.00e-03 $<$ p $\leq$ 1.00e-02, ***: 1.00e-04 $<$ p
$\leq$ 1.00e-03, ****: p $\leq$ 1.00e-04).

|
---|---
|

(a)

(b)

(c)

(d)

Figure 4: Quantitative results for all the methods on each structure with
multiple points (polygon DSC, higher values indicate better registration). The
upper part of each boxplot figure indicates the results of pairwise
statistical Wilcoxon tests: (ns: 5.00e-02 $<$ p $\leq$ 1.00e+00, *: 1.00e-02
$<$ p $\leq$ 5.00e-02, **: 1.00e-03 $<$ p $\leq$ 1.00e-02, ***: 1.00e-04 $<$ p
$\leq$ 1.00e-03, ****: p $\leq$ 1.00e-04). Figure 5: Quantitative results for
all the methods according to the SSIM metric (higher values indicate better
registration). The upper part of each boxplot figure indicates the results of
pairwise statistical Wilcoxon tests: (ns: 5.00e-02 $<$ p $\leq$ 1.00e+00, *:
1.00e-02 $<$ p $\leq$ 5.00e-02, **: 1.00e-03 $<$ p $\leq$ 1.00e-02, ***:
1.00e-04 $<$ p $\leq$ 1.00e-03, ****: p $\leq$ 1.00e-04).

|
---|---
|
|

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Qualitative example of the alignment between the reference image and
a randomly selected subject (11). The background image corresponds to the
warped image after registration, circles denote the reference points for each
structure, crosses represent the registered points and lines are used to
illustrate misalignment between the reference and the registration. Figure 7:
Probabilistic map for the cerebellum structure based on the averaged concave
hull of the landmarks for each image. Landmarks for the reference image
(background) are also provided.