# Transmission of natural scene images through a multimode fibre

Piergiorgio Caramazza1, Oisín Moran2, Roderick Murray-Smith2, Daniele Faccio1
1School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK
2School of Computing Science, University of Glasgow, Glasgow, G12 8QQ, UK
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

Abstract
The optical transport of images through a multimode fibre remains an
outstanding challenge with applications ranging from optical communications to
neuro-imaging. State of the art approaches either involve measurement and
control of the full complex field transmitted through the fibre or, more
recently, training of artificial neural networks that however, are typically
limited to image classes belong to the same class as the training data set.
Here we implement a method that statistically reconstructs the inverse
transformation matrix for the fibre. We demonstrate imaging at high frame
rates, high resolutions and in full colour of natural scenes, thus
demonstrating general-purpose imaging capability. Real-time imaging over long
fibre lengths opens alternative routes to exploitation for example for secure
communication systems, novel remote imaging devices, quantum state control
processing and endoscopy.

Introduction
Optical fibres form the backbone of the internet and are a key technology in
modern society Hecht2 . The vast majority of these fibres are ‘single mode’
i.e. they can transmit only one single, roughly Gaussian-shaped beam profile,
corresponding to the so-called fundamental mode of the fibre Hecht . It is
therefore impossible to directly transmit images through an optical fibre: any
attempt to do so simply results in transmission of the one single allowed mode
and therefore the detection at the output of this (Gaussian-shaped) mode with
all other information of the image completely lost at the fibre input. One
possibility to circumvent this limitation is to resort to an array or bundle
of single mode optical fibres, each one transmitting the information of a
single pixel in the output image. However, this quickly leads to fibre-bundled
cables that are relatively thick and not optimal for applications such as
endoscopic or neurological imaging, where the fibre bundle is inserted inside
a body fibre_review .
Another option is to resort to multimode fibres, i.e. fibres that due to a
larger core diameter can carry many optical modes that will have more complex
shapes than the fundamental mode and may encode image information cizmar3 .
For example, a typical 100 $\mu$m core diameter fibre might carry around
10,000 modes and could in principle transmit an image with roughly the same
number of pixels. However, in these fibres each of these individual modes
propagates at a slightly different velocity, thus leading to an amplitude and
phase mixing of the image as this propagates along the fibre fibre_book . The
image at the fibre output therefore appears as a random array of bright and
dark spots, referred to as a speckle pattern. This effect cannot be avoided
and completely scrambles and destroys the input image. Full a priori knowledge
of the input image and fibre details could allow to numerically model the
optical propagation theoryMMF , reconstruct the transmission matrix and then
unscramble the output data, but in practice this can be extremely hard.
Methods have been developed that allow to shape the input beam profile so as
to focus the output field into a single spot that can then be scanned MMF_holo
; MMF_holo2 ; MMF_holo3 ; MMF_holo4 ; ploschner2014gpu  with an emphasis on
endoscopy choi2012scanner ; cizmar1 ; cizmar ; cizmar_brain . Notwithstanding
this notable progress, the development of a viable method that allows to
unscramble the speckle patterns and thus retrieve high resolution, general
image information in real time is an open challenge.

Figure 1:  Experimental layout. (a) An SLM combined with a polarising beam
splitter (PBS) and a half-wave plate ($\lambda/2$) is used to imprint
intensity images onto a laser beam that is then coupled into a multimode
fibre. The fibre output is collected with a lens and recorded on a CMOS
camera. (b) A schematic overview of the computational processing steps of the
inversion process: output speckle data, $x$ from a series of images (in our
experiments, 50,000 images from the ImageNet database imagenet_cvpr09 ) are
fully connected to a complex matrix, $W$ which provides an output image
$I=|Wx|^{2}$. This image is compared to the actual original image (ground
truth) through a cost function: the total cost $\zeta$ is then back-propagated
to W and the process is repeated for a fixed number of loops (epochs),
ensuring minimisation of $\zeta$.

One promising route in this direction is based on the complete
characterisation of the optical fibre in the form of a measurement of its
transmission matrix choi2012scanner ; cizmar3 ; cizmar . This matrix connects
certain orthogonal modes at the fibre input to the fibre output and can
therefore, once known, be used to invert the speckle pattern back into the
original image. This approach requires measurements of the full complex
(amplitude and phase) profile of a large subset of modes and has been shown to
work over fibre lengths of 0.3-1 m.
Other approaches pioneered by Takagi et al. Hori , have recently been proposed
that used artificial neural networks (ANNs) using deep learning encoders to
infer images from the speckles patterns without any need for an a priori
mathematical model of the fibre psaltis1 ; psaltis2 ; Su_fibre . These have
used multi-layer convolutional ANN’s and have shown that it is possible to
reconstruct hand-written digits from the MNIST database lecun98  that consists
of patterns with $28\times 28$ pixel resolution: training of the network and
testing of its reconstruction abilities are both performed on digits from the
same database. ANNs have also been shown to allow for example to focus a beam
through a thin scattering medium (as well as through a multimode fiber) in
alex , where both single layer and multi-layer real-valued neural networks
have been implemented. When used for imaging, these approaches are mostly
expected to work for classes of objects that belong to the same class used for
the training as explicitly pointed out by Psaltis et al. psaltis1 . First
steps towards generic imaging have been made in Ref. psaltis2 : these
approaches present a very promising route forward if the suitability for
general purpose imaging applications can be addressed.

Figure 2: Reconstruction of Muybridge videos for 1 m long fibre. Individual
frames of video data taken from the Muybridge collection: a) running horse, b)
jumping cat, c) flying parrot, d) punching boxer. All videos are grayscale,
scaled to $92\times 92$ pixels and transmitted at 4 frames/second. The first
column shows the original input video. The second column shows the speckle
patterns ($x$) measured at the fibre output. SSIM and PCC indicate the
structural similarity index and Pearson correlation coefficient that quantify
the quality of the reconstruction (see Supplementary Note 2).

We have developed an approach that allows us to transmit and reconstruct
images of natural scenes at high resolution and frame rates. Our method
resorts to building an approximate model of the inverse of a complex-valued,
intensity transmission matrix of the optical fibre. This approach does not
require the use of deep (multi-layer) ANNs and enables the full reconstruction
of detailed images. Full colour images and videos can be recorded at 20 fps
and could be scaled up to kfps.

Results
Experiments. The experimental layout is shown in Fig. 1(a). We use an SLM
(maximum frame rate of 20 Hz) to impart grayscale (100 grayscale levels)
intensity images on to a continuous-wave laser beam (532 nm wavelength). This
image is then coupled into a multimode fibre (step index core, core diameter
105 $\mu$m, fibre lengths of 1 m and 10 m, $\sim 9000$ propagating optical
modes, image spot size at fibre input is $\sim 2$ $\mu$m) and then coupled out
using identical objectives for the fibre input and output (focal length,
$f=34$ mm, NA$=0.26$). The speckle pattern at the fibre output (near field) is
imaged on to a CMOS camera at $350\times 350$ pixel resolution.
The goal is to transmit ‘natural scenes’, i.e. photographs of everyday-life
scenes. The importance of this choice lies in the significant additional
complexity of natural scenes when compared to e.g. MNIST-database digits or
other simple geometric features. As sample images we use a selection of 50,000
images from the Imagenet database imagenet_cvpr09 , sized at $92\times 92$
pixels so as to have less pixels than optical fibre modes. These images have
been selected randomly from the Imagenet database while looking for almost
square images so as to facilitate projection into the fibre. The output
speckle patterns with amplitude distribution, $x$ (i.e. $x$ is the square root
of the measure speckle intensity patterns), together with the knowledge of the
image (intensity distribution, $I$) that generated each speckle pattern, is
used in the algorithm described below to approximate the inverse of a complex
transmission matrix, $W$. This matrix is then used to retrieve images that
were not part of the sample data set from intensity measurements of their
output speckle patterns, $I=|Wx|^{2}$. Examples are images and videos from
from the Muybridge collection such as a running horse, a jumping cat and a
flying parrot. We also tested the imaging on videos of a rotating Earth and
Jupiter. Both of these are in full colour, obtained by projecting and then
recombining the R, G and B channels independently.
Image reconstruction. There are two possible approaches to reconstructing an
image from a speckle pattern. The first is to attempt to build a forward model
that describes how the images or optical modes propagate down the fibre and
then invert this or, the approach followed here, one can try to directly
construct an approximation of the inverse model. This is achieved
statistically (i.e. by employing data from many images) through a single,
fully-connected complex-valued transformation matrix.

Figure 3: Full colour results for 1 m long fibre. Individual frames of video
data in full colour of a rotating Jupiter (credit: Damian Peach) (a) and Earth
(b).

Complex inversion. A schematic overview of the approach used to reconstruct
$W$ is shown in Fig. 1b). The full transmission matrix of an optical fibre is
complex-valued. This motivates the assumption that $W$ is also complex valued,
connecting the input and output of the fibre $I=|Wx|^{2}$ (Hir03a, ;
TygBruChi16, ; Gub16, ; Tra18, ). This also motivates idea that a deep-
learning ANN approach is not required here. We only measure the intensity of
the speckle pattern from which we take the amplitude $x$ (square root of the
intensity) and therefore represent $x$ with amplitude only, and zero phase.
The values of $x$ are passed to a fully connected (‘dense’) complex matrix –
equivalent to multiplication by the complex-valued matrix $W$ [arrows in Fig.
1b) show some of these connections as an example]. An image is then obtained
as $I=|W{x}|^{2}$. We calculate the derivatives ${d\zeta}/{dw_{ij}}$ of the
cost function $\zeta$, with respect to the $i,j$th element of $W$. We then
apply a stochastic gradient descent approach to make small changes to $W$ that
reduce the cost-function and the process is repeated for a fixed number of
loops (epochs), ensuring convergence of $\zeta$ to a minimum value. The
complex-weighted inversion was implemented as a novel layer with Keras
(chollet2015keras, ) and TensorFlow (tensorflow2015-whitepaper, ), (the code
is provided in the Supplementary Note 4). The final $W$, constructed from a
database of 50,000 images can then be used to obtain an estimate of the ground
truth image for all future transmitted data corresponding to images not used
as part of the training and indeed, even transmitted at a completely different
time (e.g. several days) after the training is completed.

