# Near-wall approximations to speed up simulations for atmosphere boundary
layers in the presence of forests using lattice Boltzmann method on GPU

Xinyuan Shao  Department of Mechanics and Maritime Sciences, Division of Fluid
Dynamics, Chalmers University of Technology, 41296, Göteborg, Sweden  Marta
Camps Santasmasas  School of Mechanical, Aerospace and Civil Engineering, The
University of Manchester, Sackville Street, M1 3BB, United Kingdom  Xiao Xue
Corresponding author
His email<EMAIL_ADDRESS>Department of Mechanics and Maritime Sciences,
Division of Fluid Dynamics, Chalmers University of Technology, 41296,
Göteborg, Sweden  Jiqiang Niu  School of Mechanical Engineering, Southwest
Jiaotong University, Chengdu, 610031, Sichuan, China  State Key Laboratory of
Automotive Simulation and Control, Jilin University, Changchun, 130025, Jilin,
China  Lars Davidson  Department of Mechanics and Maritime Sciences, Division
of Fluid Dynamics, Chalmers University of Technology, 41296, Göteborg, Sweden
Alistair J. Revell  School of Mechanical, Aerospace and Civil Engineering, The
University of Manchester, Sackville Street, M1 3BB, United Kingdom  Hua-Dong
Yao  Department of Mechanics and Maritime Sciences, Division of Fluid
Dynamics, Chalmers University of Technology, 41296, Göteborg, Sweden

###### Abstract

Forests play an important role in influencing the wind resource in atmospheric
boundary layers and the fatigue life of wind turbines. Due to turbulence, a
difficulty in the simulation of the forest effects is that flow statistical
and fluctuating content should be accurately resolved using a turbulence-
resolved CFD method, which requires a large amount of computing time and
resources. In this paper, we demonstrate a fast but accurate simulation
platform that uses a lattice Boltzmann method with large eddy simulation on
Graphic Processing Units (GPU). The simulation tool is the open-source
program, GASCANS, developed at the University of Manchester. The simulation
platform is validated based on canonical wall-bounded turbulent flows. A
forest is modelled in the form of body forces injected near the wall. Since a
uniform cell size is applied throughout the computational domain, the averaged
first-layer cell height over the wall reaches to $\langle\Delta
y^{+}\rangle=165$. Simulation results agree well with previous experiments and
numerical data obtained from finite volume methods. We demonstrate that good
results are possible without the use of a wall-function, since the forest
forces overwhelm wall friction. This is shown to hold as long as the forest
region is resolved with several cells. In addition to the GPU speedup, the
approximations also significantly benefit the computation efficiency. (This
paper is currently submitted to peer-review journal)

Wind energy, Atmospheric boundary layers, Forest modelling, Lattice Boltzmann
method, Large-eddy simulation, GPU

##  I Introduction

When wind farms are located within forest terrain, the presence of trees are
known to generate strong atmospheric turbulence, which can in turn lead to an
adverse impact on the performance and fatigue life of wind turbines [1]. In
recent years, there has been increased focus on the siting of wind turbines
within forests. This is in part motivated by the opportunity to reduce the
visual and environmental noise footprint of turbines, since the forest canopy
and background noise provide a natural screen. It is also driven by need,
since around 30% of the Earth’s land surface is covered by forest. It is thus
increasingly relevant to find an effective method to simulate wind field under
the influence of forests to help inform the siting at the early design phase
of a wind farm.

The effect of forests may be considered as a kind of surface roughness which
is modelled by introducing a roughness length[2]. However, this method is only
applicable to short vegetation. The permeability of the foliage should also be
taken into account for relatively high forests. Shaw and Schumann [3]
presented the effect of forests in their simulations as a drag force applied
by the trees on the flow and the drag force varies according to the varying
leaf area density (LAD). Their method has been accepted by many researchers
and has been used to investigate many aspects of canopy flows, for example,
coherent structures in canopy flow [4, 5, 6], the influence of LAD on canopy
flow [7] and edge flow due to heterogeneity in forest structure[8].

Large-eddy simulation (LES) has proven to be an useful tool to unravel the
engineering problems such as wind farms flow prediction[9, 10, 11]. The
atmospheric boundary layer (ABL) can be simulated using conventional finite
difference or finite volume methods solvers, for instance, PALM [12], ICON
[13], and MicroHH [14]. There are many studies which have endeavored to
combine LES with forest models. Shaw and Schumann [3] successfully represented
a forest as a drag layer and heat sources in their LES. Shaw and Patton [15]
further refined the forest model by including skin friction effect under the
LES framework.

LES is commonly adopted in the Navier-Stokes (NS) equation framework. However,
studies have been conducted to involve LES into the framework of lattice
Boltzmann models (LBM)[16]. In the past decades, LBM has become more
widespread on account of its broad applicability across a range of complex
fluid dynamics problems ranging from micro-nano scales[17, 18, 19] to
macroscopic scales[20, 21] under low Mach numbers. Moreover, its intrinsic
locality, whereby only information from directly adjacent points is required,
renders the algorithm efficient for parallelization [22]. LBM is thus seen as
a good alternative to more common methods based on NS solvers. A recent study
combined LBM and LES to simulate ABL flows [23]. They simulated neutral,
stable and convective ABLs over a forest in reference to the cases in [24],
where LES implemented with a finite volume method was used. Their simulations
included the forest effects using a forest model. Furthermore, since the mesh
resolution near the ground wall was not sufficient to resolve small-scale
turbulent fluctuations, the effect was recovered using the Monin-Obukhov wall
model. Their simulations results showed good agreement with the reference data
in [24].

The present study aims to explore the efficiency of an LBM solver based on a
GPU architecture, applied to the simulation of ABL flows with forests. Like
other works in the literature, LES is used to ensure the accuracy. Various
mesh refinement techniques have been proposed for LBM on GPUs, but they suffer
from reduced computational efficiency. Another concern is that as lattice
points remain in a ‘cubic’ structure, the mesh refinement is generally limited
to octree based methods with a scale of 1-2 in each direction. It means that
high aspect ratio cells, which are often used in wall-bounded turbulent flows,
are generally not possible. To avoid loss of the computational efficiency,
meshes with uniform cells can be adopted for ABL simulations on GPUs. In
principle, however, a wall function should be used because the near-wall mesh
resolution is not sufficient to resolve small scale flow fluctuations. When
the forest effects near ground walls are modelled, questions arise about
whether the forest drag force is dominant and whether the wall function can be
neglected. These are the motivations for the present study.

##  II Methodology

###  II.1 The lattice Boltzmann method

The flow is described based on the mesoscale lattice Boltzmann equations. The
classical Bhatnagar–Gross–Kroog (BGK) LBM collision operator and D3Q19 model
are adopted for the LBM. This model is a three-dimensional lattice model with
19 discretized velocity directions $\mathbf{c}_{i}$ ($i=0...Q-1$). We consider
$f_{i}(\mathbf{x},t)$ as the discretized particle’s probability distribution
function on the $i$-th direction of a lattice cell. The lattice cell is
located at position $\mathbf{x}$ at time $t$. The LBM governing equation for
the fluid, considering collision and forcing, can be written as:

$f_{i}(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta
t)=f_{i}(\mathbf{x},t)+\Omega\left[f_{i}(\mathbf{x},t)-f_{i}^{\mbox{\tiny
eq}}(\mathbf{x},t)\right]+\Delta tF_{i}(\mathbf{x},t),\\\ $ (1)

where $\Omega$ is a collision kernel [25, 26], $\Delta t$ is the marching time
step, $F_{i}$ is the volume force acting on the fluid following the approach
from Guo [27].

The collision kernel is the classical BGK collision kernel which has been
widely adopted to various applications [25]. The BGK collision operator fixes
the single relaxation time $\tau$ for the colliding process. Thus, the
collision operator can be represented as:

$\Omega\left[f_{i}(\mathbf{x},t)-f_{i}^{\mbox{\tiny
eq}}(\mathbf{x},t)\right]=\frac{\Delta
t}{\tau}\left[f_{i}(\mathbf{x},t)-f_{i}^{\mbox{\tiny
eq}}(\mathbf{x},t)\right]$ (2)

The collision kernel relaxes the distribution function towards the local
equilibrium $f_{i}^{\mbox{\tiny eq}}$:

$f_{i}^{\mbox{\tiny
eq}}\left(\textbf{x},t\right)=\omega_{i}\rho\left(\textbf{x},t\right)\bigg{[}1+\frac{\textbf{c}_{i}\cdot\textbf{u}\left(\textbf{x},t\right)}{c_{s}^{2}}+\frac{\left[\textbf{c}_{i}\cdot\textbf{u}\left(\textbf{x},t\right)\right]^{2}}{2c_{s}^{4}}-\frac{\left[\textbf{u}\left(\textbf{x},t\right)\cdot\textbf{u}\left(\textbf{x},t\right)\right]}{2c_{s}^{2}}\bigg{]}\;,$
(3)

where $\omega_{i}$ is a suitable weight needed to impose the isotropy in the
interaction, $\rho\left(\textbf{x},t\right)$ and
$\textbf{u}\left(\textbf{x},t\right)$ are the macroscopic hydrodynamic
quantities for density and velocity, respectively. The volume force
$F_{i}(\mathbf{x},t)$ in Eq. 1 can be obtained by:

$F_{i}(\mathbf{x},t)=(1-\frac{1}{2\tau})\omega_{i}\left[\frac{\textbf{c}_{i}-\textbf{u}\left(\textbf{x},t\right)}{c_{s}^{2}}+\frac{\textbf{c}_{i}\cdot\textbf{u}\left(\textbf{x},t\right)}{c_{s}^{4}}\textbf{c}_{i}\right]\mathbf{F},$
(4)

where $\mathbf{F}$ is the volume acceleration.

The macro-scale quantities for the density, momentum, and momentum flux
tensors can be calculated from the distribution function, the discrete
velocity, and the volume force:

$\rho(\mathbf{x},t)=\sum_{i=0}^{Q-1}f_{i}(\mathbf{x},t),\\\ $ (5)

$\rho(\mathbf{x},t)\mathbf{u}(\mathbf{x},t)=\sum_{i=0}^{Q-1}f_{i}(\mathbf{x},t)\mathbf{c}_{i}+\frac{1}{2}\mathbf{F}\Delta
t,$ (6)

