# Continuous Coupling of Ultracold Atoms to an Ionic Plasma via Rydberg
Excitation

T.M. Weber<EMAIL_ADDRESS>T. Niederprüm  T. Manthey  P. Langer  V.
Guarrera  G. Barontini  H. Ott  Research Center OPTIMAS, Technische
Universität Kaiserslautern, 67663 Kaiserslautern, Germany

###### Abstract

We characterize the two-photon excitation of an ultracold gas of Rubidium
atoms to Rydberg states analysing the induced atomic losses from an optical
dipole trap. Extending the duration of the Rydberg excitation to several ms,
the ground state atoms are continuously coupled to the formed positively
charged plasma. In this regime we measure the $n$-dependence of the blockade
effect and we characterise the interaction of the excited states and the
ground state with the plasma. We also investigate the influence of the quasi-
electrostatic trapping potential on the system, confirming the validity of the
ponderomotive model for states with $20\leq n\leq 120$.

Several proposals have demonstrated that dressing ultracold atoms with highly
excited Rydberg states can be an extremely powerful tool to tune the
interactions among them 1 ; 2 ; 3 . In particular long-range interactions
would lead to novel quantum phases of matter as super-solid phases or dipolar
crystals, once the dressed atoms are loaded into optical lattices 4 ; 5 . With
respect to the regime of frozen Rydberg gases, where all the relevant physics
can take place within a few $\mu$s 6 ; 7 ; 8 , the onset of such exotic phases
would require a much longer timescale, during which the coherent excitation
must be preserved. A limitation in this case comes from the tendency of
Rydberg gases to spontaneously evolve into a plasma 9 ; 10 . Hence a complete
understanding of this process is fundamental in order to find the appropriate
strategies to preserve the dressed states for sufficiently long times. In this
paper we report on the characterisation of the dynamics originating from the
two-photon excitation of an optically trapped ultracold gas of neutral atoms
to Rydberg states. We implement excitation pulses of at least a few tens of
$\mu$s, i.e., sufficiently long to ensure that the Rydberg gas spontaneously
evolves into a plasma. We exploit the excitation itself to continuously couple
the neutral ultracold atoms to the plasma. We measure the strength and the
shape of the resulting resonance lines varying the final state $nl$, with
$20\leq n\leq 120$ and with $l=0,2$. Moreover, we study the dependence of the
atom-plasma dynamics on different excitation times and trapping potentials.
The comparison of our results with those obtained with a rate equation model
allows us to fully characterise the system and to highlight the role of the
plasma-induced blockade and of the trapping potential.

Our experimental sequence starts by cooling room-temperature vapours of 87Rb
atoms in a 2D-MOT. The atoms are subsequently transferred to the main chamber
by a resonant push beam, loading a 3D-MOT with a rate of $2\times
10^{7}\,$atoms/s, thus allowing to collect $10^{8}$ atoms after $5\,$s. A
single 15 W beam at 10.6$\,\mu$m, produced by a commercial cw CO2 laser (model
C55, Coherent Inc.), is focussed at the center of the 3D-MOT to a waist of
30$\,\mu$m, creating an optical dipole trap. After a dark-MOT phase of 70 ms,
$1.2\times 10^{6}$ atoms at 60$\,\mu$K are left in the trap. By ramping down
the power of the trapping beam to 180 mW in about 6 s, we drive forced
evaporation, ending up with $4\times 10^{4}$ atoms in a cigar shaped spinorial
Bose-Einstein condensate in the $\left|5S_{1/2},F=1\right\rangle$ ground
state. The final trapping frequencies are $\nu=(205,205,17)\,$Hz. Unless
explicitly stated we stop the evaporation at a power of 220 mW obtaining a
thermal cloud of $10^{5}$ atoms at 250 nK with a density of $\cong
10^{14}\,$cm-3.

Figure 1: (Color online) Typical resonance lineshape for a thermal cloud at
250 nK. The atoms are excited to the $27S_{1/2}$ state for 1 ms. The solid
line is a Gaussian fit to the data. The inset shows the two photon Rydberg
excitation scheme that we employ.

After the end of the evaporation we drive the transition from the
$\left|5S_{1/2},F=1\right\rangle$ state to the selected $n$ Rydberg state,
using the two-photon scheme depicted in the inset of Fig. 1. The infra-red
(IR) light at 780 nm is locked to the
$\left|5S_{1/2},F=2\right\rangle\rightarrow\left|5P_{3/2},F=3\right\rangle$
transition. In this way it is always red detuned by $\Delta\simeq 2\pi\times
6.8\,$GHz from the $|5S_{1/2},F=1\rangle\rightarrow|5P_{3/2}\rangle$
transition, thus reducing the single photon scattering rate and allowing for a
longer lifetime of the trapped atoms. The IR beam is collimated into the cold
cloud perpendicularly to the CO2 laser beam with a waist of 1 mm and a typical
power of $\simeq 15\,$mW. The blue light for the second excitation step is
generated by frequency doubling the light produced by a diode laser based MOPA
system using a LBO crystal in a bow-tie ring cavity. The master laser can be
locked in a range of about 20 nm exploiting a combined transfer cavity and
offset locking scheme, generating light from 479 nm to 488 nm. Once combined
with the IR light, it provides the possibility to excite Rydberg states from
$n=20$ to $n=150$. The actual frequency of the beam is determined with a
standard EIT spectroscopy technique in a glass cell, with an accuracy of $\pm$
5 MHz. The blue beam is sent on the atoms collinearly with the IR one and
focussed at the center of the atom cloud with a waist of $40\,\mu$m and a
typical power of $\simeq 120\,$mW. The applied intensities correspond to Rabi
frequencies of $\Omega_{1}=2\pi\times 145\,$MHz and $\Omega_{2}=2\pi\times
44\,$MHz for the first and second excitation steps to the $27S_{1/2}$ state,
which result in an effective overall Rabi frequency of
$\Omega=\Omega_{1}\Omega_{2}/2\Delta=2\pi\times 470\,$kHz. The blue and the IR
light are switched on together with variable pulse lengths, typically on the
order of 1 ms. The excitation of Rydberg atoms is revealed by the reduction of
the number of trapped atoms at the end of the two-photon pulse, when scanning
the frequency of the blue beam across the atomic transition. A typical line
shape is shown in Fig. 1.

Figure 2: (Color online) Measured resonance widths for different $ns$ states.
Data (squares) are compared with the solution of the rate equation model
explained in the text (circles). In the inset the same quantities are shown
for $nd$ states. Note, that for $n<27$ the power of the blue excitation laser
significantly drops.

In Fig. 2 we show the measured widths of the resonances as a function of the
principal quantum number $n$. Notably, they are up to six orders of magnitude
larger than the natural linewidths ($\simeq$ 80 kHz for $n=27$ and $\simeq$
300 Hz for $n=120$ 11 ). Moreover they show a non-trivial dependence on $n$.
These features cannot be attributed to Doppler or saturation effects since,
for the given experimental parameters, they are in the order of 100 kHz.
Furthermore it has been repeatedly demonstrated that the bulk excitation of
Rydberg atoms in the ultracold regime rapidly leads to a series of secondary
effects like fast ionization by collisions or radiation 9 ; 14 , production of
plasmas 21 ; 22 ; 23  or blockade effects 24 ; 33 . In our case, on the
timescale of the excitation this secondary effects are certainly present and a
complete understanding of the observed features necessarily requires to take
into account all of them. For these reasons we analyse our data starting from
the following rate equation model:

$\displaystyle\dot{N}_{g}$ $\displaystyle=$ $\displaystyle
W_{gr}\xi(N_{r}-N_{g})+W_{ge}(N_{e}-N_{g})$ $\displaystyle-$ $\displaystyle
K_{HM}N_{g}N_{r}+K_{PI}N_{e}N_{r}+\Gamma_{e}N_{e}-\gamma_{g}N_{g}$
$\displaystyle\dot{N}_{e}$ $\displaystyle=$ $\displaystyle
W_{ge}(N_{g}-N_{e})+W_{er}(N_{r}-N_{e})$ $\displaystyle-$ $\displaystyle
K_{PI}N_{e}N_{r}-\Gamma_{e}N_{e}$ $\displaystyle\dot{N}_{r}$ $\displaystyle=$
$\displaystyle W_{gr}\xi(N_{g}-N_{r})+W_{er}(N_{e}-N_{r})$ $\displaystyle-$
$\displaystyle K_{BB}N_{r}-K_{RR}N_{r}-K_{HM}N_{g}N_{r}$ $\displaystyle-$
$\displaystyle K_{PI}N_{e}N_{r}-\Gamma_{r}N_{r}$ $\displaystyle\dot{N}_{i}$
$\displaystyle=$ $\displaystyle
K_{BB}N_{r}+\frac{1}{2}K_{RR}N_{r}+\frac{1}{2}K_{HM}N_{g}N_{r}$ (1)
$\displaystyle+$ $\displaystyle K_{PI}N_{e}N_{r}-\chi N_{i},$

where $g$ labels the ground state atoms, $e$ the 5p state atoms, $r$ the
Rydberg atoms and $i$ the ions. With $\Gamma_{e}$ and $\Gamma_{r}$ we indicate
the natural linewidths of the corresponding levels. The loss rate
$\gamma_{g}$, which is due to the high scattering rate from the IR laser, is
measured from the decay of the atom number in dependence of the excitation
pulse time when the blue laser is far detuned from the transition
($\delta\simeq 100\,$MHz), as shown in Fig. 3a. The coefficients $W_{ij}$
represent the excitation rates and, in our experimental regime, they are
$\simeq 2\pi\times$300 Hz for the $g\rightarrow e$ transition while they range
between a few Hz to a few mHz for the transitions $e\rightarrow r$, with
$20<n<120$. We calculate the two-photon coupling between the ground state and
the Rydberg states using the standard formula from second-order perturbation
theory 26 :

$\displaystyle W_{gr}$ $\displaystyle=$
$\displaystyle\frac{8I_{1}I_{2}}{\hbar^{4}c^{2}{\varepsilon_{0}}^{2}}{\left|\sum_{k}\frac{\widehat{d_{3k}}\widehat{d_{kg}}}{\omega_{kg}-\omega_{1}}+\frac{\widehat{d_{ek}}\widehat{d_{kg}}}{\omega_{ek}-\omega_{2}}\right|}^{2}$
(2) $\displaystyle\frac{\Gamma_{r}/2}{(\Gamma_{r}/2)^{2}+\Delta\omega^{2}},$