Figure 4: Image reconstruction of an ostrich taken from the Muybridge
collection after transmission through a 1 m and 10 m long fibre. The inversion
matrix was constructed using only random grayscale patterns.

Experimental Results. In Fig. 2 we show a first set of results obtained from
data transmitted through a 1 m long fibre (wrapped in a loose coil on the
table) by applying the estimated $W$ to a series of videos that do not form
part of the ImageNet database and are significantly diverse, in order to
demonstrate the robustness of our approach. These videos are taken from the
Muybridge recordings from the 1870’s that marked the historically important
breakthrough of the first ever high-speed photography images. The running
horse in Fig. 2a) is probably the most iconic of the Muybridge videos, but the
others, a jumping cat in Fig. 2b), a flying parrot Fig. 2c) and punching boxer
Fig. 2d) provide a broad scene variability and all show good image
reconstruction over the full gray-scale spectrum as opposed to the binary
black/white MNIST images often used in previous work. We also note that the
MNIST database, given the limited set of symbols, will tend to create what is
essentially a simple classification system, with a decoder to generate the
associated image. True imaging capability should therefore demonstrate
functionality beyond this dataset, as shown in Fig. 2 (see Supplementary
Figure 3 for futher examples). Under each reconstructed image in Fig. 2 we
also give a quantitive measure of the reconstruction quality based on the
‘structural similarity index’ (SSIM) and the ‘Pearson correlation coefficient’
(see Supplementary Note 2). These coefficients measure the similarity between
the ground truth and retrieved images with a maximum value of 1 (indicating
image identity).
We also observed no degradation in the video quality even when the data was
transmitted, recorded and reconstructed more than 48 hours after transmission
of the original first set of ‘training’ data had been completed (see
Supplementary Note 3), thus indicating robustness to subsequent environmental
changes such as temperature fluctuations (of the order of a few degrees) and
vibrations (the setup is not placed on a vibration-isolated table).
In Fig. 3 we show examples of full colour video transmission of a rotating
Jupiter and a rotating Earth. Each individual R, G and B channel was
transmitted and reconstructed separately and then recombined. We note that the
same matrix $W$ obtained for grayscale images is used for all three R, G, B
channels when imaging in full colour mode. Subtle features such as the Red
Spot on Jupiter or slightly lighter areas in the Northern region of Africa
(roughly corresponding to the Nile delta region in Egypt) can be observed in
the reconstructed images (other examples are shown in the Supplementary Fig.
2).
Tests were also performed to investigate the role of the class of images used
for retrieving $W$. The images in Fig. 4 were down-sampled to $28\times 28$
pixels in order to simplify the problem and demonstrate that, if desired, one
may also reconstruct the inversion matrix $W$ with a completely “agnostic”
approach, i.e. with no prior assumption on the images. This is obtained by
using 50,000 completely random grayscale images. As can be seen, this
completely agnostic approach is still able to correctly reconstruct the images
although with a clear loss of quality. We noticed a good insensitivity to
fibre length (as already pointed out by Psaltis et al. psaltis1 ) and an
improvement of image quality with increasing number of random images used for
the $W$ matrix retrieval, although GPU RAM limitations did not allow us to
investigate this further.
We also noted that changing the size of the focused image at the fibre input
significantly impacts the final reconstruction. By placing a telescope after
the SLM so as to rescale the image at the focusing objective input, we noticed
a significant increase of the final image quality with increasing size at the
focusing objective input pupil (corresponding to an increasing effective NA,
i.e. to an increasing angular spread at the fiber input). This was also
accompanied by a clear decrease in the average speckle spot size at the fibre
output, indeed indicating the excitation of higher spatial frequency modes
(see Supplementary Note 1 and Supplementary Figure 1).
A matter of concern in many studies is the robustness to changes in the fibre
configuration. A change in the fibre geometry (e.g. by bending the fibre) will
lead to a different propagation of the individual modes, which ultimately
leads to a different output speckle pattern. Without precise knowledge of how
the fibre has been changed, it is not possible therefore to reconstruct the
image using the inversion matrix $W$ from a different configuration cizmar . A
recent solution has been proposed by using specially designed fibres that have
a parabolic refractive index profile in the core cizmar2018 . For a long-range
transmission system, such a solution might be appropriate under the assumption
that a long-haul fibre would remain in a relatively fixed position over time.
In future realisations, we expect that a combination of fibre design, position
classification and/or extensive training over fibre configurations, will allow
to efficiently remove this last obstacle that for the time being, is beyond
the scope of this work that is aimed at demonstrating that high pixel-density,
colour images can be transmitted efficiently through a static fibre and at
video frame rates.

Discussion
Imaging through a single multimode fibre is an outstanding challenge and has
so far been limited to relatively low frame rates, small image sizes or short
fibre lengths. Imaging of natural scenes increases this challenge further as
it ideally entails video frame rates, colour detail and sufficient image
resolution to allow identification of the scene details. The technique
developed here is based on a physically informed model of the imaging system
that retrieves an approximation to the full transmission matrix. In this
sense, our approach sits somewhere between the techniques devised for
reconstructing the actual transmission matrix and deep learning approaches
that do not make any explicit assumptions on the system. This allows for
efficient video and data transmission through fibres and promises applications
beyond endoscopic imaging, such as direct video or multimode data transmission
over long fibres for communication systems and fibre sensing by, for example,
exploiting the image sensitivity to changes along the fibre length. Moreover,
one could also exploit the intrinsic random nature of the multimode coupling
and output speckle patterns to securely encode and authenticate data as
proposed recently Goorden:14 ; random  where the transmission fibre itself
would play the role of the encoding medium, with possible extensions also to
the control of quantum states for quantum sensing and simulation quantum .

Data availability
An example of the code is provided along with data for 1m at
http://dx.doi.org/10.5525/gla.researchdata.751.

Code availability
The code has been presented and explained in the Supplementary Note 4:
Software. Furthermore, a code example is available at the DOI link.

Acknowledgements
DF and RMS acknowledge financial support from EPSRC (UK, grant no.
EP/M01326X/1).

## References

  * (1) Hecht, J.  _City of Light: the story of fibre optics._ (Oxford University Press, 2004).
  * (2) Hecht, J.  _Understanding Fiber Optics_ (Pearson, 2005).
  * (3) Keiser, G., Xiong, F., Cui, Y. & Shum, P. P.  Review of diverse optical fibers used in biomedical research and clinical practice.  _J. biomed. opt._ 19, 080902 (2014).
  * (4) Čižmár, T. & Dholakia, K.  Exploiting multimode waveguides for pure fibre-based imaging.  _Nat. Commun._ 3, 1027 (2012).
  * (5) Mitschke, F.  _Fiber Optics._ (Springer, 2010).
  * (6) Liu, C., Deng, L., Liu, D. & Su, L.  Modeling of a single multimode fiber imaging system.  _arXiv:1607.07905_ (2018).
  * (7) Di Leonardo, R. & Bianchi, S.  Hologram transmission through multi-mode optical fibers.  _Opt. Express_ 19, 247–254 (2011).
  * (8) Papadopoulos, I. N., Farahi, S., Moser, C. & Psaltis, D.  Focusing and scanning light through a multimode optical fiber using digital phase conjugation.  _Opt. Express_ 20, 10583–10590 (2012).
  * (9) Papadopoulos, I. N., Farahi, S., Moser, C. & Psaltis, D.  High-resolution, lensless endoscope based on digital scanning through a multimode optical fiber.  _Opt. Express_ 40, 260–270 (2013).
  * (10) Morales-Delgado, E. E., D., Psaltis & Moser, C.  Two-photon imaging through a multimode fiber.  _Opt. Express_ 23, 32158–32170 (2015).
  * (11) Plöschner, M., Straka, B., Dholakia, K. & Čižmár, T.  Gpu accelerated toolbox for real-time beam-shaping in multimode fibres.  _Optics express_ 22, 2933–2947 (2014).
  * (12) Choi, Y. _et al._ Scanner-free and wide-field endoscopic imaging by using a single multimode optical fiber.  _Physical review letters_ 109, 203901 (2012).
  * (13) Čižmár, T. & Dholakia, K.  Shaping the light transmission through a multimode optical fibre: complex transformation analysis and applications in biophotonics.  _Optics Express_ 19, 18871–18884 (2011).
  * (14) Plöschner, M., Tyc, T. & Čižmár, T.  Seeing through chaos in multimode fibres.  _Nat. Photon._ 9, 529 (2015).
  * (15) Turtaev, S. _et al._ High-fidelity multimode fibre-based endoscopy for deep-brain in vivo imaging.  _Light:Sci. Appl._ 7, 92 (2018).
  * (16) Deng, J. _et al._ ImageNet: A Large-Scale Hierarchical Image Database.  In _CVPR09_ (2009).
  * (17) Takagi, R., Horisaki, R. & Tanida, J.  Object recognition through a multimode fiber.  _Opt. Rev._ 24, 117–120 (2017).
  * (18) Borhani, N., Kakkava, E., Moser, C. & Psaltis, D.  Learning to see through multimode fibers.  _Optica_ 5, 960 (2018).
  * (19) Rahmani, B., Loterie, D., Konstantinou, G., Psaltis, D. & Moser, C.  Multimode optical fiber transmission with a deep learning network.  _Light:Sci. Appl._ 7, 69 (2018).
  * (20) Fan, P., Zhao, T. & Su, L.  Deep learning the high variability and randomness inside multimode fibres.  _arXiv:1807.09351_ (2018).
  * (21) Lecun, Y., Bottou, L., Bengio, Y. & Haffner, P.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ 86, 2278–2324 (1998).
  * (22) Turpin, A., Vishnaikou, I. & Seelig, J.  Light scattering control in transmission and reflection with neural networks.  _Optics express_ 26, 30911 (2018).
  * (23) Hirose, A.  Complex-valued neural networks: An introduction.  In _Complex-Valued Neural Networks: Theories and Applications_ , 1–6 (World Scientific, 2003).
  * (24) Tygert, M. _et al._ A mathematical motivation for complex-valued convolutional networks.  _Neural computation_ 28, 815–825 (2016).
  * (25) Guberman, N.  On complex valued convolutional neural networks.  _arXiv preprint arXiv:1602.09046_ (2016).
  * (26) Trabelsi, C. _et al._ Deep complex networks.  In _International Conference on Learning Representations_ (2018).  URL https://openreview.net/forum?id=H1T2hmZAb.
  * (27) Chollet, F. _et al._ Keras.  https://keras.io (2015).
  * (28) Abadi, M. _et al._ TensorFlow: Large-scale machine learning on heterogeneous systems (2015).  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * (29) Boonzajer Flaes, D. E. _et al._ Robustness of light-transport processes to bending deformations in graded-index multimode waveguides.  _Phys. Rev. Lett._ 120, 233901 (2018).  URL https://link.aps.org/doi/10.1103/PhysRevLett.120.233901.
  * (30) Goorden, S. A., Horstmann, M., Mosk, A. P., Škorić, B. & Pinkse, P. W. H.  Quantum-secure authentication of a physical unclonable key.  _Optica_ 1, 421–424 (2014).  URL http://www.osapublishing.org/optica/abstract.cfm?URI=optica-1-6-421.
  * (31) Carnicer, A. & Javidi, B.  Optical security and authentication using nanoscale and thin-film structures.  _Adv. Opt. Photon._ 9, 218–256 (2017).  URL http://aop.osa.org/abstract.cfm?URI=aop-9-2-218.
  * (32) Defienne, H., Barbieri, M., Walmsley, I. A., Smith, B. J. & Gigan, S.  wo-photon quantum walk in a multimode fiber.  _Science Advances_ 2, 1e1501054 (2016).

