# Multiscale simulation of spatially correlated microstructure via a latent
space representation

Reese E. Jones
Sandia National Laboratories,
Livermore, CA USA &Craig M. Hamel
Sandia National Laboratories,
Albuquerque, NM USA &Dan Bolintineanu
Sandia National Laboratories,
Albuquerque, NM USA &Kyle Johnson
Sandia National Laboratories,
Albuquerque, NM USA &Robert Buarque de Macedo
Sandia National Laboratories,
Albuquerque, NM USA &Jan Fuhg
The University of Texas at Austin,
Austin, TX USA &Nikolaos Bouklas
Cornell University,
Ithaca, NY USA &Sharlotte Kramer,
Sandia National Laboratories,
Albuquerque, NM USA corresponding<EMAIL_ADDRESS>
###### Abstract

When deformation gradients act on the scale of the microstructure of a part
due to geometry and loading, spatial correlations and finite-size effects in
simulation cells cannot be neglected. We propose a multiscale method that
accounts for these effects using a variational autoencoder to encode the
structure-property map of the stochastic volume elements making up the
statistical description of the part. In this paradigm the autoencoder can be
used to directly encode the microstructure or, alternatively, its latent space
can be sampled to provide likely realizations. We demonstrate the method on
three examples using the common additively manufactured material AlSi10Mg in:
(a) a comparison with direct numerical simulation of the part microstructure,
(b) a push forward of microstructural uncertainty to performance quantities of
interest, and (c) a simulation of functional gradation of a part with
stochastic microstructure.

##  1 Introduction

Materials with strongly heterogeneous microstructure are becoming more common
with the rise of additively manufactured (AM) metals [1] and the advent of
exotic materials such as metamaterials [2] and microarchitectured materials
[3]. When the scales of stress and deformation are comparable to that of the
microstructure through the interplay of loading and gross geometry, common
simplifying assumptions such as _scale separation_ of the deformation field
and the underlying microstructure fail to hold. Hence many multiscale
simulation methods [4, 5] that account for microstructure but rely on scale
separation become ineffective. Likewise, classical micromechanics theory is
mostly concerned with the limits of averages/mean properties [6], although
finite-size effects and spatial correlations need to be accounted for. This
poses a challenge for the efficient design of structures that accounts for
microstructure [7] since direct numerical simulation (DNS) of the
microstructure of an entire part is typically infeasible. Furthermore, only
the spatial statistics of the microstructure are usually known, not the exact
microstructural configuration, which requires a design process based on
sampling multiple realizations of the microstructure to account for this
uncertainty.

To address this challenge we propose a multiscale simulation method that
incorporates the spatial correlations of the underlying microstructure and,
thereby, some of the finite-size effects that complicate homogenization in
this scenario. The goal is to enable the design of components with interacting
gross and microstructural scales under microstructural uncertainty by
accelerating the forward propagation of this uncertainty. The method uses a
variational autoencoder (VAE)-like neural network (NN) [8] to encode
microstructure into a compact latent vector and thereby to effective
properties of the microstructural sample. These properties, in turn, determine
the response of finite elements representing the neighboring, correlated high-
dimensional microstructural samples. In effect the NN provides the structure-
property relationship for the material and, through an intermediate structure-
to-latent encoding, a compact representation of the microstructural features
relevant to property prediction. By sampling the latent distribution of
statistically representative samples, the method can be used generatively for
uncertainty quantification (UQ); alternatively, if the full microstructure of
a component is available, the VAE can be used to encode it into fields of
corresponding latent descriptors and properties. Our use of interpretable
properties is expedient since it allows the incorporation of the method into
large-scale simulators with minimal modification.

In the next section, we provide background to give context for the proposed
method. Then, in Sec. 3, we provide the problem setting. Sec. 4 describes
the methods we use to create statistically consistent realizations. In Sec. 5
we give details of the training data generation, and in Sec. 6 we describe
the NN architecture that provides the engine for the proposed multiscale
method. Then in Sec. 7 we provide three demonstrations motivated by laser-
powder-bed-fusion additive manufacturing (L-PBF AM) with AlSi10Mg [9, 10]: (a)
a comparison of the multiscale response to that of a corresponding DNS, (b)
forward propagation of microstructural uncertainty for a case where
deformation gradients interact directly with microstructure, and (c) an
illustration of the use of the method to simulate random materials with
functional gradation. Lastly, Sec. 8 provides a summary of conclusions and
avenues for future work.

##  2 Related work

This work builds upon developments in microstructural characterization,
synthetic realization generation, homogenization, and related fields to
deliver an NN-powered multiscale simulation method that embeds correlated
microstructural features. Each of these fields has a long history of
development. Here we give a brief topical review of developments relevant to
the present work.

Classical microstructural characterization can be seen as a form of encoding
high-dimensional structures into a few interpretable statistics. For instance,
Kalidindi and coworkers have been particularly productive in the use of
2-point statistics to describe the spatial correlations of experimentally
observed microstructures [11, 12, 13, 14, 15]. These spatial statistics
provide the basis for numerous structure generation techniques such as
Gaussian random fields (GRFs) [16, 17, 18] and Karhunen-Loève expansions
(KLEs) [19, 20, 21, 22, 23] through spectral connections via the Wiener-
Khinchin and related theorems. General latent encoding of high dimensional
inputs is also a common task in machine learning. The Variational auto-encoder
(VAE) [24, 25, 26] has seen wide-spread use in this task. As an autoencoder,
the VAE reduces a high dimensional input, such as an image, into a compact
latent space, and then decodes this latent encoding to reproduce the input
using an hourglass-like stack of NN layers. A VAE is variational in the sense
that parameters of approximate latent distribution are optimized to match the
training data. Typically a multivariate Gaussian is used to represent the
latent distribution but others, such as the so-called _concrete_ distribution
[27], may be more appropriate for encoding binary images/structures [28, 29].
Chen and co-workers [8] invented a property-VAE (pVAE) which augments a
standard VAE with a regressor devoted to the secondary task of predicting the
properties of the input structure. They showed that the VAE has some facility
in organizing the latent space so that there is a smooth gradient of
properties. The secondary task also emphasized the fact that the VAE is
seeking a multi-objective solution [30] in terms of input reconstruction loss,
the divergence between the surrogate latent distribution and the data
distribution, and (c) the regression error of the predicted properties, which
need to be balanced in a single objective.

The homogenization of the properties and response of heterogeneous random
media is a mature field of research [6, 31, 32] and provides the foundation
for multiscale simulation. The Hill-Mandel [33, 34] condition, postulating the
equivalence of power expended at the micro and macro scale plays a central
role in many multiscale methods. Although generally concerned with ensemble
averages and large-sample limits, the finite-size effects/variability inherent
in small samples [35, 36] as well as the correlation structures [37] are also
themes in this discipline. Computational multiscale methods [38, 5, 39]
typically ignore these effects and appeal to scale separation of the
microstructure and the other length scales in the problem of interest. In
fact, the well-cited FE2 method [5, 4] is predicated on scales being well
separated. However, recently Chen, Guilleminot, and Soize [40] published a
treatment of random elastic media where spatial coherence of the realizations
was incorporated. The earlier work of Girolami et al. [41] also provided a
probabilistic learning framework applied to material modeling. A number of
recent reviews by Peng et al. [42], Fish et al. [43], and Bishara et al. [44]
provide a more comprehensive view of the application of machine learning to
the multiscale problem.

The treatment of the material uncertainty due to microstructure variability
[45, 46] and design under this uncertainty [47] are closely related tasks to
multiscale modeling [45, 46], and they have recently seen the rise of purely
data-driven approaches [48, 49, 50, 51], As a counterpoint, machine-learning
discovery of models of material response [52, 53, 54, 55] can provide data
reduction of a complex material response as a function of microstructure [56,
57, 58]. The data reduction of the response of microstructural samples can be
taken to the point of obtaining material properties from the samples (aka
parameters of chosen models) [59]. Each of these approaches involving encoding
the structure-response map have clear advantages and disadvantages in terms of
generality, computational expense and ease of implementation.

##  3 Problem

In pursuit of an efficient, acceptably accurate solution to a boundary value
problem with pronounced heterogeneity, concurrent multiscale (MS) methods
typically condition a field of small-scale problems on a coarse kinematic
field and use these responses to drive the solution of the overall problem.
Although there may be many (ill-defined) scales of heterogeneity, the scale-
separation assumption of traditional multiscale methods allows these small-
scale boundary value problems to represent independent samples of the
microstructure. In this work, we provide an approximation for the case where
the range of influence of the scales are overlapping by accounting for the
spatial correlation of the underlying structure. This case is prevalent when
there is interplay between the deformation-gradient field and characteristic
structural scales, as when the loading and/or the nominal geometry produce
deformation gradients that change over a scale comparable to that of the
microstructure [60]. In these scenarios, traditional MS will fail to be
predictive, while the proposed approach should bring the multiscale solution
closer to that of the DNS treatment that resolves the microstructure.

Since microstructure cannot be controlled precisely, it is typically
characterized and modeled statistically. Although many microstructures appear
random at some scales, such as the micrograph of AlSi10Mg shown in Fig. 1,
they typically have a spatial correlation that can be measured with 2-point
statistics [14]. These spatial statistics are a form of latent encoding that
we will exploit in synthetic structure generation described in Sec. 4.

We will assume that a well-defined set of properties þ exists at the smallest
relevant structural scale. This is plausible for observable structural
properties, such as phase concentration, while mechanical properties are
harder to probe directly at small scales but are clearly dependent on these
observables. For phase-separated mixtures and alloys, simple mixture rules
give an expectation of derived quantities. If the contributions are
independent, such as for mass density, these rules are exact. For other, more
complex, deformation-dependent quantities we need a model to predict the
properties that control the response of the microstructure.

Unlike in an homogeneous material, for a body $\Omega$ with microstructure,
the property  þ

${\text{ \th}}_{V}\equiv\frac{1}{V}\int_{\square}{\text{
\th}}(\mathbf{x})\,\mathrm{d}V$ (1)

of a finite-sized sample $\square$ with volume $V$ can vary by location. A
_stochastic volume element_ (SVE) is the common term for a microstructural
sample that is smaller than the limit

$\bar{{\text{ \th}}}\equiv\mathbb{E}[{\text{
\th}}_{V}]=\lim_{V\to\infty}{\text{ \th}}_{V}$ (2)

at which the particular property reaches its mean value $\bar{{\text{ \th}}}$.
The smallest sample that attains this limit, to a given accuracy, is called a
_representative volume element_ (RVE). For SVEs, there can be significant
variance in the properties:

$\mathbb{V}[{\text{ \th}}_{V}]=\frac{1}{V_{\Omega}}\int_{\Omega}({\text{
\th}}_{V}(\mathbf{x})-\bar{{\text{
\th}}})^{2}\,\mathrm{d}V\propto\frac{1}{V}\gg 0$ (3)

We expect the variance in properties with sample size to follow a central
limit trend, in the sense that it is roughly inversely proportional to sample
size since each sample property is effectively an average of contributions
throughout the sample. Fig. 2 shows the distribution of binary phase fraction
and spatial correlations for a sequence of sample sizes associated with the
material shown in Fig. 1. Clearly, the distribution narrows and the
correlation length increases with increasing sample size.