where $k$ labels every intermediate state, $\widehat{d_{ji}}$ are the dipole
matrix elements, $I_{i}$ the laser intensities and $\Delta\omega$ the overall
laser linewidth. We calculate that $W_{gr}$ ranges between a few MHz ($n$=20)
and a few hundred Hz ($n$=120). All the remaining processes are related to the
production of ions or to the interactions with them. When exciting Rydberg
atoms ions can be produced in several ways: in our experimental circumstances
the main channels are blackbody radiation (BB), Hornbeck-Molnar ionisation
(HM), Penning ionisation (PI) and Rydberg-Rydberg collisions (RR). At 300 K
the rate $K_{BB}$ never exceeds 500 s-1 for all the states that we excite 16 .
The HM ionisation is due to collisions between Rydberg atoms and ground state
atoms that produce one Rb${}_{2}^{+}$ ion and one electron. More complicated
is the PI channel, in which one Rydberg atom collides with one atom in the 5p
state producing one ion, one electron and one atom in the ground state. For
atomic densities in the order of $10^{14}-10^{15}\,$cm-3, the corresponding
rates are $K_{HM}\simeq 10^{5}-10^{6}\,$s-1 and $K_{PI}\simeq
10^{6}-10^{7}\,$s-1 12 ; 13 . The last ionisation channel is characterized by
the well-known $n^{4}$ dependence of the cross-section through the rate
$K_{RR}=\rho_{r}v(\pi a_{0}^{2}n^{4})$, where $\rho_{r}$ is the density of the
Rydberg atoms and $v$ their relative velocity 17 . Every time one ion is
produced, the corresponding electron leaves the trapping region extremely fast
leaving an excess positive charge around the atoms. The positively charged
plasma that originates from the continuous production of ions is then
subjected to the so-called Coulomb explosion: the ions repel each other via
the electrostatic force. The complete description of such a complicated
process lies beyond the purposes of this work and we model the expansion of
the plasma with an effective ion loss rate from the trapping volume
$\chi=v/\sigma+\sqrt{2e^{2}\rho_{i}/(4\pi\varepsilon_{0}m)}/\sigma$, where
$v=\sqrt{3k_{B}T/m}$ is the thermal velocity which we suppose to be the atomic
one, $\rho_{i}$ the plasma density and $\sigma$ the radius of the atomic
distribution. Moreover, the plasma that forms and expands produces an electric
field in the region of the trapped atoms that significantly shifts the Rydberg
levels. The most striking consequence is the blockade effect, i.e., an
effective reduction of the coupling between the ground state and the Rydberg
state. We model this process introducing the coefficient
$\xi=\Delta\omega/(\Delta\omega+B)$ in eqs. (1). Indeed, the spacially varying
electric field produced by the space charge of the plasma induces a broadening
of the Rydberg line that is
$B=(2\sigma\rho_{i}/(4\pi\varepsilon_{0}))^{2}\alpha/(2\hbar)$, where $\alpha$
is the $n$-dependent atomic polarisability given by
$\alpha=h\times(2.202\times 10^{-7}n^{6}+5.53\times 10^{-9}n^{7})\,$Cm2V-1 32
. We find the linewidths calculated solving the rate equation model to be in
excellent agreement with the observed ones (Fig. 2). The simulated dynamics of
the number of trapped atoms, Rydberg atoms and ions are shown in Fig. 3b for
the $27s$ state. Through the competition between the Coulomb expansion and the
blockade effect the atom-plasma system rapidly evolves into a self-balanced
situation where the feeding rate is continuously adjusted in order to
compensate the losses of ions. In practice, for a given value of the principal
quantum number $n$, the number of Rydberg atoms that are on average present in
the volume is kept almost constant. Unfortunately this self-balancing effect
is strongly altered by the huge losses induced by the scattering of photons
from the IR laser (Fig. 3a) which can not be balanced. Due to this effect the
self-balanced phase has only a limited lifetime. From the solution of the rate
equation model we can determine the number of Rydberg atoms $N_{r}$ that are
present in the volume during this reduced lifetime for different $n$ states,
as reported in Fig. 3c. For the given trapping volume the maximum number of
Rydberg atoms excited corresponds to an average distance
$r=1/\sqrt[3]{\rho_{r}}$ that ranges from 1$\,\mu$m to 13$\,\mu$m for $27\leq
n\leq 120$. In Fig. 3d the $n$-dependence of the maximum number of ions
present at a time is shown.

Figure 3: (Color online) a) Typical decay measurements in (downpointing
triangles) and out of resonance (uppointing triangles) together with the
curves obtained as solutions of the rate equation system (2) for the $55s$
state. b) Time evolution of total number of trapped atoms (solid line),
Rydberg atoms (dashed) and ions (dashed-dotted) for the $27s$ state. c)
Calculated average number of Rydberg atoms present as a function of $n$. d)
Calculated maximum number of ions present as a function of $n$.

Finally we investigate the influence of the trapping potential on the Rydberg
atoms and on the plasma dynamics. The CO2 laser is expected to create a
repulsive ponderomotive potential for any Rydberg atom 15  while for the atoms
in the $5S_{1/2}$ and $5P_{3/2}$ states it creates the attractive trapping
potential. We first verify the reliability of the ponderomotive assumption for
the dipole potential, recording the atom losses for different trapping powers
and measuring the relative shifts of the center of the resonances, as reported
in Fig.4a. We observe that the AC Stark shift is effectively the same for
every Rydberg state and that it is compatible with the expected theoretical
value of 179 MHz/W, as shown in Fig.4b.

Figure 4: (Color online) a) Typical resonance lineshapes for different powers
of the CO2 laser for the $27S_{1/2}$ state b) Measured AC Stark shifts for
different values of $n$. The solid line is the expected state independent
theoretical value of 179 MHz/W. c) Resonance lineshapes for $n=27$ for a BEC
and a thermal cloud. d) Measurement of the resonance width as a function of
the CO2 power for $n=43$.

It has been demonstrated that the plasma dynamics is not directly affected by
the density of the ground state atoms since the interaction is mediated mainly
by the Rydberg atoms 9 . We have verified this, measuring the linewidths for a
BEC and a thermal cloud when the density differs by one order of magnitude
while the trapping potential remains almost the same. As can be seen in
Fig.4c, there is only a minimal difference between them. However, as reported
in Fig.4d and as shown in Fig.4a we do observe a dependence of the linewidth
on the power of the trapping laser. This effect is due to the fact that an
increase of the power of the CO2 laser compresses the atoms in the ground
state reducing the excitation volume. In a smaller volume the ionic blockade
effects are even stronger producing a more pronounced broadening and a further
detriment of the transition probability.

In summary, we have reported the two-photon excitation and subsequent fast
ionization of Rydberg atoms in a unprecedented wide range of $n$ states.
Elongating the excitation time to a few ms we have systematically investigated
the interplay between the ultracold atoms and the formed plasma. We have
carried out a detailed yet simple analysis of the system, highlighting the
important role of the plasma on hampering the Rydberg excitation. We have
finally characterized the influence of a quasi-electrostatic trap on Rydberg
atoms and plasma dynamics.

Our results have a direct impact on the schemes that aim at the dressing of
ground state atoms with Rydberg states. Indeed the possible evolution of
Rydberg excited samples into plasma must be taken into account, since it leads
to a huge loss in the transition probability. A possible way to avoid this
dynamics is the application of electric fields to remove the produced ions
fastly. However, in the presence of electric fields Rydberg states are
strongly mixed thus leading to population distribution over the neighbouring
states 31 , and hence to loss of coherence. This, together with our findings,
suggests that the best way to obtain long-living plasma-free Rydberg-dressed
samples would require the use of high-intensity lasers in order to reach the
condition where $\delta$ is larger than the plasma-induced broadening with a
decent amount of admixture.

###### Acknowledgements.

We acknowledge financial support by the DFG within the SFB/TRR 49 and GRK 792.
V. G. and G. B. are supported by Marie Curie Intra-European Fellowships. It is
a pleasure to thank P. Pillet. T. Pohl and H. Hotop for enlightening
discussions. We are grateful to A. Widera for technical support.

## References

  * (1) J. E. Johnson and S. L. Rolston, Phys. Rev. A 82 033412 (2010)
  * (2) G. Pupillo, A. Micheli, M. Boninsegni, I. Lesanovsky, P. Zoller, Phys Rev Lett. 104, 223002 (2010).
  * (3) J. Honer, H. Weimer, T. Pfau and H. P. Büchler, Phys. Rev. Lett. 105, 160404 (2010).
  * (4) N. Henkel, R. Nath and T. Pohl, Phys. Rev. Lett. 104 195302 (2010)
  * (5) F. Cinti, P. Jain, M. Boninsegni, A. Micheli, P. Zoller and G. Pupillo, Phys. Rev. Lett. 105 135301 (2010)
  * (6) M. Saffman, T.G. Walker, K. Mølmer, Rev. Mod. Phys. 82 2313 (2010).
  * (7) I. Mourachko et al, Phys. Rev. Lett. 80 253 (1998).
  * (8) W. R. Anderson, J. R. Veale and T. F. Gallagher, Phys. Rev. Lett. 80 2 (1997).
  * (9) G. Vitrant, J. M. Raimond, M. Gross and S. Haroche, J. Phys. B: At. Mol. Phys. 15 L49-L55 (1982).
  * (10) M. P. Robinson, B. Laburthe-Tolra, M. W. Noel, T. F. Gallagher and P. Pillet, Phys Rev Lett. 85, 4466 (2000).
  * (11) Einstein-A-coefficients and dipole matrix elements are calculated by solving the time independent Schrödinger equation with wavefunctions derived from Numerov method and energy levels given by quantum defect theory.
  * (12) L. Barbier and M. Chéret, J. Phys. B: At. Mol. Phys. 20 1229-1248 (1987)
  * (13) A. Kumar, B. C. Saha, C. A. Weatherford, S. K. Verma, J. Mol. Structure (Theochem) 487 (1999)
  * (14) R. M. Potvliege and C. S. Adams, New J. Phys. 8 163 (2006).
  * (15) F. Markert, P. Würtz, A. Koglbauer, T. Gericke, A. Vogler and H. Ott, New J. Phys. 12 113003 (2010).
  * (16) I. I. Beterov, D. B. Tretyakov and I. I. Ryabtsev and A. Ekers and N. N. Bezuglov, Phys. Rev. A 75 052720 (2007)
  * (17) T. F. Gallagher, Rydberg Atoms (Cambridge University Press, 1994).
  * (18) T. C. Gallagher, P. Pillet, M. P. Robinson, B. Laburthe-Tolra and M. W. Noel, JOSA B 20 5 (2003).
  * (19) T. C. Killian, M. J. Lim, S. Kulin, R. Dumke, S. D. Bergeson and S. D. Rolston, Phys. Rev. Lett. 86 17 (2000).
  * (20) S. D. Bergeson and F. Robicheaux, Phys. Rev. Lett. 101 073202 (2008).
  * (21) G. Bannasch and T. Pohl, Phys. Rev. A 84 052710 (2011).
  * (22) T. C. Killian, T. Pattard, T. Pohl and J. M. Rost, Phys. Rep. 449 77-130 (2007)
  * (23) G. Grynberg and B. Cagnac, Rep. Prog. Phys. 40 791-841 (1977).
  * (24) A. Reinhard, T. Cubel Liebisch, B. Kuffman and G. Raithel, Phys. Rev. A 75 032712 (2007).
  * (25) M. S. O’Sullivan and B. P. Stoicheff, Phys. Rev. A 31 4 (1984)


# Automatic Video Colorization using 3D Conditional Generative Adversarial
Networks

###### Abstract

In this work, we present a method for automatic colorization of grayscale
videos. The core of the method is a Generative Adversarial Network that is
trained and tested on sequences of frames in a sliding window manner. Network
convolutional and deconvolutional layers are three-dimensional, with frame
height, width and time as the dimensions taken into account. Multiple
chrominance estimates per frame are aggregated and combined with available
luminance information to recreate a colored sequence. Colorization trials are
run succesfully on a dataset of old black-and-white films. The usefulness of
our method is also validated with numerical results, computed with a newly
proposed metric that measures colorization consistency over a frame sequence.