Transmission of natural scene images through a multimode fibre: Supplementary
Material

Piergiorgio Caramazza1, Oisín Moran2, Roderick Murray-Smith2,*, Daniele
Faccio1,*
1School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK
2School of Computing Science, University of Glasgow, Glasgow, G12 8QQ, UK
*Email<EMAIL_ADDRESS><EMAIL_ADDRESS>

Supplementary information to “Transmission of natural scene images through a
multimode fibre”
providing additional details regarding the code and measurements.

SUPPLEMENTARY NOTE 1: EFFECT OF IMAGE SIZE AT FIBRE INPUT

Supplementary Figure 1:  Comparison between same images focused with different
angular spread reported by an estimation of the effective NA (this was
realized by putting a telescope before the objective with different
magnification). Each image (a)-(d) is reported along with its relative speckle
pattern at the fiber output and ANN reconstruction. An increase in the
effective NA leads to averaged smaller speckles (indicating that higher
spatial frequency modes have been excited).

Whilst performing experiments we investigated the effect of different focusing
(i.e. image size) configurations at the fibre input. We performed three
different experiments in which we used exactly the same data/images for the
retrieval of the $W$ matrix and also the same test images. However, by
inserting an additional telescope straight after the SLM, we could control the
images sizes before the objective coupling into the fibre. By demagnifying the
input images by factors 1x, 1.4x and 4x, we expect the angular spread with
which the images are being focused to decrease. In these three cases, an
effective numerical aperture can be estimated by considering the objective
focal length ($f=34$ mm) and the image size at the input of the fiber.
Respectively, the sides of the images are: $10.5$ mm, $7.5$ mm and $2.6$ mm,
corresponding to NA${}_{eff}=0.22$, NA${}_{eff}=0.16$ and NA${}_{eff}=0.05$.
In all three case, the images were placed in a similar position on the fibre,
i.e. slightly off to the centre. Indeed, we noticed that placing the image
directly in the centre of the fibre led to very clear accumulation of the
speckles patterns towards the centre of the fibre, surrounded by broken ring-
like structures in the outer region. These ring-like structures are a clear
indication of higher modes, in keeping with previous reports also from other
groups [1]. However, the fact that these modes maintain their ring-like
structure as opposed to a more random speckle structure indicates that only a
few of the higher order modes have been excited. Conversely, the full speckle
pattern, distributed across the full fibre output, is obtained only by
exciting many modes. This condition was observed by displacing the beam
slightly to either side with respect to the central position. This is the
desired configuration as the objective here is to image with the highest
resolution possible: the fine features of any image are carried by the higher
spatial frequencies, which in turn correspond to the higher order modes in the
fibre.
The effect of changing the inout focusing condition can be clearly seen in
Fig. 1. When focusing with the smallest NAeff at the fibre input facet, the
speckles at the output are largest, corresponding to fewer modes and the final
retrieval is significantly worse when compared to the largest NAeff. The
intermediate NAeff shows slightly worse results with respect to the largest
NAeff. A consistent trend is therefore found between the effective NA, the
size and number of speckles at the fibre output and final image quality.

SUPPLEMENTARY NOTE 2: METHODS FOR IMAGE COMPARISON

Two methods have been considered in order to quantify the quality of our
predicted images: the structural similarity index (SSIM) [2] and the Parson
correlation coefficient (PCC). In both cases, a perfect match would correspond
to the maximum value 1. Considering two images $X$ and $Y$, we use the
definitions:

$SSIM(X,Y)=\frac{(2\mu_{X}\mu_{Y}+C_{1})(2\sigma_{XY}+C_{2})}{(\mu_{X}^{2}\mu_{Y}^{2}+C_{1})(\sigma_{X}^{2}\sigma_{Y}^{2}+C_{2})}$
(1)

where $\mu_{X}$ represent the average of $X$, $\mu_{Y}$ the average of $Y$,
$\sigma_{XY}$ the covariance of $X$ and $Y$, $\sigma_{X}^{2}$ the variance of
$X$ and $\sigma_{Y}^{2}$ the variance of $Y$. Whereas, $C_{1}$ and $C_{2}$ are
two parament defined as $C_{1}=(K_{1}L)**2$ and $C_{1}=(K_{2}L)**2$ where
$K_{1}$ was set to $0.01$ and $K_{2}$ to $0.03$ and $L$ is the dynamic range
of the image pixels. Instead, the Parson correlation coefficent is defined as:

$PCC(X,Y)=\frac{\sum_{i}(x_{i}-\bar{X})(y_{i}-\bar{Y})}{\sqrt{\sum_{i}(x_{i}-\bar{X})^{2}\sum_{i}(y_{i}-\bar{Y})^{2}}}$
(2)

where $x_{i}$ and $y_{i}$ indicate the pixels with index $i$ respectively of
the images $X$ and $Y$, and $\bar{X}$ (similarly with $Y$) the average of $X$.

SUPPLEMENTARY NOTE 3: DATA

For optimisation and testing of the model parameters, we use 50,000 images
from the ImageNet collection [3]. The training and experimental datasets are
supplied as additional material [4].

Supplementary Figure 2:  Additional examples of natural scenes imaged through
a mulitmode fiber. (a)-(b) Colour photographs of Glasgow University (imaged
through a 1 meter fibre). (c) Black and white photograph of a panda (imaged
through 10 m fibre). (d) Satellite photograph of Earth imaged at two different
times after the transmission of the initial training dataset (1 hour and 52
hours, imaged through 10 m fibre). (e) Comparison between green-colour-channel
speckles of the images in (d), at different times. The speckle patterns at
time 1 hour, 16 hours, 40 hours and 52 hours are shown on the first column
together with the SSIM parameters calculated respect to the time “1 hour”
speckle pattern. As can be seen, the 10 m length of the fibre leads to a slow
decorrelation of the speckle pattern, quantified here by the decreasing SSIM
with time. Similarly, the reconstructed green-channel-images show a slow
degradation of the SSIM with a linear dependence with respect to the speckle
SSIM degradation.

Validation image examples are images and videos from the Muybridge collection
such as a running horse, a jumping cat and a flying parrot. As explained in
the main text, the behaviour of the image retrieval seems to be largely
independent of the actual test images that was chosen. We focused mainly on
the Muybridge images in the main text but here show other examples in Fig. 2.
Figures 2(a) and (b) show colour photographs of the University of Glasgow,
imaged through a 1 m fibre. (c) is a grayscale image of a panda. (d) is a
satellite image of the Earth imaged through a 10 m fibre at a two different
times after the initial training an inversion process is completed showing
that the retrieved inversion matrix and setup can still be used at a later
time, as would be expected (yet still need to be verified) for a robust
inversion system. In (e) we report the speckle patterns relative to (d) at
different times (1 hour, 16 hours, 40 hours and 52 hours) of a single channel
of the RGB image. In order to allow a quantitative comparison, we used the
SSIM defined in the previous section. We note that the setup was not placed in
a specifically engineered or stabilised environment and was thus subject to
standard night-day temperature fluctuations (2-3 degrees) and environmental
vibrations. On the other hand, as we can see in Fig. 2-(e), the correlation
between the speckle pattern in time is still high, allowing a good
reconstruction. We judged out the scope of the present work to introduce
models able to deal with deep changes in the system transmission matrix, such
as in presence of relevant bending or temperature variations. Indeed, in
future work it would be interesting to explore the possibilities given by a
physics inspired artificial neural network also respect to these challenges.
In Fig. 3 we show a collection of images taken from the ImageNet database
together with their respective output speckle patterns (at the output of 1 m
long fibre) and final reconstructions. These images provide further evidence
for the image variability and robustness of the imaging reconstruction.

SUPPLEMENTARY NOTE 4: SOFTWARE

The code was developed in Python 3.6.5 with a standard Anaconda
http://www.anaconda.com/download configuration, including Keras [5] and
TensorFlow [6]. The code is supplied as additional supplementary material that
can be downloaded together with all of the training and experimental
data/images [4].

A. Model specification

The model is implemented as a simple complex, densely connected layer. The
individual weights are regularised with an $L_{2}$ minimising term, weighted
by $\lambda=0.03$. weights are initialised randomly, uniformly between $\pm
0.002$. Recorded images are collected into a training set of $N=45,000$ and
validation set of $5,000$.

Listing 1: Model specification

⬇

speckle_dim = 120

out_dim = 92

lamb = 0.03

inp = Input(shape=(speckle_dim**2,2))

comp = ComplexDense(out_dim**2, use_bias=False,

kernel_initializer=RandomUniform(-.002, .002),

kernel_regularizer=regularizers.l2(lamb))(inp)