Furthermore, the integral over the region $\square$ in Eq. (1) can be
considered a filter of an underlying random process ${\text{
\th}}(\mathbf{x})$ in the sense that ${\text{ \th}}_{V}$ is the convolution of
the field ${\text{ \th}}(\mathbf{x})$ at the smallest scale with top hat
kernel with extent given by the SVE size $\square$. Hence for a coherent
structure, such as that shown in Fig. 1, the SVE properties at any location
are clearly correlated with the neighboring values. The 2-point correlation

$\mathbb{C}[{\text{ \th}}(\mathbf{x})]=\mathbb{E}[{\text{
\th}}_{V}(\mathbf{0}){\text{
\th}}_{V}(\mathbf{\mathbf{x}})]=\frac{1}{V_{\Omega}}\int_{\Omega}{\text{
\th}}_{V}(\mathbf{y}){\text{
\th}}_{V}(\mathbf{y}+\mathbf{x})\mathrm{d}\mathbf{y}$ (4)

is a common statistic used to characterize this dependence. Typically
correlations decay with distance, with the notable exception of periodic
structures, and hence have a characteristic length. Also, it stands to reason
that, as the SVE property estimate approaches the RVE limit, the variance of
the estimate will decrease, and the correlation length will increase. In this
work, we assume $\mathbb{C}[{\text{ \th}}(\mathbf{x})]=\mathbb{C}[{\text{
\th}}(\|\mathbf{x}\|)]$ to be isotropic for simplicity in the following
developments. Note that the variance in Eq. (3) $\mathbb{V}[{\text{
\th}}]=\mathbb{C}[{\text{ \th}}(\mathbf{0})]$ i.e. the expected variance of a
property in a sample is given by the self-correlation. Fig. 2 shows the
spatial correlation of the average phase of SVEs of increasing size
corresponding to the microstructure shown in Fig. 1. The rapidly decaying
primary correlation is associated with the average feature size, and the dip
in correlation is likely due to one phase being anti-correlated with the
other, while the longer decaying trend indicates a larger scale order. More
generally this figure illustrates how: (a) the variance of a mean property
decreases with increasing sample size and (b) the spatial correlation of this
mean property increases with increasing sample size. The limit being the RVE
value of the property that is equivalent everywhere and hence has a
correlation length as large as the body.

If the mean (2) and covariance (4) do not depend on location in the larger
body $\Omega$, these statistics are sufficient to describe a stationary
second-order process. More exposition of spatial statistics can be found in
Ref. [61] for instance. These two statistics are generally insufficient to
determine the mechanical properties, but all aspects of the high-dimensional
structural configuration space are not needed to determine the necessary
properties. In Sec. 6 we will describe how we can reduce the structure-
property relationship to a compact latent space.

Figure 1: AlSi10Mg microstructure black: Al, white: Si courtesy of Sara
Dickens and Paul Specht, Sandia National Laboratories [9].

Figure 2: Volume fraction distribution and correlation of AlSi10Mg
microstructure

##  4 Microstructure generation

Since, in applications ranging from component design to part forensics, we
typically only have statistical characterization available, we need structure
generation methods. For our purposes, we need to generate synthetic
realizations with the measured spatial statistics in order to push forward the
variability to performance metrics. There are a multitude of methods to
generate microstructures that obey given spatial statistics. We employ two
methods, in different capacities, which generate fields with a given spatial
correlation and point-wise expectation: (a) a Gaussian random field (GRF) and
(b) a Karhunen-Loève expansion (KLE). We use a GRF to generate SVEs and
structures for direct numerical simulation (DNS), and we use KLE to impart
spatial correlation of latent vectors on meshes used to perform multiscale
simulation. Each method relies on a correlation function that we assume, for
simplicity, is isotropic. For structure generation we employ a correlation
function representation with multiple decay / correlation lengths

$\kappa(\|\mathbf{x}-\mathbf{x}^{\prime}\|)=\sum_{i}w_{i}\exp\left(-\frac{\|\mathbf{x}-\mathbf{x}^{\prime}\|^{p}}{\ell_{i}^{p}}\right)$
(5)

with $\sum_{i}w_{i}=1$ and $p>0$.

A Gaussian random field [17, 18] is a random field with given spatial
correlation, and can be readily constructed on a structured grid using fast
Fourier transforms. In essence this involves the convolution of a noise
process $\eta$ without spatial dependence with a kernel $\kappa$.

$\varphi(\mathbf{x})=\mathcal{F}^{-1}_{\mathbf{k}\to\mathbf{x}}\left[\mathcal{F}_{\mathbf{x}\to\mathbf{k}}[\eta]\sqrt{\left|\mathcal{F}_{\mathbf{x}\to\mathbf{k}}[\kappa]\right|}\right]$
(6)

where $\mathcal{F}$ is a spatial Fourier transform, and
$\left|\mathcal{F}_{\mathbf{x}\to\mathbf{k}}[\kappa]\right|$ is the spectral
density of the spatial correlation. Typically $\eta$ is Gaussian white noise,
and $\kappa$ is a correlation kernel resembling Eq. (5). From the form of Eq.
(6) it is clear that a GRF is defined by the power spectral density of the
spatial correlation, and the statistics of the point process are given by
$\eta$. With this Fourier transform construction, a GRF is particularly suited
for generating periodic microstructures. Although more sophisticated methods
exist [62, 63, 64], Eq. (5) is sufficient for this exposition.

The Karhunen-Loève expansion (KLE) is another method that uses kernel
$\kappa$, which is more adaptable to meshes but also more computationally
expensive than GRFs. We will use a KLE to sample the latent space generated by
the pVAE described in Sec. 6. A concise exposition of the KLE applied to this
task will be given in Sec. 7.2.

Since these methods generate continuous fields, to create binary structures, a
threshold $\varepsilon$ that depends on a target volume fraction $\upsilon$ is
typically employed

$\phi(\mathbf{x})=\begin{cases}0&\text{if}\ \varphi(\mathbf{x})>\varepsilon\\\
1&\text{else}\end{cases}$ (7)

Adding multiple levels provides a direct extension to produce multiphase
alloys.

Fig. 3 illustrates the overall procedure in one dimension. The black line is
the GRF realization on the 1D domain. The red dashed line is the threshold
$\varepsilon$, which results in the binary field in blue. The orange steps
indicate the volume fraction per cell, which have the distribution shown in
the side panel. The volume fraction distribution is erratic for this small
sample but has the most mass at the expected value.

Fig. 4 shows representative SVEs constructed with correlation

$\kappa(r)=0.95\exp(-r/1.5)+0.05\exp(-r/15.0)\ ,$ (8)

which resembles that of Fig. 2, where $r$ is in units of voxels, and the mean
of the process is $\bar{\varphi}=0.1$ In both the SVE and the larger samples
shown in Fig. 5 coherent structures and structural scales can be observed.

Figure 3: 1D illustration of the microstructure generation process. A
realization of the underlying random process is depicted by the black curve,
the binarization threshold is indicted by the horizontal dashed red line, and
the resulting binary phase is shown as a blue histogram, while the volume
fractions of SVEs are shown as an orange histogram. The histograms on the
right show the distribution of the random process (black) and the volume
fractions (orange) for the selected SVE size.  Figure 4: Representative SVEs
Dividing lines are at voxel size.  Figure 5: Representative larger scale
samples. Dividing lines are at the SVE size shown in Fig. 4.

##  5 Dataset

Our exemplar material is nominally AlSi10Mg, which we simplify to a binary
mixture of 10% Si and 90% Al. We assume the two phases are segregated as in
Fig. 1. Over the range of interest, Si is modeled as elastic, and Al is
modeled as elastic-plastic with properties given in Table 1. Specifically, we
used a standard J2 elastic-plastic model with isotropic hardening for Al. The
Al hardening follows the often-employed Swift-Voce (SV) rule [65, 66] :

$h(\epsilon_{p})=Y+H\epsilon_{p}+A(1-\exp(-B\epsilon_{p})),$ (9)

which is essentially linear hardening augmented with an exponential trend. In
Eq. (9)

$\epsilon_{p}=\sqrt{\left(\frac{2}{3}\|{\operatorname{dev}}{\boldsymbol{\epsilon}}_{p}\|\right)}$
(10)

is the equivalent plastic strain, and ${\boldsymbol{\epsilon}}_{p}$ is the
plastic strain.

We used tensile loading to extract the relevant properties of the SVE samples
and ignore the anisotropy induced by the microstructural configurations for
simplicity although it is likely a source of error. We used periodic
displacement boundary conditions to reduce the potential for boundary effects
overwhelming the bulk response.

The volume fraction of Si $\upsilon$ of each SVE was determined directly from
the structure. The mechanical properties were extracted from the stress-strain
response using standard methods. Young’s modulus $E$ was fit to the mean slope
$\sigma_{11}/\epsilon_{11}$, and Poisson’s ratio $\nu$ to the mean lateral
contractions $(\epsilon_{22}+\epsilon_{33})/2\epsilon_{11}$. Using the fitted
$E$ the plastic regime was determined by the 0.2 % offset in strain rule. All
the data in the plastic regime was used to fit the Voce rule Eq. (9) to
determine yield $Y$ and the three hardening parameters $H$,$A$,$B$ for every
SVE. The similarity of the SVE response to the SV rule enabled a high accuracy
representation, although this generally may not be valid. Note that the linear
and exponential components of the SV hardening rule can confound each other;
and so care in fitting is needed, e.g. good starting guesses to optimizer in
region of attraction or regularization/penalization of irregular fits.
Although generally not feasible, for the present case regressing to a known
functional form to extract property features was an accurate reduction of the
response complexity. The range of effective properties and the mean standard
errors of the fits are given in Table 2.

Prior to sampling a large ensemble of SVEs and their properties for a training
dataset, we explored the size dependence of the properties of SVEs with 43,
83, 163, 243, 323 elements. These SVEs were extracted from a DNS structure at
least 10 times as large as an SVE with the selected volume fraction, Fig. 6
shows a sampling of the stress-strain response curves on a log-log scale to
show the range of responses. For the smallest SVE size, 43, there is a
significant propensity to sample high Si volume fractions, and for these
samples, the stress-strain response is nearly linear. Fig. 6 also illustrates
the decreasing variance in response with increasing SVE size.

To illustrate a case where there is significant variance, and spatial
correlation, but not the extremes illustrated by the smaller SVEs we selected
an SVE size of 243 elements since this was also a good trade-off between
structure resolution and computational expense. As mentioned, the SVEs are
taken at random from multiple 103 larger realizations with 0.1 volume fraction
of Si. Fig. 7 shows representative SVEs for a range of volume fractions and
their corresponding deformations and equivalent plastic strain fields. As can
be seen, there is a complex interplay of the structure and the imposed
deformation. Properties extracted from the 2000 SVEs used for training are
summarized in Table 2. Fig. 8 shows property correlations with volume
fraction. Clearly, the response depends on structural aspects beyond just
volume fraction, and a simple mixture model of the properties for Si and Al
would be insufficient. Also apparent are the effects of fitting to a model
that has multiple capacities to match nearly linear behavior, particularly in
the $B$ parameter.