Index Terms—  Video Colorization, Generative Adversarial Networks, three-
dimensional convolution, black-and-white films

††We gratefully acknowledge the support of NVIDIA Corporation with the
donation of the Titan XP GPU used for this research.
---
Fig. 1:  Architecture of the proposed model. Sequences of luminance frames
(’L’ channels) are fed to the generator network, which creates chrominance
channels for the corresponding frames (’a’, ’b’ channels). The discriminator
network quantifies to what extend the generated chrominance sequence
corresponds to a plausible colorization. As inputs and outputs are frame
sequences, all convolutions are three-dimensional. During training, model
weights are estimated through optimization of a GAN loss that effectively
combines these constraints: how plausible the generated colorization is per
frame _and as a sequence_ , and how close the chrominance estimate is to the
ground truth. In test time, the video is processed in a sliding window manner,
producing $C$ colorization estimates (the size of the sliding window) for each
frame, which are then aggregated to obtain a single estimate per frame. The
output chrominance is finally combined with the luminance input in order to
recreate a sequence of colorized frames.

##  1 Introduction

In this paper, we address the problem of automatic colorization of monochrome
digitized videos [1, 2, 3, 4, 5]. Perhaps the most straightforward practical
application is to colorizing black-and-white footage from old films or
documentaries. Video compression is another possible application of note [5].

Video colorization methods can be categorized according to the level of user
interaction required. A group of methods assume that a partially colored frame
exists in the video, where color has been manually annotated in the form of
color seeds [5, 6, 7]. The method must then propagate color from these seeds
to the rest of the frame, then to other frames in the video. Other methods
assume instead that a reference colored image exists that is similar in
content and structure to the target monochrome video frames [3, 4, 8]. These
methods may or may not require user intervention; for example, in [9] the user
can specify matching areas between the reference and the target frames. In
reference image-based methods, the problem of video colorization is hence
converted to the problem of how to propagate color from the reference frame to
other frames and/or from frame to frame. Optical flow estimation has been used
to guide frame-to-frame color propagation [2, 3]. In [7], Gabor feature flow
is used as alternative to standard optical flow as a more robust guide to
color propagation. Naturally, methods of this vein work best for coloring
short videos or frames coming from the same scene [3].

In the present work, we propose a learning-based method for video
colorization. As such, we assume that a collection of colored frames exist,
that will be used to train the model. In particular, the proposed method is
based on an appropriately designed Generative Adversarial Network (GAN) [10].
GANs have gained a fair amount of traction in the last few years. Despite
their being harder to train even more than standard neural networks, requiring
the employment of various heuristics and careful choosing of hyperparameters
to attain convergence to a Nash equilibrium [11, 12], they have proven to be
excellent generative models. The proposed model employs a conditional GAN
(cGAN) architecture, popularized by the pix2pix model [13]. In the current
work, convolutional and deconvolutional layers are 3D (height, width, time
dimensions) to accomodate for the sequential nature of video data.

The main novel points of the current paper are as follows: (a) we present a
model for learning-based automatic video colorization that can take advantage
of the sequential nature of video, while avoiding the use of frame-by-frame
color propagation techniques that come with their own inherent limitations
(typically they require existing colored key frames and/or are practically
applicable within a single shot). Other recent works use learning methods to
color video via propagation [1], or via frame-by-frame image colorization,
with each frame processed separately [14]; (b) we elaborate on the issue of
video colorization evaluation and propose a quantitative colorization metric
specifically for video; (c) we show that the proposed method creates
colorization models that are transferable, in the sense that learning over a
particular frame sequence produces a plausible output usable on a sequence of
different content.

The rest of the paper is organized as follows. In section 2, we briefly
discuss preliminaries on adversarial nets and present the architecture and
processing pipeline of the proposed video colorization method. In section 3,
we elaborate on existing numerical evaluation methods and propose a new metric
to evaluate video colorization. In section 4, we show numerical and
qualitative results of our method, tested on a collection of old films. We
close the paper with section 5, where we discuss conclusions and future work.

Grayscale | Ground Truth | Proposed
---|---|---
|  |
|  |
|  |
|  |
Fig. 2:  Colorization results using our method. Depicted frames are samples from films: “Et Dieu..créa la femme” and “Tzéni, Tzéni” (2 top and 2 bottom rows respectively).  2D cGAN | Proposed
---|---
|
|
Fig. 3:  Comparison of proposed model vs non-sequential 2D cGAN model. The
proposed model produces better results than the non-sequential variant, as the
former can take advantage of optical flow information, with its 3D
convolution/deconvolution layers and estimate aggregation scheme. (Note for
example how each method colorizes the hand of the standing actor on the top
frame, or the color of the suit on the bottom frame). Depicted frames are
samples from the film “Dial M for Murder”.

##  2 Proposed method

The proposed method assumes the existence of a training set consisting of a
sequence of $N_{train}$ colored frames, and a test set consisting of a
sequence of $N_{test}$ monochrome frames that are to be colorized. During the
training phase, a cGAN model is used to learn how to color batches of $C$
ordered frames. Hyperparameter $C$ is fixed beforehand with
$C<<min\\{N_{train},N_{test}\\}$. During the testing phase, the model is run
on the input monochrome video in a sliding window manner. Windows are
overlapping and move by a single frame at a time, thereby producing a set of
$C$ colorization estimates _for each_ monochrome frame, hence $C$ video
colorization proposals. These $C$ estimates are then combined to produce a
single colorized output. In what follows we discuss the details of this
process.

A GAN is a generative, neural network-based model that consists of two
components, the generator network and the discriminator network. The cGAN
architecture [13] that is employed as part of the proposed method, is a
supervised variant of the original unsupervised GAN [10]. A cGAN learns a
mapping from observed input $x$ to target output $y$ 111Other variants of a
cGAN are possible; for example, a noise variable $z$ could be added to produce
a non-deterministic output [13]. We employ a deterministic cGAN variant in
this work..

Formally, the objective to be optimized is:

$\arg\min_{G}\max_{D}[\operatorname{\mathbb{E}}_{y}[logD(y)]+$
$+\operatorname{\mathbb{E}}_{x}[log(1-D(G(x)))]]+\lambda\operatorname{\mathbb{E}}_{x,y}[\left\lVert
y-G(x)\right\rVert_{1}],$ (1)

with hyperparameter $\lambda$ controlling the trade-off between the GAN
(discriminator) loss and the $L_{1}$ loss. $D(\cdot)$ and $G(\cdot)$
correspond to the discriminator and generator respectively. The GAN loss
quantifies how plausible the colorization output, while the $L_{1}$ loss
forces the colorization to be close to the ground truth. We use
representations in the CIE Lab color space (following e.g.[15]). For a
monochrome frame sequence, only luminance is known beforehand. Input is a
sequence of luminance channels (channel $L$) of C consecutive frames $x\in$
${\rm I\\!R}^{H\times W\times C\times 1}$, and the objective is to learn a
mapping from luminance to chrominance (channels $a$,$b$) $y\in$ ${\rm
I\\!R}^{H\times W\times C\times 2}$ where H, W are frame dimensions.

The generator network is comprised of a series of convolutional and
deconvolutional layers. Skip connections are added in the manner introduced by
UNet [16]. As inputs and outputs are sequences of fixed-size frames, all
convolutions and deconvolutions are three-dimensional (frame height, width and
time dimensions). The encoder and decoder stacks comprise $8$ strided
convolutional/deconvolutional layers each (stride=$2$), followed iteratively
by batch normalization (BN) and rectified linear unit (ReLU) activation
layers. Following [17], outputs are forced to lie in the $(-1,1)$ range with a
tanh activation layer at the end of the generator network, and only later
renormalized to valid $a,b$ chrominance values. The discriminator network is a
3D convolutional network comprising $5$ convolutional layers iteratively
followed by BN and ReLU layers. The discriminator is topped by a fully
connected (“dense”) layer and a sigmoid activation unit in order to map the
image to a real/fake probability figure.

At test time, we use the generator in a sliding window fashion over the
footage to be colorized. Hence, each frame is given as input to the generator
at a total of $C$ times, since $C$ is the size of the sliding window. The
produced $C$ chrominance estimates $\chi_{1},\chi_{2},\cdots,\chi_{C}$
222$\chi_{i}$ denotes the $i^{th}$ colorization estimate for a frame. $y$
denotes a colorization estimate for a sequence of $C$ frames. then need to be
used to produce a single estimate $\hat{\chi}$. We can write $\hat{\chi}$ as a
maximum-a-posteriori (MAP) estimate as:

$\hat{\chi}=\arg\max_{\chi}p(\chi|\chi_{1},\cdots,\chi_{C})$ (2)

where a prior distribution $p(\chi)$ can be assumed over possible $a,b$ values
in order to favor a particular chrominance setup. If identical distributions
centered around each $\chi_{i}$ and an uninformative prior is used, the above
formula simplifies as an average over all chrominance values per frame pixel:
$\hat{\chi}=\nicefrac{{1}}{{C}}\sum_{i=1}^{C}\chi_{i}$. Finally, the
chrominance estimate is recombined with input luminance to recreate colored
RGB frames for the input video. The architecture of the proposed model is
summarized in figure 1.

##  3 Metrics for numerical evaluation of video colorization

In this section we describe the metrics we use for numerical evaluation of
video colorization. We use two metrics that measure per-frame colorization
quality, also usable in single-image colorization. Furthermore, we propose a
new metric suitable for video colorization in particular.

_Peak Signal-to-Noise Ratio (PSNR):_ PSNR is calculated per each test frame in
the RGB colorspace, and their mean is reported as a benchmark over the whole
video.

_Raw Accuracy (RA):_ Raw Accuracy, used in [15] to evaluate image
colorization, is defined in terms of accuracy of predicted colors over a
varying threshold. Colors are classified as correctly predicted if their
Euclidean distance in the $ab$ space is lower than a threshold. Accuracy is
computed over color values for every pixel position and frame. Integrating
over the curve that is produced by taking into account varying threshold
yields the RA metric. We integrated from $0$ to $150$ distance units as in
[15].

_Color Consistency (CC):_ The aforementioned metrics measure strictly the
quality of colorization of each frame separately. We propose and use a metric
to measure both per-frame quality and also the consistency of the choice of
colors between consecutive frames. Such a metric can, for example, penalize
erratic differences in colorization from frame to frame, that would otherwise
be “invisible” to the other metrics, borrowed from single image
restoration/colorization. We define color consistency over sets of two
consecutive colorization predictions $\hat{\chi}^{(t)},\hat{\chi}^{(t+1)}$ and
corresponding ground truth values $\chi^{(t)},\chi^{(t+1)}$ as

$CC^{(t,t+1)}=\nicefrac{{1}}{{HW}}\sum_{i=1}^{H}\sum_{j=1}^{W}\nicefrac{{1}}{{2}}[A^{(t)}_{ij}+A^{(t+1)}_{ij}]A^{(t\times
t+1)}_{ij}$ (3)

where affinity matrices $A^{(t)}$ and $A^{(t\times t+1)}$ are defined as

$A_{ij}^{(t)}=\phi(||\chi_{ij}^{(t)}-\hat{\chi}_{ij}^{(t)}||),$
$A_{ij}^{(t\times
t+1)}=\phi(\left\lVert||\chi_{ij}^{(t)}-\chi_{ij}^{(t+1)}||-||\hat{\chi}_{ij}^{(t)}-\hat{\chi}_{ij}^{(t+1)}||\right\rVert),$