amp = Amplitude()(comp)

model = Model(inputs=inp, outputs=amp)

The ComplexDense Layer is a custom layer we developed for Keras. It is a
straightforward Dense layer, but with complex-valued weights. The complex
weights are represented with Complex64 64-bit data types. Its only task is to
implement the complex-valued multiplication.

Listing 2: Custom ComplexDense Layer

⬇

class ComplexDense(Layer):

def __init__(self, output_dim,

activation=None,

use_bias=True,

kernel_initializer='glorot_uniform',

bias_initializer='zeros',

kernel_regularizer=None,

**kwargs):

super(ComplexDense, self).__init__(**kwargs)

self.output_dim = output_dim

self.activation = activations.get(activation)

self.use_bias = use_bias

self.kernel_initializer = initializers.get(kernel_initializer)

self.bias_initializer = initializers.get(bias_initializer)

self.kernel_regularizer = regularizers.get(kernel_regularizer)

def build(self, input_shape):

self.kernel = self.add_weight(name='kernel',

shape=(input_shape[1], self.output_dim, 2),

initializer=self.kernel_initializer,

regularizer=self.kernel_regularizer,

trainable=True)

if self.use_bias:

self.bias = self.add_weight(name='bias',

shape=(self.output_dim, 2),

initializer=self.bias_initializer,

trainable=True)

else:

self.bias = None

super(ComplexDense, self).build(input_shape)

def call(self, X):

# True Complex Multiplication (by channel combination)

complex_X = channels_to_complex(X)

complex_W = channels_to_complex(self.kernel)

complex_res = complex_X @ complex_W

if self.use_bias:

complex_b = channels_to_complex(self.bias)

= K.bias_add(complex_res, complex_b)

output = complex_to_channels(complex_res)

if self.activation is not None:

output = self.activation(output)

return output

def compute_output_shape(self, input_shape):

return (input_shape[0], self.output_dim, 2)

def get_config(self):

config = {'output_dim': self.output_dim,

'use_bias': self.use_bias,

'kernel_initializer': initializers.serialize(self.kernel_initializer),

'bias_initializer': initializers.serialize(self.bias_initializer),

'kernel_regularizer': regularizers.serialize(self.kernel_regularizer)

}

base_config = super(ComplexDense, self).get_config()

return dict(list(base_config.items()) + list(config.items()))

B. Parameter optimisation

The model fitting uses the standard Keras routines, involving a choice of
stochastic gradient descent and mean square error for the cost function. Here
the data variable x_train, y_train, x_validation, y_validation, x_text and
y_test refer to the amplitudes of, respectively, speckle patterns (x) and
original images (y).

Listing 3: Compile and optimise parameters

⬇

model.compile(optimizer=SGD(lr=1e-5), loss='mse', metrics=['mse'])

model_chk = ModelCheckpoint(weights_filepath, monitor='mse', verbose=0,

save_best_only=False,

save_weights_only=False, mode='auto', period=1)

reduce_lr = ReduceLROnPlateau(monitor='loss', factor=0.1, patience=2,

min_lr=lr/1e3, verbose=1,)

early_stop = EarlyStopping(monitor='loss', min_delta=0.0001, patience=8)

model.fit(x_train, y_train, validation_data = (x_validation, y_validation),

epochs = 850, batch_size = 32,

callbacks = [model_chk, reduce_lr, early_stop], shuffle = True)