Fig. 9 shows the correlation lengths of various structural-property aspects of
the particular microstructural process we are using. The underlying continuous
structure-generating process has a significantly shorter correlation length
than the binary voxel phase assignments. The induced structural property
volume fraction $\upsilon$ and the mechanical properties of the SVEs have even
longer correlation lengths. The latent vectors generated by the NN encoder of
the structure $\mathsf{N}\\!\mathsf{N}_{\boldsymbol{\ell}}$, described in
Sec. 6, have a correlation length between that of the structure process and
the induced property process.

Al |  |  |
---|---|---|---
Young’s modulus | $E$ | 70.0 | GPa
Poisson’s ratio | $\nu$ | 0.32 |
Yield strength | $Y$ | 120 | MPa
Voce linear hardening coefficient | $H$ | 290.8 | MPa
Voce exponential hardening coefficient | $A$ | 1.0 | MPa
Voce exponential hardening prefactor | $B$ | 434.9 |
Si |  |  |
Young’s modulus | $E$ | 140 | GPa
Poisson’s ratio | $\nu$ | 0.30 |
Table 1: Direct numerical simulations (DNS) properties for AL and Si assumed
elastic over range of interest.  Figure 6:  Stress-strain response for a
sampling of different SVE sizes. Note the log scale on both axes.

Figure 7: Realizations of 243 cells for volume fractions $\upsilon$ = 0.1(left), 0.2 (middle), 0.3(right). Phase (upper row) red:Al, gray:Si. Equivalent plastic strain (lower row) at failure 0.0 blue, red: 1.5  Property |  | min | max | std.error |
---|---|---|---|---|---
Volume fraction | $\upsilon$ | 0.0 | 0.86 | 0.0000 |
Young’s moduli | $E$ | 69.9 | 115.1 | 0.0299 | GPa
Poisson’s ratio | $\nu$ | 0.27 | 0.32 | 0.0000 |
Yield strength | $Y$ | 120 | 185 | 0.0005 | MPa
Linear hardening coefficient | $H$ | 0.0 | 39.0 | 0.0937 | GPa
Exponential hardening coefficient | $A$ | 0.6 | 8226 | 0.0300 | MPa
Exponential hardening exponent | $B$ | 0.84 | 1152 | 34.1342 |
Table 2: Stochastic volume element (SVE) property ranges for nominal volume
fraction $\upsilon=0.10$. Overall model fit RMSE = 0.0023 MPa.  Figure 8:
Property correlations with volume fraction $\upsilon$. Note stress-like
quantities are in MPa.  Figure 9: Correlation of underlying process, resulting
phase field, volume fraction of SVEs, as well as SVE latent vectors, and
properties. Correlation lengths $\ell$ and decay exponents $p$ were extracted
by fitting to $\exp(-(r/\ell)^{p})$. Note the minor negative correlations are
likely due to numerical noise from undersampling at long correlation lengths.

##  6 Structure-property encoding network

The neural network (NN) we use to encode the microstructures into a compact
latent space is a variational autoencoder (VAE) [26] augmented with a
regressor. Since our ultimate goal is to encode the structure-to-property map
we use an architecture akin to that of Wang et al. [8], which we will call a
property variational autoencoder (pVAE). This NN has three parts: (a) an
encoder of the structure to latent representation
$\phi\to{\boldsymbol{\ell}}$, (b) decoder that takes latent representations
back to structures ${\boldsymbol{\ell}}\to\phi$, and (c) a regressor that
takes latents to properties of the structure
${\boldsymbol{\ell}}\to{\text{\bf\th}}$. The architecture is illustrated in
the schematic in Fig. 10. Here $\phi(\mathbf{x})$ is the field describing the
microstructure, in the present context the Al/Si phase,
$\phi^{\prime}(\mathbf{x})$ is its reconstruction, ${\boldsymbol{\ell}}$ is
the latent vector that encodes $\phi$, and þ is the collection of properties
associated with $\phi$.

The loss function

$\mathcal{L}(\theta_{e},\theta_{d},\vartheta;\phi,{\text{\bf\th}})=\mathbb{E}_{q_{\theta_{e}}({\boldsymbol{\ell}}|\phi)}\left[\log
p_{\theta_{d}}(\phi)\right]-\beta\,\mathbb{E}_{q_{\theta_{e}}({\boldsymbol{\ell}}|\phi)}\left[\log\frac{q_{\theta_{e}}({\boldsymbol{\ell}}|\phi)}{p_{\theta_{d}}({\boldsymbol{\ell}}|\phi)}\right]+\lambda\,\|{\text{\bf\th}}-{\text{\bf\th}}_{\vartheta}({\boldsymbol{\ell}})\|^{2}$
(11)

reflects the scalarization of three competing objectives: (a) the data fit of
the structure reconstruction $\phi\to\phi^{\prime}$ from encoder to latent
space to decoder, (b) the Kullback-Leibler (KL) divergence between the latent
distribution of the data $p_{\theta_{d}}$ and its representation by a
surrogate distribution $q_{\vartheta}$, and (c) the regression error of the
latent to property map ${\text{\bf\th}}_{\vartheta}({\boldsymbol{\ell}})$.
Here $q_{\theta_{e}}({\boldsymbol{\ell}}|\phi)$ is the distribution of latent
vectors given microstructures produced by the encoder, and
$p_{\theta_{d}}({\boldsymbol{\ell}}|\phi)$, conversely, is the distribution
produced by the decoder. The KL term promotes similiarity of the trained
latent distribution and the distribution implicit in the data. As in the
original VAE we assume that the distribution
$q_{\theta_{e}}({\boldsymbol{\ell}}|\phi)=\mathcal{N}({\boldsymbol{\mu}},{\boldsymbol{\sigma}})$
is well approximated by a multivariate, independent Gaussian distribution so
that
${\boldsymbol{\ell}}={\boldsymbol{\mu}}+{\boldsymbol{\sigma}}\odot{\boldsymbol{\varepsilon}}$
where $\odot$ is the component-wise Hadamard product and
${\boldsymbol{\varepsilon}}\sim\mathcal{N}({\boldsymbol{\mu}},\mathbf{1})$.
The vectors ${\boldsymbol{\ell}}$, ${\boldsymbol{\mu}}$,
${\boldsymbol{\sigma}}$, and ${\boldsymbol{\varepsilon}}$ are all of length
$N_{\boldsymbol{\ell}}$, which is the size of the latent space and a
hyperparameter. The parameters ${\boldsymbol{\mu}}$ and
${\boldsymbol{\sigma}}$ are trained at the same time as $\theta_{e}$,
$\theta_{d}$, and $\vartheta$. This assumption also allows for a particularly
simple form of the KL divergence. As in the $\beta$-VAE [30] we also use
hyperparameters $\beta$ and $\lambda$ to balance the relative scale and
importance of the separate objectives. In addition, these allow annealing, in
the sense of a schedule of increasing and decreasing $\beta$ or $\lambda$, to
aid training to the available data. Reconstruction loss is mainly a
regularizer of the regression task and is said to organize the latent space
[8] by promoting it to be continuous with respect to the regression task.

The trainable parameters are summarized in Table 3. We exploited the relative
computational cost of obtaining the data in training the pVAE network. Since
the encoder-decoder needed to be highly parameterized to encode the wide
variety of structures, we first trained the NN to a large dataset of 100,000
samples of the Al-Si microstructure. For this first stage the only property we
expect the regressor to predict is the volume fraction since this data is
obtainable with negligible cost. Then we warm-started the network with a much
smaller training set of 2000 samples where we added the 6 mechanical
properties that were relatively expensive to compute. Note that only the
output layer of the regressor changed size in this stage. After allowing the
latent space to adjust to the new tasks, we increased $\lambda$ to make the
training focus on the accuracy of the regressor. (This could also be
accomplished by freezing the parameters of the encoder-decoder.)

Fig. 11 shows how the latent space adapts to the new task of predicting the
mechanical properties. The latent samples are projected into the first two
principal component analysis (PCA) components. Clearly the latent space
transitions from a generally Gaussian distribution to a distribution that is
less normal. Fig. 12 gives projected views in the first 10 PCA components. The
distribution in the first four components appear to be non-Gaussian, while the
remaining, less significant components appear to be reasonably Gaussian.

Fig. 13 shows that the property regression errors for the final pVAE are
generally between 2 and 5 percent. The properties with more complex
distributions and less simple correlations with the SVE volume fraction,
namely the hardening parameters $H$ and $B$, have the highest errors. We
attribute the high relative error in the regression of the Poisson’s ratio to
the narrow range of the true values (Al and Si have similar $\nu$) and hence
the low sensitivity of the Poisson’s ratio to the structure.

Figure 10: Property variational autoencoder (pVAE) schematic. The encoder (red) encodes the microstructure $\mathsf{X}=\phi(\mathbf{x})$ to a compact latent description ${\boldsymbol{\ell}}$ which is then passed to the decoder (blue) to output a reconstruction $\mathsf{X}^{\prime}$ and the property regressor (green) which outputs a property prediction þ. The latent vector ${\boldsymbol{\ell}}$ has a parametrized distribution $\mathcal{N}({\boldsymbol{\mu}},{\boldsymbol{\sigma}})$.  Layer | output | parameters
---|---|---
Encoder |  |
convolution | 243 $\times$ 16 | 448
max pooling | 123 $\times$ 16 | 0
convolution | 123 $\times$ 32 | 13856
max pooling | 63 $\times$ 32 | 0
convolution | 63 $\times$ 64 | 55360
Latent |  |
dense | 2 $\times$ $N_{\boldsymbol{\ell}}$ | 840
Decoder |  |
convolution | 123 $\times$ 32 | 55328
convolution | 243 $\times$ 16 | 13840
convolution | 243 $\times$ 1 | 433
Regressor |  |
dense | $N_{\boldsymbol{\ell}}$ | 420
dense | $N_{\boldsymbol{\ell}}$ | 420
dense | $N_{\boldsymbol{\ell}}$ | 420
dense | $N_{\boldsymbol{\ell}}$ | 420
dense | $N_{\text{\bf\th}}$ | 147
Table 3: Property variational autoencoder (pVAE) configuration. GELUs were
used for all activations including the property output layer since all
properties are non-negative. Dependence on batch size suppressed for clarity.

Figure 11: Latent sample: first two PCA components colored by volume fraction:
(left) pVAE trained to volume fraction only, (right) pVAE trained to all
properties.  Figure 12: Latent PCA components. Note only 10 of the 20
components are shown for clarity.

Figure 13: Property prediction correlations with held out data (left) and CDF
of relative/normalized errors (right). Note the CDFs terminate at the highest
sample in the data.

Latent encoding of properties serves the purpose of giving a compact
representation of the correlation of properties at any point, and hence the
structure-property correlations through the latent features. Fig. 14 shows the
variation of each property with the first and third PCA components of the
latent space, which are the most aligned with variations in properties.
Clearly, there is some correlation between all properties. We conjecture that
the second PCA component does not correlate significantly with properties
since it is used to facilitate the structure reconstruction of the pVAE. Also,
note that more diffuse lobes in the scatter plots are low density and likely
due to fitting the hardening with two competing functions, see Eq. (9).

We associate the components of the _semantic direction_ of a property ${\text{
\th}}({\boldsymbol{\ell}})$ with the correlation coefficient of the
(orthogonal) PCA directions ${\boldsymbol{\ell}}_{I}$