with function $\phi(\cdot)$ a positive, strictly decreasing function that is
used to convert distances to similarities. We use $\phi_{ij}(X)=\lfloor
60X_{ij}/(\max(X)+\epsilon)+1\rfloor^{-1}$. Total CC over a video sequence is
calculated as the average CC over all consecutive frames. Higher values
correspond to better results.

##  4 Experiments

We have tested our method over a collection of old films: $(a)$ “Dial M for
Murder” (USA, 1954; 63,243 frames) [18] $(b)$ “Et Dieu..créa la femme”
(France, 1956; 54,922 frames) [19] $(c)$ “Tzéni, Tzéni” (Greece, 1965; 58,932
frames) [20] $(d)$ “A streetcar named desire” (USA, 1951; 18,002 frames) [21]
$(e)$ “Twelve angry men” (USA, 1957; 12,000 frames) [22]. Frames were sampled
off these films at $10$ fps. Films $(a),(b),(c)$ are colored, while $(d)$ and
$(e)$ are originally black-and-white. Consequently, only the colored films
could be used for training, while the black-and-white ones could be used only
for testing with a colorizer trained on another film.

We have first experimented with training and testing on different parts of the
same (colored) film. For training/testing we have used the first 75%/last 25%
from each of the colored films. The proposed 3D cGAN model was used, with
model parameters set to $C=3$ (sliding window size), $\lambda=100$
(GAN-$L_{1}$ loss tradeoff), and compared against a 2D cGAN model that learned
to colorize each frame separately. We have also use data augmentation on our
training set, with random horizontal flips ($50\%$ chance to use a flipped
input during training) and gaussian additive noise
($\sim\mathcal{N}(0,1.2e-3)$). For estimate aggregation (eq. 2) we present
results with an uninformative prior (preliminary tests with priors learned
over data statistics did not give any definite improvement). We also compare
with a greyscale baseline, i.e. the case where the “colorized” video estimate
uses only luminance information. Numerical results can be examined in table 1.
Qualititative results can be examined in figures 2 and 3. While in general
both models fare satisfactorily, the proposed model can avoid erroneous
colorizations in several cases (cf. fig. 3). This point is validated by our
numerical results, where we calculate the metrics presented in section 3.
While w.r.t. to PSNR and RA the proposed model still is better, it could be
argued that the difference in the result is statistically insignificant. This
is not the case with the proposed CC metric however, where the performance of
the proposed model is markedly better. These results validate our expectation,
as the 3D structure of the proposed model can take into account the sequential
structure of the video, in contrast to its 2D counterpart.

We have also run tests for training and testing on different films. The case
that is perhaps closest to a practical application of the current model is
using trained models on one of the colored films to color black-and-white
footage, i.e. in our case films $(d)$ and $(e)$. Results for this case can be
examined at fig. 4 (training performed on film $(a)$). Video colorization
demos are available online 333http://www.cs.uoi.gr/~sfikas/video_colorization
.

| PSNR | RA | CC
---|---|---|---
(a) “Dial M for Murder”
Grayscale | 32.69 | 96.55 | 73.09
2D cGAN | 34.97 | 96.67 | 82.07
Proposed | 35.66 | 96.73 | 85.59
(b) “Et Dieu..créa la femme”
Grayscale | 30.23 | 94.07 | 47.82
2D cGAN | 32.08 | 95.17 | 56.67
Proposed | 32.32 | 95.31 | 58.80
(c) “Tzéni, Tzéni”
Grayscale | 29.83 | 92.85 | 39.17
2D cGAN | 31.44 | 93.87 | 50.81
Proposed | 31.77 | 94.14 | 55.16
Table 1:  Numerical results for colorization evaluation. Training and testing
is performed on different clips of the same film. PSNR is measured in dB; RA
and CC values are percentages. Higher values are better. The proposed model
performs best, in all cases.

##  5 Conclusion and Future work

We have presented a method for automatic video colorization, based on a novel
cGAN-based model with 3D convolutional and deconvolutional layers and an
estimate aggregation scheme. The usefulness of our model has been validated
with tests on colorizing old black-and-white film footage. Model performance
has also been evaluated with single-image based metrics as well as a newly
proposed metric that measures sequential color consistency. As future work, we
envisage exploring the uses of the color prior in our aggregation scheme.

|  |  |
---|---|---|---
|  |  |
Fig. 4:  Colorization results where color ground-truth is unavailable.
Depicted are samples from “A streetcar named desire” and “Twelve angry men” (2
leftmost, 2 rightmost columns respectively), colorized with the proposed model
trained on “Dial M for Murder”.

## References

  * [1] Simone Meyer, Victor Cornillère, Abdelaziz Djelouah, Christopher Schroers, and Markus Gross,  “Deep video color propagation,”  arXiv preprint arXiv:1808.03232, 2018.
  * [2] Mayu Otani and Hirohisa Hioki,  “Video colorization based on optical flow and edge-oriented color propagation,”  in Computational Imaging XII. International Society for Optics and Photonics, 2014, vol. 9020, p. 902002.
  * [3] VS Rao Veeravasarapu and Jayanthi Sivaswamy,  “Fast and fully automated video colorization,”  in Signal Processing and Communications (SPCOM), 2012 International Conference on. IEEE, 2012, pp. 1–5.
  * [4] Sifeng Xia, Jiaying Liu, Yuming Fang, Wenhan Yang, and Zongming Guo,  “Robust and automatic video colorization via multiframe reordering refinement,”  in IEEE International Conference on Image Processing. IEEE, 2016, pp. 4017–4021.
  * [5] Liron Yatziv and Guillermo Sapiro,  “Fast image and video colorization using chrominance blending,”  IEEE Transactions on Image Processing, vol. 15, no. 5, pp. 1120–1129, 2006.
  * [6] Anat Levin, Dani Lischinski, and Yair Weiss,  “Colorization using optimization,”  in ACM transactions on graphics (TOG). ACM, 2004, vol. 23, pp. 689–694.
  * [7] Bin Sheng, Hanqiu Sun, Marcus Magnor, and Ping Li,  “Video colorization using parallel optimization in feature space,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 24, no. 3, pp. 407–417, 2014.
  * [8] Nir Ben-Zrihem and Lihi Zelnik-Manor,  “Approximate nearest neighbor fields in video,”  in IEEE International Conference on Computer Vision and Pattern Recognition (CVPR), 2015, pp. 5233–5242.
  * [9] Tomihisa Welsh, Michael Ashikhmin, and Klaus Mueller,  “Transferring color to greyscale images,”  in ACM Transactions on Graphics (TOG). ACM, 2002, vol. 21, pp. 277–280.
  * [10] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio,  “Generative adversarial nets,”  in Advances in neural information processing systems (NIPS), 2014, pp. 2672–2680.
  * [11] Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen,  “Improved techniques for training GANs,”  in Advances in neural information processing systems (NIPS), 2016, pp. 2234–2242.
  * [12] Constantinos Daskalakis, Andrew Ilyas, Vasilis Syrgkanis, and Haoyang Zeng,  “Training GANs with optimism,”  CoRR, vol. abs/1711.00141, 2017.
  * [13] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros,  “Image-to-image translation with conditional adversarial networks,”  arXiv preprint arXiv:1611.07004, 2016.
  * [14] A.W. Juliani,  “Pix2Pix-Film,” http://github.com/awjuliani/Pix2Pix-Film, 2017,  [Online; accessed 2-January-2018].
  * [15] Richard Zhang, Phillip Isola, and Alexei A Efros,  “Colorful image colorization,”  in IEEE European Conference in Computer Vision (ECCV). Springer, 2016, pp. 649–666.
  * [16] Olaf Ronneberger, Philipp Fischer, and Thomas Brox,  “U-Net: Convolutional networks for biomedical image segmentation,”  in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2015, pp. 234–241.
  * [17] S. Chintala, E. Denton, M. Arjovsky, and M. Mathieu,  “How to train a GAN? Tips and tricks to make GANs work,” http://github.com/soumith/ganhacks, 2016,  [Online; accessed 25-January-2018].
  * [18] “Dial M for murder,” https://www.imdb.com/ title/tt0046912/, 1954.
  * [19] “Et Dieu..créa la femme,” https://www.imdb.com/ title/tt0049189/, 1956.
  * [20] “Tzéni, tzéni,” https://www.imdb.com/ title/tt0145006/, 1966.
  * [21] “A streetcar named desire,” https://www.imdb.com/ title/tt0044081/, 1951.
  * [22] “Twelve angry men,” https://www.imdb.com/ title/tt0050083/, 1957.


# Towards Microstructural State Variables in Materials Systems

Veera Sundararaghavan
Department of Aerospace Engineering
University of Michigan
Ann Arbor, MI
<EMAIL_ADDRESS>
& Megna Shah and Jeff Simmons
Materials and Manufacturing Directorate
Air Force Research Laboratory
Wirght Patterson Air Force Base, OH
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
Corresponding author: Prof. Sundararaghavan, Email<EMAIL_ADDRESS>Tel:
734-615-7242

###### Abstract

The vast combination of material properties seen in nature are achieved by the
complexity of the material microstructure. Advanced characterization and
physics based simulation techniques have led to generation of extremely large
microstructural datasets. There is a need for machine learning techniques that
can manage data complexity by capturing the maximal amount of information
about the microstructure using the least number of variables. This paper aims
to formulate dimensionality and state variable estimation techniques focused
on reducing microstructural image data. It is shown that local dimensionality
estimation based on nearest neighbors tend to give consistent dimension
estimates for natural images for all p-Minkowski distances. However, it is
found that dimensionality estimates have a systematic error for low-bit depth
microstructural images. The use of Manhattan distance to alleviate this issue
is demonstrated. It is also shown that stacked autoencoders can reconstruct
the generator space of high dimensional microstructural data and provide a
sparse set of state variables to fully describe the variability in material
microstructures.

_Keywords_ Intrinsic dimensionality $\cdot$ maximum likelihood estimation
$\cdot$ Minkowski distances $\cdot$ Microstructures $\cdot$ autoencoders

##  1 Introduction

It is widely accepted that controlling the microstructure of a material will
enable control of its properties. But it is less clear which, or even how
many, of the features of the microstructure represent its variability.
Recently, Chen, et al. [8] identified _intrinsic dimensions_ in complex and
chaotic dynamical systems, using only short videos of their behavior and
proposed that _state variables_ of complex systems may be identified in this
way. This suggests the tantalizing prospect of identification of a minimal set
of microstructural state variables that would govern the material’s behavior.
This minimum number of features would encode all of the dimensions in the
microstructure necessary to make design decisions, much like when the Wright
brothers ‘invented the airplane’ by discovering how to control _all_
dimensions of motion. Finding and controlling all dimensions of the
microstructure could enable a completely new way of exploiting design spaces.