Once the network parameters have converged (we ran the network for 850
iterations which takes ca 2 days on a PC with Nvidia TitanXp GPU card, you can
generate predictions of outputs using

⬇

pred_test = model.predict(x_test)**2

SUPPLEMENTARY REFERENCES

  1. ]

  2. [1

Papadopoulos, I. N., Farahi, S., Moser, C. & Psaltis, D Focusing and scanning
light through a multimode optical fiber using digital phase conjugation.
Optics Express 20, 10583–10590 (2012).

  3. [2

Wang, Z., et al. Image quality assessment: from error visibility to structural
similarity. IEEE transactions on image processing 13, 600–612 (2004).

  4. [3

Deng, J., et al. ImageNet: A Large-Scale Hierarchical Image Database. In
CVPR09 (2009).

  5. [4

Caramazza, P., Moran, O., Murray-Smith, R. & Faccio, D. Data can be downloaded
from the University of Glasgow repository: URL
DOI:http://dx.doi.org/10.5525/gla.researchdata.751

  6. [5

Chollet, F., et al. Keras. URL https://keras.io (2015).

  7. [6

Abadi, M. et al. TensorFlow: Large-scale machine learning on heterogeneous
systems (2015). URL https://www.tensorflow.org/. Software available from
tensorflow.org.

Supplementary Figure 3: A sample collection of images from the ImageNet
database, as projected into the fibre together with their relative output
speckle patterns and retrieved images. These testing images were not present
in the training dataset.


# Physics-guided deep reinforcement learning for flow field denoising

Mustafa Z. Yousif, Meng Zhang, Yifan Yang, Haifeng Zhou, Linqi Yu and HeeChang
Lim
School of Mechanical Engineering, Pusan National University
2, Busandaehak-ro 63beon-gil,Geumjeong-gu, Busan 46241, Republic of Korea
<EMAIL_ADDRESS>

###### Abstract

A multi-agent deep reinforcement learning (DRL)-based model is presented in
this study to reconstruct flow fields from noisy data. A combination of
reinforcement learning with pixel-wise rewards, physical constraints
represented by the momentum equation and the pressure Poisson equation and the
known boundary conditions is utilised to build a physics-constrained deep
reinforcement learning (PCDRL) model that can be trained without the target
training data. In thePCDRL model, each agent corresponds to a point in the
flow field and it learns an optimal strategy for choosing pre-defined actions.
The proposed model is efficient considering the visualisation of the action
map and the interpretation of the model operation. The performance of the
model is tested by utilising direct numerical simulation-based synthetic noisy
data and experimental data obtained by particle image velocimetry. Qualitative
and quantitative results show that the model can reconstruct the flow fields
and reproduce the statistics and the spectral content with commendable
accuracy. Furthermore, the dominant coherent structures of the flow fields can
be recovered by the flow fields obtained from the model when they are analyzed
using Proper orthogonal decomposition and dynamic mode decomposition. This
study demonstrates that the combination of DRL-based models and the known
physics of the flow fields can potentially help solve complex flow
reconstruction problems, which can result in a remarkable reduction in the
experimental and computational costs.

_K_ eywords Flow field denoising $\cdot$ PIV $\cdot$ Physics-guided DRL
$\cdot$ Machine learning

##  1 Introduction

The understanding of fluid flows plays a crucial role in life (for instance,
in medicine, construction, transportation, aerospace and astronomy). However,
fluid flow problems are usually complex with high non-linear behaviour,
especially turbulent flows, which occur at generally high Reynolds numbers. In
most cases, data from experiments and simulations are utilised to understand
and describe the behaviour of fluids with various accuracy levels that are
related to the experimental and numerical setups. Numerous methods have been
developed to improve the accuracy and practicality of the obtained flow
fields. However, several limitations still exist. One of the most notable
limitations of the experimental approach is the noise of the obtained flow
fields due to the experimental setup. Herein, obtaining measurements with an
acceptable signal-to-noise ratio (SNR) is practically impossible in some
cases. Therefore, several methods for the reconstruction of flow fields have
been introduced. Methods based on Linear data-driven approaches, such as
proper orthogonal decomposition [1] and dynamic mode decomposition [2], have
shown their capability to enhance the resolution of the flow data and filter
noisy flow data [3, 4, 5, 6, 7]. Additionally, various denoising methods for
particle image velocimetry (PIV) measurements, such as convolution filters,
wavelet methods and the Wiener filters, have revealed various levels of
success [8]. All the aforementioned methods showed limited success in terms of
denoising of flow fields because they are based on linear mapping or
handcrafted filtering processes, which are mostly incapable of dealing with
highly non-linear fluid problems [9].

With the recent rapid development in machine learning (ML) and graphic
processing unit, new data-driven methods have been introduced to provide
efficient solutions for problems in various fields, such as image processing,
natural language processing, robotics and weather forecasting. Several ML
algorithms have been recently utilised to address problems in fluid dynamics
and have shown promising results [9, 10, 11]. In contrast to linear methods,
ML-based techniques can deal with complex non-linear problems. This feature
has paved the way to explore the feasibility of applying ML to various
problems in complex turbulent flows [12, 13]. Several supervised and
unsupervised ML-based methods have been introduced considering flow
reconstruction from spatially limited or corrupted data [14]. Recently
promising results have been reported using deep learning (DL) by applying end-
to-end trained convolutional neural network (CNN)-based models [15, 16] and
generative adversarial network (GAN)-based models [17, 18, 19], where deep
learning is a subset of machine learning, in which neural networks with
multiple layers are utilised in the model [20]. Herein, the GAN-based models
have shown better performance than the traditional CNN-based models.
Nonetheless, the drawback of such methods lies in the need for the target
(high-resolution or uncorrupted) flow data to train the model, which is
difficult or impossible to obtain in most cases. Therefore, attempts have been
recently made to address this issue under certain conditions; for instance, in
the case of super-resolution reconstruction of randomly seeded flow fields
[21] or applying physical constraints in the loss function of the model to
reconstruct high-resolution steady flows from low-resolution noisy data [22].
However, the insufficient explainability and interpretability are the main
concerns of using ML-based methods, where no concrete explanation nor control
of the model performance is available.

On the other hand, Reinforcement Learning (RL), which is an ML method where an
agent learns to make decisions by interacting with an environment has shown
remarkable results in areas like robotics, game playing, and optimization
problems [23]. In RL The agent takes actions and receives feedback in the form
of rewards or penalties. Over time, it aims to learn the optimal actions to
maximize cumulative rewards and achieve its objectives through trial and
error. This approach to learning makes deep reinforcement learning (DRL) a
good candidate to be applied to several problems in fluid dynamics, such as
flow control [24], design optimisation [25], computational fluid dynamics [26]
and others [27, 28].

This paper presents a DRL-based approach that can be utilised for
reconstructing flow fields from noisy data. The main advantages of the
presented model lie in overcoming the necessity of the target data in the
training process and the explainable filtering process of the noisy data.

The remainder of this paper is organised as follows. Section 2 explains the
reconstruction methodology of denoised flow fields using the proposed DRL
model. Section 3 describes the generation and preprocessing of the data used
for training and testing the model. Section 4 discusses the results of testing
the proposed model. Finally, the conclusions of this study are presented in
Section 5.

##  2 Methodology

Different from supervised and unsupervised learning, reinforcement learning is
based on the Markov decision process, which is an iterative process where an
agent interacts with an environment. This process comprises four elements:
state $s$, action $a$, policy $\pi(a|s)$ and reward $r$. The action is an
operation that is applied by the agent. The policy represents the action
selection strategy of the agent. In other words, at each iteration step, the
agent obtains a state and chooses an action according to the policy. Owing to
the taken action, the state in the environment is then changed and the agent
receives an immediate reward, which is feedback showing the usefulness of the
taken action. The agent learns experience from the collected states, actions
and rewards after several iterations to find an optimal policy $\pi^{*}(a|s)$
that maximises the long-term reward. In DRL, a deep neural network is utilised
to obtain the optimal policy.

This study presents a physics-constrained deep reinforcement learning (PCDRL)
model that is built on the reinforcement learning with pixel-wise rewards
(PixelRL) [29], which is a CNN-based multi-agent DRL method for image
processing [30, 31, 32]. In PixelRL, the asynchronous advantage actor–critic
(A3C) algorithm [33] is applied for learning policies, which determine the
actions that are represented by the choice of basic filters for each pixel. In
other words, each pixel has one agent in PixelRL. A model that applies optimal
policies to change the velocity values is investigated in this study by
choosing the suitable actions for each point in the flow field in an instant,
as shown in Figure 1. In contrast to image processing problems that require
the target data in the training process [29], the physics of the flow
represented by the governing equations and the known boundary conditions are
utilised to train the model.

Figure 1: Learning process in the PCDRL model. Each agent at each iteration
step in the episode obtains a state from a point in the flow, calculates the
reward and applies an action according to the policy.

Let $\chi^{n}_{i,j}$ be the value of an instantaneous velocity component at
iteration step $n$ and in the location $(i,j)$ of the field. Herein, each
location has its agent with a policy $\pi_{i,j}(a^{n}_{i,j}|s^{n}_{i,j})$,
where $a^{n}_{i,j}\in\mathcal{A}$, which is a pre-defined set of actions
(Appendix D). Each agent obtains the next state, that is, $s^{n+1}_{i,j}$ and
reward $r^{n+1}_{i,j}$ from the environment by taking the action
$a^{n}_{i,j}$.

Physical constraints represented by the momentum equation and the pressure
Poisson equation and the known boundary conditions are embedded in the reward
function, which enables the model to follow an optimal denoising strategy that
results in changing the noisy data to the true flow field distribution. Hence,
the objective of the model is to learn the policy that maximises the expected
long-term rewards:

$\pi^{*}_{i,j}=\mathop{\mathrm{argmax}}\limits_{\pi_{i,j}}E_{\pi_{i,j}}(\displaystyle\sum_{n=1}^{N}\gamma^{(n-1)}r^{n}_{i,j}),$
(1)

where $\gamma^{(n-1)}$ is the $(n-1)$-th power of the discount factor
$\gamma$, which determines the weights of the immediate rewards in the
iteration steps. In this study, the value of $\gamma$ is set to 0.95.

The combination of the momentum equation,

$\frac{\partial\textbf{{u}}}{\partial
t}+(\textbf{{u}}\cdot\bm{\nabla})\textbf{{u}}=-\bm{\nabla}p+\nu\bm{\nabla}^{2}\textbf{{u}},$
(2)

and the pressure Poisson equation,

$\bm{\nabla}\cdot(\textbf{{u}}\cdot\bm{\nabla})\textbf{{u}}=-\bm{\nabla}^{2}p,$
(3)

is utilised to build the physics-based immediate reward,
$(r^{n}_{i,j})_{Physics}$, where u, $p$, $t$ and $\nu$ are the velocity
vector, pressure (divided by density), time and kinematic viscosity,
respectively.

At each iteration step, the pressure field is obtained by numerically solving
equation (3). Herein, the pressure gradient calculated from the pressure field
(($\bm{\nabla}{p}^{n}_{i,j})_{Poisson}$) is utilised in equation (2) such that

$(r^{n}_{i,j})_{Physics}=-|(\frac{\partial\textbf{{u}}}{\partial
t}+(\textbf{{u}}\cdot\bm{\nabla})\textbf{{u}}-\nu\bm{\nabla}^{2}\textbf{{u}})^{n}_{i,j}+(\bm{\nabla}{p}^{n}_{i,j})_{Poisson}|.$
(4)

The pressure integration in equation (3) is done by utlising a Poisson solver
that applies a standard 5-point scheme (second-order central difference
method)[34] with the initial pressure field being estimated from numerically
integrating the pressure gradient obtained from the initial noisy data in
equation (2) [35]. Notably, the central difference method is applied for all
the spatial discretisations. Regarding the temporal discretisation, for the
first and the last time steps in each training mini-batch, the forward
difference and the backward difference are used, respectively, and the central
difference is applied for the other time steps. Furthermore, Neumann and
Dirichlet boundary conditions according to each case used in this study are
enforced in the calculations.

Additionally, the velocity values obtained after each action $a^{n}_{i,j}$ are
directly made divergence-free by applying Helmholtz–Hodge decomposition [36]
using Fourier transformation. Furthermore, the known boundary conditions are
utilised to obtain the boundary conditions-based immediate reward
$(r^{n}_{i,j})_{BC}$ for the velocity by considering the absolute error of the
reconstructed data at the boundaries of the domain.

Thus, the combined immediate reward function can be expressed as

$r^{n}_{i,j}=(r^{n}_{i,j})_{Physics}+\beta(r^{n}_{i,j})_{BC},$ (5)

where $\beta$ is a weight coefficient and its value is empirically set to 20.

This approach considers the convergence of the model output to satisfy the
governing equations and boundary conditions as a measure of the model
performance without the need for the target training data. Furthermore, the
reward function is designed to mimic the denoising process of PIV velocity
field data without the need for measured pressure field data in the model.
Nine iteration steps for each episode, that is, $N=9$, are used in this study.
In addition, the size of the training mini-batch is set to 4. The model is
applied to direct numerical simulation (DNS)-based data (corrupted by
different levels of additive zero-mean Gaussian noise) and real noisy PIV data
of two-dimensional flow around a square cylinder at Reynolds number, $Re_{D}$
= 100 and 200, respectively. Herein, $Re_{D}=u_{\infty}D/\nu$, where
$u_{\infty}$ and $D$ are the free-stream velocity and the cylinder width,
respectively. Details regarding the source code of the proposed model, A3C,
PixleRL and the selected pre-defined denoising action set can be found in
appendices A, B, C and D, respectively.

##  3 Data description and preprocessing

###  3.1 Synthetic data

DNS data of a two-dimensional flow around a square cylinder at a Reynolds
number, $Re_{D}=100$, are considered as an example of synthetic data. The
open-source computational fluid dynamics finite-volume code OpenFOAM-5.0x is
used to perform the DNS. The domain size is set to be $(x_{D}\times
y_{D})=(20\times 15)$, where $x$ and $y$ are the streamwise and spanwise
directions, respectively. The corresponding grid size = $(381\times 221)$.
Local mesh refinement is applied using the stretching mesh technique near the
cylinder walls. Uniform inlet velocity and pressure outlet boundary conditions
are applied to the inlet and the outlet of the domain, respectively. No-slip
boundary condition is applied to the cylinder walls and the symmetry plane to
the sides of the domain. The dimensionless time step of the simulation, that
is, $u_{\infty}\Delta t/D$, is set to $10^{-2}$. The DNS data are corrupted by
additive zero-mean Gaussian noise, that is,
$\mathcal{S}\sim\mathcal{N}(0,\sigma^{2})$, where $\mathcal{S}$, $\mathcal{N}$
and $\sigma^{2}$ represent the noise, the normal distribution and the
variance, respectively. The signal-to-noise ratio (SNR), wherein a large SNR
yields a low noise level, is used to evaluate the noise level. Herein,
SNR$=\sigma^{2}_{DNS}/\sigma^{2}_{noise}$, where $\sigma^{2}_{DNS}$ and
$\sigma^{2}_{noise}$ denote the variance of the DNS and the noise data,
respectively. Three levels of noise are applied, $1/{\rm SNR}=$ 0.01, 0.1 and
1. The interval between the collected snapshots of the flow fields is set to
be 10 times the simulation time step; 1000 snapshots are used for training the
model, whereas 200 snapshots are used for testing the performance of the
model.

###  3.2 Experimental data

Two PIV experiments are performed to generate noisy and clear (uncorrupted)
data (for comparison) of flow over a square cylinder to investigate the
performance of the proposed PCDRL model on real experimental data. The noisy
data are generated by utilising a return-type water channel. The test section
size of the water channel is 1 m (length) $\times$ 0.35 m (height) $\times$
0.3 m (width). The free-stream velocity is set to 0.02 m/s, with the
corresponding $Re_{D}$ of 200. The background noise is generated at relatively
high levels due to the external noise and the sparse honeycomb of the water
channel. The channel was seeded by polyamide12 seed particles from INTECH
SYSTEMS with 50 $\mu$m diameter. The high-speed camera (FASTCAM Mini UX 50)
and a continuous laser with a 532 nm wavelength are utilised to build the
complete PIV system. The snapshot frequency is set to 24 Hz. Herein, 2,000 and
500 instantaneous flow fields are used for the model training and testing of
its performance, respectively. Meanwhile, clear data of the flow are generated
by utilising a return-type wind tunnel. The test section size of the wind
tunnel is 1 m (length) $\times$ 0.25 m (height) $\times$ 0.25 m (width). The
free-stream velocity is set to 0.29 m/s, with the corresponding $Re_{D}$ of
200. The turbulence intensity of the free-stream is less than $0.8\%$. The
wind tunnel is seeded by the olive oil droplets generated by TSI 9307 particle
generator. The PIV system used in the wind tunnel comprises a two-pulsed laser
(Evergreen, EVG00070) and a CCD camera (VC-12MX) with $4,096\times 3,072$
pixel resolution. Herein, the snapshot frequency is set to 15 Hz. In the water
channel and wind tunnel experiments, the square cylinder model comprised an
acrylic board and the cross-section of the model is set to 1 cm $\times$ 1 cm.
The model is not entirely transparent. Thus, a shadow region is generated in
the area below the bluff body when the laser goes through the model.

##  4 Results and discussion

###  4.1 Performance of the model

The capability of the PCDRL model to denoise flow fields is investigated in
this study qualitatively and quantitatively by utilising the DNS and PIV data.
The model is primarily applied to DNS-based data. Figure 2 shows the progress
of the mean reward during the training process, that is,

$\bar{r}=\frac{1}{IJN}\sum_{i=1}^{I}\sum_{j=1}^{J}\sum_{n=1}^{N}r^{n}_{i,j}.$
(6)

The solid line and light area indicate $\bar{r}$ and the standard deviation of
the reward on nine iteration steps, respectively. As shown in the figure, the
reward for the three different noise levels rapidly increases and approaches
its optimal level after a few number of episodes. This finding indicates that
the agents in PixelRL learn the policy in a few episodes in the training
process compared with the other multi-agent networks because they share the
information represented by the network parameters and also due to the averaged
gradients [29]. Thus, this approach can remarkedly reduce the computational
cost of the model. Furthermore, as expected, the magnitude of the optimal
$\bar{r}$ decreases with the increase in the noise level.

Figure 2: Progress of the mean reward during the training process. Cases 1, 2
and 3 represent the noisy DNS data at noise levels $1/{\rm SNR}=$ 0.01, 0.1
and 1, respectively.

Figure 3 shows a visual overview of the prediction process of the PCDRL model.
The figure reveals that the choice of filters changes with the spatial
distribution of the velocity data and also with each iteration step in the
episode. The visualisation of the action map is one of the model features,
providing additional access to the model considering the action strategy.
Furthermore, it can be seen that the action map is strongly correlated with
the physics of the flow, which is represented in this case by the vortex
shedding behind the square cylinder.

Figure 3: Action map of the prediction process for an instantaneous streamwise
velocity field. The top part of the figure shows the types of filters used in
the process and the action map in each iteration step and the bottom part of
the figure shows the corresponding velocity field. Result for DNS noisy data
at noise level $1/{\rm SNR}=$ 0.1.

The instantaneous denoised flow data are presented in Figure 4(a) by employing
the vorticity field ($\omega$). The figure reveals that the model shows a
remarkable capability to reconstruct the flow field even when using an extreme
level of noise in the input data of the model.

The general reconstruction accuracy of the model is examined via the relative
$L_{2}$-norm error of the reconstructed velocity fields,

$\epsilon(\chi)=\frac{1}{K}\sum_{k=1}^{K}\frac{||\chi^{PCDRL}_{k}-\chi^{DNS}_{k}||_{2}}{||\chi^{DNS}_{k}||_{2}},$
(7)

where $\chi^{PCDRL}_{k}$ and $\chi^{DNS}_{k}$ represent the predicted velocity
component and the ground truth (DNS) one, respectively, and $K$ is the number
of test snapshots. Figure 4(b) shows that the values of the error are
relatively small for the velocity components and are proportional to the
increase in the noise level.

Figure 5 shows the probability density function (p.d.f.) plots of the
streamwise ($u$) and spanwise ($v$) velocity components. Herein, the p.d.f.
plots obtained from the reconstructed velocity fields are generally consistent
with those obtained from DNS, indicating that the proposed model could
successfully recover the actual distribution of flow data. Furthermore, the
scatter plots of the maximum instantaneous velocity values in all the test
data set are presented in Figure 6. The figure reveals that the predicted data
are generally in commendable agreement with the DNS data for the entire range
of each velocity component with a slight reduction in the consistency as the
noise level increases.

Figure 4: (a) Instantaneous vorticity field; (b) Relative $L_{2}$-norm error
of the reconstructed velocity fields. Cases 1, 2 and 3 represent the noisy DNS
data at noise levels $1/{\rm SNR}=$ 0.01, 0.1 and 1, respectively.

Figure 5: Probability density function plots of the (a) streamwise and (b)
spanwise velocity components. Cases 1, 2 and 3 represent the noisy DNS data at
noise levels $1/{\rm SNR}=$ 0.01, 0.1 and 1, respectively.

Figure 6: Scatter plots of the maximum instantaneous values of the (a)
streamwise and (b) spanwise velocity components. Cases 1, 2 and 3 represent
the results from the PCDRL model using noisy DNS data at noise levels $1/{\rm
SNR}=$ 0.01, 0.1 and 1, respectively. The contour colours (from blue to red)
are proportional to the density of points in the scatter plot.

The power spectral density (PSD) of the streamwise velocity fluctuations at
two different locations is plotted in Figure 7 to examine the capability of
the model to reproduce the spectral content of the flow. Herein, a commendable
agreement with the DNS results can be observed with a slight deviation in the
high frequencies for the noise level, $1/{\rm SNR}=1$.

Figure 7: Power spectral density plots of the streamwise velocity fluctuations
at two different locations: (a) ($x/D,y/D$) = (1,1) and (b) ($x/D,y/D$) =
(6,1). The dimensionless frequency is represented by the Strouhal number,
$St=fD/u_{\infty}$, where $f$ is the frequency. Cases 1, 2 and 3 represent the
noisy DNS data at noise levels $1/{\rm SNR}=$ 0.01, 0.1 and 1, respectively.

The statistics of the velocity fields, represented by the spanwise profiles of
the root mean square of the velocity ($u_{rms},v_{rms}$) and Reynolds shear
stress ($\overline{u^{\prime}v^{\prime}}$), are presented in Figure 8. The
figure reveals an accurate reconstruction of the statistics at two different
streamwise locations in the domain, indicating that the model could
successfully reproduce the statistics of the flow despite of the extreme noise
level.

Figure 8: Spanwise profiles of flow statistics ($u_{rms}$ (left), $v_{rms}$
(middle) and $\overline{u^{\prime}v^{\prime}}$ (right)) at two different
streamwise locations. (a) $x/D=3$; (b) $x/D=6$. Cases 1, 2 and 3 represent the
results from the PCDRL model using noisy DNS data at noise levels $1/{\rm
SNR}=$ 0.01, 0.1 and 1, respectively.

The model performance is further examined by utilising actual noisy PIV data.
The reconstructed instantaneous vorticity field is shown in Figure 9(a). The
figure reveals that the model could successfully denoise the velocity fields
with commendable accuracy considering the input noisy data to the model. In
addition, the model shows the capability of recovering the corrupted regions
in the flow due to the experimental setup. Furthermore, the relative
difference of the spanwise profile of the vorticity root mean square
($\omega_{rms}$) between the reconstructed data and the clear PIV data
($\varepsilon(\omega_{rms})$) presented in Figure 9(b) shows that the results
from the model reveal a smooth behaviour that is generally consistent with
that from the clear PIV data. These results indicate that the PCDRL model can
be practically applied to noisy PIV data.

Figure 9: (a) Instantaneous vorticity field of the noisy (left) and denoised
PIV data by the PCDRL (right); (b) Relative difference of the spanwise profile
of the vorticity root mean square at two different streamwise locations.

###  4.2 POD and DMD results

In this section, the accuracy of the results from the PCDRL model is examined
in terms of flow decomposition. Firstly, the results of applying POD to the
denoised data are compared with the POD results of the ground truth data.
Figure 10 shows the contour plots of the leading POD modes for the vorticity
field obtained from the DNS data. As can be observed from the figure, for the
case of the highest noise level, i.e 1/SNR = 1, all the seven true leading
modes can be recovered using the denoised data, while only three modes can be
recovered using the noisy data and no distinguishable features can be seen for
the other modes. Furthermore, the energy plots in Figure 11 represented by the
normalised POD eigenvalues show that even for the case of the flow with the
highest noise level, the energy contribution values of the POD modes are
consistent with those obtained from the ground truth DNS data. As expected,
the results from the noisy data reveal a different behaviour, especially for
the cases of noise levels 1/SNR = 0.1 and 1. Figure 12 shows a reconstructed
instantaneous vorticity field of the DNS data using the first ten POD modes.
As shown in the figure, the result from PCDRL model reveals a commendable
reconstruction accuracy as compared with the ground truth DNS results, whereas
the obtained result from the noisy data indicates the limitation of POD to
recover the flow with the right physics.

Figure 10: Leading POD modes obtained from the DNS data. Results from ground
truth DNS (left), PCDRL (middle), and noisy data with 1/SNR = 1 (right).

Figure 11: Normlised energy (left) and cumulative energy (right) of the POD
modes obtained from the DNS data; (a) Noisy data, where Cases 1, 2 and 3
represent the noisy DNS data at noise levels $1/{\rm SNR}=$ 0.01, 0.1 and 1,
respectively. (b) Results from the PCDRL model.

Figure 12:  Reconstructed instantaneous vorticity field obtained from the DNS
data using first ten POD modes. Cases 1, 2 and 3 represent the results of
using noisy DNS data at noise levels $1/{\rm SNR}=$ 0.01, 0.1 and 1,
respectively.

Similar results can be obtained by applying the POD to the PIV data. As can be
observed from Figure 13, the seven leading POD modes obtained from the
denoised PIV data are relatively consistent with the modes obtained from the
clear PIV data, considering that the clean PIV data are obtained using a
different experimental setup, whereas the noisy PIV data fails to recover the
modes after the third mode. Notably, the shadow region is clearly shown in
some of the modes obtained from the clear PIV data, whereas no such region can
be seen in the modes obtained from the denoised data. This is consistent with
results from Figure 9(a). The results from Figure 14 further indicate the
ability of the model to reconstruct the flow data with POD modes that
generally have a behaviour similar to that of the clear PIV data. Furthermore,
as shown in Figure 15, the reconstructed instantaneous vorticity field using
the first ten modes of the denoised data shows a realistic flow behaviour that
is expected from the case of flow around a cylinder.

Figure 13: Leading POD modes obtained from the PIV data. Results from clear
PIV (left), PCDRL (middle), and noisy PIV data (right).

Figure 14: Normlised energy (left) and cumulative energy (right) of the POD
modes obtained from the PIV data.

Figure 15: Reconstructed instantaneous vorticity field obtained from the
results of the PCDRL model (left) and noisy PIV data (right) using first ten
POD modes.

To further investigate the dynamics of the denoised flow data, DMD is then
applied to the flow data. As shown in Figure 16, even in the case of the DNS
data corrupted with the level of noise 1/SNR = 1, the DMD eigenvalues of the
vorticity field show a behaviour close to that of the ground truth DNS data,
whereas for the noisy data, the eigenvalues scatter inside the unit circle
plot indicating a non-realistic behaviour of the system.

As for the denoised PIV data, Figure 17 reveals that the eigenvalues also show
a well-agreement with those from the clear PIV data. Notably, the leading DMD
eigenvalues of the clear PIV data are not exactly located on the circumference
of the unit circle as in the case of DNS data. This behaviour can be
attributed to the fact that DMD is known to be sensitive to noise [37, 38, 39,
7], and unlike the DNS data, the clear PIV contains a relatively small level
of noise, which can affect the flow decomposition.

Figure 16: DMD eigenvalues of the noisy DNS data at noise level 1/SNR = 1
(left) and the results from the PCDRL model (right) visualised on the unit
circle.

Figure 17: DMD eigenvalues of the noisy PIV data (left) and the results from
the PCDRL model (right) visualised on the unit circle.

##  5 Conclusions

This study proposed a DRL-based method to reconstruct flow fields from noisy
data. The PixelRL method is utilised to build the proposed PCDRL model,
wherein an agent that applies actions represented by basic filters according
to a local policy is assigned for each point in the flow. Hence, the proposed
model is a multi-agent model. The physical constraints represented by the
momentum equation and the pressure Poisson equation and the boundary
conditions are utilised to build the reward function. Hence, the PCDRL model
is label-training data-free, that is, target data are not required for the
model training. Furthermore, the visualisation and interpretation of the model
performance can be easily achieved due to the model setup.

The model performance was firstly investigated through DNS-based noisy data
with three different noise levels. The instantaneous results and the flow
statistics revealed a commendable reconstruction accuracy of the model.
Furthermore, the spectral content of the flow was favourably recovered by the
model with reduced accuracy with the increase in the noise level.
Additionally, the reconstruction error showed relatively low values,
indicating the general reconstruction accuracy of the model.

Real noisy and clear PIV data are utilised to examine the model performance.
Herein, the model demonstrates its capability to recover the flow fields with
the appropriate behaviour.

Furthermore, the accuracy of the denoised flow data from both DNS and PIV was
investigated in terms of flow decomposition by means of POD and DMD. Herein,
most of the leading POD modes that describe the main features (coherent
structures) of the flow were successfully recovered with commendable accuracy
and outperformed the results of directly applying POD to the noisy data.
Additionally, the DMD eigenvalues obtained from the denoised flow data
revealed a behaviour that is similar to that of true DMD modes. These results
further indicate the model’s ability to recover the flow data with most of the
flow physics.

This study demonstrates that the combination of DRL, the physics of the flow,
which is represented by the governing equations, and prior knowledge of the
flow boundary conditions, can be effectively utilised to recover high-fidelity
flow fields from noisy data. This approach can be further extended to be
utlised in the reconstruction of three-dimensional turbulent flow fields.
Herein, more sophisticated DRL models with more complex spatial filters are
needed. Applying such models to flow reconstruction problems can result in a
considerable reduction in the experimental and computational costs.

## Acknowledgments

This work was supported by ‘Human Resources Program in Energy Technology’ of
the Korea Institute of Energy Technology Evaluation and Planning (KETEP),
granted financial resource from the Ministry of Trade, Industry & Energy,
Republic of Korea (no. 20214000000140). In addition, this work was supported
by the National Research Foundation of Korea (NRF) grant funded by the Korea
government (MSIP) (no. 2019R1I1A3A01058576)

##  Appendix A Open-source code

The open-source library Pytorch 1.4.0 [40] is utilised for the implementation
of the model. The source code of the proposed PGDRL model is available at
click here.

##  Appendix B Asynchronous Advantage Actor–Critic (A3C)

The asynchronous advantage actor–critic (A3C) [33] algorithm is applied in
this paper. A3C is a variant of the actor–critic algorithm, which combines the
policy- and value-based networks to improve its performance. Figure B1 shows
that the actor generates an action $a^{n}$ for the given state $s^{n}$ based
on the current policy, whilst the critic provides the value function
$V(s^{n})$ to evaluate the effectiveness of the action.

Figure B1: Architecture of the actor–critic algorithm.

Based on the understanding of the actor–critic algorithm, A3C also has two
sub-networks: policy and value networks. Herein, $\theta_{p}$ and $\theta_{v}$
are used to represent the parameters of each network. The gradients of
$\theta_{p}$ and $\theta_{v}$ can be calculated as follows:

$R^{n}=r^{n}+\gamma r^{n+1}+\cdots+\gamma^{(N-1)}r^{N}+\gamma^{(N)}V(s^{N}),$
(B1) $d\theta_{v}=\bigtriangledown_{\theta_{v}}(R^{n}-V(s^{n}))^{2},$ (B2)
$A(a^{n},s^{n})=R^{n}-V(s^{n}),$ (B3)
$d\theta_{p}=-\bigtriangledown_{\theta_{p}}\log\pi(a^{n},s^{n})A(a^{n},s^{n}),$
(B4)

where $A(a^{n},s^{n})$ is the advantage.

##  Appendix C PixelRL

A3C is modified in this study to a fully convolutional form [29], and its
architecture can be found in Figure C1. Through this approach, all the agents
share the same parameters, which saves computational cost and trains the model
more efficiently compared to the case where agents need to train their models
individually. The size of the receptive field can also affect the performance
of the CNN network, and a large receptive field can result in superior capture
connections between points. Therefore, a receptive field ($3\times 3$) is used
in the architecture; that is, the outputs of the policy and value networks at
a specific pixel will be affected by the pixel and its surrounding neighbor
pixels. Figure C1 shows that the input flow field data firstly pass through
four convolutional and Leaky rectified linear unit (ReLU) [41] layers and are
then inputted to the policy and value networks, respectively. The policy
network comprises three convolutional layers with a ReLU activation function,
a ConvGRU layer and a convolutional layer with a SoftMax activation function
[41], and its output is the policy. The first three layers of the value
network are the same as the first three layers of the policy network, and the
value function is finally obtained through the convolutional layer with a
linear function.

Figure C1: Architecture of the fully convolutional A3C.

The gradient of the parameters $\theta_{p}$ and $\theta_{v}$ is then defined
on the basis of the architecture of the fully convolutional A3C.

$\textbf{{R}}^{n}=\textbf{{r}}^{n}+\gamma\textbf{{W}}\ast\textbf{{r}}^{n+1}+\cdots+\gamma^{(N-1)}\textbf{{W}}^{(N-1)}\ast\textbf{{r}}^{N}+\gamma^{(N)}\textbf{{W}}^{(N)}\ast\textbf{{V}}(\textbf{{s}}^{N}),$
(C1) $d\theta_{v}=\bigtriangledown_{\theta_{v}}\frac{1}{I\times
J}\textbf{1}^{\top}\\{(\textbf{{R}}^{n}-\textbf{{V}}(\textbf{{s}}^{n}))\odot(\textbf{{R}}^{n}-\textbf{{V}}(\textbf{{s}}^{n}))\\}\textbf{1},$
(C2)
$\textbf{{A}}(\textbf{{a}}^{n},\textbf{{s}}^{n})=\textbf{{R}}^{n}-\textbf{{V}}(\textbf{{s}}^{n}),$
(C3) $d\theta_{p}=-\bigtriangledown_{\theta_{p}}\frac{1}{I\times
J}\textbf{1}^{\top}\\{\log\bm{\pi}(\textbf{{a}}^{n},\textbf{{s}}^{n})\odot\textbf{{A}}(\textbf{{a}}^{n},\textbf{{s}}^{n})\\}\textbf{1}.$
(C4)

where $\textbf{{R}}^{n}$, $\textbf{{r}}^{n}$,
$\textbf{{V}}(\textbf{{s}}^{n})$,
$\textbf{{A}}(\textbf{{a}}^{n},\textbf{{s}}^{n})$ and
$\bm{\pi}(\textbf{{a}}^{n},\textbf{{s}}^{n})$ are the matrices whose
$(i,j)$-th elements are $R^{n}_{i,j}$, $r^{n}_{i,j}$, $V(s^{n}_{i,j})$,
$A(a^{n}_{i,j},s^{n}_{i,j})$ and $\pi(a^{n}_{i,j},s^{n}_{i,j})$, respectively.
$\ast$ is the convolution operator, 1 is the all-ones vector and $\odot$ is
the element-wise multiplication. W is a convolution filter weight, which is
updated simultaneously with the $\theta_{p}$ and $\theta_{v}$ such that

$d\textbf{{W}}=-\bigtriangledown_{\textbf{{W}}}\frac{1}{I\times
J}\textbf{1}^{\top}\\{\log\bm{\pi}(\textbf{{a}}^{n},\textbf{{s}}^{n})\odot\textbf{{A}}(\textbf{{a}}^{n},\textbf{{s}}^{n})\\}\textbf{1}+\bigtriangledown_{\textbf{{W}}}\frac{1}{I\times
J}\textbf{1}^{\top}\\{(\textbf{{R}}^{n}-\textbf{{V}}(\textbf{{s}}^{n}))\odot(\textbf{{R}}^{n}-\textbf{{V}}(\textbf{{s}}^{n}))\\}\textbf{1}.$
(C5)

Notably, after the agents complete their interaction with the environment, the
gradients are acquired simultaneously, which means that the number of
asynchronous threads is one; that is, A3C is equivalent to advantage
actor–critic (A2C) in the current study [42].

##  Appendix D Denoising action set

The action set for removing the noise of the flow fields is shown in Table D1.
The agent can apply the following nine possible actions: do nothing, six
classical image filters and plus/minus a Scalar. The actions in this study are
discrete and determined empirically. The table shows that the parameters
$\sigma_{c}$, $\sigma_{s}$ and $\sigma$ represent the filter standard
deviation in the color space, the coordinate space and the Gaussian kernel,
respectively. The Scalar in the $8^{th}$ and $9^{th}$ actions is determined on
the basis of the difference between the variance of the clear and noisy data.

Table D1: Action set for the denoising process Number | Action type | Parameters
---|---|---
1 | Do nothing | filter size = $5\times 5$
2 | Box filter | filter size = $5\times 5$
3 | Bilateral filter1 | filter size = $5\times 5$, $\sigma_{c}=1.0,\sigma_{s}=5.0$
4 | Bilateral filter2 | filter size = $5\times 5$, $\sigma_{c}=0.1,\sigma_{s}=5.0$
5 | Median filter | filter size = $5\times 5$
6 | Gaussian filter1 | filter size = $5\times 5$, $\sigma=1.5$
7 | Gaussian filter2 | filter size = $5\times 5$, $\sigma=0.5$
8 | Value += Scalar |
9 | Value -= Scalar |

## References

  * [1] J. L. Lumley.  The structure of inhomogeneous turbulent flows.  Atmospheric Turbulence and Radio Wave Propagation., pages 166–177, 1967.
  * [2] P. J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics, 656:5–28, 2010.
  * [3] H. Gunes and U. Rist.  Spatial resolution enhancement/smoothing of stereo–particle-image-velocimetry data using proper-orthogonal-decomposition–based and kriging interpolation methods.  Physics of Fluids., 19(6):064101, 2007.
  * [4] C. He and Y. Liu.  Proper orthogonal decomposition-based spatial refinement of tr-piv realizations using high-resolution non-tr-piv measurements.  Experiments in Fluids., 58(7):86, 2017.
  * [5] Mojtaba F. Fathi, Ali Bakhshinejad, Ahmadreza Baghaie, and Roshan M. D’Souza.  Dynamic denoising and gappy data reconstruction based on dynamic mode decomposition and discrete cosine transform.  Applied Sciences, 8(9), 2018.
  * [6] Taku Nonomura, Hisaichi Shibata, and Ryoji Takaki.  Extended-kalman-filter-based dynamic mode decomposition for simultaneous system identification and denoising.  PLOS ONE, 14:1–46, 2019.
  * [7] Isabel Scherl, Benjamin Strom, Jessica K. Shang, Owen Williams, Brian L. Polagye, and Steven L Brunton.  Robust principal component analysis for modal decomposition of corrupt fluid flows.  Physical Review Fluids., 5(5):054401, 2020.
  * [8] J. Vétel, A. Garon, and D. Pelletier.  Denoising methods for time-resolved piv measurements.  Experiments in Fluids., 51(4):893–916, 2011.
  * [9] S. L. Brunton, B. R. Noack, and P. Koumoutsakos.  Machine learning for fluid mechanics.  Annual Review of Fluid Mechanics., 52(1):477–508, 2020.
  * [10] K. Duraisamy, G. Iaccarino, and H. Xiao.  Turbulence modeling in the age of data.  Annual Review of Fluid Mechanics., 51(1):357–377, 2019.
  * [11] R. Vinuesa and S. L. Brunton.  Enhancing computational fluid dynamics with machine learning.  Nature Computational Science., 2(6):358–366, 2022.
  * [12] L. Guastoni, A. Güemes, A. Ianiro, S. Discetti, P. Schlatter, H. Azizpour, and R. Vinuesa.  Convolutional-network models to predict wall-bounded turbulence from wall quantities.  Journal of Fluid Mechanics., 928:A27, 2021.
  * [13] Mustafa Z. Yousif, Meng Zhang, Linqi Yu, Ricardo Vinuesa, and HeeChang Lim.  A transformer-based synthetic-inflow generator for spatially developing turbulent boundary layers.  Journal of Fluid Mechanics., 957:A6, 2023.
  * [14] S. Discetti and Y. Liu.  Machine learning for flow field measurements: A perspective.  Measurement Science and Technology., 34(2):021001, 2022.
  * [15] K. Fukami, K. Fukagata, and K. Taira.  Super-resolution reconstruction of turbulent flows with machine learning.  Journal of Fluid Mechanics., 870:106–120, 2019.
  * [16] B. Liu, J. Tang, H. Huang, and X.-Y. Lu.  Deep learning methods for super-resolution reconstruction of turbulent flows.  Physics of Fluids., 32(2):025105, 2020.
  * [17] M. Z. Yousif, L. Yu, S. Hoyas, R. Vinuesa, and H. Lim.  A deep-learning approach for reconstructing 3D turbulent flows from 2D observation data.  Sci Rep, 13:2529, 2023.
  * [18] H. Kim, J. Kim, S. Won, and C. Lee.  Unsupervised deep learning for super-resolution reconstruction of turbulence.  Journal of Fluid Mechanics., 910:A29, 2021.
  * [19] L. Yu, M. Z. G. Yousif, M. Zhang, S. Hoyas, R. Vinuesa, and H. Lim.  Three-dimensional enhanced super-resolution generative adversarial network for super-resolution reconstruction of turbulent flows with tricubic interpolation-based transfer learning.  Physics of Fluids., 34:125126, 2022.
  * [20] Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  Nature., 521:436–444, 2015.
  * [21] A. Güemes, C. S. Vila, and S. Discetti.  Super-resolution gans of randomly-seeded fields.  Nat Mach Intell., 4:1165–1173, 2022.
  * [22] H. Gao, L. Sun, and J.-X. Wang.  Super-resolution and denoising of fluid flow using physics-informed convolutional neural networks without high-resolution labels.  Physics of Fluids., 33(7):073603, 2021.
  * [23] Thomas Hickling, Abdelhafid Zenati, Nabil Aouf, and Phillippa Spencer.  Explainability in deep reinforcement learning, a review into current methods and applications.  CoRR, abs/2207.01911, 2022.
  * [24] J. Rabault, M. Kuchta, A. Jensen, U. Réglade, and N. Cerardi.  Artificial neural networks trained through deep reinforcement learning discover control strategies for active flow control.  Journal of Fluid Mechanics., 865:281–302, 2019.
  * [25] Jonathan Viquerat, Jean Rabault, Alexander Kuhnle, Hassan Ghraieb, Aurélien Larcher, and Elie Hachem.  Direct shape optimization through deep reinforcement learning.  Journal of Computational Physics, 428:110080, 2021.
  * [26] G. Novati, H.L. de Laroussilhe, and P. Koumoutsakos.  Automating turbulence modelling by multi-agent reinforcement learning.  Nature Machine Intelligence, 3:87–96, 2021.
  * [27] Paul Garnier, Jonathan Viquerat, Jean Rabault, Aurélien Larcher, Alexander Kuhnle, and Elie Hachem.  A review on deep reinforcement learning for fluid mechanics.  Computers and Fluids, 225:104973, 2021.
  * [28] J. Viquerat, P. Meliga, A. Larcher, and E. Hachem.  A review on deep reinforcement learning for fluid mechanics: An update.  Physics of Fluids, 34:111301, 2022.
  * [29] R. Furuta, N. Inoue, and T. Yamasaki.  Pixelrl: Fully convolutional network with reinforcement learning for image processing.  IEEE Transactions on Multimedia., 22(7):1704–1719, 2020.
  * [30] P. Jarosik, M. Lewandowski, Z. Klimonda, and M. Byra.  Pixel-wise deep reinforcement learning approach for ultrasound image denoising.  2021 IEEE International Ultrasonics Symposium (IUS)., pages 1–4, 2021.
  * [31] W. Li, X. Feng, H. An, X. Y. Ng, and Y.-J. Zhang.  Mri reconstruction with interpretable pixel-wise operations using reinforcement learning.  Proceedings of the AAAI Conference on Artificial Intelligence., 34(01):792–799, 2020.
  * [32] K. Vassilo, C. Heatwole, T. Taha, and A. Mehmood.  Multi-step reinforcement learning for single image super-resolution.  2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)., pages 2160–2168, 2020.
  * [33] V. Mnih, A. P. Badia, M. Mirza, A. Graves, T. Lillicrap, T. Harley, D. Silver, and K. Kavukcuoglu.  Asynchronous methods for deep reinforcement learning.  Proceedings of The 33rd International Conference on Machine Learning., pages 1928–1937, 2016.
  * [34] J. W. Van der Kindere, A. Laskari, B. Ganapathisubramani, and R. de Kat.  Pressure from 2D snapshot PIV.  Experiments in Fluids., 60:32, 2019.
  * [35] Bas W. van Oudheusden, Francesco Scarano, E. W. M. Roosenboom, E. W. F. Casimiri, and Louis Souverein.  Evaluation of integral forces and pressure fields from planar velocimetry data for incompressible and compressible flows.  Experiments in Fluids, 43:153–162, 2007.
  * [36] H. Bhatia, G. Norgard, V. Pascucci, and P.-T Bremer.  The helmholtz-hodge decomposition—a survey.  IEEE Transactions on Visualization and Computer Graphics., 19(8):1386–1404, 2013.
  * [37] Shervin Bagheri.  Effects of weak noise on oscillating flows: Linking quality factor, Floquet modes, and Koopman spectrum.  Physics of Fluids, 26(9):094104, 09 2014.
  * [38] Scott T. M. Dawson, Maziar S. Hemati, Matthew O. Williams, and Clarence W. Rowley.  Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition.  Experiments in Fluids., 57:42, 2016.
  * [39] Maziar S. Hemati, Clarence W. Rowley, Eric A. Deem, and Louis N. Cattafesta.  De-biasing the dynamic mode decomposition for applied Koopman spectral analysis of noisy datasets.  Comput. Fluid Dyn., 31:349–368, 2017.
  * [40] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In Proceedings of the 33rd International Conference on Neural Information Processing Systems., 32, 2019.
  * [41] I. Goodfellow, Y. Bengio, and A. Courville.  Deep Learning. MIT Press, 2016.
  * [42] A. V. Clemente, H. N. Castejón, and A. Chandra.  Efficient parallel methods for deep reinforcement learning.  arXiv:1705.04862., 2017.