$\mathbf{\Pi}(\mathbf{x},t)=\sum_{i=0}^{Q-1}f_{i}(\mathbf{x},t)\mathbf{c}_{i}\mathbf{c}_{i},$
(7)

where the momentum flux, $\mathbf{\Pi}$, can be presented by the sum of the
equilibrium and non-equilibrium parts,
$\mathbf{\Pi}(\mathbf{x},t)=\mathbf{\Pi}_{\text{eq}}(\mathbf{x},t)+\mathbf{\Pi}_{\text{neq}}(\mathbf{x},t)$.

###  II.2 Smagorinsky subgrid-scale model

The basic idea of LES is not only to resolve relatively large scales, but also
to use a subgrid-scale (SGS) model to estimate unresolved small scales. The
detailed implementation of the Smagorinsky SGS model can be found in [16, 28].
Here, we briefly summarize the essential part of this model. The key of the
LES modeling of Smagorinsky is to model the effective viscosity
$\nu_{\text{total}}$, which can be seen as the sum of the molecular viscosity
$\nu_{0}$ and the turbulent viscosity $\nu_{t}$:

$\nu_{\text{total}}=\nu_{0}+\nu_{t},$ (8)

where $\nu_{t}$ is represented by:

$\nu_{t}=C_{smag}\Delta^{2}\left|\bar{\mathbf{S}}\right|,$ (9)

where smag is the Smagorinsky constant, $\Delta$ represents the filter size,
and $\left|\mathbf{\bar{S}}\right|$ is the filtered strain rate tensor, which
can be obtained by:

$\left|\bar{\mathbf{S}}\right|=\frac{-\tau_{0}\rho c\Delta
x+\sqrt{(\tau_{0}\rho c\Delta x)^{2}+18\sqrt{2}\rho
C_{smag}\delta^{2}Q^{1/2}}}{6\rho C\delta^{2}},$ (10)

where $\tau_{0}$ is the original relaxation time from the input, $c=\Delta
x/\Delta t$, and $Q^{1/2}$ can be written as:

$Q^{1/2}=\sqrt{\mathbf{\Pi}_{\text{neq}}:\mathbf{\Pi}_{\text{neq}}},$ (11)

where $\mathbf{\Pi}_{\text{neq}}$ is the non-equilibrium part of the momentum
flux tensor $\mathbf{\Pi}$ shown in Eq. 7. Following [28], we can obtain the
total relaxation time $\tau_{\text{total}}$, which is written as:

$\tau_{\text{total}}=\frac{\tau_{0}}{2}+\frac{\sqrt{(\tau_{0}\rho
c)^{2}+18\sqrt{2}C_{smag}Q^{1/2}}}{2\rho c}.$ (12)

Finally, we replace $\tau$ in Eq. 2 with $\tau_{\text{total}}$ and obtain the
Smagorinsky SGS collision kernel.

###  II.3 Forest model near the wall

We model the effects of the forest by introducing an additional drag force
near the wall. The drag force is modelled with the help of the forest drag
coefficient $C_{D}$, the LAD $a_{f}$ and the local wind velocity
$\bm{u}=(u,v,w)$ [3]. Thus, the drag force of the forest, which is added in
the streamwise direction, can be written as:

$F_{f}(x,y,z,t)=-C_{D}a_{f}(y)\left|\bm{u}(x,y,z,t)\right|u(x,y,z,t),$ (13)

where $y$ refers to the vertical distance to the wall and $u(x,y,z)$ is the
streamwise velocity. The forest drag coefficient $C_{D}$ is usually taken
between $0.15$ and $0.2$ [7]. In this study, $C_{D}=0.15$.

According to [30], the LAD distribution $a_{f}(y)$ can be described by the
function:

$a_{f}(y)={a_{f}}_{m}\left(\frac{h-y_{m}}{h-y}\right)^{n}exp{\left[n\left(1-\frac{h-y_{m}}{h-y}\right)\right]}$
(14)

where