Recent advances in characterization techniques and computing have led to
generation and analysis of large datasets, enabling improved understanding of
microstructure. Much work has been done to quantify various aspects of the
microstructure, such as particle size and shape distributions, orientation
distributions, n-point statistics among others [39, 28, 5, 35], enabling
significant advancements in processing-structure-properties understanding. But
this has relied on domain experts manually identifying which features should
be characterized. While advancing the understanding, this still leaves the
uncertainty as to the degree to which the important variation in the structure
was actually quantified. Although each microstructure can be represented as a
vector of size $n$, the actual dimensionality of an entire database of
microstructures is expected to be much lower. In formal terms, a data set
containing points of dimensionality $n$ is said to have intrinsic
dimensionality (ID) equal to $\mu$ < $n$ if every point lies entirely within
an $\mu$-dimensional manifold of $\Re^{n}$. The methods of dimensionality
estimation can be categorized as local and global approaches. Global methods
for ID estimation rely on the spread of the entire dataset, as exemplified by
projection methods such as the principal component analysis (PCA). Linear
methods such as PCA and multidimensional scaling were explored for
microstructural data in Refs [38, 34, 36]. However, it is known that such
methods tend to fail on non–linear manifolds [15]. Other global approaches to
dimension reduction such as ISOMAP and its variants treat non–linear manifolds
using geodesic distances [37] and have been used to reduce dimensionality of
microstructures [15, 24]. Local approaches use the local geometry of the high
dimensional space to estimate the intrinsic dimension and tend to be more
computationally efficient [10]. Levina and Bickel [23] developed such an
estimate by choosing an optimal dimension in which the local neighborhood of
points would be uniformly spaced. Pope et. al. [31] applied this methodology
to estimate the dimensionality of some well known benchmark datasets such as
the MNIST [22] and CIFAR [20] datasets and found that the information in those
had a surprisingly low number of dimensions, i.e. degrees of freedom.
Estimates ranged from 10 to 25 dimensions from the simplest to most complex
dataset.

Much of the work cited above relied on human judgement as to the
reasonableness of the dimensionality estimates and did not have any ground
truths by which to evaluate such reasonableness. Consequently, assessing the
validity of the methods becomes problematic. This paper addresses itself to
the problem of developing a self-consistent methodology for estimating the
dimensionality of random media such as microstructures through validation
against datasets with known dimensionality and by employing additional
distance measures, all of which should yield the same estimated
dimensionality. The classic work of Levina and Bickel [23] and subsequent
papers use the $L_{2}$ norm (Euclidean distance) to estimate the intrinsic
dimension. We find that this approach is inaccurate for low bit depth images,
due to the sparsity of the data. This is a systematic error that persists in
recent papers, for example in Ref. [8], where MLE estimates are higher than
the intrinsic dimension for binary images. This is especially concerning for
microstructure images that have a significantly reduced bit depth representing
a handful of material phases. Ability to obtain consistent dimensionality
estimates for generalized Minkowski distances is shown, so long as the
histogram of pixel values covers a wide range, but that the estimates become
inconsistent when this range is significantly reduced (known as _sparsity_ in
the imaging literature). In this case, its is found that the $L_{1}$ norm,
Minkowski distance for $p=1$ is the most accurate.

ID estimators provide only the true dimensionality leaving other questions
such as what state variables are actually encoded in these dimensions. More
recently, deep learning generative methods have created representations that
automatically capture the key variables accounting for most of the variation
in image based datasets [21, 19], and such models have been trained on
microstructural data [33, 11, 9, 4, 13, 18]. Machine learned representations
are expected to parsimoniously capture the maximal amount of information about
the microstructure, as was demonstrated in Ref [25] by combining neural
network representation of images with manifold learning. In Ref. [8], a
stacked autoencoder was employed to reduce physical dynamics data to the
intrinsic dimensional space. The minimum number of variables (matching the
intrinsic dimension) found from the autoencoder network are referred to as
‘state variables’ in [8], a terminology that is adopted in this work. To test
the technique, microstructure image datasets were upsampled from a synthetic
low dimensional space and passed to the stacked autoencoder. The results show
a successful reduction of the images back to space describing the state
variables providing a promising route to capture useful information in
microstructures.

##  2 Methodology

###  2.1 Microstructures as Random Variables

In this paper, microstructures are modelled as images whose contents are
outcomes of observations of random variables [27]. More formally, a random
variable $\mathcal{M}$ is defined to describe the _Microstructure_. In this
context, ‘Microstructure’ is that used by, say, a process engineer who wants a
certain microstructure because of its desirable properties. The outcomes
($\mathbf{m}\in\Re^{n}$) from sampling $\mathcal{M}$ represent the images that
would be observed, say, by a microscopist investigating the _microstructure_ ,
where the lower case ‘$\mathbf{m}$’ is used to distinguish an instance from
the class. Here, $n$ is the number of pixels in the image. This way, one can
make use of the considerable results from sampling theory, particularly point
processes[32], in the analyses.

###  2.2 Microstructures on a Manifold

Modeling _microstructure_ observations as images, an image is an outcome of
sampling $\mathcal{M}$ to give $\mathbf{m}$. If this image is, say $256\times
256$ in dimension, $\mathbf{m}\in\Re^{256\times 256}$. This is a huge space,
from which all images of this spatial resolution may be sampled. The vast
majority of these images simply represent random noise. By hypothesis, natural
(or microstructural) images occupy a very small subset of this space. That is,
the valid images that would plausibly represent a Microstructure occupy a
_manifold_ in $\Re^{256\times 256}$.

Speaking loosely, a manifold is a lower-dimensional space that is contained in
our $\Re^{256\times 256}$ space, but having fewer total number of dimensions.
A plane embedded in a 3-D space is an example of a _linear manifold_ , having
only 2 dimensions. More generally, the term ‘manifold’ means some non-linear
subspace that can be distorted within the embedding space. Figure 1(a) shows
an example of a manifold in $\Re^{3}$, which is known as the ‘swiss roll’
manifold. Essentially, this is a plane that contains all of the data, but has
been ‘rolled up’ into a spiral, so that it exists in $\Re^{3}$, but the
points, themselves only occupy $\Re^{2}$. In this work, the manifold is
referred to as a _latent space_ and the high–dimensional embedding space as
the _ambient space_. This is motivated by the fact that one would observe the
images in the ambient space ($\Re^{256\times 256}$ in the current example).

Figure 1: (i) A ‘swiss roll’ manifold containing image data represented as
points. Near–neighbor shells around a data point are illustrated which can be
used to estimate the intrinsic dimensionality. (ii) A mainifold representation
of binary images with n pixels, which exist on vertices of a cube of dimension
n. The space of 3 pixel images are shown with a 2D domain representing images
(marked a,b,c) with pixel values that sum to two.

Visualizing images on a latent manifold becomes problematic for greater than
three dimensions: very simple images must be used for illustration, with the
extension to higher dimensions being made in a more abstract sense. Using a
very simple image, consisting of only 3 pixels, one can illustrate a latent
manifold in an embedding space in Figure 1(b). The actual ambient space is the
closed set

$A=\\{(x,y,z)\in\Re^{3}|x\in[0,1],y\in[0,1],z\in[0,1]\\}$ (1)

The ‘corners’ of $A$ represent _binary images_ , i.e. 1-bit images, where the
pixels can only have values of 0 or 1.

Within this space, a (linear) manifold is embedded as:

$B=\\{(x,y,z)\in A|x+y+z=2\\}$ (2)

which represents binary images in which one of the pixels has a value of 0 and
two have a value of 1, as well as all _convex combinations_ [6] of these
images to form a constrained set of grayscale images.

By hypothesis, microstructure images occupy some latent manifold in an
enormous ambient space. Obtaining this manifold is the subject of _manifold
learning_[30, 26]. Our hypothesis is that the Microstructure may be controlled
by identifying state variables for its description and that these may be
enumerated if one knows the dimensionality of the latent manifold on which the
microstructure images lie. It is the subject of _disentanglement_ , an active
area of research in machine learning [17, 14], to make these dimensions
interpretable.

###  2.3 Nearest Neighbor Approach to Dimensionality Estimation

The nearest neighbor method [29] is a geometric estimator of the intrinsic
dimensionality of the manifold on which the data lies. The assumptions behind
this approach are (1) that the samples are independent and identically
distributed (iid) from some distribution, (2) that, in a space of proper
dimension, they will be uniformly distributed, (3) that the mapping between
the latent space and the ambient space is continuous, and (4) that the
distance between two points in the ambient space is the same as that in the
latent space.

The intuitive meaning of ‘random placement,’ where there is no bias towards
one area in space or another. This is a common one made with modeling, say,
trees in a forest. The unique point process that will assure such a random
placement is the Poisson process [2]. The intuitive meaning of ‘continuous’ is
that neighboring points in the latent space correspond to neighboring points
in the ambient space. Topology[7] provides a more precise statement of this,
but the essential intuitive interpretation is this.

There is one subtle complication that arises because data is generally not on
a linear manifold, but on one that is curved and twisted. The distance between
two points on a curved manifold would be measured as its geodesic distance,
whereas, in the ambient space, it would be measured as a Euclidean distance or
similar. Since differentiable manifolds are approximately Euclidean for small
distances, this amounts to a requirement that the distance between points be
made small.

With these assumptions, the dimensionality of a dataset may be made, knowing
only a distance between the points. Levina and Bickel used the Euclidean
distance, but we use the generalized $p$-norm approach of Minkowski, which
reduces to the Euclidean distance for $p=2$. All $p$-norms are required to
estimate the same dimensionality, which yields a ‘best practice’ for intrinsic
dimensionality estimation.

The nearest neighbor (NN) method aims to find the intrinsic dimensionality
$\mu\leq n$ using the number of nearest neighbors $k$ of each data point [23].
The data are modeled as being iid samples from a probability density in the
low dimensional latent space $\Re^{\mu}$. By hypothesis, there is a locally
homogeneous Poisson process, of dimension $\mu$, such that the density is
constant within a neighborhood of $\mathbf{m}$ [32], which will uniformly (at
least, locally) distribute the data points in this space.

Let $\mathbf{m}_{1},\mathbf{m}_{2},..,\mathbf{m}_{s}\in\Re^{n}$ be the
instances of $s$ microstructures. Under these assumptions, the average number
of data points ($\bar{k}$) that fall into a hypersphere in $\Re^{\mu}$ around
a point $\textbf{m}_{i}$ will be proportional to the volume of the
hypersphere:

$\bar{k}=f(\textbf{m})\mathcal{V}(\mu,p)$ (3)

where the proportionality constant, $f(\mathbf{m})$, defines the uniform
probability density defining number of points per unit volume in $\Re^{\mu}$
and $\mathcal{V}(\mu,p)$ refers to the volume of the hypersphere of
dimensionality $\mu$ that has an expected number $\bar{k}$ nearest neighbors
with distances represented using a $L_{p}$ norm.

The volume of the hypersphere is given by the particular choice of the
distance measure. Levina and Bickel[23] used the Euclidean distance measure
($p=2$) where the volume is given by the formula:

$\mathcal{V}(\mu,2)=V(\mu,2){[T_{k}(2)]}^{\mu}$ (4)

Where $V(\mu,2)$ is the volume of a hypersphere of unit radius in $\Re^{\mu}$,
$T_{k}(2)$ is the distance from a fixed point m to its $k^{th}$ nearest
neighbor in the ambient space, and the constant $2$ within brackets in Eq. 4
indicates the Minkowski 2-norm, which is the Euclidean distance measure, is
being used. By the locally isometric hypothesis, this is the same as the
distance would be measured in the latent space.

####  2.3.1 Generalized Distance Measures

For the Euclidean distance, the volume of a unit hypersphere is
$\frac{\pi^{\mu/2}}{\Gamma(\mu/2+1)}$[23]. We extend this analysis to apply to
the general Minkowski distances of order $p$, ($d_{p}$).