$\boldsymbol{\zeta}_{I}=\frac{\langle{\text{
\th}},{\boldsymbol{\ell}}_{I}\rangle}{\sqrt{\langle{\text{ \th}},{\text{
\th}}\rangle\langle{\boldsymbol{\ell}}_{I},{\boldsymbol{\ell}}_{I}\rangle}}$
(12)

where $\langle{\text{ \th}},{\text{ \th}}\rangle$ is the variance of  þ. The
semantic direction, in the sense of the direction of highest variance (and
gradient), may not be constant in latent space with a complex locus of
representative samples. Fig. 15 shows that generally all the properties are
strongly correlated with the first and third PCA component of the latent space
although the correlation is lower for the hardening parameters that have a
more complex dependence. There are also significant correlations of all
properties with less significant PCA directions. Clearly, there is no perfect,
simple alignment of the latent space variance with the property variations.

Figure 14: Latent sample: first and third PCA components colored by
corresponding properties indicating semantic directions.  Figure 15: Latent-
property correlations in terms of PCA components.

##  7 Demonstrations

In this section, we provide: (a) a comparison of the proposed method with a
DNS, (b) a demonstration of its use in propagating microstructural uncertainty
to mechanical quantities of interest, and (c) an illustration of how it can be
used to model materials with microstructural functional gradation. Each
demonstration utilizes a component geometry with three competing holes, which
is meant to be representative of a topological optimization scenario with an
uncertain material. Each is loaded in tension for ease of interpretation by
applying fixed $x_{1}$ boundary conditions on one end, $x_{2}$ and $x_{3}$
boundary conditions on edges of this end and prescribed boundary conditions on
the opposing end. For the DNS comparison, we used a 84$\times$48$\times$4
$\mu$m component, and for the generative demonstrations we used
240$\times$40$\times$4 $\mu$m configurations. Each multiscale simulation was
performed with a standard finite element simulator [67] that was modified to
take the property regressor component of the trained pVAE and a field of
latent vectors. Since the imprinting of the mesh with a latent vector field is
done as a preprocessing step and the complexity of the constitutive model is
comparable to that of a traditional model, the speed up is proportional to the
reduction in the number of elements, i.e. for a DNS boundary value problem
encoded into $n^{3}$-element SVEs the speed up is $n^{3}$.

###  7.1 DNS comparison

First, we show how well the proposed method approximates a corresponding DNS.
Using the GRF structure generation algorithm we created a large
504$\times$288$\times$24 voxel structure. A structured rectangular grid was
employed to facilitate partitioning the overall structure into disjoint SVEs.
Using the trained pVAE we encoded each of the SVEs into latent vectors that
will be subsequently used in a multiscale simulation and mapped to properties.
This DNS was compared to a corresponding multiscale simulation with a selected
SVE size.

Fig. 16 shows that if there are no holes, as in Fig. 17, there are weak
gradients and the multiscale method approximates the response of the DNS. Fig.
17 also shows that the coarse field resembles that of the DNS in this case.
Holes, as in the configuration shown in Fig. 18, create gradients that
interact with the microstructure. There is apparently a bias in the hardening
response that leads the multiscale methods to overestimate the nominal stress
in this scenario. In this case, the multiscale response with a coarse 243 SVEs
is moderately far from the DNS response; however, employing 83 SVEs provides
better resolution of the gradients and brings the response closer to that of
the DNS.

Figure 16: DNS comparison of sample norminal stress-strain response. Nominal
strain was obtained from the end displacement of the sample and nominal strain
was calculated from the reaction force on the end and its initial area.

Figure 17: DNS comparison of localization for the part without holes.
Equivalent plastic strain at nominal strain 0.007, blue: 0.0, red DNS 0.028
SVE 243 0.009,

Figure 18: DNS comparison of localization for part with 3 holes. Equivalent
plastic strain at nominal strain 0.007, blue: 0.0, red DNS 0.069, SVE 243
0.020, SVE 83 0.028.

###  7.2 Uncertainty quantification

To propagate uncertainty in performance quantities of interest (QoI) due to
uncertainty in the material microstructure, we construct an ensemble of
realizations of the latent fields on a mesh of the component and then evaluate
the resulting ensemble of responses. To construct realizations of the latent
field ${\boldsymbol{\ell}}(\mathbf{x})$ consistent with the chosen
microstructure we need a method that embeds the spatial correlations of the
latent space and the point-wise latent component correlations that embed the
correlations between properties. A KLE is adaptable to this purpose. A KLE
uses a kernel to construct a (deterministic) correlation matrix for the mesh
by encoding spatial interactions via

${\boldsymbol{\Sigma}}_{IJ}=\kappa(\mathbf{x}_{I},\mathbf{x}_{J})=\kappa(\|\mathbf{x}_{I}-\mathbf{x}_{J}\|)$
(13)

by taking $\mathbf{x}_{I}$ to be cell/element centroids of the mesh employed
for the realization. The corresponding eigenvectors $\mathbf{v}_{a}$ and
eigenvalues $\lambda_{a}$ of ${\boldsymbol{\Sigma}}$ are then used in the
expansion:

$\varphi(\mathbf{x})=\sum_{a}\eta_{a}\sqrt{\lambda_{a}}\mathbf{v}_{a}$ (14)

which is truncated so that the largest included eigenvalues represent a given
fraction of the total variance (we used 0.99). Here, $\eta$ is a uncorrelated
noise process, e.g. $\eta\sim\mathcal{N}(0,1)$, with realizations $\eta_{a}$.
This construction is sufficient to represent the spatial correlation of a
scalar random field. To account for the point-wise correlations of
${\boldsymbol{\ell}}$ we estimated the correlation
$C_{ab}=\langle{{\ell}}_{a},{{\ell}}_{b}\rangle$ from encoding the large
training sample and then construct:

${\boldsymbol{\ell}}(\mathbf{x})=\sum_{a}\eta_{a}\sqrt{\lambda_{a}}\mathsf{L}\mathbf{v}_{a}$
(15)

using the Cholesky decomposition $\mathsf{L}$ of
$\mathsf{C}=\mathsf{L}\mathsf{L}^{T}$ [23]. Modifications of this
representation exist for non-Gaussian fields [19], but we assumed this
representation was sufficient for the current context. Also, the samples used
in training are equivalent to an empirical distribution that could have been
used for this purpose, but the route we chose is extensible to as many
independent realizations as needed and provides some filtering of outliers.

To illustrate the effect of the spatial correlation across the realizations we
simulated a scenario where the structure correlation in the latent vectors is
comparable to deformation gradients and compare it to a corresponding one
where the spatial correlation is too short to be effectively uncorrelated over
the mesh. Fig. 19 shows a realization of the correlated structure where the
latents have been mapped to their corresponding properties. Fig. 20 compares
the two response ensembles of 64 realizations each, interestingly the
spatially correlated ensemble clearly has the higher variance apparently due
to the fact that the spatial correlation reinforces high or low properties at
the necks near the holes. In contrast, the uncorrelated ensemble tends to have
elements with moderate properties near elements with outlying properties. Fig.
21 corroborates this conjecture. The necks and near the interior holes have
the highest equivalent plastic strain and stress, respectively, and also the
highest variance in those critical fields.

Figure 19: Uncertainty quantification: One realization of the bar with 3
holes. Variation of (upper) modulus $E$ and (lower) yield strength $Y$. The
blue-red range of $E\in[69.9,115.1]$ and $Y\in[120,185]$ are given in Table 2.

Figure 20: Uncertainty quantification: stress-strain response (left)
uncorrelated and (right) spatially correlated realizations.

Figure 21: Uncertainty quantification: equivalent plastic strain (top) at
nominal strain 0.007 (blue: 0.00, red:0.08, white:0.000, black:0.056) and
stress (bottom, blue: -85 MPa red:60 MPa, white:0 MPa, black:66 MPa) fields.
Mean field is in color and the standard deviation field is in gray.

###  7.3 Functional gradation

Functional gradation is achieved by modulating the manufacturing process to
create gradations of properties over a component. To mimic this process, we
use the semantic direction for volume fraction $\upsilon$ to create a non-
stationary distribution

$\pi({\boldsymbol{\ell}};\xi)=\mathcal{N}({\boldsymbol{\mu}}+\xi\,\mathbf{s},{\boldsymbol{\Sigma}})$
(16)

which is a function of $\xi$ which we associate with the normalized distance
to the surface, as shown in Fig. 22. The semantic direction of the volume
fraction was chosen since this is directly accessible to a manufacturing
process, however the mechanical properties also vary with it. We chose $\mu$
to be the center of the training latent space and ${\boldsymbol{\Sigma}}$ to
be 0.5 of the covariance of the training latent space to represent a more
tightly controlled process.

We compare two scenarios: (a) where the Si volume fraction increases toward
the surface (outward gradient), and (b) where the Si volume fraction increases
toward the interior of the part (inward gradient). Fig. 22 illustrates the
functional gradation in properties for the inward scenario, where the trend is
clear as is the stochastic variation. Fig. 23 shows that these two scenarios
produce distinct ensembles of (32) responses and corresponding distributions
performance QoIs. Clearly, having higher concentration of Si near the surface
leads to tougher parts. Fig. 24 and Fig. 25 illustrate that the local changes
in response are fairly subtle. The outward variation appears to spread the
concentration of equivalent plastic strain and stress near the holes.

Figure 22: Functional gradation: (top) distance field. (middle) yield stress
$Y$ inward gradation, (bottom) yield stress $Y$ outward gradation.

Figure 23: Functional gradation: stress-strain response inward (top) and
outward (bottom) gradation.

Figure 24: Functional gradation: equivalent plastic strain field at nominal
strain 0.007. Inward (top) blue:0, red:0.075, white: 0.0, black:0.029. Outward
(bottom) blue:0, red:0.050, white: 0.0, black:0.022. Mean (upper, color),
standard deviation (lower, gray).

Figure 25: Functional gradation: 11-stress field at nominal strain 0.007.
Inward (top) blue:-87 MPa, red: 63 MPa, white: 0.0, black:67 MPa. Outward
(bottom) blue:-11, red: 70 MPa, white: 0.0, black:11 MPa. Mean (upper, color),
standard deviation (lower, gray).

##  8 Conclusion

Using a pVAE structure-property representation, KLE realization generation,
and a standard finite element simulator we developed a multiscale method that
accounts for the finite-size effects and spatial correlations that complicate
simulation scenarios where there is no clear separation of the scales of
deformation and microstructure. The method produces a latent space that
compactly and extensibly encodes the microstructure-response map. This latent
space and the semantic directions implicit in it are applicable to a variety
of applications. We demonstrated its ability to replicate DNS QoIs, push
forward microstructure uncertainty, and simulate functional gradation of
stochastic media. The demonstrations employed a common AM material, AlSi10Mg;
however the proposed method is amenable to other microstructures such as those
with multiple phases and/or with pores (i.e. a null material). In materials
such as those generated by AM, defects can strongly affect the local stress
gradients, and therefore scale separation.

The expediency of mapping structures to properities allowed us to interpret
the maps and focus on the manipulation of the latent spaces to practical
applications. Encoding the structure-property map had straightforward
implementation in a standard finite element code and was computationally
efficient because the imprinting of latents on a mesh is done prior to
simulation. More general applications will likely require a more expressive
structure-to-response map. In future work we would like to extend the
application of this method to microstructure ranging from random media with
statistical characterization to highly structured lattice/metamaterials [68,
69, 70] and also to more expressive latent distributions. The method also has
applications beyond forward propagation of material variability/uncertainty
[45], it could be indispensable in topology optimization [71] and material and
structural optimization [72, 73].