$n=\left\\{\begin{array}[]{ll}6&0\leq y<y_{m}\\\ \frac{1}{2}&y_{m}\leq y\leq
h\end{array}\right.$ (15)

The parameter $h$ is the total height of the forest, and $y_{m}$ is the
location where $a_{f}(y)$ gets its maximum value ${a_{f}}_{m}$.

##  III Numerical Implementation

The simulations presented in this paper are performed using the open source
CFD software GASCANS, which implements LBM on multiple graphic processing
units (GPUs) using CUDA and C++. The main aim of GASCANS is to provide fast
and accurate simulation of turbulent flows over complex geometries. To further
increase its usability, it incorporates a synthetic eddy method (SEM) to
generate turbulent velocity fluctuations from mean flow data and GASCANS can
be coupled at run time with another CFD software [31]. It also implements
immersed boundary method (IBM) and a structural solver for the simulation of
filaments submerged in fluid.

The flow solver implements the BGK collision scheme, the LES Smagorinsky
turbulence model and forcing scheme described in section II. The algorithm
follows a merged stream-collision time loop and stores the particle
distribution functions using a double-population scheme [32]. The simulation
data is stored in a mesh formed by equally sized cubic cells using a structure
of arrays (SOA) pattern.

Fig. 1 shows the GASCANS architecture and interactions between its components;
the features not used for this paper are greyed out. The code is divided into
two sections: the Application contains the user interface and inputs, while
the Library contains the core software. This structure allows for the
Application to be reworked or substituted to suit the user’s needs without
modifying the core of GASCANS. The user inputs are through the Parameters
class, and the terrain and forest geometry is read from point cloud files. See
[31] for a detailed explanation of the GASCANS functionality and structure.

Figure 1: Illustration of GASCANS architecture and relationship between its
components. The solid arrows indicate aggregation (i e. the classes at the
beginning of the arrow are part of the class that the arrow points to) and the
dashed arrows indicate primary communication. The isolated boxes represent
definitions and utilities accessible by the core LB solver.

The point cloud file contains a list of coordinates covering the forest’s
volume. GASCANS reads the file and assigns a code to each cell that contains
one or more forest coordinates. The code for each forest cell contains the
following information:

  * •

$h$: height of the forest, calculated as the bounding box of the forest object
in the vertical direction.

  * •

$y$: distance in the vertical direction between the lower bounding box of the
forest object and the current cell.

  * •

$y_{m}$: distance in the vertical direction between the lower bounding box of
the forest and the maximum value of the LAD distribution.

This information is stored into the forest cell’s code using a bit mask to
reduce the GPU memory required to store the information. Algorithm 1 shows the
forest implementation in GASCANS.

Algorithm 1 Implementation of forest forcing with varying leaf area density in
GASCANS.

Read the user input data: $C_{D}$, $a_{f_{m}}$, $y_{m}$ in the vertical
direction.

Read the forest geometry and generate the forest code for each forest cell.

Set the body force for all cells to 0.

for time = 0 to time = last time do

for all cells do

Stream

if cell $\mathbf{x}$ is forest then

Calculate local cell velocity $\mathbf{u}(\mathbf{x},t)$

Calculate forest force Eq. 13

Correct the local cell velocity using the forest force

end if

Collide using the calculated body force Eq. 1

end for

end for

##  IV Validation of simulation platform

Channel flow is a typical case of wall-bounded turbulence. The accuracy of the
simulation platform, which uses the LES with the Smagorinsky SGS model, is
validated based on the database of direct numerical simulation (DNS) at
$Re_{\tau}=180$ in [33].

The computational domain is shown in Fig. 2. The domain dimensions are non-
dimensionalized based on the characteristic length $H$. That is, the
dimensions are $L_{x}\times L_{y}\times L_{z}=12\times 2\times 4$ in the
streamwise, vertical and spanwise directions. The mesh consists of cubic
lattice cells. There are $31$ cells per length of $H$. Correspondingly, the
mesh contains $372\times 124\times 62$ lattice cells.

Figure 2: The computational domain of the turbulent channel flow. The domain
dimensions, normalized with the length scale $H$, are $L_{x}\times L_{y}\times
L_{z}=12\times 2\times 4$ along the $x$, $y$ and $z$ axes. The walls are
marked in gray, and the arrows indicate the flow direction.

The upper and lower walls are set with the no-slip boundary condition, and the
remaining boundaries of the computational domain with the periodic boundary
condition. The flow is driven by a volume body force from Eq. 4, which is
defined as:

$\mathbf{F}=[\frac{\rho u_{\tau}^{2}}{H},0,0].$ (16)

The body force acts equivalently as the pressure gradient that balances the
wall shear stress in the fully developed channel flow.

The parameters set for the simulation are listed in Table 1. Here $C_{smag}$
is the Smagorinsky model constant, as shown in Eq. 9. The dimensionless
density $\rho$ is set to unity in the lattice Boltzmann unit (LBU). The fluid
is incompressible.

Table 1: The simulation parameters for the turbulent channel flow Parameter | Value | Parameter | Value
---|---|---|---
$\Delta x$ [m] | $1/31$ |  $L_{x}$ [m] | 12
$\Delta t$ [s] | $0.001$ |  $L_{y}$ [m] | 2
$\nu$ [$m^{2}$/s] | $1/3250$ |  $L_{z}$ [m] | 4
$F$ [N] | 0.002986 | $C_{smag}$ | 0.01

Figure 3 shows the instantaneous streamwise velocity distribution and the
Q-criterion of $Q=0.0002$. A plane is placed at the vertical position of
$y^{+}=9$ to visualize streamwise velocity contours. Streaks featuring the
boundary layer are observed. This is consistent with the previous DNS [33].
Plentiful turbulent structures are identified within the flow field based on
the Q-criterion. Hairpin vortices exist near the walls, leading to streak
contours of the streamwise velocity.

Figure 3: (a) The instantaneous non-dimensional streamwise velocity,
$u/u_{\tau}$, of the turbulent channel flow at $Re_{\tau}=180$. The bottom
plane shown here is located at $y^{+}=9$. (b) Isosurfaces of the Q-criteron of
$Q=0.0002$, which are colored with the streamwise velocity.

As shown in Fig. 4, the statistics of the flow quantities from the present LBM
LES method are compared to the data computed using the LBM LES reported by
Koda et al. [28] and the finite volume method of DNS by Moser et al. [33]. The
non-dimensional time-averaged streamwise velocity is defined as
$\langle{u}^{+}\rangle=\left<u\right>/u_{\tau}$. Here $\left<\cdot\right>$ is
the averaging operator, and $u_{\tau}$ denotes the friction velocity. The
averaging is done in time and then in x-z planes that are parallel to the top
and bottom walls. The non-dimensional room mean square (RMS) velocity
components ($u_{rms}^{+}$, $v_{rms}^{+}$ and $w_{rms}^{+}$) are also
normalized with the friction velocity $u_{\tau}$. The present results are
consistent with Koda’s LBM LES data, especially in the near wall region. But
these results exhibit discrepancies in comparison with the Moser’s DNS data.
Both LBM LES methods underestimate $\langle{u}^{+}\rangle$. The reason is that
the mass flow rate computed with the coarse meshes in LES is lower than that
from DNS. The RMS velocity components from LES and DNS are mainly different
below $y^{+}=100$. It suggests that the LES method does not completely
reproduce the subgrid-scale flow structures near the wall. Nonetheless, the
largest discrepancy of $\langle{u}^{+}\rangle$ is less than $5\%$, and the
largest RMS velocity difference less than $10\%$. Therefore, the present
simulation platform, which is developed based on GASCANS, show overall correct
prediction in comparison with DNS for wall-bounded turbulence.

Figure 4: The comparison of the present results and the DNS data computed
using LBM of LES [28] and a finite volume method of DNS [33]. (a) The non-
dimensional averaged streamwise velocity $\langle{u}^{+}\rangle$ as a function
of $y^{+}$. (b-d) The non-dimensional RMS velocity components in the
streamwise, wall normal and spanwise directions.

##  V Atmospheric boundary layers over forest

###  V.1 Case description

A classical ABL case reported in [34, 24] is studied using the present
simulation platform. The computational domain is shown in Fig. 5. The physical
dimensions are $L_{x}\times L_{y}\times L_{z}=1000\times 400\times
800\;m^{3}$. To produce the representative forest effects, a forest with a
homogeneous distribution over the ground is assumed [34]. The height of the
forest is $h=20\;m$. The freestream flow direction is set along the $x$ axis.
The forest friction velocity is defined as
$u_{*}=\left[\left<u^{\prime}v^{\prime}\right>^{2}+\left<v^{\prime}w^{\prime}\right>^{2}\right]^{1/4}$
at the location of $y/h=2$ [34, 24].

Figure 5: The computational domain of the ABL. The green color marks out the
zone where a homogeneous forest with a height of $h=20\;m$ is modelled above
the ground. The freestream flow moves along the x-axis.

The flow is driven by a constant non-dimensional body force of $2.38\times
10^{3}$ per unit volume on the basis of $u_{*}$ and $L_{x}$. The force is
imposed in the whole domain. At the bottom wall below the forest, a no-slip
boundary condition is applied. The inlet and outlet are set with the periodic
boundary condition, and the spanwise side boundaries are also set with the
same boundary condition. The upper boundary is specified with a forced
equilibrium boundary condition. This boundary condition specifies the upper
wall to move with the same velocity as the free-stream flow. Because the upper
wall is sufficiently far away from the bottom, and because only the boundary
layer region near the bottom wall where wind turbines are embedded are of
interest for the wind energy harvesting, this forced equilibrium boundary
condition is suitable for the present ABL case.

In reality, the LAD $a_{f}$ varies with the vertical height of a tree, as the
tree crown has dense leaves and the trunk is less obstructing. In the current
work, we investigate three distributions of LAD, as listed in Table
LABEL:tab:forest_cases. Case 1 has constant LAD in the vertical direction,
indicating that the leaf and obstructions are uniform along a tree. Case 2 has
varying LAD, which is more similar to the real situation [34, 24]. Case 3 has
no forest force imposed (i.e., constant LAD of $a_{f}=0$). However, the LAD
distributions of cases 1 and 2 have the same total obstruction, which is
defined as $A_{f}=\int_{0}^{h}a_{f}(z)dz$, whereas case 3 has no obstruction.
The LAD distributions of cases 1 and 2 are illustrated in Fig. 6. In case 1,
$a_{f}=0.215$. In the varying LAD case, $a_{f}$ is set on the basis of Eq.(14)
with the coefficients $z_{m}=0.7h$ and ${a_{f}}_{m}=0.38$.

Table 2: The simulated ABL cases in terms of the LAD. | Case 1 | Case 2 | Case 3
---|---|---|---
$a_{f}(y)$ | $0.215$ | Eq. 14 ($z_{m}=0.7h$, ${a_{f}}_{m}=0.38$) |  $0.0$ (i.e., no forest modelling)

To develop the turbulence, a set of small boxes are set near the inlet in the
initial stage of the simulations. The non-dimensional box size is $20\times
20\times 20$. Once the turbulence is developed, the boxes are removed from the
computational domain. The flow remains in the turbulent status afterwards.

Figure 6: The distributions of LAD in two situations such as constant LAD of
$a_{f}=0.215$, indicated by the dash line, and varying LAD with $a_{f}$
following Eq.(14) with $z_{m}=0.7h$ and ${a_{f}}_{m}=0.38$, indicated by the
solid curve. The green colors the area of $a_{f}=0.215$, which is equal to the
area of the varying LAD case colored in red.

###  V.2 Mesh independence study

The mesh quality for the simulation accuracy is examined with the meshes with
coarse, medium and fine cell size. The constant LAD of $a_{f}=0.215$ is
adopted for the forest modelling. The information of the meshes are presented
in Table LABEL:tab:mesh_cost. The coarse mesh has the number of cells in the
three dimensions as $N_{x}\times N_{y}\times N_{z}=200\times 80\times 160$,
resulting in $2.6\times 10^{6}$ cells in total. We define the refinement ratio
as the ratio of a refined cell size to the reference cell size. In reference
to the coarse mesh, the refinement ratios of the medium and fine meshes are
$0.8$ and $1/1.75\approx 0.57$.

The computational time and memory costs relevant to the meshes are also
reported in Table LABEL:tab:mesh_cost. The GPU, GeForce RTX 2070 SUPER, with a
total memory of $8$ GB is used for the computations. The coarse mesh
simulation consumes the computing time of approximately $8.5$ hours per
$1\times 10^{6}$ time steps and the memory of nearly $1.44$ GB. For the medium
mesh simulation, the computing time is $1.65$ times larger than the coarse
mesh, and the memory cost is $1.5$ times larger. The fine mesh simulation is
$4.59$ times faster and $4.5$ times more memory consumption. The increasing
ratios of the computing time are comparable with those of the memory costs. On
the other hand, the total number of cells in the medium and fine meshes
increase with ratios of $1.95$ and $5.36$, which are larger than the ratios of
the computing time and memory costs. It suggests that the increased mesh size
improves the computational efficiency. The reason is that there are modules
responsible for processing and storing solution data in the present simulation
platform. The additional computational and memory costs from these modules are
independent on the mesh size, so they are more weighted in the coarse mesh
case, as compared with the other refined mesh cases. Moreover, the GPU solves
the LBM equations on all cells at the same time. An increased number of cells
has less impact on the computing time if the former coarse mesh does not fill
all capacity of the GPU.

Table 3: The meshes and respective computational costs in the mesh
independence study.

Meshes |  $N_{x}\times N_{y}\times N_{z}$ § |  $N_{tot}$ ¶ |  $\Delta t$ † |  $T_{c}$ (hours) ‡ | Memory (GB)
---|---|---|---|---|---
Coarse | $200\times 80\times 160$ | $2.6\times 10^{6}$ | $0.3$ | $8.5$ | $1.44$
Medium | $250\times 100\times 200$ | $2.6\times 10^{6}\cdot 1.95$ | $0.2$ | $8.5\cdot 1.65$ | $1.44\cdot 1.5$
Fine | $350\times 140\times 280$ | $2.6\times 10^{6}\cdot 5.36$ | $0.1$ | $8.5\cdot 4.59$ | $1.44\cdot 4.5$

  * •

$\S$ The number of cells in the three dimensions.
$\P$ The total number of cells.
$\dagger$ The non-dimensional time interval of each time step.
$\ddagger$ The computing time per $1\times 10^{6}$ time steps.

The normalized averaged streamwise velocity along the wall normal direction
for the simulations of the three meshes is shown in Fig. 7. The averaging is
done in time and then in $x-z$ planes that are parallel to the bottom wall.
Hereafter, the averaged quantities reported in this study are obtained using
the same averaging method. As compared to the experimental data from Bergström
et al. [35], the fine mesh simulation shows the best agreement. In the range
of the experimental data, the coarse mesh simulation gives better prediction
than the medium mesh simulation. But these two simulations become similar
above $y/h=14$, which is far from the bottom wall and forest. The profile of
the coarse mesh simulation rises more slowly for $y/h>2$. It implies that
there are less turbulent vortices resolved by the coarse mesh.

Figure 7: The averaged streamwise velocity obtained from the present
simulations with the constant LAD and from the experiments by Bergström et al.
[35]

The Reynolds stress tensor components, $\left<u^{\prime}u^{\prime}\right>$,
$\left<v^{\prime}v^{\prime}\right>$, $\left<w^{\prime}w^{\prime}\right>$ and
$\left<u^{\prime}v^{\prime}\right>$, normalized with the averaged wall
friction velocity, $u_{*}$, are shown in Fig. 8. It is observed that the
coarse mesh leads to the worst prediction of all components in reference to
the experimental data from Bergström et al. [35], whereas the prediction of
the fine mesh is most consistent. Therefore, a general trend is that the mesh
refinement improves the simulation accuracy.

For the streamwise component $\left<u^{\prime}u^{\prime}\right>/u_{*}^{2}$ in
Fig. 8, the slope of the profile from the coarse mesh simulation is more
straight than the other meshes. This is a typical phenomenon when the
turbulent flow is less developed. As a consequence, velocity fluctuations
owing to turbulence vortices are insufficient to form a steeper profile. This
effect is also found for the wall normal component
$\left<v^{\prime}v^{\prime}\right>/v_{*}^{2}$ and the spanwise component
$\left<w^{\prime}w^{\prime}\right>/w_{*}^{2}$, but it is not so significant as
$\left<u^{\prime}u^{\prime}\right>/u_{*}^{2}$. Since streamwise velocity
fluctuations have larger magnitudes than spanwise and wall normal velocity
fluctuations, the streamwise velocity ones are more influential in the
formation of the Reynolds stress component
$\left<u^{\prime}v^{\prime}\right>$. As seen in subfigure 8d, the profiles of
$\left<u^{\prime}v^{\prime}\right>/u_{*}^{2}$ show similar trends to
$\left<u^{\prime}u^{\prime}\right>/u_{*}^{2}$.

Figure 8: The normalized Reynolds stress tensor components computed with the
three meshes and reported in the experiments by Bergström et al. [35] Here the
constant LAD is set in the forest modelling.

The fine mesh results in overall the best prediction, as indicated in Figs. 7
and 8. Based on Tabl LABEL:tab:mesh_cost, the averaged first-layer cell height
over the wall is $\left<\Delta y+\right>=165$. Moreover, the results of the
medium and fine meshes are nearly the same. This converged prediction means
that the fine mesh resolution is sufficient to simulate the current ABL case.
In addition, the computational time and memory costs of this mesh are lowest,
which is found in Table LABEL:tab:mesh_cost. The fine mesh is therefore
adopted in the following simulations and analysis.

Based on Tabl LABEL:tab:mesh_cost, the averaged first-layer cell height over
the wall is $\left<\Delta y+\right>=165$. This infers that the wall modelling
should be used in principle. However, the practice in the current simulations
show that good agreement is achieved without using a wall function. The reason
is that the forest forces are much larger than the wall friction. This effect
will be discussed in the subsequent section on the forest modelling.

###  V.3 Effects of forest modelling

To investigate the effects of the forest modelling, the simulation cases with
different LAD are defined in Table LABEL:tab:forest_cases. Case 1 (the
modelled forest with the constant LAD) has been discussed in the preceding
section, and it is found producing consistent prediction in comparison to the
experiments.

For case 2 (the modelled forest with the varying LAD) and case 3 (zero LAD,
i.e., no modelled forest), the instantaneous normalized streamwise velocity
$u/u_{*}$ are shown in Fig. 9. Case 2 exhibits generally much smaller velocity
fluctuations than case 3 in the whole computational domain. Fluctuations are
visualized at a $x-z$ plane of $y/h=0.5$, which is the center position of the
modelled forest zone. Velocity magnitudes in this plane in both cases are
smaller than those in the region far from the bottom wall. As case 2 includes
the drag forces from the forest modelling and from the wall friction, the
small fluctuations within the forest zone are explainable for this case. There
are no forest drag forces imposed in case 3. In other words, case 3 does not
have a forest zone. Its near-wall flow is only dragged by the wall friction.
Therefore, the forest dragging effect is dominant in reducing the velocity of
the ABL flow, as compared with the wall friction. It is observed in case 2
that the velocity quickly increases to larger values outside the forest zone.
This also qualitatively indicates that the significant dragging effect from
the modelled forest. The same finding was reported in [34, 24].

Figure 9: A snapshot of the instantaneous normalized streamwise velocity
$u/u_{*}$. (a) The simulation with the forest modelling of the varying LAD;
(b) the simulation with no modelled forest. The bottom surface visualized here
is located at $y/h=0.5$, which is the middle of the forest zone.

The vertical profiles of the averaged streamwise velocity
$\left<u\right>/u_{*}$ in all simulated cases listed in
LABEL:tab:forest_cases, as well as the FVM LES with the same varying LAD
setting as case 2 in [34, 24] and the experiment data in [35], are plotted in
Fig. 10. In the near ground zone of $y/h<6$ where large part of the boundary
layer is located, cases 1 and 2 exhibit good agreement with the experimental
data. The FVM LES also agrees with the experiment, while discrepancies are
observed for $y/h$ above approximately $5$. The The mesh for the FVM LES
contains $10$ cells in the forest zone. It has been demonstrated that the mesh
resolution is sufficient to resolve the profile of the varying LAD illustrated
in Fig. 6 [34, 24]. There are $7$ cells in the fine mesh for the current
simulations, which also resolves the same varying LAD profile. The comparison
of cases 2 and 3 indicates that the constant and varying LADs lead to slight
differences in the prediction of the mean velocity in the forest zone and most
of the zone above the forest.

In case 3 that has no forest, the boundary layer development is completely
different from cases 1 and 2 considering the forest effects. The velocity
gradient in the forest zone in case 3 is much larger than the other cases. A
similar phenomenon was also reported in [34, 24]. The reason is that the
forest model introduces significant drag forces into the flow. Based on Eq.
13, the forest drag force is independent from the wall friction. Its vertical
distribution is not governed by the classical boundary layer statistics, which
can be resumed using a wall function. Although the wall function is ignored,
only considering the forest model is still able to predict the mean velocity
profile matching the experiments. This is explainable because the forest force
overwhelms wall friction.

Figure 10: The averaged streamwise velocity in the present simulations with
different forest modelling settings, as well as from the simulations using the
FVM LES by Nebenführ et al. [34, 24] and the experiments by Bergström et al.
[35]

The vertical profiles of the normalized stress tensor components are displayed
in Fig. 11. In the simulations set with the forest model (cases 1 and 2, as
well as the FVM LES in [34, 24]), the predicted trends of
$\left<u^{\prime}u^{\prime}\right>$ are consistent. But a mismatch between
these methods and the experimental data is observed in the near forest zone of
$2<y/h<4$. Cases 1 and 2 also show $\left<v^{\prime}v^{\prime}\right>$ and
$\left<u^{\prime}v^{\prime}\right>$ comparable with the experimental data. The
different LAD settings of cases 1 and 2 lead to discrepancies between their
results. But the results are approximately identical at the first cell height.
Recalling Eq. 13, the forest drag force modelled is proportional to the local
kinetic energy. Since the velocity at the wall is zero based on the no-slip
boundary condition, the local forest drag force diminishes there. The small
velocity near the wall leads to the negligible force, although the constant
LAD is much larger than the varying LAD near the wall (see Fig. 6). As the
majority difference of the two LAD settings are defined near the wall for
$y/h$ less than $0.4$, their influence on the drag force production is limited
owing to the small local velocity. As a consequence, the profiles of the
tensor statistics simulated with the two LAD settings are similar.
Nevertheless, the present study only demonstrates the forest effect modelled
with the LADs shown Fig. 6. If the LAD definition is changed, the forest
effects could change and, thus, should be thoroughly examined.

As the forest drag force is not modelled, case 3 gives different prediction in
comparison with cases 1 and 2, in particular, for
$\left<u^{\prime}u^{\prime}\right>/u_{*}^{2}$ displayed in Fig. 11. Obvious
differences are identified for all stress tensor components within the forest
zone, and $\left<u^{\prime}u^{\prime}\right>/u_{*}^{2}$ also shows
inconsistency at vertical locations far from the forest zone. Therefore, the
mismatch conveys that the modelled forest drag force play a critical role in
predicting the statistics of flow fluctuations in the boundary layer, in
addition to its importance for the mean flow statistics that is shown in Fig.
10.

Figure 11: The stress tensor components in the present simulations with
different forest modelling settings, as well as from the simulations using the
FVM LES by Nebenführ et al. [34, 24] and the experiments by Bergström et al.
[35].

##  VI Conclusions

A GPU-based simulation platform using an LBM method for LES has been proven to
be suitable for computing ABL flows with the forest presence, which are of
interest in the site selection of wind turbine farms. The forest effects are
modelled by defining a volume drag force in relation to the local velocity and
LAD. The way of implementing the force for LBM is presented in detail.

To refrain from influencing the GPU computational efficiency, the meshes are
not locally refined. Uniform-sized lattices are generated in the whole
computational domain. This leads to large $\Delta y+$ near the wall, which in
principle requires a wall function. However, we introduce the approximations
near the ground wall that the wall function is not taken into account, and
that the forest effects are considered. These approximations are made because
the forest drag force overwhelms the the wall friction. This is validated
based on the qualitative and quantitative analysis of the results. The
statistics of the mean flow velocity and stress tensor components are well
consistent with the previous FVM-LES simulations and experiments. The mesh
independence study illustrates that the modelled forest zone is resolved with
$7$ cells.

One advantage of the near-wall approximations is that the LBM program is
compatible with the GPU architecture, so the GPU speedup is preserved.
Moreover, the wall treatment is ignored to reduce the computational cost. The
platform extends the capabilities of the open-source simulation tool GASCANS
developed at the University of Manchester.

The forest effects are examined for a classic vertically varying LAD
distribution, which was simulated in previous studies using the FVM-LES, and a
constant distribution possessing the same total area as the varying one. The
results between these two LADs differ only slightly. Nevertheless, it would be
interesting to investigate the effects of the LAD vertical profile on the
characteristics of ABL flows. In the present study, the standard Smagorinsky
model is used to model the SGS turbulence. A future work is to implement other
SGS models.

##  VII ACKNOWLEDGEMENTS

The authors would like to acknowledge the funding from Chalmers Area of
Advance Transport. The computations and data handling were enabled by
resources provided by the Swedish National Infrastructure for Computing
(SNIC), partially funded by the Swedish Research Council through grant
agreement no. 2018-05973.

## References

  * Nebenführ and Davidson [2017] B. Nebenführ and L. Davidson, Prediction of wind-turbine fatigue loads in forest regions based on turbulent les inflow fields, Wind Energy 20, 1003 (2017).
  * Lo [1990] A. K. Lo, On the determination of zero-plane displacement and roughness length for flow over forest canopies, Boundary-layer meteorology 51, 255 (1990).
  * Shaw and Schumann [1992] R. H. Shaw and U. Schumann, Large-eddy simulation of turbulent flow above and within a forest, Boundary-Layer Meteorology 61, 47 (1992).
  * Nezu and Sanjou [2008] I. Nezu and M. Sanjou, Turburence structure and coherent motion in vegetated canopy open-channel flows, Journal of hydro-environment research 2, 62 (2008).
  * Dupont and Brunet [2009] S. Dupont and Y. Brunet, Coherent structures in canopy edge flow: a large-eddy simulation study, Journal of Fluid Mechanics 630, 93 (2009).
  * Finnigan  _et al._ [2009] J. J. Finnigan, R. H. Shaw, and E. G. Patton, Turbulence structure above a vegetation canopy, Journal of Fluid Mechanics 637, 387 (2009).
  * Dupont and Brunet [2008a] S. Dupont and Y. Brunet, Influence of foliar density profile on canopy flow: a large-eddy simulation study, Agricultural and forest meteorology 148, 976 (2008a).
  * Dupont and Brunet [2008b] S. Dupont and Y. Brunet, Edge flow and canopy structure: a large-eddy simulation study, Boundary-Layer Meteorology 126, 51 (2008b).
  * Stevens  _et al._ [2014] R. J. Stevens, D. F. Gayme, and C. Meneveau, Large eddy simulation studies of the effects of alignment and wind farm length, Journal of Renewable and Sustainable Energy 6, 023105 (2014).
  * Stevens  _et al._ [2015] R. J. Stevens, D. F. Gayme, and C. Meneveau, Coupled wake boundary layer model of wind-farms, Journal of renewable and sustainable energy 7, 023115 (2015).
  * Stevens and Meneveau [2017] R. J. Stevens and C. Meneveau, Flow structure and turbulence in wind farms, Annual review of fluid mechanics 49, 311 (2017).
  * Maronga  _et al._ [2015] B. Maronga, M. Gryschka, R. Heinze, F. Hoffmann, F. Kanani-Sühring, M. Keck, K. Ketelsen, M. O. Letzel, M. Sühring, and S. Raasch, The parallelized large-eddy simulation model (palm) version 4.0 for atmospheric and oceanic flows: model formulation, recent developments, and future perspectives, Geoscientific Model Development 8, 2515 (2015).
  * Dipankar  _et al._ [2015] A. Dipankar, B. Stevens, R. Heinze, C. Moseley, G. Zängl, M. Giorgetta, and S. Brdar, Large eddy simulation using the general circulation model icon, Journal of Advances in Modeling Earth Systems 7, 963 (2015).
  * Van Heerwaarden  _et al._ [2017] C. C. Van Heerwaarden, B. J. Van Stratum, T. Heus, J. A. Gibbs, E. Fedorovich, and J. P. Mellado, Microhh 1.0: A computational fluid dynamics code for direct numerical simulation and large-eddy simulation of atmospheric boundary layer flows, Geoscientific Model Development 10, 3145 (2017).
  * Shaw and Patton [2003] R. H. Shaw and E. G. Patton, Canopy element influences on resolved-and subgrid-scale energy within a large-eddy simulation, Agricultural and forest meteorology 115, 5 (2003).
  * Hou  _et al._ [1994] S. Hou, J. Sterling, S. Chen, and G. Doolen, A lattice boltzmann subgrid model for high reynolds number flows, arXiv preprint comp-gas/9401004 (1994).
  * Belardinelli  _et al._ [2015] D. Belardinelli, M. Sbragaglia, L. Biferale, M. Gross, and F. Varnik, Fluctuating multicomponent lattice boltzmann model, Phys. Rev. E 91, 023313 (2015).
  * Xue  _et al._ [2018] X. Xue, M. Sbragaglia, L. Biferale, and F. Toschi, Effects of thermal fluctuations in the fragmentation of a nanoligament, Physical Review E 98, 012802 (2018).
  * Xue  _et al._ [2020] X. Xue, L. Biferale, M. Sbragaglia, and F. Toschi, A lattice boltzmann study on brownian diffusion and friction of a particle in a confined multicomponent fluid, Journal of Computational Science 47, 101113 (2020).
  * Filippova  _et al._ [2001] O. Filippova, S. Succi, F. Mazzocco, C. Arrighetti, G. Bella, and D. Hänel, Multiscale Lattice Boltzmann Schemes with Turbulence Modeling, Journal of Computational Physics 170, 812 (2001).
  * Toschi and Bodenschatz [2009] F. Toschi and E. Bodenschatz, Lagrangian properties of particles in turbulence, Annual review of fluid mechanics 41, 375 (2009).
  * Krüger  _et al._ [2017] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen,  _The Lattice Boltzmann Method_, Graduate Texts in Physics (Springer International Publishing, Cham, 2017).
  * Feng  _et al._ [2021] Y. Feng, J. Miranda-Fuentes, S. Guo, J. Jacob, and P. Sagaut, Prolb: A lattice boltzmann solver of large-eddy simulation for atmospheric boundary layer flows, Journal of Advances in Modeling Earth Systems 13, e2020MS002107 (2021).
  * Nebenführ and Davidson [2015] B. Nebenführ and L. Davidson, Large-eddy simulation study of thermally stratified canopy flow, Boundary-layer meteorology 156, 253 (2015).
  * Succi [2001] S. Succi,  _The lattice Boltzmann equation: for fluid dynamics and beyond_ (Oxford university press, 2001).
  * Krüger  _et al._ [2017] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, The lattice boltzmann method, Springer International Publishing 10, 978 (2017).
  * Guo  _et al._ [2002] Z. Guo, C. Zheng, and B. Shi, Discrete lattice effects on the forcing term in the lattice boltzmann method, Physical review E 65, 046308 (2002).
  * Koda and Lien [2014] Y. Koda and F.-S. Lien, The lattice boltzmann method implemented on the gpu to simulate the turbulent flow over a square cylinder confined in a channel, Flow, Turbulence and Combustion 94, 495 (2014).
  * [29] Large-eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms, Journal of Wind Engineering and Industrial Aerodynamics 99, 154 (2011).
  * Lalic and Mihailovic [2004] B. Lalic and D. T. Mihailovic, An empirical relation describing leaf-area density inside the forest for environmental modeling, Journal of Applied Meteorology 43, 641 (2004).
  * Santasmasas [2021] M. C. Santasmasas, Hybrid gpu / cpu navier-stokes lattice boltzmann method for urban wind flow (2021).
  * Latt  _et al._ [2021] J. Latt, C. Coreixas, and J. Beny, Cross-platform programming model for many-core lattice boltzmann simulations, PLoS One, Public Library of Science 16, 1 (2021).
  * Moser  _et al._ [1999] R. D. Moser, J. Kim, and N. N. Mansour, Direct numerical simulation of turbulent channel flow up to $re_{\tau}=590$, Physics of Fluids 11, 943 (1999).
  * Nebenführ and Davidson [2014] B. Nebenführ and L. Davidson, Influence of a forest canopy on the neutral atmospheric boudary layer-a les study, in  _ETMM10: 10th International ERCOFTAC Symposium on Turbulence Modelling and Measurements_ (2014).
  * Bergström  _et al._ [2013] H. Bergström, H. Alfredsson, J. Arnqvist, I. Carlén, E. Dellwik, J. Fransson, H. Ganander, M. Mohr, A. Segalini, and S. Söderberg, Wind power in forests: wind and effects on loads (2013).


# Towards Autonomous Eye Surgery by Combining Deep Imitation Learning with
Optimal Control

Ji Woong Kim  Department of Mechanical Engineering, Johns Hopkins University
Peiyao Zhang  Department of Mechanical Engineering, Johns Hopkins University
Peter Gehlbach  Wilmer Eye Institute, Johns Hopkins University School of
Medicine  Iulian Iordachita  Department of Mechanical Engineering, Johns
Hopkins University  Marin Kobilarov  Department of Mechanical Engineering,
Johns Hopkins University

###### Abstract

During retinal microsurgery, precise manipulation of the delicate retinal
tissue is required for positive surgical outcome. However, accurate
manipulation and navigation of surgical tools remain difficult due to a
constrained workspace and the top-down view during the surgery, which limits
the surgeon’s ability to estimate depth. To alleviate such difficulty, we
propose to automate the tool-navigation task by learning to predict relative
goal position on the retinal surface from the current tool-tip position. Given
an estimated target on the retina, we generate an optimal trajectory leading
to the predicted goal while imposing safety-related physical constraints aimed
to minimize tissue damage. As an extended task, we generate goal predictions
to various points across the retina to localize eye geometry and further
generate safe trajectories within the estimated confines. Through experiments
in both simulation and with several eye phantoms, we demonstrate that our
framework can permit navigation to various points on the retina within 0.089mm
and 0.118mm in xy error which is less than the human’s surgeon mean tremor at
the tool-tip of 0.180mm. All safety constraints were fulfilled and the
algorithm was robust to previously unseen eyes as well as unseen objects in
the scene. Live video demonstration is available here:
https://youtu.be/n5j5jCCelXk

> Keywords: deep imitation learning, optimal control, eye surgery,
> ophthalmology

Figure 1: (Left) the surgical tool is inserted into the human eye through a
sclera port and the light source projects a shadow which can be used as cues
to estimate depth (right) demonstration of tool-shadow dynamics

##  1 Introduction

Retinal surgery is among the most challenging of microsurgical procedures
requiring precise manipulation of the delicate retinal tissue. The surgery
proceeds with the surgeon inserting a surgical tool through a sclera port and
visualizing the retina and the surgical tools from a top-down view using an
operating microscope. During surgery, one of the most challenging tasks is
precisely navigating the surgical tool-tip to the desired position on the
retinal tissue. For example, while performing retinal vein cannulation or
retinal membrane peeling, the surgeon must intuitively localize the tool-tip
with respect to the retina in order to land the surgical tool-tip precisely at
the desired tissue location. Several factors contribute to its difficulty,
such as viewing the surgery from top-down view which can hinder the surgeon’s
depth-estimation accuracy, visual distortions caused by the human lens, and
the limited depth of field of the operating microscope. Precise navigation of
surgical tools is required for most retinal procedures, yet it remains among
the more challenging tasks, requiring years of experience to master.

Figure 2: Data collection, training and inference process for achieving
autonomous navigation setup

In this work, we propose to automate the tool navigation task by learning to
predict relative goal positions from video and by generating smooth
trajectories to these goals while satisfying known physical kinematic
constraints. The learned goal-prediction network estimates distance and
direction from the tool-tip to various target locations on the retinal
surface, in a manner similar to how surgeons utilize visual cues to safely
navigate the tool-tip with respect to the retina. Specifically, the input to
our proposed network is a top-down monocular view of the surgery and user-
input defining the 2-d goal location to be reached i.e. as simple as clicking
the desired location using a mouse Fig.2. Given these inputs, the output of
the network is a 3-d vector (direction and distance) to the clicked position.
The advantage of this approach is that the user specifies the goal in 2-d, and
the network outputs a prediction in 3-d space; since estimating depth is the
challenging task for humans, the network takes the burden by predicting
distance along the depth dimension based on its training experience. The
further advantage of this learning setup is that with goal predictions across
various points on the retina, the retinal geometry can be reconstructed and
employed for safe tool navigation. We include a demonstration of this
approach, and as an extended implementation, we demonstrate a blood-vessel
following task where the surgical tool-tip hovers over the localized retinal
geometry while staying as close as 0.200$\mu m$ to the retinal surface. This
enables safe planning of various trajectories inside the eye.

Another important consideration in achieving autonomous navigation is the
physical constraint posed by the sclera point, i.e. the hole through which the
surgical tool is inserted into the eye. The sclera point can be considered as
a point constraint which should not be violated in lateral directions; large
deviations form the sclera point could cause eye rotation and thus eye-muscle
damage. To address this issue, we integrate an optimal control framework that
enables path-planning of trajectories to the retina while respecting the
sclera constraint. Specifically, we utilize differential dynamic programming
(DDP) [1], which is a numerical optimization approach to path-planning based
on specified costs. To satisfy the point constraint of the sclera, we simply
add a cost to minimize the distance between the sclera point and the tool-axis
when searching for an optimal trajectory. In extended tasks such as blood-
vessel following, we utilize this framework to constrain the trajectories to
remain close to and without colliding with the retinal surface. Ultimately,
our system combines both a learning-based approach and optimal control
framework, where the goal-prediction network generates a goal waypoint on the
retina and the DDP algorithm plans an optimal trajectory to that goal. During
execution of the perceptual inference (including goal prediction and eye
geometry reconstruction), the optimal trajectory is repeatedly re-computed at
some fixed frequency, e.g. 5 Hz, until the tool-tip reaches its target.

We also note that our approach is based on the assumption that learning
autonomous navigation is possible primarily using vision. In fact, surgeons
rely on their visual perception to localize objects and spatially estimate
their proximity to perform retinal surgery. Furthermore, the tool-shadow
dynamics provide an important cue that helps with detecting proximity between
the tool-tip and the retina. Specifically, when the tool approaches the
retina, the tool and its shadow converge, which can be used as cues to train
the network (Fig.1).

To assess the accuracy of our system, we employ a benchmark task where the
surgical tool-tip must reach a predefined array of positions on the retina.
The objective of this experiment is to assess how well the network can
navigate to various positions on the retina given a specific goal image. As an
extension of this task, we also consider localization of the retina by
predicting distances to various points across the retinal surface and
reconstructing the eye-geometry. Furthermore, using the localized eye, we
demonstrate the usefulness of the localization by demonstrating a blood-vessel
following task, which can be challenging even for the surgeons when achieving
constant micron-scale off-set is the goal.

Figure 3: Phantom experiment setup

##  2 Related Work

###  2.1 Robotic Eye Surgery

Past works in robotic eye surgery have focused primarily on state estimation
or detection system that assist the surgeon by providing extra information.
Tayama et al. [2] modeled a proximity detection system that forecast detection
when the tool-tip and the shadow-tip converged by a pre-defined pixel
distance. Richa et al. [3] applied stereo vision setup that estimated the
depth of the tool and the retina respectively thereby demonstrating a
proximity detection system. Optical coherence tomography (OCT) has also been
integrated into both the surgical microscope and robotic tools. OCT is a laser
imaging modality that gives distance readings from the optical fiber tip to
the retina with micron-scale accuracy [4] [5]. However, due to its slow
processing speed and high-cost, its use in the operating room has been
limited. Contrary to laser-imaging based depth predictions, some works have
attempted to use stereo imaging to reason about depth including [6] and [3].
However, their studies have been limited to open-sky phantoms (i.e. no
refractive lens between the retina and the microscope) and dealing with
unknown distortions and out-of-focus images pose a significant challenge for
stereo-based methods. Ultrasound imaging, which is another alternative, has
generally lacked the necessary precision for retinal tool navigation but is
useful for certain anterior segment applications.

###  2.2 Learning

Many works have demonstrated the effectiveness of deep learning in
sensorimotor problems, such as in playing computer games [7, 8] and vision-
based navigation [9, 10, 11, 12, 13, 14]. Our network design borrows from an
approach proposed in [12] and [14], which uses a user-defined goal represented
as topographical binary image. Similarly, in our approach for designing the
network, we include a topographical binary image defining the mouse-click
position as input to a network. Directly relevant to eye surgery, recent works
have utilized OCT imaging and reinforcement learning to demonstrate needle
insertion using tool-pose states as input [15], although vision is not used.
More recently, a vision-based autonomous navigation system was developed for
retinal surgery, though without the inclusion of optimal control to impose
appropriate physical constraints [16]. Furthermore, our approach considers
predicting goal point to a surrounding object (i.e. the retina) rather than
predicting an incremental waypoint to the goal, which enables further
capabilities such as localization of the retina itself.

##  3 Problem Formulation

Consider a robotic manipulator with end-effector (i.e. surgical needle) state
defined by $x=(p,R,v,\omega)$ where $p\in\mathbb{R}^{3}$ defines the tool-tip
position, $R\in SO(3)$ the orientation matrix, $v\in\mathbb{R}^{3}$ the base-
frame translational velocity and $\omega\in\mathbb{R}^{3}$ the end-effector-
frame angular velocity. The end-effector pose is defined as $q=(p,R)\in
SE(3)$. The robot is fully actuated and for our purposes the control inputs
are the forces $u=(u_{v},u_{\omega})$ defining the translational forces
$u_{v}\in\mathbb{R}^{3}$ and torques $u_{\omega}\in\mathbb{R}^{3}$ in end-
effector frame (the robot internal controller can map those to appropriate
joint-angle torques). Let the robot end-effector occupy a region
$A(q)\subset\mathcal{W}$ in the workspace $\mathcal{W}\subset\mathbb{R}^{3}$.

The end-effector state $x$ is fully observable using high-precision joint
encoders and precise knowledge of the robot forward kinematics. In addition,
the system is equipped with workspace-fixed microscope video camera generating
observations $o(t)\in\mathcal{I}$ from the space of images $\mathcal{I}$
capturing the surgical scene.

The procedure begins when the surgeon inserts the tool into the sclera. The
sclera entry point $p_{s}\in\mathbb{R}^{3}$ is recorded by the robot and must
remain fixed after each entry during surgery. 111this point is often called
_remote center of rotation_ (RCM) as it acts as a fixed point through which
the instrument can rotate but only translate ”longitudinally” thus avoiding
lateral forces that could damage the sclera. The surgeon then selects a
sequence of $N_{g}$ pixels (i.e. 2-d points in image space) denoted by
$g_{i}\in o$ for $i=1,...,N_{g}$, that correspond to a set of euclidean points
(i.e. 3-d points in robot frame) on the surface of the retina inside the eye
to which the robot must navigate, denoted by $p_{i}\in E$, where
$E\subset\mathcal{W}$ defines the surface of the eye retina, i.e. the boundary
of the eye internal geometry.

The objective is to control the robot to safely navigate to the specified goal
locations. More formally, we seek to generate an end-effector trajectory
$x([t_{0},t_{f}])$ over some time-interval $[t_{0},t_{f}]$ defined as:

$\displaystyle\arg\\!\min_{\\!\\!\\!\\!\\!\\!\\!\\!x(\cdot),u(\cdot)}\int_{t_{0}}^{t_{f}}C(x(t),u(t))dt,$
:minimize cost, subject to: (1) $\displaystyle\dot{x}(t)=f(x(t),u(t)),$
:robot-tissue dynamics (2)
$\displaystyle(I+r_{x}(t)r_{x}(t)^{T})(p_{s}-p(t))=0,$ :sclera constraint (3)
$\displaystyle\mathcal{A}(q(t))\cap E=\emptyset,$ :collision constraint (4)
$\displaystyle p(t_{i})=p(t_{0})+R(t_{0})F(o(t_{0}),g_{i}),\text{ for
}i=1,\dots,N_{g},$ :visit goal locations (5) $\displaystyle x(t_{f})\in
X_{f},$ :terminal constraint (6)

where $C(x,u)$ is a given cost function e.g. ensuring smooth and safe motion,
and $r_{x}=Re_{x}$ with $e_{x}=(1,0,0)$ being the first basis vector. The
function $F(o,g)$ is the mapping from an image $o$ and a selected goal pixel
$g\in o$ at time $t_{0}$ to a 3-d position in end-effector-frame corresponding
to the exact goal location on the retinal surface. Note that neither the eye
geometry $E$ nor the mapping $F$ are available a priori and must be
constructed from available sensor data.

We consider the following two problems:

  1. 1.

_Reaching a given goal point_ : safely reach a single goal location, i.e.
$N_{g}=1$, corresponding to e.g. a desired retinal vein puncture point that we
reach with zero velocity, i.e. $x(t_{f})\in X_{f}$ simply corresponds to
$v(t_{f})=0$ and $\omega(t_{f})=0$.

  2. 2.

_Vessel following_ : generate a smooth and safe navigation through a sequence
of given goal image locations, where the robot does not need to stop at each
point, i.e. $N_{g}>1$ while the terminal zero-velocity constraint is applied
to the last goal point only.

Figure 4: Simulation setup; (left) side view (right) camera view navigation
setup

##  4 Technical Approach

To solve (1)–(6) we first address the computation of the relative goal point
map $F(o,g)$ which is accomplished through machine learning using prior
demonstrations with labeled contact points. The eye geometry $E$ is then
constructed by fusing multiple pixel position predictions (from the neural
network) through a least-squares approach to estimate the local euclidean
shape assuming locally spherical geometry. These two components are then used
in a discrete-time numerical optimal control framework that additionally
enforces sclera and collision constraints.

###  4.1 Learning to Predict Goal Point

We formulate the problem as a goal-based imitation learning setup, where the
objective is to estimate the required motion to place the surgical tool-tip
precisely at the user-specified goal position on the retina, similar to how
surgeons reason about moving the tool-tip to retina. We denote this relative
motion by the vector $d\in\mathbb{R}^{3}$, which encodes both the _distance_
to the goal $\|d\|$ and also the unit vector $\frac{d}{\|d\|}$ encoding the
_direction_ of motion. Specifically, given a data-set of tool trajectories
$D=\\{(o_{i},g_{i}),d_{i}\\}_{i=1}^{N}$, the objective is to construct a
function approximator $d=F((o,g);\theta)$ with parameters $\theta$, that maps
observation-goal pair to vector-to-goal generated by the expert. The objective
function can be expressed as the following:

$\displaystyle\arg\min_{\theta}\sum_{i=1}^{N}L(d_{i},F(o_{i},g_{i};\theta)),$
(7)

where $L$ is a given loss function. In our work, we choose the observation to
be a monocular image $o\in\mathcal{I}$ of the surgical scene viewed top-down
via a microscope, the goal input to be $g_{i}=(x_{i},y_{i})\in\mathbb{R}^{2}$
which specifies the desired projected 2-d position on the retinal surface,
corresponding to vector-to-goal $d$. The ground truth goal-prediction is
generated via robot kinematics.

###  4.2 Localizing the Retina

Using goal predictions across various points on the retina as described in
4.1, we fit a sphere to these points by fitting a radius and center point
using least-squares method. Specifically, we algebraically rearrange the terms
of the sphere equation, $\|p-p_{0}\|^{2}=r^{2}$, to obtain as the linear
relationship $\vec{f}=A\vec{c}$ with:

$\displaystyle\vec{f}=\left[\begin{array}[]{c}\|p_{i}\|^{2}\\\
\|p_{i+1}\|^{2}\\\ \vdots\\\ \|p_{n}\|^{2}\end{array}\right],\quad
A=\left[\begin{array}[]{cccc}2p_{i}^{T}&1\\\ 2p_{i+1}^{T}&1\\\
\vdots&\vdots\\\
2p_{n}^{T}&1\end{array}\right],\quad\vec{c}=\left[\begin{array}[]{c}p_{0}\\\
r^{2}-\|p_{0}\|^{2}\end{array}\right],$ (18)

where we solve for $\vec{c}$. In the equation,
$\\{p_{i},...,p_{n}\\}\in\mathcal{W}\subset\mathbb{R}^{3}$, are the estimated
points on the retinal surface obtained using the goal-prediction network, and
$p_{0}$ and $r$ is the estimated center position and radius of the fitted
sphere. The least-squares method is combined with random sample consensus
(RANSAC) to search for the solution with the highest number of inliers, based
on a prescribed inlier threshold.

###  4.3 Optimal Control Formulation

In order to reach a waypoint generated by the network, we must navigate within
physical constraints to avoid damage to the eye. To accomplish this reliably
in real-time we re-formulate the optimal control problem to include the
constraints (3)–(4) as least-squared penalties, according to:

$\displaystyle\begin{split}C=\frac{1}{2}x\left(t_{f}\right)^{T}P_{f}x\left(t_{f}\right)+\int_{t_{0}}^{t_{f}}&\frac{1}{2}\left[u(t)^{T}R(t)u(t)\right]\\\
&+w_{s}\cdot\|(I+r_{x}(t)r_{x}(t)^{T})(p_{s}-p(t))\|^{2}\\\
&+w_{e}\cdot\delta_{\\{\|d_{c}(t)\|>r-\epsilon\\}}\left\|d_{c}(t)-r\frac{d_{c}(t)}{\|d_{c}(t)\|}\right\|^{2}\mathrm{d}t,\end{split}$
(19)

where $p_{c}\in\mathbb{R}^{3}$ is the center, $d_{c}(t)=x(t)-p_{c}$, and
$r\in\mathbb{R}^{3}$ is the radius of the estimated sphere obtained using the
least-squares approach described in 4.2. In general the cost aims to minimize
the error in reaching the goal (encoded using $P_{f}\geq 0$ gain matrix),
control effort (using $R>0$ gain matrix), penalize deviation from sclera
points (using weight $w_{s}$), and penalize collision with the locally
spherical retinal surface if the tool is within some $\epsilon>0$ from it
(using weight $w_{e}$), encoded by the delta function
$\delta_{\\{\|d_{c}(t)\|>r-\epsilon\\}}$. The optimal control problem is
solved numerically using differential dynamic programming (DDP) based on a
discrete-time quantization of the robot motion using some fixed time-step
$dt$. For autonomous navigation to a single goal point, we omit the non-
penetration constraint to the retina. For vessel following we include both
sclera the non-penetration constraint in order navigate along the blood-
vessels while staying close to the retina. We note that with this re-
formulation, the hard constraints are now enforced as soft constraints. In
practice, the constraints can still be ensured by setting the gains of the
safety constraints to be significantly higher than the gains of the smoothness
and the goal-error terms. Furthermore, this re-formulation enables increased
numerical stability and efficiency.

Figure 5: Results on autonomous navigation, eye localization, and vessel
following tasks

##  5 Method

###  5.1 Eye Phantom and Simulation Experimental Setup

Our setup consists of a surgical robot, rubber eye model, microscope, surgical
needle tool, and light source (chandelier) as shown in Fig.3. The rubber eye
is mounted on a piezo-electric XYZ stage in order to move the eye to simulate
involuntary movements of the patient. We also mount a refractive lens between
the eye model and the microscope to simulate the distortions encountered
during real surgery e.g. the human lens and the vitreous. The dimension of the
eye model we used is 25.4mm, which is slightly larger than the average human
eye size, which ranges from 20 - 22.4mm [17]. There are however 25.4 mm human
eyes in myopic individuals. For sense of scale, the width of the surgical
tool-tip measures 0.32mm, and the operation space is 12x16mm (Fig.5).

We also demonstrate our setup in a simulated environment using Unity3D. We
created an environment similar to our phantom setup by including refractive
lens in to the scene. Furthermore, we explored a more involved domain
randomization by changing the orientation and the position of the light
source, and changing the orientation and position of the refractive lens. The
benefit of using the simulator is that it gives an approximate upper-bound
performance we can expect from phantom experiments, since the environment
parameters can be perfectly controlled without errors.

###  5.2 Data Collection

The data collection procedure consists of navigating the surgical tool-tip to
various locations on the retinal surface. We first estimate the retinal
geometry by sampling its surface with the surgical tool and fitting a sphere
to the sampled points using the least squares method (4.2). Once the retinal
surface is localized, we demonstrate many tool navigation trajectories across
the retinal surface starting from various initial configurations of the tool.
For simplicity, for each collected trajectory, we navigate to a single
targeted location in a straight-line trajectory. Since the retinal surface had
been localized, the tool-tip makes very gentle contact with the retina.
Throughout the trajectory, images from the microscope and tool-tip positions
based on robot kinematics are saved. After each trajectory demonstration, the
piezo-stage was moved $\pm 0.30$mm in XYZ directions for domain randomization.
We also used two different eyes at varying orientation to generalize to
changing eye textures. The pseudocode for data collection is summarized in
Appendix A.

###  5.3 Network Training

To train the goal-prediction network, we chose Resnet-18. The input of the
network are current surgical image (640x480x3) and the goal image (640x480x1)
stacked along the channel dimension. The goal image marks the desired landing
position of the tool-tip on the retina. Specifically, given $n$ frames
$I_{1},...I_{n}\in\mathcal{I}$, $n$ vector-to-goal values
$d_{1},\dots,d_{n}\in\mathbb{R}^{3}$, and a goal image coordinate
$g\in\mathbb{R}^{2}$ of this trajectory, the input-output relationship can be
expressed as $(\text{input},\text{output})=\left((I_{t},g),d_{t}\right)$, for
$t=1,\dots,n$. We also took the approach of discretizing the output space of
the network, as it proved to learn better than a continuous output case.
Specifically, we discretized the continuous x, y, z coordinates into 580,
1345, and 320 bins. We set the optimizer to be the Adam optimizer [18]. The
initial learning rate and batch size were set as 0.0003 and 128. We express
the loss functions as

$\displaystyle
L(b,\hat{p})=\sum_{j\in\\{x,y,z\\}}\sum_{c=1}^{M_{j}}-b_{j,c}\log(\hat{p}_{j,c}),$
(20)

where $b_{j,c}$ are binary indicators for the true class label $c$, and
$\hat{p}_{j,c}$ are the predicted probability that the coordinate $j$ is of
class $c$. Errors from all three dimensions $j\in\\{x,y,z\\}$ have been
combined into one cost function. As specified above, we employed $M_{x}=580,\
M_{y}=1345,\ M_{z}=320$ bins.

##  6 Results and Discussions

To assess the accuracy of our results, we implemented a benchmark task where
the network must navigate the tool to fifty predefined locations on the retina
as shown in Fig. 5. The objective of this task is to assess how well and the
network can navigate to various points on the retina. For our second task, we
use the goal-prediction network to sample an array of points across the
retinal surface to obtain many distance readings to the retina. These points
were then used to reconstruct the eye surface using the least squares
approach. To confirm our solution, we used the estimated retinal geometry to
navigate to the phantom retina at random and checked from side view whether it
was making gentle contact. Also, the estimated geometry was compared to a
solution achieved by manually sampling the eye with the tool-tip, which we
regarded as the ground-truth. For our last task, we used the estimated
geometry of the eye to follow a specific segment of a blood vessel while
hovering above to avoid collision. The result of these tasks are shown in Fig.
5.

Our numeric results are summarized in Table 1. In the first column of Table 1,
we demonstrate that we can navigate to various positions on the retina with an
average xy accuracy of 0.072mm and 0.130mm respectively. The visual result is
shown in the leftmost figure of Fig 5. The white squares denote the fifty
array of positions to visit and the blue squares denote the final landing
position of the tool-tip navigated by the network. For sense of scale, the
tool-tip width measures 0.32mm. In the second column, we demonstrate that our
navigation system is robust to the presence of an unseen object. The unseen
distractor object is a forcep, which is a commonly used surgical tool.
Throughout testing, We did not observe a visual difference in performance, and
thus the average error was similar to without-distraction scenario at 0.089mm
and 0.118mm in xy error. The visual result of this test is shown in the second
image in Fig.5. For our second task of localizing the eye model, we
demonstrate an error of 0.53mm, 0.084mm and 0.292mm in estimating the center
location. Specifically, the ground-truth center of the eye was located at
(39.26, -89.58, -14.32)mm and our estimated center was estimated at (38.73,
-89.66, -14.61)mm location. To confirm whether the estimated center was
reasonable, we sampled the eye 40 times at random using the estimated center
and radius and visually confirmed that the tool-tip gently touched phantom
retina every time. Using this estimated retinal geometry, we applied it to a
vessel-following task to follow a trajectory along a chosen blood-vessel. To
ensure that no collision occurs, we decreased radius of the estimated sphere
by 0.2mm, and an appropriate trajectory was generated by the optimal control
framework above the retinal surface at 0.2mm height. The trajectory-tracking
error of vessel following is reported in the last column of Table 1.
Furthermore, the sclera error of all task is shown on the fourth row of Table
1. No trajectories exceeded a mean sclera error of 0.704mm across all tasks,
satisfying required constraints for safe surgery. We also note that eye
geometry estimation and vessel following tasks require short and easy
procedures. Specifically, eye geometry estimation may take less than a second
as it simply requires multiple surgical images. Vessel-following may require
2-3 seconds as it requires the surgeon to click the desired waypoints of
interest using a mouse.

For our simulation results, we only demonstrated the autonomous navigation
task. The overall combined xy error in reaching 50 predefined locations was
0.119mm, which is significantly less than the errors encountered in the
phantom experiments. This is expected since simulation environment states and
collision detection are perfect. Still, given the fact the simulation involved
more substantial randomization in lens movement and light source movement, it
implies similar sophistication can be learned in the real-world and
potentially in real surgery (Appendix D).

Table 1:  Error on Various Tasks Error category |  Navigation to retina (50 positions) (mm, pixels) |  Navigation to retina w/ forcep distraction (50 positions) (mm, pixels) |  Eye localization center, radius (mm) |  Vessel-following trajectory tracking (mm)
---|---|---|---|---
x error | 0.072, 1.8 | 0.089, 2.24 | 0.53 | 0.002
y error | 0.130, 3.26 | 0.118, 2.94 | 0.084 | 0.003
z error | NA | NA | 0.292 | 0.006
sclera error | 0.704mm | 0.613mm | NA | 0.059
radius | NA | NA | 0.152 | NA

##  7 Conclusion

In this work, we demonstrate that autonomous navigation is possible by
learning to predict goal positions on the retina from user-specified image
pixel locations. Our work has been tested in simulated eyes as well as eye
phantom models where eye movement, different textures, refractive view, and
various lighting conditions are encountered. We demonstrate that through a
benchmark experiment, our method is capable of achieving higher accuracy in
navigating to a target tissue than a surgeon, whose error is characterized by
mean-amplitude of hand-tremor at the tool-tip. Still, more work is necessary
to adopt our work to real human surgeries. In addition, human eyes have a
slightly elliptical shape and the proposed spherical regression must be
generalized accordingly. Finally, rigorous experiments with real physical
tissue are critical to enable to application of the proposed techniques in a
clinical setting. To this end, we are beginning to employ porcine eyes, which
share a similar anatomy as human eyes.

## References

  * MAYNE [1966] B. D. MAYNE.  A second-order gradient method for determining optimal trajectories of non-linear discrete-time systems.  _International Journal of Control_ , 3(1):85–95, 1966.  doi:10.1080/00207176608921369.  URL https://doi.org/10.1080/00207176608921369.
  * Tayama et al. [2018] T. Tayama, Y. Kurose, M. M. Marinho, Y. Koyama, K. Harada, S. Omata, F. Arai, K. Sugimoto, F. Araki, K. Totsuka, M. Takao, M. Aihara, and M. Mitsuishi.  Autonomous positioning of eye surgical robot using the tool shadow and kalman filtering.  In _2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , pages 1723–1726, July 2018\.  doi:10.1109/EMBC.2018.8512633.
  * Richa et al. [2012] R. Richa, M. Balicki, R. Sznitman, E. Meisner, R. Taylor, and G. Hager.  Vision-based proximity detection in retinal surgery.  _IEEE Transactions on Biomedical Engineering_ , 59(8):2291–2301, Aug 2012.  ISSN 0018-9294.  doi:10.1109/TBME.2012.2202903.
  * Carrasco-Zevallos et al. [2016] O. M. Carrasco-Zevallos, B. Keller, C. Viehland, L. Shen, G. Waterman, B. Todorich, C. Shieh, P. Hahn, S. Farsiu, A. N. Kuo, C. A. Toth, and J. A. Izatt.  Live volumetric (4d) visualization and guidance of in vivo human ophthalmic surgery with intraoperative optical coherence tomography.  _Scientific Reports_ , 6(1):31689, Aug 2016.  ISSN 2045-2322.  doi:10.1038/srep31689.  URL https://doi.org/10.1038/srep31689.
  * Smits et al. [2018] J. Smits, M. Ourak, A. Gijbels, L. Esteveny, G. Borghesan, L. Schoevaerdts, K. Willekens, P. Stalmans, E. Lankenau, H. Schulz-Hildebrandt, G. Hüttmann, D. Reynaerts, and E. B. Vander Poorten.  Development and experimental validation of a combined fbg force and oct distance sensing needle for robot-assisted retinal vein cannulation.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 129–134, 2018.
  * Probst et al. [2018] T. Probst, K. Maninis, A. Chhatkuli, M. Ourak, E. V. Poorten, and L. Van Gool.  Automatic tool landmark detection for stereo vision in robot-assisted retinal surgery.  _IEEE Robotics and Automation Letters_ , 3(1):612–619, Jan 2018.  ISSN 2377-3774.  doi:10.1109/LRA.2017.2778020.
  * Mnih et al. [2015] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, et al.  Human-level control through deep reinforcement learning.  _Nature_ , 518(7540):529, 2015.
  * Dosovitskiy and Koltun [2016] A. Dosovitskiy and V. Koltun.  Learning to act by predicting the future, 2016.
  * Wang et al. [2018] X. Wang, Q. Huang, A. Çelikyilmaz, J. Gao, D. Shen, Y. fang Wang, W. Y. Wang, and L. Zhang.  Reinforced cross-modal matching and self-supervised imitation learning for vision-language navigation.  _ArXiv_ , abs/1811.10092, 2018.
  * Bojarski et al. [2016] M. Bojarski, D. D. Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D. Jackel, M. Monfort, U. Muller, J. Zhang, X. Zhang, J. Zhao, and K. Zieba.  End to end learning for self-driving cars, 2016.
  * Xu et al. [2017] H. Xu, Y. Gao, F. Yu, and T. Darrell.  End-to-end learning of driving models from large-scale video datasets.  In _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 3530–3538, July 2017.  doi:10.1109/CVPR.2017.376.
  * Amini et al. [2018] A. Amini, W. Schwarting, G. Rosman, B. Araki, S. Karaman, and D. Rus.  Variational autoencoder for end-to-end control of autonomous driving with novelty detection and training de-biasing.  In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pages 568–575, Oct 2018.  doi:10.1109/IROS.2018.8594386.
  * Codevilla et al. [2018] F. Codevilla, M. Miiller, A. López, V. Koltun, and A. Dosovitskiy.  End-to-end driving via conditional imitation learning.  _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1–9, 2018.
  * Bansal et al. [2018] M. Bansal, A. Krizhevsky, and A. Ogale.  Chauffeurnet: Learning to drive by imitating the best and synthesizing the worst, 2018.
  * Keller et al. [2020] B. Keller, M. Draelos, K. Zhou, R. Qian, A. N. Kuo, G. Konidaris, K. Hauser, and J. A. Izatt.  Optical coherence tomography-guided robotic ophthalmic microsurgery via reinforcement learning from demonstration.  _IEEE Transactions on Robotics_ , 2020.
  * Kim et al. [2020] J. W. Kim, C. He, M. Urias, P. Gehlbach, G. D. Hager, I. Iordachita, and M. Kobilarov.  Autonomously navigating a surgical tool inside the eye by learning from demonstration.  In _2020 International Conference on Robotics and Automation (ICRA)_. IEEE, 2020.
  * Bekerman et al. [2014] Bekerman, Gottlieb, Paul, and Michael.  Variations in eyeball diameters of the healthy adults, Nov 2014.  URL https://www.hindawi.com/journals/joph/2014/503645/.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.

##  A Pseudocode for Data Collection and Inference

sample 30 points on the eye phantom;

center, radius = leastSquares(sampled points);

while  _iter $<$ n_ do

choose random initial tool-tip point and random sclera point;

goal point = generateRandomGoalOnRetina(center, radius);

randomize xyz position of the eye;

while  _goal point not reached_ do

moveToward(goal point);

record image and robot kinematics;

end while

end while

Algorithm 1 Data Collection

for  _each [initial point, goal point]_ do

while  _numIter $<$ maxIter_ do

waypoint = CNN(goal point, current image);

trajPoints, trajVelocity = optimalControl(waypoint, current tool-tip
position);

if  _z-distance to goal $>$ 0.1mm_ then

closedLoopControlToReachGoal(trajPoints, trajVelocity);

end if

end while

end for

Algorithm 2 Inference

##  B Error Calculation for the Benchmark Task

For the benchmark task of navigating to fifty points across the retina, the
error in reaching each goal point was calculated using the following formula:
$x_{error}=\|{(x-x^{{}^{\prime}})}\|$, $y_{error}=\|{(y-y^{{}^{\prime}})}\|$
where $(x,y)$ is the desired point to reach in image coordinates and
$(x^{{}^{\prime}},y^{{}^{\prime}})$ is the final landing position of the tool-
tip in image coordinates. The image coordinates were converted into mm by
finding the number of pixels of a known dimension. Specifically, the tool-tip
measured 0.320mm and 8pixels in image coordinates. Thus, the pixel error
values were converted to mm by multiplying the ratio of 0.320mm/8pixels.

##  C Eye Models

Figure 6: The two eyes on the left were used for training and the right eye
was used for testing.

As shown in Fig. 6, the two left eyes were used for training and the right eye
was used for testing. We fabricated all the eyes by molding the appropriate
polymer and drawing the blood-vessels using water-based paint. All eyes were
molded using a 25.4 $\pm$ 0.05mm diameter sphere.

##  D Simulation Details

Figure 7: (Left) Benchmark task of visiting 100 array of positions in
simulation; the center of each block represents each point to visit for
visualization purpose. (Right) Raycast is used to find the ground-truth xyz
points of the goal points to calculate errors accurately.

As shown in Fig.7, the benchmark task of visiting 100 array of points is
demonstrated. In the left picture of Fig.7, each goal-point-to-visit is
represented as a block for visualization purpose. In the right image of Fig.7,
100 raycasts are visualized, which are used to find the exact xyz locations of
the retina at the desired goal positions for accurate error analysis. We also
randomized the positions of the objects in the scene for robust generalization
to changing objects in the scene: the eye was moved $\pm$0.5mm in xyz-axes.
The lens position was moved $\pm$0.5mm in xyz-axes and rotated $\pm$3degrees
in xyz-axes. The light source was also moved $\pm$1mm in xy-axes and
$\pm$0.5mm in z-axis and its orientation was moved $\pm$1.5degrees along
x-axis and $\pm$2.5 degrees along yz-axes.