Between points $\textbf{m}_{q}$ and $\textbf{m}_{l}$, $d_{p}$ is defined as:

$d_{p}(\textbf{m}_{q},\textbf{m}_{l})\triangleq\left(\sum_{i=1}^{n}|\textbf{m}_{q,i}-\textbf{m}_{l,i}|^{p}\right)^{1/p}$
(5)

Particular cases of the Minkowski distance family are $d_{1}$, commonly known
as the Manhattan distance or the $L_{1}$ norm and $d_{2}$, commonly known as
the Euclidean distances or the $L_{2}$ norm. A geometric representation of a
2D _circle_ for $p=1,2,4$, and $\infty$ is shown in Fig. 2(a), where the
surface describes all points equidistant from the origin under the respective
$d_{p}$.

Figure 2: (a) Minkowski circles (b) Our estimation of the ratio of average
distances to $j^{th}$ and $k^{th}$ nearest neighbor shells for the swiss roll
dataset as a function of the Minkowski parameter.

The volume of a hypersphere of dimensionality $\mu$ when using a $d_{p}$
distance measure (see below) is:

$\mathcal{V}(\mu,p)=V(\mu,p){[T_{k}(p)]}^{\mu}$ (6)

where, $V(\mu,p)=\frac{2^{\mu}{[\Gamma(1/p+1)]}^{\mu}}{\Gamma(\mu/p+1)}$ is
the volume of a hypersphere of unit radius in $\Re^{\mu}$ and $T_{k}(p)$ is
the distance to the $k^{th}$ nearest neighbor, both being measured in terms of
the $d_{p}$ distance.

For a choice of the Minkowski parameter $p$, the relationship in Eq. 3 can be
used to estimate the dimension by regressing $\log{\bar{T}_{k}(p)}$ on $\log
k$ over a suitable range of $k$ (eg. from $k=k_{a}$ to $k=k_{b}$), where
$\bar{T}_{k}(p)$ denotes the mean $d_{p}$ distance of points to their $k^{th}$
nearest neighbor. The intrinsic dimension is obtained as the slope:

$\mu=\frac{\log k_{b}-\log k_{a}}{\log T_{k_{b}}(p)-\log
T_{k_{a}}(p)}=\log\left(\frac{k_{b}}{k_{a}}\right)\left(\log\frac{T_{k_{b}}(p)}{T_{k_{a}}(p)}\right)^{-1}$
(7)

Since $\mu$ is a unique intrinsic dimension, the above equation implies that
the ratio $\frac{T_{k_{b}}(p)}{T_{k_{a}}(p)}$ is independent of $p$.

This can be seen as follows, based on a Poisson point process. Using Eq. 3 and
Eq. 6, the expected number of points within a distance $r_{p}$ from a point m
can be written as:

$\bar{k}=c(p)r_{p}^{\mu}$ (8)

where $c(p)=f(\textbf{m})V(\mu,p)$. The hypersphere defined by the points
between $\mathbf{m}$ and the $k^{th}$ neighbor contains $k-1$ points in its
interior, the $k^{th}$ being on the boundary, itself. The Poisson distribution
($P$) for finding $k-1$ points within a distance of $r_{p}$ from point m is
given by:

$P(k-1)=\frac{(c(p)r_{p}^{\mu})^{k-1}}{\Gamma(k)}\exp(-c(p)r_{p}^{\mu})$ (9)

From which one can infer that the rate of the Poisson process is
$\lambda(p)=\frac{d}{dr}(c(p)r_{p}^{\mu})$.

The density function ($F_{k}$) of a distance $r_{p}$ from m to its $k^{th}$
neighbor can be written as (using Eq. 8 and 9, and performing change of
variables for the probability density), [29]:

$F_{k}(r_{p})=\left(\frac{(c(p)r_{p}^{\mu})^{k-1}}{\Gamma(k)}\exp(-c(p)r_{p}^{\mu})\right)c(p)\mu
r_{p}^{\mu-1}$ (10)

From this expression, the expectation of a distance $r_{p}$ from m to its
$k^{th}$ neighbor can be found as (see appendix 2):

$E_{k}(r_{p})=\int_{0}^{\infty}r_{p}F_{k}(r_{p})dr_{p}=c(p)^{-\frac{1}{\mu}}\frac{\Gamma(k+\frac{1}{\mu})}{\Gamma(k)}$
(11)

The leading term $c(p)^{-\frac{1}{\mu}}$, which is a function of Minkowski
parameter p, is independent of $k$. This implies that the ratio of average
distances for different values of $k$ (eg. in Eq. 7) will be independent of
the Minkowski parameter. This is, indeed, correct, as our estimates of the
ratio of average distances to $j^{th}$ and $k^{th}$ nearest neighbor shells
($\frac{T_{j}(p)}{T_{k}(p)}$) for different Minkowski parameters for the swiss
roll shows, Fig. 2(b).

###  2.4 MLE estimation using the p–norm

The intrinsic dimensionality estimator in Ref. [23] is a variant of the
nearest neighbor theory which seeks a maximum likelihood estimate (MLE)
instead of a mean estimate. The difference is subtle: The nearest neighbor
approach estimates the dimension as a statistic that can be computed from data
(as in Eq. 7), while the MLE approach seeks the optimum parameter in the
Poisson distribution (eq. 9), which in practise yields a more robust estimate
of dimensionality.

The log likelihood of the Poisson process can be written as:

$L(\mu,\theta,p)=\int_{0}^{R}\log(\lambda(p))~{}dN(r_{p})-\int_{0}^{R}\lambda(p)~{}dr_{p}$
(12)

where $N(r_{p})$ is the number of points within a distance $r_{p}$ from
$\mathbf{m}$ and $\theta=\log f(\mathbf{m})$. Maximizing the likelihood using
$\frac{\partial L}{\partial\theta}=0$ and $\frac{\partial L}{\partial\mu}=0$,
an optimal value of $\mu$ is obtained, also containing ratios of distances
[23]:

$\mu_{k}(\mathbf{m_{i}},p)=\frac{1}{k-1}\left[\sum_{j=1}^{k-1}\log\left(\frac{T_{k}(p)}{T_{j}(p)}\right)\right]^{-1}$
(13)

As described in Levina and Bickel, a denominator of $k-2$ instead of $k-1$
gives an unbiased estimate and is employed in this work. Fig. 3 shows the
variation of computed intrinsic dimension by this approach against the choice
of Minkowski parameter for a Helix and a broken swiss roll dataset. The
correct intrinsic dimension is found for all Minkowski parameters tested: 1
for the helix and 2 for the broken swiss roll.

Figure 3: MLE based intrinsic dimension for different Minkowski parameters:
(left) Helix, $\mu=1$ (right) Broken swiss roll, $\mu=2$.

###  2.5 Numerical considerations

Datasets of interest in this work are material microstructures that typically
consist of two phases (eg. a fiber composite, containing a fiber and a matrix)
or a finite number of phases (eg. steels containing austenite, martensite,
ferrite etc.). In two phase images, pixels are labeled such that the
precipitate is one and the matrix is zero (binary image). Note that the
independence of intrinsic dimensionality to Minkowski parameter is true only
if the distance $r$ is non-discrete (eg. cases in Fig. 3). Here, we show an
example where this breaks down using the case of binary images. Here, one has
a discrete manifold where images occupy vertices of a cube as shown in Fig.
1(ii). Consider the distance between any two images under the Minkowski
distance family for the case of binary image. Since the absolute difference
between any two pixels is either zero or one, the distance between any two
image instances $\textbf{m}_{q}$ and $\textbf{m}_{l}$ can be written as:

$d_{p}(\textbf{m}_{q},\textbf{m}_{l})=(||\textbf{m}_{q}-\textbf{m}_{l}||_{1})^{1/p}$
(14)

Here, $||..||_{1}$ refers to the $L_{1}$ norm (Manhattan distance) which
derives from using $\textbf{m}_{q,i}-\textbf{m}_{l,i}$ in eq. 14 is either
zero or one for binary images. For Minkowski parameter $p=1$, lets say that
fitting $\log{\bar{T}_{k}(1)}$ against $\log k$ would give the slope as the
intrinsic dimension $\mu^{*}$. When using $p=2$, substituting
$T_{k}(2)=(T_{k}(1))^{1/2}$ (eq. 14) would give an intrinsic dimension (slope)
of $2\mu^{*}$ instead. The same behavior is also obtained with the MLE
equation (Eq. 13). In general, an intrinsic dimension of $p\mu^{*}$ would be
obtained for the p-Minkowski measure. This inconsistency in our computed value
of intrinsic dimension for binary images is related to a discretization error.
One of the objectives of this paper is to identify a p-Minkowski measure that
mitigates this issue in binary microstructural images (or in general, images
with low bit depth) via numerical examples.

##  3 Results

The results employ the MLE estimate in Eq. 13. In this equation, every data
point $\textbf{m}_{i}$ gives a dimension estimate for every neighbor count
$k$. A dimension estimate matrix of size $k\times n$ is obtained where $n$ is
the number of images in the dataset. The intrinsic dimension is estimated as
the mean of the values in the dimension estimate matrix. An important aspect
in dimensionality estimation is removing duplicates from the dataset so that
the same datapoint is not double–counted as its own neighbor. This was done on
all datasets before computing the intrinsic dimension.

###  3.1 Binary Datasets

In rest of the results, the datasets are split into four categories:

  * •

Dataset 1. Rectangles and squares in a matrix. Ten different synthetic
datasets were tested containing rectangular and square shapes in a matrix
following different size and positional constraints which dictate the
intrinsic dimensionality (shown in Fig. 5). Images are of size $128\times
128$. In cases A, B, F, G, the shapes were randomly placed leading to two free
dimensions of x and y coordinates of the center of the shape. In addition, for
cases A and B varying sizes were used adding one more dimension in the case of
squares (the width) and two more in the case of rectangles (width and height).
In cases C, D, H and I, the shapes were placed linearly along a horizontal
axis at the center, eliminating one intrinsic dimension (the y-coordinate of
center) from cases A, B, F, G, respectively. The last two cases, E and J,
include centered shapes with a variation only in the size of the object. In
all cases, the shapes in the images are complete, ie, they are not cutoff by
the image boundary.

  * •

Dataset 2. Randomly centered circles of a constant radius in a matrix. The
coordinates of the center $(x,y)$ are independently sampled using uniform
random variables. Images are of size $256\times 256$. Four datasets 2A, 2B,
2C, 2D containing radii of $r=$ 24, 36, 48 and 58 pixels, respectively were
generated with circle centers ($(x,y)$) selected within a range such that the
circles do not intersect the boundaries.

$\\{(x,y)\in\Re^{2}|r<x<256-r,r<y<256-r\\}$ (15)

  * •

Dataset 3. A known 3D point cloud is used as the generator for microstructural
images, where each point in the 3D cloud $(x_{1},x_{2},x_{3})$ is mapped to
$(x,y,R)$ in a binary image of a circle in a matrix, where $R$ is the circle
radius, $(x,y)$ is the center of the circle. In this way, the radius is no
longer a free variable and is related to the center of the circle via the
topology of the point cloud (termed the ‘generator space’ to differentiate
from the latent space, which could be lower in dimension). The mapping for the
circle from the generator space to the ambient space would be
$f:\Re^{3}\rightarrow\Re^{128\times 128}$. Two cases were used
(i) Dataset 3A (swiss roll circles). The generator space is a 3D swiss roll
given as:

$x=t\cos t+c_{1},y=30\eta_{2}+c_{2},z=t\sin t+c_{3}$ (16)

where $t=\frac{3\pi}{2}(1+2\eta_{1})$, $\eta_{1}$ and $\eta_{2}$ are uniform
random variables in the range of 0 to 1 and $(c_{1},c_{2},c_{3})$ are
translation factors chosen as $(62,50,20)$ respectively. The latent space will
be two dimensional, governed by the choice of the two random variables
$\eta_{1}$ and $\eta_{2}$.

(ii) Dataset 3B (helix circles). Generator space is a 3D helix, given by the
equations:

$x=5(13+(2+\cos 8t)(\cos t)),y=5(13+(2+\cos 8t)(\sin t)),z=9(4+\sin 8t)$ (17)

where $t=2\pi\eta$, $\eta$ being a uniform random variable in the range of 0
to 1. The latent space will be one dimensional, with the location in the
manifold governed by the choice of $\eta$. Note that the coordinates of points
in the generator space of Dataset 3A and 3B are rounded before mapping to
images because of the integer (pixel) representation of the centers and radii.

  * •

Dataset 4: Contains results from a phase field simulation of grain growth
based on the Allen-Cahn equation following the numerical formulation of Fan
and Chen [12]. The data is in the form of binary images ($128\times 128$)
containing grain boundaries at different time steps of a single simulation.
Since all the model parameters are fixed at the start of the simulation and
images are only a function of time, the intrinsic dimensionality of all images
from a single simulation run is expected to be one. Dataset contains results
from three different simulation runs.

###  3.2 Intrinsic dimension estimates

Fig. 4(a) shows the variation of intrinsic dimension with the choice of
Minkowski parameter for binary image dataset–2D. The intrinsic dimension is
two, and corresponds to the $(x,y)$ coordinate of the center of the circle. In
the binary case, a linear increase in the intrinsic dimension with Minkowski
parameter is obtained as explained previously. While it was expected that the
$L_{1}$ norm gives the minimum intrinsic dimension estimate among these cases,
it is also seen that the $L_{1}$ distance measure matches the expected
intrinsic dimensionality. To further confirm this, the circle is replaced with
a Gaussian distribution centered at the circle origin and with a constant
standard deviation (of 20) for all datasets. Since each of the four cases in
dataset 2 contains a circle of different radii, each case spans different
number of grayscale levels. An example is shown in Fig. 4(e) with the blue
line spanning a part of the Gaussian curve, distance between the blue lines is
the circle diameter). Cases with radii of 24, 36, 48 and 58 pixels contain
201, 419, 706 and 1001 grayscale levels, respectively. As the number of
grayscale levels increase, the intrinsic dimension obtained from the use of
higher p-Minkowski distance measures converge toward the true estimate as
given by the $L_{1}$ norm.

Figure 4: a. Variation of intrinsic dimension with the choice of Minkowski
parameter for a circle with varying position, represented as binary and
grayscale images. The true dimensionality is 2. b. Shows a binary
microstructure c. Intensity across the dotted line in (b) is shown. (d) A
grayscale microstructure based on a Gaussian intensity profile. As the number
of grayscale levels increase, the answers for higher p-Minkowski distance
measures converge toward the true estimate.

To further test the use of $L_{1}$ norm for MLE estimation of binary images,
all ten cases in dataset 1 were tested (results shown in Fig. 5). The
intrinsic dimension expected for each case is indicated by numbers in green
near the cases in Fig. 5. A histogram of values in the dimension estimate
matrix is also plotted and the standard deviation of the histogram is reported
in addition to the mean. As seen from these results, the MLE approach with the
L1 norm gives a sound estimate of the intrinsic dimension in all cases.

Figure 5: Variation of MLE based intrinsic dimension for a variety of
synthetic datasets using $L_{1}$ norm. The predicted dimension for various
data points are shown as a histogram. The numbers in red indicate mean
predictions for each case. The numbers in green are the expected
dimensionality.

###  3.3 Retrieving state variables using an autoencoder

While the MLE algorithm recovers the intrinsic dimensionality, it is of
interest to identify the geometry of the latent space and to correlate the
dimensions to microstructural features. A variety of applications can benefit
from such analysis, including identification of novel processing paths and
inverse design of microstructures for a given property as shown in Ref. [33].
To generate a proof–of–concept, synthetic dataset C and D are used where the
generator space is known. Our objective was to check if the generator space
can be retrieved solely from the image data.

Figure 6:  Use of autoencoders for building a generator space is shown using a
synthetic dataset. (a) The swiss roll is the generator space, each point has a
$(x_{1},x_{2},x_{3})$ coordinate equal to $(R,x,y)$ in an image where $R$ is
the circle radius, $(x,y)$ is the center of the circle. (b) Shows images from
the database. (c) Stacked autoencoder reads in the images into a network with
a bottleneck equal to the generator space dimension (=3) (d) Decoded images
from the last layer of the network. (e) The 1000 images are reduced to three
variables in the bottleneck. These variables when plotted show the generator
space.

The state variables were identified using an autoencoder architecture. An
autoencoder (AE) [3] is a multi-layer neural network that learns the identity
function, such that the output $\hat{\mathbf{x}}$ approximates the input
$\mathbf{x}$. In the architecture, the hidden layers have fewer nodes than the
input dimension and act as a bottleneck. In the first few layers, the
autoencoder compresses the input to a compressed (latent space) representation
in a process called ‘encoding’. At its simplest, a single hidden layer
operates on the input $\mathbf{x}\in\Re^{n}$ and generates an encoding
$\mathbf{y}_{1}\in\Re^{j},j<n$ such that:

$\mathbf{y}_{1}=\sigma(\mathbf{W}_{1}\mathbf{x}+\mathbf{b}_{1})$ (18)

where $\mathbf{W}_{1}$ represents the $j\times n$ weight matrix,
$\mathbf{b}_{1}$ is the $j\times 1$ bias vector for the first layer. The
function $\sigma$ is typically a non-linear activation function and a logistic
sigmoid function is used in this work.

Later layers reconstruct the output from this latent space representation in a
process called ‘decoding’. An example is another layer that maps the latent
vector $\mathbf{y}_{1}$ in the previous step to output
$\hat{\mathbf{x}}\in\Re^{n}$ such that

$\hat{\mathbf{x}}=\sigma(\mathbf{W}_{2}\mathbf{y}_{1}+\mathbf{b}_{2})$ (19)

where $\mathbf{W}_{2}$ represents the $n\times j$ weight matrix,
$\mathbf{b}_{2}$ is the $n\times 1$ bias vector of the second layer. Multiple
layers can be used to develop a deep network. The parameters in $W$ and $b$
are found by minimizing the cost,
$\frac{1}{2}(||\mathbf{x}-\hat{\mathbf{x}}||_{2})^{2}$, by training via
backpropagation.

In this work, a stacked autoencoder configuration comprised of four layers in
total as shown in Fig. 6(c) is employed. The first autoencoder comprised of
two layers (encoder and decoder) was trained to reduce the dimensions to 100
first. This was followed by a second autoencoder with two layers (encoder and
decoder) that uses the 100 dimensional feature from the first autoencoder as
input and reduces it to the intrinsic dimension identified by the MLE
algorithm or the generator dimension. The two autoencoders were sequentially
trained first, followed by re–training a combined four–layer autoencoder.

Fig. 6 shows a schematic of the approach using dataset 3A (circles sampled
from a swiss roll, Fig. 6(b)). The stacked autoencoder reads in the images
into a network with a bottleneck equal to the generator space dimension (=3).
The decoded images from the last layer are shown in Fig. 6(d). The three
variables corresponding to each image are plotted in Fig. 6(e) which shows
that the generated topology is similar to the actual generator space in Fig.
6(a). The points shown are colored according to the circle radius in the
images. Note that the autoencoder, by default, restricts the range of values
to between 0 and 1. However, the topology of the space is generally well
reconstructed demonstrating a proof–of–concept that variables that define the
microstructural state can be identified using stacked autoencoders.

Figure 7:  Helix circles dataset 3B: (a) Generator space where images are
sampled (b) Generator space computed from the images using the stacked
autoencoder (c) Intrinsic dimensional (1D) space computed by the stacked
autoencoder.

Figure 7(b,c) shows both the computed generator space and the 1D latent space
for the helix circles dataset 3B (Figure 7(a) shows the actual generator space
from which the images were sampled with $(x_{1},x_{2},x_{3})$ coordinate equal
to $(r,x,y)$ in the image). As in the swiss roll case, the generator space for
the helix data set also looks similar to the generator space considering that
the points are mapped in the range $[0,1]$ by the autoencoder. The computed
latent spaces are colored according to a microstructural feature, the radius
of the circle. Clustering of this feature in the latent space demonstrates
promise towards the use of proximity analysis to find new microstructures with
interesting properties through interpolation. Such an application will form a
part of our future efforts.

The last example is from dataset 4 (phase field data). In the results from
this dataset in Fig. 8, one simulation trajectory was used in two stacked
autoencoders of bottleneck 3 and 2. Fig. 8(a) compares the reconstructions
from the autoencoder against the original images at five randomly chosen time
steps from this trajectory. Each consecutive time step in the phase field data
results in incremental changes in topology of grain boundaries, hence it is
expected that data points are arranged in the order of time steps in the
latent space. This is indeed seen in the topologies of the space constructed
by the autoencoder as shown in Fig. 8(b,c). Similar to the helix case seen
earlier, the intrinsic dimensionality is expected to be one which is confirmed
by the results of the MLE algorithm in Fig. 8(c).

Figure 8: (a) Decoder outputs for a single phase field trajectory (b)
Comparison of 3D versus 2D representation of the autoencoder (c) MLE estimate
of dimensionality

In the last example, the complete dataset containing three different phase
field simulations is employed with an autoencoder bottleneck of 3. The
resulting latent space is shown in Fig. 9. In all three simulations, the
initial image was the same represented by the central point in the latent
space. The trajectories from the three simulations emerge in different
directions from the initial point, resulting in different final
microstructures. This latent space is an example of a microstructural space
for a grain coarsening process. To estimate the true dimensionality of the
processing space, a large number of trajectories need to be superposed.
Although a simplistic set of three trajectories are shown, this example shows
how the framework can be used to visualize a multitude of complex processes
within a single graph. Past work in [36, 1] employed similar visualization of
linear PCA components to perform process design, the use of non–linear
manifolds as demonstrated here is expected to significantly improve
state–of–the–art and will form part of our future work.

Figure 9: (a) Decoder outputs for combined dataset containing all three phase
field trajectories (b) 3D representation from autoencoder showing all three
trajectories.

##  4 Conclusions

A methodology for reliably estimating the intrinsic dimensionality of random
media is developed. The method resolves the ambiguity in results for images
with low bit depth when using state–of–the–art techniques that employ the
Euclidean ($L_{2}$) norm. Particular novel contributions of this work are
listed below:

  * •

It is shown that the NN and MLE formulae for intrinsic dimensionality
estimation can work with all Minkowski distance measures. Further, the
examples show that all Minkowski measures give the same intrinsic
dimensionality for non-discrete datasets.

  * •