## Acknowledgments

The authors would like to thank Sara Dickens and Paul Specht, Sandia National
Laboratories, for the microphotograph of AlSi10Mg. The simulations were
performed with the SIERRA finite element simulator developed by Sandia
National Laboratories. Sandia National Laboratories is a multimission
laboratory managed and operated by National Technology and Engineering
Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell
International, Inc., for the U.S. Department of Energy’s National Nuclear
Security Administration under contract DE-NA-0003525. This paper describes
objective technical results and analysis. Any subjective views or opinions
that might be expressed in the paper do not necessarily represent the views of
the U.S. Department of Energy or the United States Government.

## References

  * [1] Yihong Kok, Xi Peng Tan, Pan Wang, MLS Nai, Ngiap Hiang Loh, Erjia Liu, and Shu Beng Tor.  Anisotropy and heterogeneity of microstructure and mechanical properties in metal additive manufacturing: A critical review.  Materials & Design, 139:565–586, 2018.
  * [2] Pengcheng Jiao, Jochen Mueller, Jordan R Raney, Xiaoyu Zheng, and Amir H Alavi.  Mechanical metamaterials and beyond.  Nature Communications, 14(1):6004, 2023.
  * [3] Norman A Fleck, Vikram S Deshpande, and Michael F Ashby.  Micro-architectured materials: past, present and future.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 466(2121):2495–2516, 2010.
  * [4] Jörg Schröder.  A numerical two-scale homogenization scheme: the FE2-method.  In Plasticity and beyond: microstructures, crystal-plasticity and phase transitions, pages 1–64. Springer, 2014.
  * [5] Frédéric Feyel.  A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua.  Computer Methods in applied Mechanics and engineering, 192(28-30):3233–3244, 2003.
  * [6] Siavouche Nemat-Nasser and Muneo Hori.  Micromechanics: overall properties of heterogeneous materials.  Elsevier, 2013.
  * [7] David T Fullwood, Stephen R Niezgoda, Brent L Adams, and Surya R Kalidindi.  Microstructure sensitive design for performance optimization.  Progress in Materials Science, 55(6):477–562, 2010.
  * [8] Liwei Wang, Yu-Chin Chan, Faez Ahmed, Zhao Liu, Ping Zhu, and Wei Chen.  Deep generative modeling for mechanistic-based learning and design of metamaterial systems.  Computer Methods in Applied Mechanics and Engineering, 372:113377, 2020.
  * [9] Paul E Specht and Nathan P Brown.  Shock compression response of additively manufactured alsi10mg.  Journal of Applied Physics, 130(24), 2021.
  * [10] Naoki Takata, Hirohisa Kodaira, Asuka Suzuki, and Makoto Kobashi.  Size dependence of microstructure of AlSi10Mg alloy fabricated by selective laser melting.  Materials Characterization, 143:18–26, 2018.
  * [11] Stephen R Niezgoda, David M Turner, David T Fullwood, and Surya R Kalidindi.  Optimized structure based representative volume element sets reflecting the ensemble-averaged 2-point statistics.  Acta Materialia, 58(13):4432–4445, 2010.
  * [12] Surya R Kalidindi, Stephen R Niezgoda, and Ayman A Salem.  Microstructure informatics using higher-order statistics and efficient data-mining protocols.  Jom, 63:34–41, 2011.
  * [13] Stephen R Niezgoda, Anand K Kanjarla, and Surya R Kalidindi.  Novel microstructure quantification framework for databasing, visualization, and analysis of microstructure data.  Integrating Materials and Manufacturing Innovation, 2:54–80, 2013\.
  * [14] Surya R Kalidindi.  Hierarchical materials informatics: novel analytics for materials data.  Elsevier, 2015.
  * [15] Ahmet Cecen, Tony Fast, and Surya R Kalidindi.  Versatile algorithms for the computation of 2-point spatial correlations in quantifying material structure.  Integrating Materials and Manufacturing Innovation, 5(1):1–15, 2016\.
  * [16] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [17] FS Gabrielli, A Labini and L Joyce, M Pietronero.  Statistical Physics for Cosmic Structures.  Springer, 2006.
  * [18] Dirk P Kroese and Zdravko I Botev.  Spatial process simulation.  In Stochastic geometry, spatial statistics and random fields: Models and algorithms, pages 369–404. Springer, 2014.
  * [19] Kok-Kwang Phoon, HW Huang, and Ser Tong Quek.  Simulation of strongly non-Gaussian processes using karhunen–loeve expansion.  Probabilistic engineering mechanics, 20(2):188–198, 2005.
  * [20] Diego Lorenzo Allaix and Vincenzo Ilario Carbone.  Karhunen–loève decomposition of random fields based on a hierarchical matrix approach.  International journal for numerical methods in engineering, 94(11):1015–1036, 2013.
  * [21] Heyrim Cho, D Venturi, and GE Karniadakis.  Karhunen–loève expansion for multi-correlated stochastic processes.  Probabilistic Engineering Mechanics, 34:157–167, 2013.
  * [22] Alfonso M Panunzio, Régis Cottereau, and Guillaume Puel.  Large scale random fields generation using localized karhunen–loève expansion.  Advanced Modeling and Simulation in Engineering Sciences, 5:1–29, 2018.
  * [23] Ranadeep Daw, Matthew Simpson, Christopher K Wikle, Scott H Holan, and Jonathan R Bradley.  An overview of univariate and multivariate karhunen loève expansions in statistics.  Journal of the Indian Society for Probability and Statistics, 23(2):285–326, 2022.
  * [24] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [25] Carl Doersch.  Tutorial on variational autoencoders.  arXiv preprint arXiv:1606.05908, 2016.
  * [26] Diederik P Kingma, Max Welling, et al.  An introduction to variational autoencoders.  Foundations and Trends® in Machine Learning, 12(4):307–392, 2019.
  * [27] Chris J Maddison, Andriy Mnih, and Yee Whye Teh.  The concrete distribution: A continuous relaxation of discrete random variables.  arXiv preprint arXiv:1611.00712, 2016.
  * [28] Francisco Mena and Ricardo Ñanculef.  A binary variational autoencoder for hashing.  In Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications: 24th Iberoamerican Congress, CIARP 2019, Havana, Cuba, October 28-31, 2019, Proceedings 24, pages 131–141. Springer, 2019.
  * [29] Jason Tyler Rolfe.  Discrete variational autoencoders.  arXiv preprint arXiv:1609.02200, 2016.
  * [30] Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner.  beta-VAE: Learning basic visual concepts with a constrained variational framework.  In International conference on learning representations, 2017.
  * [31] Toshio Mura.  Micromechanics of defects in solids.  Springer Science & Business Media, 2013.
  * [32] Salvatore Torquato.  Random Heterogeneous Materials: Microstructure and Macroscopic Properties, volume 16.  Springer Science & Business Media, 2013.
  * [33] Rodney Hill.  Elastic properties of reinforced solids: some theoretical principles.  Journal of the Mechanics and Physics of Solids, 11(5):357–372, 1963\.
  * [34] J Mandel.  Plasticité classique et viscoplasticité (cism, udine, 1971), 1972\.
  * [35] Shivakumar I Ranganathan and Martin Ostoja-Starzewski.  Scaling function, anisotropy and the size of RVE in elastic random polycrystals.  Journal of the Mechanics and Physics of Solids, 56(9):2773–2791, 2008.
  * [36] Martin Ostoja-Starzewski.  Microstructural randomness and scaling in mechanics of materials.  CRC Press, 2007.
  * [37] Martin Ostoja-Starzewski.  Microstructural disorder, mesoscale finite elements and macroscopic response.  Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 455(1989):3189–3199, 1999.
  * [38] Ioannis G Kevrekidis, C William Gear, James M Hyman, Panagiotis G Kevrekidis, Olof Runborg, Constantinos Theodoropoulos, et al.  Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis.  Commun. Math. Sci, 1(4):715–762, 2003.
  * [39] E Weinan, Bjorn Engquist, Xiantao Li, Weiqing Ren, and Eric Vanden-Eijnden.  Heterogeneous multiscale methods: a review.  Communications in computational physics, 2(3):367–450, 2007.
  * [40] Peiyi Chen, Johann Guilleminot, and Christian Soize.  Concurrent multiscale simulations of nonlinear random materials using probabilistic learning.  Computer Methods in Applied Mechanics and Engineering, 422:116837, 2024.
  * [41] Mark Girolami, Eky Febrianto, Ge Yin, and Fehmi Cirak.  The statistical finite element method (statfem) for coherent synthesis of observation data and model predictions.  Computer Methods in Applied Mechanics and Engineering, 375:113533, 2021.
  * [42] Grace CY Peng, Mark Alber, Adrian Buganza Tepole, William R Cannon, Suvranu De, Savador Dura-Bernal, Krishna Garikipati, George Karniadakis, William W Lytton, Paris Perdikaris, et al.  Multiscale modeling meets machine learning: What can we learn?  Archives of Computational Methods in Engineering, 28:1017–1037, 2021.
  * [43] Jacob Fish, Gregory J Wagner, and Sinan Keten.  Mesoscopic and multiscale modelling in materials.  Nature materials, 20(6):774–786, 2021.
  * [44] Dana Bishara, Yuxi Xie, Wing Kam Liu, and Shaofan Li.  A state-of-the-art review on machine learning-based multiscale modeling, simulation, homogenization and design of materials.  Archives of computational methods in engineering, 30(1):191–222, 2023.
  * [45] Francesco Rizzi, Mohammad Khalil, Reese E Jones, Jeremy Alan Templeton, Jakob T Ostien, and Brad L Boyce.  Bayesian modeling of inconsistent plastic response due to material variability.  Computer Methods in Applied Mechanics and Engineering, 353:183–200, 2019.
  * [46] Mohammad Khalil, Gregory Herlin Teichert, Coleman Alleman, NM Heckman, Reese E Jones, Krishnakumar Garikipati, and BL Boyce.  Modeling strength and failure variability due to porosity in additively manufactured metals.  Computer Methods in Applied Mechanics and Engineering, 373:113471, 2021.
  * [47] Miguel A Bessa, Ramin Bostanabad, Zeliang Liu, Anqi Hu, Daniel W Apley, Catherine Brinson, Wei Chen, and Wing Kam Liu.  A framework for data-driven analysis of materials under uncertainty: Countering the curse of dimensionality.  Computer Methods in Applied Mechanics and Engineering, 320:633–667, 2017.
  * [48] Trenton Kirchdoerfer and Michael Ortiz.  Data-driven computational mechanics.  Computer Methods in Applied Mechanics and Engineering, 304:81–101, 2016.
  * [49] Robert Eggersmann, Trenton Kirchdoerfer, Stefanie Reese, Laurent Stainier, and Michael Ortiz.  Model-free data-driven inelasticity.  Computer Methods in Applied Mechanics and Engineering, 350:81–99, 2019.
  * [50] Pietro Carrara, Laura De Lorenzis, Laurent Stainier, and Michael Ortiz.  Data-driven fracture mechanics.  Computer Methods in Applied Mechanics and Engineering, 372:113390, 2020.
  * [51] Konstantinos Karapiperis, Laurent Stainier, Michael Ortiz, and José E Andrade.  Data-driven multiscale modeling in mechanics.  Journal of the Mechanics and Physics of Solids, 147:104239, 2021\.
  * [52] Reese E Jones, Ari L Frankel, and KL Johnson.  A neural ordinary differential equation framework for modeling inelastic stress response via internal state variables.  Journal of Machine Learning for Modeling and Computing, 3(3), 2022\.
  * [53] Moritz Flaschel, Siddhant Kumar, and Laura De Lorenzis.  Discovering plasticity models without stress data.  npj Computational Materials, 8(1):91, 2022.
  * [54] Moritz Flaschel, Siddhant Kumar, and Laura De Lorenzis.  Automated discovery of generalized standard material models with euclid.  Computer Methods in Applied Mechanics and Engineering, 405:115867, 2023.
  * [55] Burigede Liu, Eric Ocegueda, Margaret Trautner, Andrew M Stuart, and Kaushik Bhattacharya.  Learning macroscopic internal variables and history dependence from microscopic models.  Journal of the Mechanics and Physics of Solids, 178:105329, 2023\.
  * [56] Ari L Frankel, Reese E Jones, Coleman Alleman, and Jeremy A Templeton.  Predicting the mechanical response of oligocrystals with deep learning.  Computational Materials Science, 169:109099, 2019.
  * [57] Ari Frankel, Kousuke Tachida, and Reese Jones.  Prediction of the evolution of the stress field of polycrystals undergoing elastic-plastic deformation with a hybrid neural network model.  Machine Learning: Science and Technology, 1(3):035005, 2020.
  * [58] Xiaolong He, Qizhi He, and Jiun-Shyan Chen.  Deep autoencoders for physics-constrained data-driven nonlinear materials modeling.  Computer Methods in Applied Mechanics and Engineering, 385:114034, 2021.
  * [59] Julian N Heidenreich, Maysam B Gorji, and Dirk Mohr.  Modeling structure-property relationships with convolutional neural networks: Yield surface prediction based on microstructure images.  International Journal of Plasticity, 163:103506, 2023.
  * [60] Gregory H Teichert, Mohammad Khalil, Coleman Alleman, Krishna Garikipati, and Reese E Jones.  Sensitivity of void mediated failure to geometric design features of porous metals.  International Journal of Solids and Structures, 236:111309, 2022\.
  * [61] Alan E Gelfand, Peter Diggle, Peter Guttorp, and Montserrat Fuentes.  Handbook of spatial statistics.  CRC press, 2010.
  * [62] Chunsheng Ma.  Construction of non-gaussian random fields with any given correlation structure.  Journal of Statistical Planning and Inference, 139(3):780–787, 2009\.
  * [63] EB Dynkin.  Gaussian and non-gaussian random fields associated with markov processes.  Journal of Functional Analysis, 55(3):344–376, 1984.
  * [64] Paolo Bocchini and George Deodatis.  Critical review and latest developments of a class of simulation algorithms for strongly non-gaussian random fields.  Probabilistic Engineering Mechanics, 23(4):393–407, 2008.
  * [65] E Voce.  The relationship between stress and strain for homogeneous deformation.  Journal of the Institute of Metals, 74:537–562, 1948.
  * [66] H. W. Swift.  Plastic instability under plane stress.  Journal of the Mechanics and Physics of Solids, 1(1):1–18, 1952\.
  * [67] Frank N Beckwith, Guy L Bergel, Gabriel J de Frias, Kevin L Manktelow, Mark T Merewether, Scott T Miller, Krishen J Parmar, Timothy R Shelton, Jesse D Thomas, Jeremy Trageser, et al.  Sierra/solidmechanics 5.6 user’s manual.  Technical report, Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 2022.
  * [68] Chuang Wang, Xiaojun Gu, Jihong Zhu, Han Zhou, Shaoying Li, and Weihong Zhang.  Concurrent design of hierarchical structures with three-dimensional parameterized lattice microstructures for additive manufacturing.  Structural and Multidisciplinary Optimization, 61:869–894, 2020\.
  * [69] Giuseppe Del Guercio, Manuela Galati, Abdollah Saboori, Paolo Fino, and Luca Iuliano.  Microstructure and mechanical performance of ti–6al–4v lattice structures manufactured via electron beam melting (ebm): a review.  Acta Metallurgica Sinica (English Letters), 33:183–203, 2020.
  * [70] Kavan Hazeli, Behzad Bahrami Babamiri, Joseph Indeck, Andrew Minor, and Hesam Askari.  Microstructure-topology relationship effects on the quasi-static and dynamic behavior of additively manufactured lattice structures.  Materials & Design, 176:107826, 2019.
  * [71] Liang Xia and Piotr Breitkopf.  Concurrent topology optimization design of material and structure within FE2 nonlinear multiscale analysis framework.  Computer Methods in Applied Mechanics and Engineering, 278:524–542, 2014.
  * [72] Mauricio Fernández, Felix Fritzen, and Oliver Weeger.  Material modeling for parametric, anisotropic finite strain hyperelasticity based on machine learning with application in optimization of metamaterials.  International Journal for Numerical Methods in Engineering, 123(2):577–609, 2022.
  * [73] Oliver Weeger, Iman Valizadeh, Yash Mistry, and Dhruv Bhate.  Inelastic finite deformation beam modeling, simulation, and validation of additively manufactured lattice structures.  Additive Manufacturing Letters, 4:100111, 2023.


11institutetext: Department of Radiology, Lausanne University Hospital and
University of Lausanne (UNIL), Lausanne, Switzerland 22institutetext: CIBM
Center for Biomedical Imaging, Switzerland 33institutetext: Department for
Early Life Imaging, School of Biomedical Engineering & Imaging Sciences,
King’s College London, London, UK 44institutetext: Smart Imaging Lab,
Diagnostic Radiology, FAU Erlangen-Nuremberg, Erlangen, Germany

# Improving cross-domain brain tissue segmentation in fetal MRI with synthetic
data††thanks: This research was funded by the Swiss National Science
Foundation (215641), ERA-NET Neuron MULTI-FACT project (SNSF 31NE30_203977),
UKRI FLF (MR/T018119/1) and DFG Heisenberg funding (502024488); We acknowledge
the Leenaards and Jeantet Foundations as well as CIBM Center for Biomedical
Imaging, a Swiss research center of excellence founded and supported by CHUV,
UNIL, EPFL, UNIGE and HUG.

Vladyslav Zalevskyi 1122 Thomas Sanchez 11 2 2 Margaux Roulet 11 2 2 Jordina
Aviles Verddera 33 Jana Hutter 33 4 4 Hamza Kebiri 11 2 2 Meritxell Bach
Cuadra
22 1 1<EMAIL_ADDRESS>
###### Abstract

Segmentation of fetal brain tissue from magnetic resonance imaging (MRI) plays
a crucial role in the study of in utero neurodevelopment. However, automated
tools face substantial domain shift challenges as they must be robust to
highly heterogeneous clinical data, often limited in numbers and lacking
annotations. Indeed, high variability of the fetal brain morphology, MRI
acquisition parameters, and super-resolution reconstruction (SR) algorithms
adversely affect the model’s performance when evaluated out-of-domain. In this
work, we introduce FetalSynthSeg, a domain randomization method to segment
fetal brain MRI, inspired by SynthSeg. Our results show that models trained
solely on synthetic data outperform models trained on real data in out-of-
domain settings, validated on a 120-subject cross-domain dataset. Furthermore,
we extend our evaluation to 40 subjects acquired using low-field (0.55T) MRI
and reconstructed with novel SR models, showcasing robustness across different
magnetic field strengths and SR algorithms. Leveraging a generative synthetic
approach, we tackle the domain shift problem in fetal brain MRI and offer
compelling prospects for applications in fields with limited and highly
heterogeneous data.

###### Keywords:

Domain shifts segmentation fetal brain MRI synthetic data low-field MRI

##  1 Introduction

Fetal brain magnetic resonance imaging (MRI) is a growing diagnostic tool for
studying neurodevelopment in fetuses [1, 2, 3, 4]. Despite its potential,
creating automated pipelines for fetal MRI faces challenges due to the limited
availability of annotated datasets and data heterogeneity. The fetal brain
undergoes significant morphological changes during gestation and can be
severely altered by pathologies, which complicates its automatic analysis [5,
6, 7]. Additionally, datasets encounter distribution shifts from variations in
acquisition sites, scanners, and imaging protocols [8, 9]. Utilizing super-
resolution-reconstructed (SR) volumes addresses issues like fetal motion
artefacts [10] and low through-plane acquisition resolution, but introduces
additional heterogeneity in texture, tissue contrast, intensity values and
other reconstruction artefacts [11] (see Figure 1). A recent study [12] on the
FeTA 2022 MICCAI challenge revealed significant performance drops in white
matter (WM), gray matter (GM), and ventricles segmentation when models were
tested on diverse clinical datasets, highlighting the impact of domain shifts
on the analysis of fetal MRI.

Figure 1: Domain shifts across data splits in fetal SR MRI. (GA in weeks,
site-SR). A & C - pathological, B & D - neurotypical.

Numerous techniques exist to mitigate domain shifts, including domain
adaptation [13, 14], transfer learning [15], meta-learning [16], style
transfer [17] and data harmonization [18]. However, these methods typically
rely on the availability of at least target domain images which is challenging
in domains with limited data. Other approaches that use multi-centre learning
[19, 20] require multiple training domains which are costly and labour-
intensive to acquire. In fetal brain imaging, where new SR algorithms and MRI
scanners introduce significant diversity [21], available domains may not
sufficiently capture the required variability for generalizable models. Recent
advancements in single-source domain generalization (SSDG) involve techniques
such as global intensity non-linear augmentation (GIN) and causal
interventions [22]. These modifications eliminate spurious correlations,
enhancing model robustness to variations in image intensities and textures. In
a related study [23], authors analyzed frequency’s effect on domain
discrepancy, using a mixed frequency spectrum for self-supervised
augmentation. While effective, limitations may arise in inducing spatial and
intensity transformations, especially with diverse SR algorithms introducing
artefacts associated with skull stripping or the inclusion of extra-cerebral
tissue.

A promising alternative approach involves synthetic data generation based on
segmentation maps [24, 25], achieving domain generalization through domain
randomization [26], only requiring labels but no images. Leveraging shape
information from segmentations, these models introduce diverse spatial and
intensity transformations, along with flexible artefact simulations,
mitigating many factors causing domain shifts in MRI.

In our study, we delve into exploring how synthetic generative models can be
used to construct a diverse fetal brain dataset for training segmentation
models. Our contributions are the following: i) We adapt the domain
randomization of Billot et al. [24] to fetal brain MRI, accommodating specific
fetal anatomical properties, acquisition artefacts and heterogeneity due to
fetal brain development and SR algorithms; ii) We show that our method,
trained only using synthetic data, performs better than models trained using
real data when evaluated out-of-domain and performs on par with state-of-the-
art SSDG algorithms; iii) We extend our evaluation to low-field (0.55T) MRI
data, showing the robustness of our approach to unseen magnetic field strength
and SR algorithms.

##  2 Methodology