Dimensionality estimates are dependent on distance measures for image data
with discrete levels as was shown for binary images. Through examples, it is
shown that the use of $L_{1}$ distance in the MLE estimate produces a reliable
estimate in such cases.

  * •

Examples show that stacked autoencoders can reconstruct the low dimensional
generator spaces of microstructures and provide a sparse set of state
variables to fully describe material microstructures.

Use of these state variables to quantify microstructure–property–process
relationships will form a part of our future work.

## Data availability

Codes and datasets will be made freely available upon publication in a
peer–reviewed journal.

## Acknowledgements

This research was supported in part by the Air Force Research Laboratory
Materials and Manufacturing Directorate, through the Air Force Office of
Scientific Research Summer Faculty Fellowship Program, Contract Numbers
FA8750-15-3-6003 and FA9550-15-0001. This research was supported in part
through computational resources and services provided by Advanced Research
Computing at the University of Michigan, Ann Arbor.

## References

  * [1] Pinar Acar and Veera Sundararaghavan.  Linear solution scheme for microstructure design with process constraints.  AIAA Journal, 54(12):4022 – 4031, 2016.
  * [2] Adrian Baddeley, Imre Bárány, and Rolf Schneider.  Spatial point processes and their applications.  Stochastic Geometry: Lectures Given at the CIME Summer School Held in Martina Franca, Italy, September 13–18, 2004, pages 1–75, 2007.
  * [3] Pierre Baldi and Kurt Hornik.  Neural networks and principal component analysis: Learning from examples without local minima.  Neural networks, 2(1):53–58, 1989.
  * [4] Ramin Bostanabad.  Reconstruction of 3d microstructures from 2d images via transfer learning.  Computer-Aided Design, 128:102906, 2020.
  * [5] Ramin Bostanabad, Yichi Zhang, Xiaolin Li, Tucker Kearney, L Catherine Brinson, Daniel W Apley, Wing Kam Liu, and Wei Chen.  Computational microstructure characterization and reconstruction: Review of the state-of-the-art techniques.  Progress in Materials Science, 95:1–41, 2018.
  * [6] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe.  Convex optimization.  Cambridge university press, 2004.
  * [7] Ronald Brown.  Topology and groupoids, 2006.
  * [8] Boyuan Chen, Kuang Huang, Sunand Raghupathi, Ishaan Chandratreya, Qiang Du, and Hod Lipson.  Automated discovery of fundamental variables hidden in experimental data.  Nature Computational Science, 2(7):433–442, 2022.
  * [9] Kamal Choudhary, Brian DeCost, Chi Chen, Anubhav Jain, Francesca Tavazza, Ryan Cohn, Cheol Woo Park, Alok Choudhary, Ankit Agrawal, Simon JL Billinge, et al.  Recent advances and applications of deep learning methods in materials science.  npj Computational Materials, 8(1):1–26, 2022.
  * [10] Jose Costa and Alfred Hero.  Manifold learning with geodesic minimal spanning trees.  arXiv preprint cs/0307038, 2003.
  * [11] Dennis M Dimiduk, Elizabeth A Holm, and Stephen R Niezgoda.  Perspectives on the impact of machine learning, deep learning, and artificial intelligence on materials, processes, and structures engineering.  Integrating Materials and Manufacturing Innovation, 7(3):157–172, 2018.
  * [12] Danan Fan and L-Q Chen.  Computer simulation of grain growth using a continuum field model.  Acta Materialia, 45(2):611–622, 1997.
  * [13] Daria Fokina, Ekaterina Muravleva, George Ovchinnikov, and Ivan Oseledets.  Microstructure synthesis using style-based generative adversarial networks.  Physical Review E, 101(4):043308, 2020.
  * [14] Marco Fumero, Luca Cosmo, Simone Melzi, and Emanuele Rodolà.  Learning disentangled representations via product manifold projection.  In International Conference on Machine Learning, pages 3530–3540. PMLR, 2021.
  * [15] Baskar Ganapathysubramanian and Nicholas Zabaras.  A non-linear dimension reduction methodology for generating data-driven stochastic input models.  Journal of Computational Physics, 227(13):6612–6637, 2008.
  * [16] Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik.  Table of integrals, series, and products.  Academic press, 2014.
  * [17] Daniella Horan, Eitan Richardson, and Yair Weiss.  When is unsupervised disentanglement possible?  Advances in Neural Information Processing Systems, 34:5150–5161, 2021.
  * [18] Yubo Huang, Zhong Xiang, and Miao Qian.  Deep-learning-based porous media microstructure quantitative characterization and reconstruction method.  Physical Review E, 105(1):015308, 2022.
  * [19] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114, 2013.
  * [20] Alex Krizhevsky.  Learning multiple layers of features from tiny images.  University of Toronto, 2009.
  * [21] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [22] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [23] Elizaveta Levina and Peter Bickel.  Maximum likelihood estimation of intrinsic dimension.  Advances in neural information processing systems, 17, 2004.
  * [24] Zheng Li, Bin Wen, and Nicholas Zabaras.  Computing mechanical response variability of polycrystalline microstructures through dimensionality reduction techniques.  Computational Materials Science, 49(3):568–581, 2010.
  * [25] Nicholas Lubbers, Turab Lookman, and Kipton Barros.  Inferring low-dimensional microstructure representations using convolutional neural networks.  Physical Review E, 96(5):052111, 2017.
  * [26] Yunqian Ma and Yun Fu.  Manifold learning theory and applications, volume 434.  CRC press Boca Raton, 2012.
  * [27] Stephen Richard Niezgoda.  Stochastic representation of microstructure via higher-order statistics: theory and application.  Drexel University, 2010.
  * [28] Joachim Ohser and Frank Mücklich.  Statistical analysis of microstructures in materials science.  Wiley, New York, 2001.
  * [29] Karl W Pettis, Thomas A Bailey, Anil K Jain, and Richard C Dubes.  An intrinsic dimensionality estimator from near-neighbor information.  IEEE Transactions on pattern analysis and machine intelligence, 1(1):25–37, 1979.
  * [30] Robert Pless and Richard Souvenir.  A survey of manifold learning for images.  IPSJ Transactions on Computer Vision and Applications, 1:83–94, 2009.
  * [31] Phillip Pope, Chen Zhu, Ahmed Abdelkader, Micah Goldblum, and Tom Goldstein.  The intrinsic dimension of images and its impact on learning.  arXiv preprint arXiv:2104.08894, 2021.
  * [32] Donald L Snyder and Michael I Miller.  Random point processes in time and space.  Springer Science & Business Media, 2012.
  * [33] Srihari Sundar and Veera Sundararaghavan.  Database development and exploration of process–microstructure relationships using variational autoencoders.  Materials Today Communications, 25:101201, 2020.
  * [34] V Sundararaghavan and N Zabaras.  A dynamic material library for the representation of single-phase polyhedral microstructures.  Acta Materialia, 52(14):4111 – 4119, 2004.
  * [35] V Sundararaghavan and N Zabaras.  Classification and reconstruction of three-dimensional microstructures using support vector machines.  Computational Materials Science, 32(2):223 – 239, 2005.
  * [36] Veera Sundararaghavan and Nicholas Zabaras.  Linear analysis of texture property relationships using process-based representations of rodrigues space.  Acta Materialia, 55(5):1573 – 1587, 2007.
  * [37] Joshua B Tenenbaum, Vin de Silva, and John C Langford.  A global geometric framework for nonlinear dimensionality reduction.  science, 290(5500):2319–2323, 2000.
  * [38] Sanket Thakre, Vishnu Harshith, and Anand K Kanjarla.  Intrinsic dimensionality of microstructure data.  Integrating Materials and Manufacturing Innovation, 10(1):44–57, 2021.
  * [39] Salvatore Torquato.  Random heterogeneous materials: Microstructure and macroscopic properties.  Springer-Verlag, New York, 2002.

## Appendix 1

In the datasets used in this paper, we avoided the intersection of the
particle shape with the image boundary to get an unbiased estimation of the
dimensionality. In general, image boundaries can play a role in biasing the
intrinsic dimension. Consider the case of a single circular shape placed in a
matrix. If the shape intersects the boundary, only a part of the shape is seen
and the intrinsic dimensionality of three as judged from the dataset assumes
that the shape always remains a circle. To test this, we have plotted results
from two datasets, one with and one without boundary intersections. We
consider sufficiently high sample size ($3000$ images) and image size
($128\times 128$) for each case.

Figure S1: Histograms of estimated dimension per data point (a) Two cases are
considered, a circle of varying radii is freely placed in one case and avoids
boundary in another case. (b) Here, the circle is linearly placed along the
centerline. A case with periodicity is included. Inset images show
superposition of a few different images in each dataset.

The histogram of estimated dimension per data point is plotted in Fig. S1a
showing that the algorithm is able to predict the correct mean dimension for
both cases. The key difference being that the histogram is broader and a
higher standard deviation is obtained when the circles intersect the boundary.
Another case is shown in Fig. S1b where the circle is linearly placed along
the centerline. Here, two distinct peaks are seen in the histogram where the
circles are freely placed, with the first peak at a lower intrinsic dimension.
the mean dimension is again correctly estimated for both cases with a higher
standard deviation for the case where circles intersect the image boundaries.
A case where the circles are periodic is also shown here, where the circles
wrap around on the opposite side when they intersect the boundary. While this
case shows a single sharp peak as in the case where boundaries are avoided,
the standard deviation is higher than that case.

## Appendix 2

Equation (10) simplifies to:

$\begin{split}F_{k}(r_{p})&=\left(\frac{(c(p)r_{p}^{\mu})^{k-1}}{\Gamma(k)}e^{-c(p)r_{p}^{\mu}}\right)c(p)\mu
r_{p}^{\mu-1}\\\
&=\left[c(p)^{k}\mu\left(\frac{1}{\Gamma(k)}\right)\right]r_{p}^{\mu
k-1}e^{-c(p)r_{p}^{\mu}}\\\ \end{split}$

Inserting this into Equation (11):

$\begin{split}E_{k}(r_{p})&=\int_{0}^{\infty}r_{p}F_{k}(r_{p})dr_{p}\\\
&=\left[c(p)^{k}\mu\left(\frac{1}{\Gamma(k)}\right)\right]\int_{0}^{\infty}r_{p}^{\mu
k}e^{-c(p)r_{p}^{\mu}}dr_{p}\\\ \end{split}$

Changing the integration variable with $\xi\triangleq c(p)r_{p}^{\mu}$
($r_{p}=\xi^{1/\mu}c(p)^{-1/\mu}$):

$\begin{split}E_{k}(r_{p})&=\left[c(p)^{k}\mu\left(\frac{1}{\Gamma(k)}\right)\right]\int_{0}^{\infty}(\xi^{1/\mu}c(p)^{-1/\mu})^{\mu
k}e^{-\xi}d(\xi^{1/\mu}c(p)^{-1/\mu})\\\
&=c(p)^{-1/\mu}\left(\frac{1}{\Gamma(k)}\right)\int_{0}^{\infty}\xi^{k+1/\mu-1}e^{-\xi}d\xi\\\
\end{split}$

The integral is, by definition, $\Gamma(k+1/\mu)$ [16]. Substituting:

$E_{k}(r_{p})=c(p)^{-1/\mu}\frac{\Gamma(k+1/\mu)}{\Gamma(k)}$