###  2.1 Data

Various datasets are used in our experiments to validate the SSDG efficacy of
the models we explore. These datasets come from multiple institutions and were
acquired using MRI scanners from various manufacturers, with different field
strengths, acquisition parameters and reconstructed with different SR
algorithms. The acquisition details are given in Table 1.

FeTA dataset. We used the publicly available data from the MICCAI 2022 FETA
challenge (KISPI) [27, 28]. It consists of 80 subjects, among which are 40
reconstructed using MIALSRTK [29] and 40 using Simple-IRTK [30]. Expert
annotators delineated seven tissue labels (external cerebrospinal fluid
(eCSF), GM, WM, ventricles, cerebellum, deep GM, and brainstem). The ethical
committee of the Canton of Zurich, Switzerland approved the prospective and
retrospective studies that collected and analysed KISPI MRI data (Decision
numbers: 2017-00885, 2016-01019, 2017-00167).

Clinical 1.5T dataset. Additionally, we include a proprietary clinical
dataset, named CHUV, containing 40 subjects reconstructed with MIALSRTK [29]
and manually annotated following the FeTA protocol [27]. Data was
retrospectively collected from acquisitions done between January 2013 to April
2021. All images were anonymized. This dataset is part of a larger research
protocol approved by the ethics committee of the Canton de Vaud (decision
number CER-VD 2021-00124) for re-use of their data for research purposes and
approval for the release of an anonymous dataset for non-medical reproducible
research and open science purposes. This private dataset was used to evaluate
methods submitted for the FETA 2022 challenge and a detailed description of it
can be found in [12].

Clinical 0.55T dataset. We also evaluate our model on 40 neurotypical cases
acquired on a low-field scanner from another centre, referred to as KCL.
Subjects were SR reconstructed twice, using two novel methods, NeSVoR [11] and
SVRTK [31]. There are no available manual segmentations for this dataset.
Fetal MRI was acquired at Kings College London and approved for sharing with
interested academic researchers around the world by the Ethics Committee
London Bromley (Ethics code 21/LO/0742). The data has been acquired during a
prospective single-center study and has been fully anonymised in line with
local procedures

###  2.2 FetalSynthSeg

Our generative model is inspired by SynthSeg [24] which is based on domain
randomization [26]. SynthSeg [24] leverages image segmentation as a structural
prior, integrating randomization across resolution, intensity, contrast, and
spatial distortions. This approach yields a diverse dataset that is well-
suited for training models capable of robustly handling these sources of
variation. We adapt this method to fetal brain segmentation by introducing
some crucial changes related to the tissue generation classes (see Figure 2).

Table 1: Dataset properties. Site | Scanner |  | Acquisition
---
Parameters
SR |  | Res.
---
($mm^{3}$)
| GA
---
$(weeks)$
$N_{n}$ | $N_{p}$
KISPI* |  GE Signa Discovery MR450/MR750 (1.5T/3T) |  ssFSE TR: 2500-3500/120 ms 0.5 x 0.5 x 3.5 $mm^{3}$ | mial | $0.5^{3}$ | 20-34 | 25 | 15
irtk | $0.5^{3}$ | 20-35 | 24 | 16
CHUV* |  | Siemens
---
MAGNETOM
Aera (1.5T)
| HASTE
---
TR/TE: 1200/90 ms
1.1 x 1.1 x 3 $mm^{3}$
mial | $1.1^{3}$ | 21-35 | 25 | 15
KCL |  Siemens MAGNETOM FREE.MAX (0.55T) |  HASTE TR/TE: 2500/106 ms 1.5 x 1.5 x 4.5 $mm^{3}$ | svrtk | $0.8^{3}$ | 21-35 | 40 | 0
nesvor
*FeTA Data [27] |  $N_{n}$: neurotypical, $N_{p}$: pathological
Figure 2: Synthetic image generation framework. Original segmentation labels
are merged to create a 4-meta-label tissue map (CSF, WM, GM, skull). EM
clustering then divides each meta-label into 1 to 4 subclasses, capturing
tissue heterogeneity. A generative model uses these split meta-labels to
produce synthetic images.

Instead of directly using target segmentations as generation classes, our
approach first introduces an intermediate seed generation step. In it, we
initiate synthetic image generation by defining four primary meta-labels: CSF,
WM, GM, and skull with surrounding tissue. Then we employ the
expectation–maximisation (EM) algorithm [32] for intensity clustering within
each meta-class, resulting in 1-4 subclasses per meta-class. In a second step,
these subclasses serve as inputs for synthetic data generators, producing
images and segmentations that faithfully reflect the observed heterogeneity in
SR and those intrinsic to fetal brains, for example, related to WM maturation
[33]. By leveraging tissue-type subclasses rather than original segmentation
labels during synthetic image generation, we mitigate reliance on artificial
contrast disparities between original labels, a critical consideration given
the heterogeneity within a single class. This strategy is similar to the one
used in [24] for cross-domain cardiac MRI segmentation, which showed the
importance of such splitting in representing the target structures with
different levels of heterogeneity. Following subclass creation, voxel
intensities are independently sampled from Gaussian distributions with
randomly selected means and standard deviations. Random artefact corruptions,
including bias field simulation, Gaussian blur and noise addition, are applied
to introduce common noise and artefacts prevalent in SR images. The generative
model applies a battery of spatial transformations (affine and elastic) to
simulate spatial distortions. Table S1 in the Supplementary Material provides
all detailed parameters of the generative model. We use offline image
generation, creating 200 synthetic images per real image to train the model
based on purely synthetic images, which is referred to as FetalSynthSeg (or
_synth as a suffix) further in the text. Code and pre-trained models will be
released upon acceptance of the paper.

###  2.3 Segmentation model

Architecture. Based on Billot et al. [24] and Valabregue et al. [34], our
study employed a 3D U-Net with five levels, featuring instance normalization,
max-pooling, and upsampling operations in the expanding path. Each level
includes a 3 × 3 × 3 kernel convolutional layer with LeakyReLU activation,
except the final layer using softmax activation. Starting with 32 feature
maps, the initial layer doubles and halves after max-pooling and upsampling
layers, respectively. Skip connections facilitate information flow between the
contracting and expanding paths. We use the same architecture and training
hyperparameters across all datasets and splits to ensure comparability.

Pre-processing. During model training, both real and synthetic images undergo
identical pre-processing steps before being fed to the model, including
resampling to 0.5 x 0.5 x 0.5 mm3 with centre crop and crop-padding to
256x256x256 (when needed), random contrast adjustment via gamma
transformation, random affine transformations (scaling, rotation, shearing,
and translation), random Gaussian noise and smoothing with detailed
hyperparameters presented in Supplementary Material Table S1. Subsequently,
image intensities are normalized between 0 and 1 via min-max normalization.

Training. Models are trained using Adam ($\text{LR}=10^{-3}$) on a combination
of Dice and Cross-Entropy losses [34] with a ReduceLROnPlateau scheduler [35]
($\text{factor}=0.1$, $\text{patience}=10$) for up to 500 epochs ($\text{batch
size}=1$). Training halts on persistent validation dice plateau in the last 10
epochs. Internal validation used 5 randomly selected cases per split, ensuring
35 real and 7,000 synthetic images ($200\times 35$) for the training of
baseline and FetalSynthSeg respectively.

###  2.4 Experimental settings

In our experiments, we compare FetalSynthSeg model to i) a baseline model
trained on real images and labels (denoted from now on as baseline), as well
as to ii) fit_nnUnet111The model is available on the FeTA challenge DockerHub
page [12]., FeTA 2022 challenge winner, an ensemble of nnUnet models trained
with all 80 images from KISPI using the GIN SSDG approach [22]. We aim to show
that our model using only synthetically generated images can outperform models
trained on real images when tackling out-of-domain generalization and reach
comparable performance to SSDG SotA domain generalization approaches.

Experiment 1 - High-field generalization. We perform the first domain
generalization experiment by using the FeTA data as well as the clinical 1.5T
dataset. This setting replicates the evaluations carried out in the FeTA
challenge 2022 [12]. We consider three data splits, defined by images acquired
at specific sites and reconstructed with specific SR, namely KISPI - mial,
KISPI - IRTK and CHUV - mial (each containing 35 real images for the baseline
training, and 7000 synthetic images for synth training). Two models are
trained on each split, one using the original FeTA data and labels, and
another using the synthetic data generated from the labels only. This yields a
total of six models (denoted as <Site>_<SR>_base/synth in figures). The models
trained on one (Site, SR) data split are tested on the two other (Site, SR)
pairs. The mean dice score (mdsc) and 95th percentile Hausdorff distance of
each model variant are reported (95th HD in Supplementary Material), and
models are compared using the non-parametric Wilcoxon rank-sum test with
Bonferroni correction. The p-value for statistical significance is set to
0.05.

Experiment 2 - Low-field generalization. We assess our trained model’s
adaptability to new, unseen data with a dataset comprising 40 neurotypical
subjects acquired at KCL on a low-field (0.55T) MRI scanner and reconstructed
either using SVRTK [30] or NeSVoR [11] (none of them used in FeTA data). We
select the top-performing model variant trained on original data and the
counterpart trained on synthetic data, both derived from the same data split.
As no ground truth is available, we evaluate model predictions by comparing
tissue volume growth through gestational age (GA) with FeTA reference data. A
second-order polynomial fit with confidence interval is evaluated.

##  3 Results

####  3.0.1 High-field generalization.

A comparison of the out-of-distribution model predictions is shown in Figure
3. Models trained on synthetic data consistently outperformed those trained on
original images across all out-of-domain testing splits. Statistical tests
confirmed the significance of the differences in mean Dice scores, with
p-values $<$ 0.00005 for all paired comparisons between corresponding baseline
and FetalSynthSeg models trained and tested on the same split. The same trend
is observed with the 95th percentile Hausdorff distances, available in
Supplementary Figure S1. Additionally, we compared the segmentation accuracy
of our model with the FeTA 2022 challenge winners on the CHUV-mial dataset,
which was not used for fit_nnUnet training. We highlight that we reached a
close albeit slightly lower performance to their solution, although our models
were trained on half the amount of data (as they rely on only 35 subjects).

The lowest Dice scores were obtained for pathological cases, particularly
evident in the KISPI-mial split (see Figure 3B). This discrepancy can be
attributed to the fact that approximately 24% of the KISPI-mial dataset
exhibits poor SR quality often in the severe pathological cases, as noted in
[27] and illustrated in Figure 4 (top row). Nonetheless, we noted qualitative
improvements in the FetalSynthSeg model compared to the baseline model, as
illustrated in Figure 4 (and in more detail in Figure S2 in the Supplement).
Differences in skull stripping by IRTK and MIAL algorithms lead to erroneous
segmentation of the skull and surrounding tissue using baseline models.
However, the synthetic model is more robust to SR-induced domain shifts and
artefacts and avoids these errors, even though it was trained on the
segmentations from the same split.

Figure 3: Comparison of the out-of-distribution performance of the
segmentation models (mdsc - mean Dice score across all tissues). (A) baseline
(light) vs FetalSynthSeg (dark). Data split: KISPI-mial (green), KISPI-irtk
(red), CHUV-mial (blue). See Figure S1 from the appendix for a comparison with
in-distribution performance as well as results split by gestational age. The
dashed horizontal line inside the boxplot corresponds to the mean value. (B)
Pathological vs neurotypical, aggregated across all models. Figure 4:  Cross-
domain model inference qualitative results. (A) Model trained on CHUV-mial and
tested on KISPI-irtk/KISPI-mial. (B) Model trained on KISPI-irtk and tested on
KCL-svrtk/KCL-nesvor.

####  3.0.2 Low-field generalization.

Figure 5: Segmented tissue volumes vs GA for KCL data reconstructed with SVRTK
(top row) and NeSVoR (bottom row). KISPI_irtk_base (blue) and KISPI_irtk_synth
(orange) model predictions are compared to FeTA reference values (green) which
are based on the ground truth segmentation of 40 healthy subjects selected
across all splits. Lines are second-order polynomial fit and a corresponding
shaded area is a confidence interval. See Figure S3 in the Supplement for all
tissues evaluation.

Segmented tissue volumes (total brain, GM, WM and CSF volumes) as a function
of GA for KCL-svrtk and KCL-nesvor segmentations are illustrated in Figure 5.
The volume growth curves obtained from SVRTK reconstructions remain within the
confidence interval for both FetalSynthSeg and baseline models. However, a
notable deviation is observed in all estimated tissue volumes for NeSVoR
reconstructions predicted by KISPI-irtk baseline model and to a lesser extent
from our proposed FetalSynthSeg model. This discrepancy highlights a
substantial domain gap within SR algorithms, resulting in an underestimation
of expected tissue volumes that do not occur on a closer domain to the KISPI-
irtk of SVRTK reconstructions. Remarkably, the model trained on synthetic data
demonstrates greater robustness to this domain shift compared to the model
trained on real data and exhibits minimal deviation from the expected values
calculated on the FeTA dataset while qualitatively showing a superior
performance as seen in Figure 4B. The baseline model struggles with correct
tissue segmentation on NeSVoR reconstructions due to out-of-domain appearance,
while performance is improved on SVRTK reconstructions, which exhibit a
smaller SR domain gap with IRTK reconstructions. While cortex topology could
still be more precise, the synthetic model shows consistent qualitative
performance across both scenarios.

##  4 Conclusion

Our study demonstrates that FetalSynthSeg allows robust fetal brain tissue
segmentation across datasets with significant domain shifts. We showed how
strong randomization of spatial and intensity properties during the synthetic
image generation helps models overcome differences caused by MRI acquisition
variations and super-resolution reconstruction. Even with half the data, our
approach achieved performance close to state-of-the-art SSDG models trained
for fetal brain segmentation. Compared to models trained solely on real data,
those trained exclusively on synthetic data showed superior performance,
particularly in cases of novel SR algorithms or images acquired at different
field strengths. The generalization of segmentation models to low-field MRI is
of utmost significance, offering an avenue to enhance fetal MRI accessibility
in underserved cohorts and low-income regions, by providing a cost-effective
diagnostic solution. Our findings suggest that synthetic data can mitigate
performance drops caused by limited data and diverse imaging conditions,
offering promising applications in fields with highly heterogeneous data, such
as fetal imaging.

## References

  * [1] O. M. Benkarim et al., “Toward the automatic quantification of in utero brain development in 3d structural MRI: A review,” Human Brain Mapping, vol. 38, p. 2772–2787, Feb. 2017.
  * [2] A. Jakab et al., “Emerging magnetic resonance imaging techniques in open spina bifida in utero,” European Radiology Experimental, vol. 5, no. 1, 2021.
  * [3] J. Aviles Verdera et al., “Reliability and feasibility of low-field-strength fetal MRI at 0.55 t during pregnancy,” Radiology, vol. 309, no. 1, 2023.
  * [4] M. C. Cortes-Albornoz et al., “MR insights into fetal brain development: what is normal and what is not,” Pediatric Radiology, pp. 1–11, 2024.
  * [5] J. Wu et al., “Age-specific structural fetal brain atlases construction and cortical development quantification for chinese population,” NeuroImage, vol. 241, 2021.
  * [6] C. M. Pfeifer et al., “MRI depiction of fetal brain abnormalities,” Acta Radiologica Open, vol. 8, Dec. 2019.
  * [7] F. Vahedifard et al., “Automatic ventriculomegaly detection in fetal brain MRI: A step-by-step deep learning model for novel 2d-3d linear measurements,” Diagnostics, vol. 13, no. 14, 2023.
  * [8] H. Guan and M. Liu, “Domain adaptation for medical image analysis: A survey,” IEEE Transactions on Biomedical Engineering, no. 3, 2022.
  * [9] R. Lin et al., “Cross-age and cross-site domain shift impacts on deep learning-based white matter fiber estimation in newborn and baby brains,” arXiv preprint arXiv:2312.14773, 2023.
  * [10] T. Sanchez et al., “FetMRQC: Automated quality control for fetal brain MRI,” in Perinatal, Preterm and Paediatric Image Analysis, (Cham), pp. 3–16, Springer Nature Switzerland, 2023.
  * [11] J. Xu et al., “NeSVoR: Implicit neural representation for slice-to-volume reconstruction in MRI,” IEEE Transactions on Medical Imaging, 2023.
  * [12] K. Payette et al., “Multi-center fetal brain tissue annotation (FeTA) challenge 2022 results,” arXiv preprint arXiv:2402.09463, 2024.
  * [13] P. de Dumast et al., “Synthetic magnetic resonance images for domain adaptation: Application to fetal brain tissue segmentation,” in 2022 IEEE 19th International Symposium on Biomedical Imaging (ISBI), IEEE, 2022.
  * [14] Z. Xu et al., “Asc: Appearance and structure consistency for unsupervised domain adaptation in fetal brain mri segmentation,” in Medical Image Computing and Computer Assisted Intervention – MICCAI, 2023.
  * [15] N. Karani et al., “A lifelong learning approach to brain mr segmentation across scanners and protocols,” in MICCAI, 2018.
  * [16] D. Li et al., “Learning to generalize: meta-learning for domain generalization,” in Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence and Thirtieth Innovative Applications of Artificial Intelligence Conference and Eighth AAAI Symposium on Educational Advances in Artificial Intelligence, AAAI Press, 2018.
  * [17] K. Zhou et al., “Domain generalization with mixstyle,” in International Conference on Learning Representations, 2021.
  * [18] F. Hu et al., “Image harmonization: A review of statistical and deep learning methods for removing batch effects and evaluation metrics for effective harmonization,” NeuroImage, vol. 274, 2023.
  * [19] D. Li et al., “Episodic training for domain generalization,” in 2019 IEEE/CVF International Conference on Computer Vision (ICCV), (Los Alamitos, CA, USA), pp. 1446–1455, IEEE Computer Society, 2019.
  * [20] Q. Liu et al., “Shape-aware meta-learning for generalizing prostate MRI segmentation to unseen domains,” in MICCAI, 2020.
  * [21] J. P. Marques et al., “Low-field MRI: An MR physics perspective,” Journal of magnetic resonance imaging, vol. 49, no. 6, pp. 1528–1542, 2019.
  * [22] C. Ouyang et al., “Causality-inspired single-source domain generalization for medical image segmentation,” IEEE Transactions on Medical Imaging, vol. 42, no. 4, p. 1095–1106, 2023.
  * [23] H. Li et al., “Frequency-mixed single-source domain generalization for medical image segmentation,” in Medical Image Computing and Computer Assisted Intervention – MICCAI 2023, Springer Nature Switzerland, 2023.
  * [24] B. Billot et al., “Synthseg: domain randomisation for segmentation of brain scans of any contrast and resolution,” arXiv:2107.09559, 2021.
  * [25] B. Billot et al., “Robust machine learning segmentation for large-scale analysis of heterogeneous clinical brain MRI datasets,” Proceedings of the National Academy of Sciences, vol. 120, no. 9, 2023.
  * [26] J. Tremblay et al., “Training deep networks with synthetic data: Bridging the reality gap by domain randomization,” in Proceedings of the IEEE conference on computer vision and pattern recognition workshops, 2018.
  * [27] K. Payette et al., “An automatic multi-tissue human fetal brain segmentation benchmark using the fetal tissue annotation dataset,” Scientific Data, vol. 8, no. 1, 2021.
  * [28] K. Payette and A. Jakab, “Fetal tissue annotation dataset feta,” 2021.
  * [29] S. Tourbier et al., “Medical-image-analysis-laboratory/mialsuperresol- utiontoolkit: MIAL super-resolution toolkit v2.0.1.” _Zenodo_ , 2020.
  * [30] M. Kuklisova-Murgasova et al., “Reconstruction of fetal brain MRI with intensity matching and complete outlier removal,” Medical image analysis, vol. 16, no. 8, pp. 1550–1564, 2012.
  * [31] A. U. Uus et al., “Automated 3d reconstruction of the fetal thorax in the standard atlas space from motion-corrupted MRI stacks for 21–36 weeks ga range,” Medical image analysis, vol. 80, p. 102484, 2022.
  * [32] A. P. Dempster et al., “Maximum likelihood from incomplete data via the em algorithm,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 39, no. 1, pp. 1–22, 1977.
  * [33] H. Lajous et al., “A fetal brain magnetic resonance acquisition numerical phantom (FaBiAN),” Scientific Reports, vol. 12, May 2022.
  * [34] R. Valabregue et al., “Comprehensive analysis of synthetic learning applied to neonatal brain MRI segmentation,” arXiv:2309.05306, 2023.
  * [35] A. Paszke et al., “PyTorch: An imperative style, high-performance deep learning library,” 2019.

## Supplementary material

Figure S1: Detailed segmentation evaluation. (A) Comparison between models
including in-domain evaluation on the 5 validation cases. (B) Aggregated
performance across all models stratified by GA of subjects. (C) Boxplot
presenting mean 95th-percentile Hausdorff-Distance scores across all tissues
for our experiments on a log scale. (D) Per-model mean dice score results of
models evaluated on KISPI-mial split with the distinction between pathological
and neurotypical cases. Table S1: Hyperparameters of the synthetic generator
and augmentations. Only parameters deviating from the default ones used by the
fully randomized generative model without any tissue priors [24] are reported.
Intensities are in $[0,255]$ range, rotations in degrees and spatial
parameters in mm. Synthetic Data Generator Hyperparameters
---
$a_{\text{sc}}$ | 0.9 | $a_{\text{tr}}$ | -10 | $r_{\text{HR}}$ | 0.5 |  |
$b_{\text{sc}}$ | 1.1 | $b_{\text{tr}}$ | 10 | $b_{\text{res}}$ | 0.5 |  |
Augmentations Hyperparameters-Probability
$\gamma_{\text{range}-0.5}$ | 0.5 – 1.5 | $\text{scale}_{\text{range}-0.5}$ | -0.1 – 0.1 | $\sigma_{\text{noise}-0.5}$ | 0.1 | $\mu_{\text{noise}-0.5}$ | 0
$\text{rotation}_{\text{range}-0.5}$ | -0.2 – 0.2 | $\text{shear}_{\text{range}-0.5}$ | -0.1 – 0.1 | $\sigma_{\text{smooth}-0.7}$ | 0.5 – 1.5 |  |
Figure S2: Cross-domain model inference qualitative results Figure S3:
Segmented tissue volumes vs GA for KCL data for all tissues.