# High-accuracy three-dimensional surface detection in smoothed particle
hydrodynamics for free-surface flows

Wen-Bin Liu  Dong-Jun Ma  Jian-Zhen Qian  Ming-Yu Zhang  An-Min He  Nan-Sheng
Liu  Pei Wang  Department of Modern Mechanics, University of Science and
Technology of China, Hefei, Anhui, 230027, China  Institute of Applied Physics
and Computational Mathematics, Beijing 100094, China  Graduate School of China
Academy of Engineering Physics, Beijing 100088, China  Center for Applied
Physics and Technology, Peking University, Beijing 100871, China

###### Abstract

In this study, we investigate high-accuracy three-dimensional surface
detection in smoothed particle hydrodynamics for free-surface flows. A new
geometrical method is first developed to enhance the accuracy of free-surface
particle detection in complex flows. This method detects free-surface
particles via continuous global scanning inside the sphere of a particle
through a cone region whose vertex corresponds to the particle position. The
particle is identified as a free-surface particle if there exists a cone
region with no neighboring particles. Next, an efficient semi-geometrical
method is proposed based on the geometrical method to reduce the computational
cost. It consists of finding particles near the free surface via position
divergence and then detecting these particles using the geometrical method to
identify free-surface particles. The accuracy and robustness of the proposed
method are demonstrated by performing tests on several model problems.

###### keywords:

Free-surface particle detection , Geometrical method , Semi-geometrical method
, Smoothed particle hydrodynamics

##  1 Introduction

Smoothed particle hydrodynamics (SPH) is a representative Lagrangian particle
method. In the SPH method, the materials are discretized into a series of
particles with independent information. Because of its Lagrangian and meshless
nature, the SPH method has the advantage of simulating free-surface and multi-
phase flows with large deformations, hence, it is widely applied to scientific
and engineering problems [1, 2, 3, 4, 5, 6].

In many aspects of the SPH method, such as the implementation of boundary
conditions [7, 8], particle shifting technology (PST) [9, 10], and continuum
surface force (CSF) model [11, 12], it is important to know exactly which
particles belong to the free surface. Generally, two methods can be applied to
detect free-surface particles: the kernel and geometrical methods [13, 14, 15,
12, 16, 17, 18, 19, 20, 21, 22]. In the kernel method, the sum of the kernel
function or its gradients over the neighboring particles is used to detect
free-surface particles [14, 15]. Although the kernel method is easy to
implement, its accuracy depends on the particle distribution. In the
geometrical method, Dilts [13] first developed a two-dimensional exposure
method. If the circle of a particle is not completely covered by the circles
of neighboring particles, the particle is identified as a free-surface
particle. Sun et al. [12] developed an improved two-dimensional geometrical
method by performing continuous global scanning inside the circle of a
particle through a sector region. The vertex of the sector corresponds to the
particle position. If there exists a sector region with no neighboring
particles, the particle is identified as a free-surface particle. Haque et al.
[17] extended the two-dimensional exposure method to be a three-dimensional
version. If the sphere of a particle is not completely covered by the spheres
of neighboring particles, the particle is identified as a free-surface
particle. Marrone et al. [18] proposed a semi-geometrical method for free-
surface particle detection. The particles near the free surface are first
found by the minimum eigenvalue of a renormalization matrix. Then, the free-
surface particles are accurately detected by a scan region along the normal
direction of the particle. Based on the method proposed by Marrone et al.
[18], Wang et al. [20] proposed an optimized semi-geometrical method for free-
surface particle detection. This method makes use of the position divergence
to find free-surface candidate particles. This considerably reduces the
computational time for three-dimensional problems. However, owing to the
estimation of the normal direction of the particle, both methods [18, 20] are
sensitive to particle distribution.

To improve the accuracy of the free-surface particle detection in three-
dimensional complex flows, we develop a new geometrical method for the three-
dimensional SPH method, extending the two-dimensional geometrical method of
Sun et al. [12]. Our method detects free-surface particles by performing
continuous global scanning within the sphere of a particle through a cone
region. The vertex of the cone corresponds to the particle position. If there
exists a cone region with no neighboring particles, the particle is identified
as a free-surface particle, otherwise, it is identified as an inner particle,
i.e., there is at least one neighboring particle in any cone region. Next,
based on the proposed geometrical method, a semi-geometrical method is
developed to reduce computational costs. Finally, the proposed method is
applied to several tests, and the results of the free-surface particle
detection are compared with those of the previous methods [17, 18, 20].

This paper is organized as follows. In Section 2, the SPH method for weakly-
compressible flows is introduced. Then, in Section 3, we present a new method
for free-surface particle detection. In Section 4, the proposed method is
validated through four tests. Finally, the conclusions are drawn in Section 5.

##  2 SPH method

###  2.1 Governing equations

The Navier$-$Stokes equations for weakly-compressible fluids in Lagrangian
description are given as follows:

$\left\\{\begin{array}[]{l}\displaystyle\frac{D\rho}{Dt}=-\rho\nabla\cdot\bm{u},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\displaystyle\frac{D\bm{u}}{Dt}=-\displaystyle\frac{{1}}{{\rho}}\nabla
p+\text{div}\mathbf{V}+\bm{g}+\frac{1}{\rho}\bm{f}^{s},\end{array}\right.$ (1)

where $\rho$ and $p$ denote density and pressure, respectively. The $\bm{u}$,
$\bm{g}$, and $\bm{f}^{s}$ represent velocity vector, gravitational
acceleration, and surface tension force, respectively. $\mathbf{V}$ denotes
the viscous stress tensor.

A linear state equation is used to close Eqs. (1) are as follows:

$p=c_{0}^{2}\left(\rho-\rho_{0}\right),$ (2)

where $\rho_{0}$ and $c_{0}$ denote the density and speed of sound under rest
condition, respectively.

For weakly-compressible fluids [23], the speed of sound needs to satisfy the
following conditions:

$c_{0}\geq 10\max\left(\sqrt{\frac{p_{\max}}{\rho_{0}}},U_{\max}\right),$ (3)

where $p_{\max}$ and $U_{\max}$ denote the expected maximum pressure and
velocity, respectively.

###  2.2 SPH equations

The $\delta$-SPH method [24] is used to discretize the governing equations, in
which a diffusive term for stabilizing the pressure field is added to the
continuity equation. The discrete governing equations with diffusive terms are
given as follows:

$\left\\{\begin{array}[]{l}\displaystyle\frac{D\rho_{i}}{Dt}=-\rho_{i}\sum\limits_{j=1}^{N}{\left(\bm{u}_{j}-\bm{u}_{i}\right)\cdot\nabla_{i}W_{ij}V_{j}}+\delta
h_{i}c_{0}D_{i},\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\displaystyle\frac{D\bm{u}_{i}}{Dt}=\displaystyle-\frac{1}{\rho_{i}}\sum\limits_{j=1}^{N}{\left(p_{i}+p_{j}\right)\nabla_{i}W_{ij}V_{j}}+\displaystyle\frac{1}{\rho_{i}}2\left(d+2\right)\mu\sum\limits_{j=1}^{N}{\pi_{ij}\nabla_{i}W_{ij}V_{j}}+\bm{g}+\frac{1}{\rho_{i}}\bm{f}_{i}^{s},\end{array}\right.$
(4)

where $N$, $d$, and $V$ denote the number of neighboring particles, spatial
dimensions, and volume of a particle, respectively. The parameter $\delta$ is
set to 0.1 [25].

The term $D_{i}$ in the continuity equation is given as follows:

$\left\\{\begin{array}[]{l}\displaystyle
D_{i}=2\sum\limits_{j=1}^{N}{\psi_{ij}\frac{\left(\bm{x}_{j}-\bm{x}_{i}\right)\cdot\nabla_{i}W_{ij}}{\left|\bm{x}_{j}-\bm{x}_{i}\right|^{2}}},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\displaystyle\psi_{ij}=\left(\rho_{j}-\rho_{i}\right)-\frac{1}{2}\left(\left<\nabla\rho\right>_{i}^{L}+\left<\nabla\rho\right>_{j}^{L}\right)\cdot\left(\bm{x}_{j}-\bm{x}_{i}\right),\end{array}\right.$
(5)

where $\bm{x}_{i}$ denotes the particle position vector.
$\left<\nabla\rho\right>^{L}$ is expressed as:

$\left\\{\begin{array}[]{l}\left<\nabla\rho\right>_{i}^{L}=\displaystyle\sum\limits_{j=1}^{N}{\left(\rho_{j}-\rho_{i}\right)\mathbf{L}_{i}\nabla_{i}W_{ij}V_{j}},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\mathbf{L}_{i}=\displaystyle\left[\sum\limits_{j=1}^{N}{\left(\bm{x}_{j}-\bm{x}_{i}\right)\otimes\nabla_{i}W_{ij}V_{j}}\right]^{-1}.\end{array}\right.$
(6)

The term $\pi_{ij}$ in the momentum equation is given as follows:

$\pi_{ij}=\frac{\left(\bm{u}_{j}-\bm{u}_{i}\right)\cdot\left(\bm{x}_{j}-\bm{x}_{i}\right)}{\left|\bm{x}_{j}-\bm{x}_{i}\right|^{2}}.$
(7)

The Wendland C2 kernel function is used, which is written as follows:

$W\left(\left|\bm{x}_{j}-\bm{x}_{i}\right|,h\right)=\left\\{\begin{array}[]{ll}\displaystyle\alpha_{d}\left(1+2q\right)\left(1-\frac{q}{2}\right)^{4},&0\leq
q\leq 2\vskip 6.0pt plus 2.0pt minus 2.0pt\\\ \displaystyle
0,&q>2\end{array}\right.$ (8)

where $q=\left|\bm{x}_{j}-\bm{x}_{i}\right|/h$ and $h$ denotes the smoothing
length. The parameters $\alpha_{d}$ is given as follows:

$\alpha_{d}=\left\\{\begin{array}[]{ll}\displaystyle\frac{3}{4h},&d=1\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ \displaystyle\frac{7}{4\pi h^{2}},&d=2\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ \displaystyle\frac{21}{16\pi
h^{3}},&d=3\end{array}\right.$ (9)

The PST developed by Wang et al. [20] is used to improve the spatial
configuration of the particles and is given as follows:

$\bm{x}_{i}^{*}=\bm{x}_{i}+\delta\bm{x}_{i},$ (10)
$\delta\bm{x}_{i}=\left\\{\begin{array}[]{lll}0&&{\rm if}\;i\in F\vskip 6.0pt
plus 2.0pt minus 2.0pt\\\ -k_{\rm CFL}\cdot{\rm Ma}\cdot
2h^{2}\cdot\displaystyle\sum\limits_{j=1}^{N}{\frac{m_{j}}{\rho_{j}}R\left(\frac{W_{ij}}{W\left(\Delta
x_{i}\right)}\right)^{n}\nabla_{i}W_{ij}}&j\in\left[\left|\bm{x}_{ij}\right|<l_{i}\right]&{\rm
if}\;i\in V\vskip 6.0pt plus 2.0pt minus 2.0pt\\\ -k_{\rm CFL}\cdot{\rm
Ma}\cdot
2h^{2}\cdot\displaystyle\sum\limits_{j=1}^{N}{\frac{m_{j}}{\rho_{j}}\left[1+R\left(\frac{W_{ij}}{W\left(\Delta
x_{i}\right)}\right)^{n}\right]\nabla_{i}W_{ij}}&j\in\left[\left|\bm{x}_{ij}\right|<kh\right]&{\rm
if}\;i\in I\end{array}\right.$ (11)

where $\delta\bm{x}_{i}$ and $\bm{x}_{i}^{*}$ denote the particle shifting
vector and new particle position vector, respectively. $F$, $V$, and $I$
denote the free-surface region, free-surface vicinity region, and inner
region, respectively. $n$, $R$, and $k_{\rm CFL}$ are set to 4, 0.2, and 0.25,
respectively.

##  3 High-accuracy method for three-dimensional free-surface particle
detection

This section demonstrates the high-accuracy method for three-dimensional free-
surface particle detection in the SPH method. First, a high-accuracy
geometrical method for three-dimensional free-surface particle detection is
developed. Continuous global scanning within the sphere of a particle is
performed through a cone region, rather than along the normal direction alone,
which is difficult in ensuring accurate computation of complex flows. The
particle is identified as a free-surface particle if there exists a cone
region with no neighboring particles. Then, based on the proposed geometrical
method, a semi-geometrical method is developed to reduce computational costs.

The symbols employed in this section are defined as follows. Let $S_{i}$
denote the sphere with radius $R$ centered at particle $i$, and $C_{j}$ denote
the cone of rotation angle $\theta$ and rotation axis $l_{ij}$ with particle
$i$ as the vertex, where $l_{ij}$ is the line passing through particle $i$ and
its neighboring particle $j$. Moreover, $C_{ij}$ is used to denote the circle
obtained by intersecting the sphere $S_{i}$ with cone $C_{j}$. See Fig. 1 for
details.

###  3.1 Geometrical method for free-surface particle detection

The key idea of the geometrical method is to perform continuous global
scanning inside the sphere $S_{i}$ through a cone region $C$. The vertex of
the cone $C$ corresponds to the position of particle $i$, and the rotation
angle of cone $C$ is $\theta$. If there exists a cone region with no
neighboring particles, the particle $i$ is identified as a free-surface
particle; otherwise, it is identified as an inner particle, i.e., there is at
least one neighboring particle in any cone region, as shown in Fig. 1a.
Because the number of such cones is infinite during the continuous global
scanning, the aforementioned idea cannot be directly implemented in the
program. However, this idea can be equivalently expressed as follows: the
sphere $S_{i}$ intersects the cone $C_{j}$ of all neighboring particles and
the circles $C_{ij}$ are drawn on the surface of sphere $S_{i}$. If the
surface of sphere $S_{i}$ is not completely covered by circles $C_{ij}$, i.e.,
some parts of the circle $C_{ij}$ are exposed, the particle $i$ is identified
as a free-surface particle; otherwise, it is identified as an inner particle,
as shown in Fig. 1b. The logical proof of the aforementioned equivalence
transformation is given as follows: when the surface of sphere $S_{i}$ is not
completely covered by the circles $C_{ij}$, for any point $m$ on the surface
of sphere $S_{i}$ not completely covered by circles $C_{ij}$, the angle of ray
$l_{im}$ and all neighboring rays $l_{ij}$ is greater than $\theta$, and there
must be no neighboring particles in the cone $C_{m}$ with rotation angle
$\theta$ and rotation axis $l_{im}$. Therefore, the particle $i$ is identified
as a free-surface particle. In contrast, when the surface of sphere $S_{i}$ is
completely covered by the circles $C_{ij}$, for any point $m$ on the surface
of sphere $S_{i}$, it must be within a circle $C_{ij}$, and the angle of ray
$l_{im}$ and $l_{ij}$ is less than $\theta$, then there is at least one
neighboring particle in the cone $C_{m}$ with the rotation angle $\theta$ and
the rotation axis $l_{im}$. Therefore, the particle $i$ is identified as an
inner particle. Algorithm 1 gives the main operations of the present
geometrical method. The remainder of this section focuses on the following
functions:

  * 1.

Intersection of sphere and cone

  * 2.

Intersection of two circles

  * 3.

Checking coverage of circle

Figure 1: Black circle denotes candidate free-surface sphere $S_{i}$ with
radius $R$ centered at particle $i$. (a) Black cone denotes scanning cone of
rotation angle $\theta$ and rotation axis $l$ with particle $i$ as the vertex.
(b) Red cone denotes cone $C_{j}$ of rotation angle $\theta$ and rotation axis
$l_{ij}$ with particle $i$ as the vertex, where $l_{ij}$ is the line passing
through particle $i$ and neighboring particle $j$. Blue cone denotes cone
$C_{k}$ of rotation angle $\theta$ and rotation axis $l_{ik}$ with particle
$i$ as the vertex, where $l_{ik}$ is the line passing through particle $i$ and
neighboring particle $k$. Red $C_{ij}$ denotes the circle given by the
intersection of sphere $S_{i}$ and cone $C_{j}$. Blue $C_{ik}$ denotes the
circle given by the intersection of sphere $S_{i}$ and cone $C_{k}$.

for  _$i$ over all particles_ do

for  _$j$ over neighboring particles of particle $i$_ do

$C_{ij}$ = _Intersection of sphere and cone_ ($S_{i},C_{j}$)

for  _$k$ over neighboring particles of particle $i$_ do

if  _angle $\varphi_{ijk}$ of rays $l_{ij}$ and $l_{ik}$ is greater than
$0^{\circ}$ and less than $2\theta$_ then

Add the particle $k$ to the set $\mathcal{L}$

end if

end for

Sort the set $\mathcal{L}$ from smallest to largest $\varphi_{ijk}$

for  _$k$ over the particles of set $\mathcal{L}$_ do

$C_{ik}$ = _Intersection of sphere and cone_ ($S_{i},C_{k}$)

_Intersection of two circles_ ($C_{ij},C_{ik}$)

_Checking coverage of circle_

if  _$C_{ij}$ is completely covered by the circles $C_{ik}$_ then

$\mathbf{continue}$ with the next $j$

end if

end for

if  _$C_{ij}$ is not completely covered by the circles $C_{ik}$_ then

Particle $i$ is identified as the free-surface particle

$\mathbf{continue}$ with the next $i$

else

$\mathbf{continue}$ with the next $j$

end if

end for

if  _for every $j$, the $C_{ij}$ is completely covered by the circles
$C_{ik}$_ then

Particle $i$ is identified as the inner particle

else

Particle $i$ is identified as the free-surface particle

end if

end for

Algorithm 1 Geometrical method for free-surface particle detection

####  3.1.1 Intersection of sphere and cone

Let $\bm{x}_{i}$ = ($x_{i},y_{i},z_{i}$) denote the position of particle $i$,
which is also the center of sphere $S_{i}$. Likewise, $\bm{x}_{j}$ =
($x_{j},y_{j},z_{j}$) denotes the position of the neighboring particle $j$.
Hence, the distance $d_{ij}$ between particles $i$ and $j$ is given as
follows:

$d_{ij}=\sqrt{\left(x_{j}-x_{i}\right)^{2}+\left(y_{j}-y_{i}\right)^{2}+\left(z_{j}-z_{i}\right)^{2}}.$
(12)

As shown in Fig. 2, if $d_{ij}\leq R$, the center $\bm{x}_{ij}$ =
($x_{ij},y_{ij},z_{ij}$) and radius $r_{ij}$ of circle $C_{ij}$ are given as
follows:

$\bm{x}_{ij}=\bm{x}_{i}+\frac{R\cos\theta}{d_{ij}}\left(\bm{x}_{j}-\bm{x}_{i}\right),$
(13) $r_{ij}=R\sin\theta.$ (14)

Figure 2: Center $\bm{x}_{ij}$, radius $r_{ij}$ and local coordinate system
for the circle $C_{ij}$.

To easily calculate the angle interval intersected by two circles, a local
coordinate system is established, in which $\bm{x}_{ij}$ is the origin, and
the unit axis vectors $\bm{\hat{x}}_{ij}$, $\bm{\hat{y}}_{ij}$, and
$\bm{\hat{z}}_{ij}$ are given as follows:

$\bm{\hat{z}}_{ij}=\frac{1}{d_{ij}}\left(\bm{x}_{j}-\bm{x}_{i}\right),$ (15)
$\bm{\hat{x}}_{ij}=\left\\{\begin{array}[]{ll}\displaystyle\frac{\bm{\hat{x}}\times\bm{\hat{z}}_{ij}}{\lVert\bm{\hat{x}}\times\bm{\hat{z}}_{ij}\rVert},&\;\text{if}\;\;\bm{\hat{y}}\times\bm{\hat{z}}_{ij}=0\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\displaystyle\frac{\bm{\hat{y}}\times\bm{\hat{z}}_{ij}}{\lVert\bm{\hat{y}}\times\bm{\hat{z}}_{ij}\rVert},&\;\text{if}\;\;\bm{\hat{y}}\times\bm{\hat{z}}_{ij}\neq
0\end{array}\right.$ (16)
$\bm{\hat{y}}_{ij}=\displaystyle\frac{\bm{\hat{z}}_{ij}\times\bm{\hat{x}}_{ij}}{\lVert\bm{\hat{z}}_{ij}\times\bm{\hat{x}}_{ij}\rVert},$
(17)

where $\bm{\hat{x}}$ and $\bm{\hat{y}}$ are the unit vectors of the original
coordinate system, respectively. The unit vector
$\bm{\hat{z}}_{ij}=\left(z_{ij}^{\left(1\right)},z_{ij}^{\left(2\right)},z_{ij}^{\left(3\right)}\right)$
is also the outward normal vector of the circle $C_{ij}$.

####  3.1.2 Intersection of two circles

Define
$\varphi_{ijk}=\cos^{-1}\left(\bm{\hat{z}}_{ij}\cdot\bm{\hat{z}}_{ik}\right)$
as the angle of the unit vectors $\bm{\hat{z}}_{ij}$ and $\bm{\hat{z}}_{ik}$.
Here, four possible scenarios can be derived:

(1) $\varphi_{ijk}=0^{\circ}$, circles $C_{ij}$ and $C_{ik}$ coincide.

(2) $0^{\circ}<\varphi_{ijk}<2\theta$ , circles $C_{ij}$ and $C_{ik}$
intersect at two points.

(3) $\varphi_{ijk}=2\theta$ , circles $C_{ij}$ and $C_{ik}$ intersect at one
point.

(4) $\varphi_{ijk}>2\theta$ , circles $C_{ij}$ and $C_{ik}$ do not intersect.

In cases 1 and 3, circles $C_{ij}$ and $C_{ik}$ are considered to have no
intersection.

Before calculating the angle intervals $A_{ijk}$ of the intersection of
circles $C_{ij}$ and $C_{ik}$, the intersection particles are sorted based on
the $\varphi_{ijk}$ from smallest to largest. The smaller the angle of
$\varphi_{ijk}$, the greater the angle intervals $A_{ijk}$, in this case, the
possibility of fully covering circle $C_{ij}$ increases rapidly. Once the
circle $C_{ij}$ is determined as covered, no more angle intervals $A_{ijk}$
need to be computed.

For case 2 where circles $C_{ij}$ and $C_{ik}$ intersect at two points, we
consider the plane containing the points $\bm{x}_{i}$, $\bm{x}_{ij}$ and
$\bm{x}_{ik}$. Its intersection with sphere $S_{i}$ is the circle $C_{ijk}$,
as shown in Fig. 3. The unit normal vector
$\bm{\hat{z}}_{ijk}=\left(z_{ijk}^{\left(1\right)},z_{ijk}^{\left(2\right)},z_{ijk}^{\left(3\right)}\right)$
of circle $C_{ijk}$ is:

$\bm{\hat{z}}_{ijk}=\displaystyle\frac{\bm{\hat{z}}_{ij}\times\bm{\hat{z}}_{ik}}{\lVert\bm{\hat{z}}_{ij}\times\bm{\hat{z}}_{ik}\rVert}.$
(18)

Figure 3: Intersection of two circles. The point $\bm{x}_{ijk}$ is the
intersection of the planes of circles $C_{ij}$, $C_{ik}$ and $C_{ijk}$.

The equations describing the planes of the three circles $C_{ij}$, $C_{ik}$,
and $C_{ijk}$ in the point-direction are:

$\left\\{\begin{array}[]{l}\bm{\hat{z}}_{ij}\cdot\left(\bm{x}-\bm{x}_{ij}\right)=0,\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\bm{\hat{z}}_{ik}\cdot\left(\bm{x}-\bm{x}_{ik}\right)=0,\vskip 3.0pt plus
1.0pt minus 1.0pt\\\
\bm{\hat{z}}_{ijk}\cdot\left(\bm{x}-\bm{x}_{i}\right)=0,\end{array}\right.$
(19)

where $\bm{x}=\left(x,y,z\right)$. The above equations can be simplified as:

$\left\\{\begin{array}[]{l}z_{ij}^{\left(1\right)}x+z_{ij}^{\left(2\right)}y+z_{ij}^{\left(3\right)}z=z_{ij}^{\left(1\right)}x_{ij}+z_{ij}^{\left(2\right)}y_{ij}+z_{ij}^{\left(3\right)}z_{ij},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
z_{ik}^{\left(1\right)}x+z_{ik}^{\left(2\right)}y+z_{ik}^{\left(3\right)}z=z_{ik}^{\left(1\right)}x_{ik}+z_{ik}^{\left(2\right)}y_{ik}+z_{ik}^{\left(3\right)}z_{ik},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
z_{ijk}^{\left(1\right)}x+z_{ijk}^{\left(2\right)}y+z_{ijk}^{\left(3\right)}z=z_{ijk}^{\left(1\right)}x_{i}+z_{ijk}^{\left(2\right)}y_{i}+z_{ijk}^{\left(3\right)}z_{i}.\end{array}\right.$
(20)

The point $\bm{x}_{ijk}=\left(x_{ijk},y_{ijk},z_{ijk}\right)$ , which is the
intersection of the planes of circles $C_{ij}$, $C_{ik}$, and $C_{ijk}$, can
be obtained using:

$\left[\begin{array}[]{l}x_{ijk}\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
y_{ijk}\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
z_{ijk}\end{array}\right]=\left[\begin{array}[]{lll}z_{ij}^{\left(1\right)}&z_{ij}^{\left(2\right)}&z_{ij}^{\left(3\right)}\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
z_{ik}^{\left(1\right)}&z_{ik}^{\left(2\right)}&z_{ik}^{\left(3\right)}\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
z_{ijk}^{\left(1\right)}&z_{ijk}^{\left(2\right)}&z_{ijk}^{\left(3\right)}\end{array}\right]^{-1}\left[\begin{array}[]{l}z_{ij}^{\left(1\right)}x_{ij}+z_{ij}^{\left(2\right)}y_{ij}+z_{ij}^{\left(3\right)}z_{ij}\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
z_{ik}^{\left(1\right)}x_{ik}+z_{ik}^{\left(2\right)}y_{ik}+z_{ik}^{\left(3\right)}z_{ik}\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
z_{ijk}^{\left(1\right)}x_{i}+z_{ijk}^{\left(2\right)}y_{i}+z_{ijk}^{\left(3\right)}z_{i}\end{array}\right].$
(21)

The local coordinates of point
$\bm{x}_{ijk}^{*}=\left(x_{ijk}^{*},y_{ijk}^{*},z_{ijk}^{*}\right)$ in the
plane of circle $C_{ij}$ are:

$\left\\{\begin{array}[]{l}x_{ijk}^{*}=\left(x_{ijk}-x_{ij},y_{ijk}-y_{ij},z_{ijk}-z_{ij}\right)\cdot\bm{\hat{x}}_{ij},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
y_{ijk}^{*}=\left(x_{ijk}-x_{ij},y_{ijk}-y_{ij},z_{ijk}-z_{ij}\right)\cdot\bm{\hat{y}}_{ij},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ z_{ijk}^{*}=0.\end{array}\right.$ (22)

Figure 4: Finding the angle interval $A_{ijk}$ of the intersection of circles
$C_{ij}$ and $C_{ik}$.

Consider the plane of circle $C_{ij}$, as shown in Fig. 4. The angles $\beta$
and $\gamma$ are:

$\beta=\cos^{-1}\left(\frac{d_{ijk}^{*}}{r_{ij}}\right),$ (23)
$\gamma=\left\\{\begin{array}[]{ll}\cos^{-1}\left(\displaystyle\frac{x_{ijk}^{*}}{d_{ijk}^{*}}\right),&\;\text{if}\;\;y_{ijk}^{*}\geq
0\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
2\pi-\cos^{-1}\left(\displaystyle\frac{x_{ijk}^{*}}{d_{ijk}^{*}}\right),&\;\text{if}\;\;y_{ijk}^{*}<0\end{array}\right.$
(24)

where $d_{ijk}^{*}=\sqrt{x_{ijk}^{*}\cdot x_{ijk}^{*}+y_{ijk}^{*}\cdot
y_{ijk}^{*}}$. The angle interval $A_{ijk}$ of the intersection of circles
$C_{ij}$ and $C_{ik}$ is:

$A_{ijk}=\left[\alpha_{1},\alpha_{2}\right]=\left[\gamma-\beta,\gamma+\beta\right].$
(25)

####  3.1.3 Checking coverage of circle

Figure 5: Example of the dynamic list update procedure.

After obtaining the angle intervals $A_{ijk}$ of the intersection between
circles $C_{ij}$ and $C_{ik}$, the checking coverage of a circle is
accomplished using the method developed by Haque et al. [17]. This method
employs a dynamic list that represents the parts of the circle $C_{ij}$ that
have not been intersected by the neighboring circles $C_{ik}$. When the
dynamic list is empty, it can be concluded that the circle $C_{ij}$ is
entirely covered by the neighboring circles $C_{ik}$ and the calculation of
the angle intervals $A_{ijk}$ can be terminated. Figure 5 shows an example of
the dynamic list update procedure, in which the dynamic list
$\mathcal{H}_{ij}$ is initialized to the entire interval
$\mathcal{B}=\left[0,2\pi\right)$ and represents the complement of the angle
interval $\tilde{A}_{ijk}$ in the current list $\mathcal{H}_{ij}$ in the
following calculations, where it is expressed as:

$\mathcal{H}_{ij}=\mathcal{H}_{ij}-\tilde{A}_{ijk},$ (26)

where $\tilde{A}_{ijk}$ is expressed as follows:

$\tilde{A}_{ijk}=\left\\{\begin{array}[]{lll}\left[0,\alpha_{2}\right]\cup\left[\alpha_{1}+2\pi,2\pi\right)&\text{if}&\,\alpha_{1}<0\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\left[0,\alpha_{2}-2\pi\right]\cup\left[\alpha_{1},2\pi\right)&\text{if}&\,\alpha_{2}\geq
2\pi\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
A_{ijk}&\text{else}&\end{array}\right.$ (27)

When the linked list is empty, the circle $C_{ij}$ is completely covered by
the angle intervals $A_{ijk}$:

$\mathcal{H}_{ij}=\oslash\Longleftrightarrow C_{ij}\;\text{is covered}$ (28)

Because the intersection particles are sorted from smallest to largest
$\varphi_{ijk}$, the corresponding angle intervals $A_{ijk}$ are obtained from
largest to smallest. In fact, for most interior circles $C_{ij}$, only a few
angle intervals $A_{ijk}$ are required to determine whether the circle
$C_{ij}$ is covered.

####  3.1.4 Analysis of the geometrical method

The performance of the proposed geometrical method is determined by the radius
of the sphere and the rotation angle of the cone. If the radius of the sphere
is too large, the small hole inside the fluid will not be detected, and a huge
detection time will be consumed. If it is too small, several inner particles
will be misidentified as free-surface particles. In this study, the radius of
the sphere is set to $2.5\Delta x$. If the rotation angle of the cone is too
large, the detection accuracy of the concave surface will decrease, and a huge
detection time will be consumed. If it is too small, several inner particles
will be misidentified as a free-surface particles. In this study, the rotation
angle of the cone is set to $45^{\circ}$.

###  3.2 Semi-geometrical method for free-surface particle detection

The geometrical method proposed in Section 3.1 can accomplish free-surface
particle detection of the three-dimensional model using only the particle
coordinates. However, if the proposed geometrical method is directly employed
in a numerical simulation, it will result in high computation costs because it
is necessary for each particle to perform continuous global scanning through a
cone region at each time step. As a result, a semi-geometrical method is
proposed to reduce computational costs based on the proposed geometrical
method. The semi-geometrical method for detecting free-surface particles
comprise the following two steps.

In the first step, the particles near the free surface are found by position
divergence [20], which considerably reduces the number of particles to be
accurately detected in the second step. The position divergence $\lambda_{i}$
is calculated as follows:

$\lambda_{i}=\nabla\cdot\bm{x}_{i}=\sum_{j=1}^{N}{\frac{m_{j}}{\rho_{j}}\left(\bm{x}_{j}-\bm{x}_{i}\right)\cdot\nabla_{i}W_{ij}}.$
(29)

Figure 6 demonstrates the position divergence of particles in three
dimensions. The position divergence of the free-surface particles is 1.874
when the particles are uniformly distributed. As a result, particles with a
position divergence of less than 1.874 are directly identified as free-surface
particles. In the proposed semi-geometrical method, we selected 1.8 as the
threshold for detecting free-surface particles. The position divergence of the
particles in the two layers near the free surface ranges from 1.874 to 2.571.
Therefore, particles with a position divergence exceeding 2.571 are directly
identified as inner particles. We selected 2.8 as the threshold for detecting
inner particles in the semi-geometrical method. Finally, for particles with a
position divergence greater than 1.8 but less than 2.8, an accurate detection
of the free-surface particle is executed in the second step. Note that the
thresholds vary with respect to the ratio $h/\Delta x$ and kernel function.

Figure 6: The position divergences of particles in three dimensions.

The second step involves accurate detection executed on free-surface candidate
particles with a position divergence greater than 1.8 but less than 2.8. In
contrast to the method used in the study [18, 20], i.e., scanning the
“umbrella-shaped” zone along the normal direction of the particle, continuous
global scanning within the sphere $S_{i}$ is accomplished via a cone region
$C$. If there exists a cone region with no neighboring particles, particle $i$
is identified as a free-surface particle, otherwise, it is identified as an
inner particle. Finally, the conditions for determining whether particle $i$
is identified as a free-surface particle are given by:

$\left\\{\begin{array}[]{lll}\lambda_{i}<1.8&\Longrightarrow&i\in F\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ 1.8\leq\lambda_{i}\leq
2.8\;\&\&\;\exists\;C=\varnothing&\Longrightarrow&i\in F\vskip 6.0pt plus
2.0pt minus 2.0pt\\\ \text{otherwise}&\Longrightarrow&i\notin
F\end{array}\right.$ (30)

where $C$ denotes the scanning cone, and $F$ denotes the free-surface zone.
The performance of the proposed semi-geometrical method is tested in section
4.1 and section 4.2.

##  4 Results and discussion

Four tests are selected in this section to illustrate the present method for
3D free-surface particle detection. The advantages of the present method are
analyzed in detail in the first two tests, i.e., surface detection of the free
surface with periodic perturbations and dam-breaking, where the detection
results and detection time are compared with those of the Haque method [17],
the Marrone method [18], and the Wang method [20]. Then, the present method is
applied to the visualization tools and the particle shifting technology.

In the Haque method, the radius of the sphere is set to $1.25\Delta x$,
therefore, the spheres of neighboring particles within $2.5\Delta x$ from the
candidate particle intersect the sphere of the candidate particle, and the
number of the neighboring particles is 81, which is the same as the present
method. In the Marrone method, the renormalized Gaussian kernel function is
used, and the smoothing length is set to $1.33\Delta x$ [18]. The QR
decomposition method is used to solve the minimum eigenvalue of the
renormalization matrix $\mathbf{L}^{-1}$. The rotation angle of the “umbrella-
shaped” zone in the second step is $45^{\circ}$. In the Wang method, the
improved Gaussian kernel function is adopted, and the smoothing length is set
to $1.5\Delta x$ [20]. The rotation angle of the “umbrella-shaped” zone in the
second step is $45^{\circ}$. In the present method, the radius of the sphere
is set to $2.5\Delta x$, and the rotation angle of the cone is set to
$45^{\circ}$.

###  4.1 Surface detection of the free surface with periodic perturbations

In this section, the present method is applied to the free surface with
periodic perturbations and the detection results are compared with those of
the Haque method [17], the Marrone method [18], and the Wang method [20]. The
distance between the free-surface particles is maintained as the initial
inter-particle distance at varied angles of surface perturbation in addition
to a uniform distribution of inner particles. The 3D detection methods are
transformed into the corresponding 2D detection methods and applied to the 2D
free surface with periodic perturbations to facilitate the analysis and
presentation of the detection results. The Haque 3D method corresponds to the
Dilts 2D method [13], i.e., if the circle of a candidate particle is not
completely covered by the circles of its neighboring particles, the candidate
particle is identified as a free-surface particle. The Marrone method is
applicable to both 2D and 3D cases, i.e., the particles near the free surface
are initially found by the minimum eigenvalue of the renormalization matrix
$\mathbf{L}^{-1}$, then, the free-surface particles are accurately detected by
a scan region along the normal direction of the candidate particle. The Wang
method is applicable to both 2D and 3D cases, i.e., the particles near the
free surface are initially found by the position divergence, then, the free-
surface particles are accurately detected by a scan region along the normal
direction of the candidate particle. The present 3D method corresponds to the
Sun 2D method [12], i.e., continuous global scanning within the circle of a
candidate particle is performed through a sector region. If there exists a
sector region with no neighboring particles, the candidate particle is
identified as a free-surface particle.

Figure 7: Detection results of the Dilts 2D (Haque 3D) method and the Sun 2D
(Present 3D) method applied to the 2D free surface with periodic
perturbations, where the angles of the surface perturbations are
$175^{\circ}$, $130^{\circ}$, and $85^{\circ}$. Figure 8: Detection results of
the Marrone 2D (3D) method, the Wang 2D (3D) method, and the Sun 2D (Present
3D) method applied to the 2D free surface with periodic perturbations for
different inner particle distributions, where the angles of the surface
perturbations are $100^{\circ}$.

Figure 7 shows the detection results of the Dilts method and the Sun method
applied to the 2D free surface with periodic perturbations, where the angles
of the surface perturbations are $175^{\circ}$, $130^{\circ}$, and
$85^{\circ}$. In the Dilts method, a portion of the candidate circle is
exposed, thus, the candidate particle $i$ at the trough is identified as a
free-surface particle when the angle of the surface perturbations is
$175^{\circ}$, however, because the entire arc of the candidate circle is not
exposed, the candidate particle $i$ at the trough is identified as an inner
particle when the angles of the surface perturbations are $130^{\circ}$ and
$85^{\circ}$. In uniform particle distributions, the minimum angle
$\theta_{\min}$ that can be detected by the Dilts method is

$\theta_{\min}=2\arccos{(\frac{\Delta x}{2R})},$ (31)

where $R$ is the radius of the circle. When the radius of the circle is
$1.25\Delta x$, the minimum angle $\theta_{\min}$ is $133^{\circ}$, where
$\theta_{\min}$ increases as the radius of the circle increases. In the Sun
method, the angle of the scanning sector is set to $90^{\circ}$. When the
angles of the surface perturbations are $175^{\circ}$ and $130^{\circ}$ which
are greater than $90^{\circ}$, because continuous global scanning inside the
circle of the candidate particle is performed through a sector region, there
must exist a scanning sector region of particle $i$ with no neighboring
particles. Therefore, the candidate particle $i$ at the trough is identified
as a free-surface particle. When the angle of the surface perturbations is
$85^{\circ}$ that is less than $90^{\circ}$, there is at least one neighboring
particle in any scanning sector region of particle $i$. Therefore, the
candidate particle $i$ at the trough is identified as an inner particle. To
summarize, compared with the Dilts 2D method, the detection accuracy of the
Sun 2D method is improved because of the enhanced capability to detect concave
curves. Similarly, compared with the Haque 3D method, the detection accuracy
of the present 3D method is improved because of the enhanced capability to
detect concave surfaces.

Figure 8 depicts the detection results of the Marrone method, the Wang method,
and the Sun method applied to the 2D free surface with periodic perturbations
at various inner particle distributions, where the angles of the surface
perturbations are $100^{\circ}$. In the Marrone method and the Wang method,
because the scan region along the accurate normal direction does not contain
neighboring particles, the candidate particle $i$ is accurately identified as
a free-surface particle when the inner particles are uniformly distributed.
However, when maintaining a constant distribution of free-surface particles,
an inaccurate estimation of the normal vector is obtained when the inner
particles are lost or aggregated, and the scan region along the inaccurate
normal direction contains a neighboring particle $j$ such that the candidate
particle $i$ at the trough is misidentified as an inner particle. In the Sun
method, the angle of the scanning sector is set to $90^{\circ}$. When the
angle of the surface perturbation is $100^{\circ}$ that is greater than
$90^{\circ}$, because continuous global scanning inside the circle of the
candidate particle is performed through a sector region, there must exist a
scanning sector region of particle $i$ with no neighboring particles inside
it. Therefore, the candidate particle $i$ at the trough is identified as a
free-surface particle at various inner particle distributions. To summarize,
compared with the Marrone 2D method and the Wang 2D method, the detection
accuracy of the Sun 2D method is improved because the detection accuracy is
independent of the estimation of the normal vector of particles. Similarly,
compared with the Marrone 3D method and the Wang 3D method, the detection
accuracy of the present 3D method is improved because the detection accuracy
is independent of the estimation of the normal vector of particles.

###  4.2 Surface detection of the dam-breaking

In this section, the present method is applied to the dam-breaking, and the
detection results and detection time are compared with those of the Haque
method [17], the Marrone method [18], and the Wang method [20]. Because the
detection capability for concave surfaces is improved and the detection
accuracy is independent of the estimation of the normal vector of the
particles, the detection accuracy of the present method is improved over the
Haque method, the Marrone method, and the Wang method, and the free-surface
particles are well detected using the present method. The free-surface
particles detected by the present method are used as the reference free-
surface particles in the subsequent comparison study. Figure 9 shows the
configuration of the 3D dam-breaking simulation. The water tank is 3.2 m long,
0.3 m wide, and 1.8 m high, and the water column is 1.2 m long, 0.3 m wide,
and 0.6 m high. The initial inter-particle distance is 0.01 m. Ghost particle
technology [26] is applied to impose boundary conditions. Furthermore, we
record the total time of the first ten time steps and obtain the mean time of
each method. All free-surface particle detections are performed on a computer
with an i5-3470 processor at 3.20 GHz.

Figure 9: Configuration of the dam-breaking. Left: $x-z$ direction view. Right: $y-z$ direction view. Figure 10: Three-dimensional dam-breaking problem. The free-surface particles are detected by the Haque method [17], the Marrone method [18], the Wang method [20], and the present method. Red, blue, and green particles denote free-surface particles, inner particles, and misidentified particles, respectively. Table 1: Quantitative comparisons of the detection results of the three-dimensional dam-breaking problem. The first number denotes the number $TP$ of reference free-surface particles, the second number denotes the number $EP$ of misidentified particles, and the third number denotes the accuracy $M$ of the method. |  | Haque method
---
| Marrone method
---
| Wang method
---
| Present method
---
| 0.600 s
---
10440/3120/0.701 |  10440/1650/0.842 |  10440/1650/0.842 |  10440/0/1.000
| 1.450 s
---
13380/4510/0.663 |  13380/3780/0.717 |  13380/3690/0.724 |  13380/0/1.000
| 1.800 s
---
19380/5400/0.721 |  19380/4260/0.780 |  19380/4530/0.766 |  19380/0/1.000

Figure 10 shows the free-surface particle detection using the present method
in dam-breaking, and the detection results are compared with those of the
Haque method, the Marrone method, and the Wang method. Red, blue, and green
particles denote free-surface particles, inner particles, and misidentified
particles, respectively. The results show that the free-surface particles are
well detected using the present method even when the particles are highly
disordered in distribution. Several free-surface particles are misidentified
as the inner particles owing to the reduced accuracy of the Haque method in
detecting concave surfaces. Because the free-surface particles must be
accurately detected by scanning neighboring particles along the normal
direction in the second step of the Marrone method and the Wang method,
several free-surface particles are misidentified as the inner particles owing
to the inaccurate estimation of the normal vector in complex flows.
Furthermore, quantitative comparisons of the detection results of the
different methods are performed, as shown in Tables 1. We measure the accuracy
$M$ of each method using the following expression:

$M=\frac{TP-EP}{TP}$ (32)

where $TP$ denotes the number of reference free-surface particles, and $EP$
denotes the number of misidentified particles. As time goes by, the free-
surface flows are gradually in a state of large deformation, and the detection
accuracy of the Haque method, the Marrone method, and the Wang method all
occur a significant decrease to just over seventy percent. Figure 11 shows the
particle distributions of single particles in the right-hand fixed area of the
mid-longitudinal section at 1.45 s and 1.80 s, and the free-surface particles
are detected by different methods. As illustrated in Fig. 11, the accuracy of
the Haque method is reduced for detecting concave surfaces, thus, several
free-surface particles on concave surfaces are misidentified as the inner
particles. Because the free-surface particles must be accurately detected by
scanning neighboring particles along the normal direction in the second step
of the Marrone method and the Wang method, several free-surface particles have
at least one neighboring particle in the scan region and are misidentified as
the inner particles owing to the inaccurate estimation of the normal vector in
complex flows. Although the Marrone method and the Wang method have the same
procedure in the second step, i.e., the free-surface particles are further
detected by a scan region along the normal direction of the candidate
particle, the difference in the way the free-surface candidate particles are
identified in the first step leads to a difference in the detection results in
regions where only a small number of particles are present. Because the
detection capability for concave surfaces influences the detection accuracy of
free-surface particles, the effect of the rotation angle of the cone on the
detection results in the present method is analyzed, as shown in Fig. 12. When
the angle of the cone is $66.5^{\circ}$, which is half of the $\theta_{\min}$
in Eq. (31), the detection results of free-surface particles are essentially
the same as those of the Haque method. Owing to the limited detection
capability for concave surfaces, only a small number of free-surface particles
are detected. When the rotation angle is reduced to $56^{\circ}$, the number
of misidentified free-surface particles decreases due to the higher detection
capability for concave surfaces. When the angle of the cone is $45^{\circ}$,
the free-surface particles are well detected using the present method. It is
critical to note that the rotation angle of the cone cannot be too small,
otherwise, several inner particles may be misidentified as the free-surface
particles.

Figure 11: Particle distributions of single particles in the right-hand fixed
area of the mid-longitudinal section at 1.45 s and 1.80 s, and the free-
surface particles are detected by the Haque method [17], the Marrone method
[18], the Wang method [20], and the present method. Red, blue, and green
particles denote free-surface particles, inner particles, and misidentified
particles, respectively. Because the detection accuracy of the present method
is improved over the Haque method, the Marrone method, and the Wang method,
the free-surface particles detected by the present method are used as the
reference free-surface particles. Figure 12: Effect of the rotation angle of
the cone on the detection results in the present method. Red, blue, and green
particles denote free-surface particles, inner particles, and misidentified
particles, respectively. The free-surface particles detected by the present
method with a rotation angle of $45^{\circ}$ are used as the reference free-
surface particles. Figure 13: The detection time of the Haque method [17], the
Marrone method [18], the Wang method [20], and the present method with the
number of particles.

Figure 13 gives the detection times of the Haque method, the Marrone method,
the Wang method, and the present method. In the Haque method, because the
circles obtained by the intersection of sphere $S_{i}$ and its neighboring
spheres $S_{j}$ are assigned an uncertain normal direction, certain
computations must be performed even if the circles do not intersect. When the
particles are uniformly distributed, the upper limit of the computational time
of the Haque method is $O(MN^{2})$, where $M$ is the particle number and $N$
is the average number of neighboring particles per particle. Therefore, it
leads to huge computational costs. In the Marrone method, the minimum
eigenvalue of the matrix $\mathbf{L}^{-1}$ is computed for all particles in
the first step, which takes most of the detection time. Meanwhile, the
scanning by the “umbrella-shaped” zone along the normal direction is performed
for a few particles in the second step, which takes only a small portion of
the detection time. In the Wang method, the position divergence of the
particle is computed for all particles in the first step, which takes only a
small amount of time. Meanwhile, the scanning by the “umbrella-shaped” zone
along the normal direction is performed for a few particles in the second
step, which also takes a small amount of time. In the present method, the
position divergence of the particle is computed for all particles in the first
step, which takes only a small portion of the detection time. Meanwhile,
continuous global scanning within the sphere is performed through the cone
region for a few particles in the second step, which takes most of the
detection time. Furthermore, because the circles given by the intersection of
a sphere $S_{i}$ and its neighboring cones $C_{j}$ have the same radius and
outward normal direction in the continuous global scanning, only when the
angle $\varphi_{ijk}$ is greater than $0^{\circ}$ and less than $2\theta$,
circles $C_{ij}$ and $C_{ik}$ intersect at two points, and the point
$\bm{x}_{ijk}$ of the intersection of the planes of the circles $C_{ij}$,
$C_{ik}$ and $C_{ijk}$ is unique and can be calculated directly. When the
particles are uniformly distributed, because the angle $\theta$ is set to
$45^{\circ}$, only half of the neighboring circles $C_{ik}$ intersect the
circle $C_{ij}$, and the upper limit of the computational time of the present
geometrical method is $O(0.5MN^{2})$. Therefore, the detection time of a
single particle in the second step of the present method is significantly
reduced compared with the Haque method. When the particle number is small, the
detection time of the Marrone method is shorter than that of the present
method. When the particle number increases, the detection time of the Marrone
method increases approximately linearly with respect to the number of
particles. However, in the present method, the proportion of particles
detected in the second step compared to the total number of particles
gradually decreases with the number of particles, hence, the increasing rate
of the detection time for the present method gradually decreases. When the
particle number is one million, the detection time is 139.79 s for the Haque
method, 4.04 s for the Marrone method, 0.65 s for the Wang method, and 3.34 s
for the present method. Compared with the Haque method and the Marrone method,
the detection time of the present method is reduced by 97.61% and 17.33%,
respectively. Compared with the Wang method, the detection time of the present
method is increased by 413.85%.

Furthermore, Figures 14 and 15 depict the free-surface particle detection
using the present method in complex flows resulting from the impact of a
double dam-breaking and the impact of a single dam-breaking against a rigid
tall obstacle, respectively. In each case, the detection results are compared
with those of the Haque method, the Marrone method, and the Wang method. The
results show that the free-surface particles are well detected using the
present method even when the particles are highly disordered in distribution.
Quantitative comparisons of the detection results of the different methods are
also performed, as shown in Tables 2 and 3. At the initial moment, all free-
surface particles are accurately detected by the different methods because of
the uniform distribution of particles. However, as time goes by, the free-
surface flows are gradually in a state of large deformation, and the detection
accuracy of the Haque method, the Marrone method, and the Wang method all
occur a significant decrease to just over sixty percent.

Figure 14: Three-dimensional dam-breaking problem of two liquid columns at opposite corners of a container with a square base of size 2.4 and height 3.5. The free-surface particles are detected by the Haque method [17], the Marrone method [18], the Wang method [20], and the present method. Red, blue, and green particles denote free-surface particles, inner particles, and misidentified particles, respectively. Figure 15: Impact against a rigid obstacle after a dam-break. The free-surface particles are detected by the Haque method [17], the Marrone method [18], the Wang method [20], and the present method. Red, blue, and green particles denote free-surface particles, inner particles, and misidentified particles, respectively. Table 2: Quantitative comparisons of the detection results of the three-dimensional dam-breaking problem of two liquid columns at opposite corners of a container. The first number denotes the number $TP$ of reference free-surface particles, the second number denotes the number $EP$ of misidentified particles, and the third number denotes the accuracy $M$ of the method. |  | Haque method
---
| Marrone method
---
| Wang method
---
| Present method
---
| 0.000 s
---
16882/0/1.000 |  16882/0/1.000 |  16882/0/1.000 |  16882/0/1.000
| 0.400 s
---
17274/904/0.948 |  17274/712/0.959 |  17274/657/0.962 |  17274/0/1.000
| 0.800 s
---
38281/14922/0.610 |  38281/12025/0.686 |  38281/12225/0.681 |  38281/0/1.000
Table 3: Quantitative comparisons of the detection results of the impact against a rigid obstacle after a dam-break. The first number denotes the number $TP$ of reference free-surface particles, the second number denotes the number $EP$ of misidentified particles, and the third number denotes the accuracy $M$ of the method. |  | Haque method
---
| Marrone method
---
| Wang method
---
| Present method
---
| 0.000 s
---
16680/0/1.000 |  16680/0/1.000 |  16680/0/1.000 |  16680/0/1.000
| 0.400 s
---
37647/12207/0.676 |  37647/10057/0.733 |  37647/9869/0.738 |  37647/0/1.000
| 0.700 s
---
71173/24405/0.657 |  71173/26106/0.633 |  71173/23495/0.670 |  71173/0/1.000

###  4.3 Visualization tool

Figure 16: Three-dimensional dam-breaking. (a) 0.52 s; (b) 1.42 s. Top:
particle distributions. Middle: free surface. Bottom: pressure contours
interpolated on the mid-longitudinal section.

The SPH method always outputs computation results in the form of discrete
points, making it difficult to obtain slices, iso-surfaces, contour maps, and
streamlines. The aforementioned analysis, however, can be conducted directly
using existing tools if the particle data can be interpolated from the
neighboring regular meshes. Marrone et al. [18] developed the following method
for analyzing free-surface flows in detail.

(1) Detect free-surface particles.

(2) Construct regular meshes with spatial resolution $dx$ containing all
computational domains. For each mesh $M$, search for the nearest free-surface
particle $F_{M}$, and obtain $d_{MF_{M}}$ as:

$d_{MF_{M}}=\left(\bm{x}_{F_{M}}-\bm{x}_{M}\right)\cdot\bm{n}_{F_{M}},$ (33)

where $\bm{n}_{F_{M}}$ is the normal vector of particle $F_{M}$. Then, a
level-set function $\phi\left(\bm{x}_{M}\right)$ on each mesh can be defined
as follows:

$\phi\left(\bm{x}_{M}\right)=\left\\{\begin{array}[]{ll}-1,&d_{MF_{M}}\leqslant-2h\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ d_{MF_{M}}/2h,&-2h<d_{MF_{M}}<2h\vskip 6.0pt
plus 2.0pt minus 2.0pt\\\ 1,&d_{MF_{M}}\geqslant 2h\end{array}\right.$ (34)

The level-set function is negative outside the fluid, zero on the free surface
where $\phi=-dx/4h$, and positive inside the fluid.

(3) For each mesh $M$, the physical quantity on the mesh can be interpolated
from the neighboring scattered point data using the moving least squares (MLS)
method [27].

Figure 16 depicts the 3D dam-breaking at 0.52 s and 1.42 s, and the free-
surface particles are detected by the present method. The top shows particle
distributions, the middle shows the free surface, and the bottom shows
pressure contours on the mid-longitudinal section. It can be seen that the
visualization tools can clearly describe the surface position and internal
flow.

###  4.4 Particle shifting technology

Figure 17: Sketch of the PST. For the inner particle $a$, all neighboring
particles participate in the calculation of the shifting vector. For the free-
surface vicinity particle $b$, only the neighboring particles whose distance
from the particle is less than $l_{b}$ participate in the calculation of the
shifting vector. For the free-surface particle $c$, the shifting distance is
zero.

PST has become a popular technology in recent years for improving the
robustness and accuracy of the SPH method [28, 29, 9, 20, 30]. Xu et al. [28]
first developed the PST. In this method, the particles are shifted slightly
across streamlines to avoid the extreme bunching and stretching of particles.
When the PST is applied to free-surface flows, special treatment of the free-
surface particles is required to maintain the kinematic boundary condition [9,
10, 20, 30]. Wang et al. [20] proposed an improved PST that enables accurate
shifting of particles near the free surface. This method is depicted in Fig.
17 and is implemented as follows:

(1) Detect free-surface particles.

(2) For free-surface vicinity particles, record the distance $l$ between them
and the nearest free-surface particle.

(3) For inner particles, all neighboring particles participate in the
calculation of the shifting vector. For free-surface vicinity particles, only
neighboring particles whose distance from the particle is less than $l$
participate in the calculation of the shifting vector. For free-surface
particles, the shifting distance is zero.

Figure 18 shows the visualization results of the single-layer particles on the
mid-longitudinal sector at 2.24 s and 2.97 s, and the free-surface particles
in PST are detected by the Wang method and the present method, respectively.
The results show that the diffusion of the free surface obtained by applying
the Wang method to PST is observed compared with the results obtained by the
present method. Because the free-surface particles must be accurately detected
by scanning neighboring particles along the normal direction in the second
step of the Wang method, several free-surface particles are misidentified as
the inner particles owing to the inaccurate estimation of the normal vector in
complex flows. This part of misidentified particles obtains a PST displacement
towards the free surface, therefore, the number of particles in the free-
surface region gradually increases and the free surface diffusion occurs.
Because the detection accuracy is independent of the estimation of the normal
vector of particles, the free-surface particles are accurately detected by the
present method, and the PST displacement of the free-surface particles is
zero, therefore, there is no particle aggregation in the free-surface region
and no diffusion on the free surface.

Figure 18: Visualization results of the single-layer particles on the mid-
longitudinal sector at 2.24 s and 2.97 s, and the free-surface particles in
PST are detected by the Wang method [20] and the present method, respectively.
The red line denotes the free surface obtained by applying the Wang method to
PST, and the black line denotes the free surface obtained by applying the
present method to PST. The arrow indicates the direction of free surface
diffusion.

##  5 Conclusions

A high-accuracy three-dimensional surface detection algorithm is studied for
the SPH method. First, a geometrical method for free-surface particle
detection is developed to improve detection accuracy in complex flows. This
method detects free-surface particles by performing continuous global scanning
within the sphere of a particle through a cone region. The vertex of the cone
corresponds to the particle position. If there exists a cone region with no
neighboring particles, the particle is identified as a free-surface particle.
Then, based on the proposed geometrical method, a two-step semi-geometrical
method is developed to reduce computational costs. In the first step, the
particles near the free surface are found by position divergence. In the
second step, the free-surface particles are accurately detected using the
proposed geometrical method.

The accuracy and robustness of the proposed method are demonstrated via four
tests. The test demonstrating the surface detection of the free surface with
periodic perturbations shows that the detection accuracy of the proposed
method is significantly improved compared with the traditional methods because
the detection capability for concave surfaces is enhanced and the detection
accuracy is independent of the estimation of the normal vector of particles.
The surface detection of 3D dam-breaking tests shows that the proposed method
can well detect the free-surface particles in complex flows but with no
increase in computational costs compared with the traditional methods. The
results of the tests on visualization tools and the particle shifting
technology show that the proposed method can be applied to various fields with
satisfactory results and has a wide range of application prospects. In
conclusion, all tests demonstrate that the proposed method is effective and
robust for complex flows.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China
under Grant Nos. U1530261, 11772065, and the Science Challenge Project, China
(Grant No. TZ2016001).

## Reference

## References

  * [1] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3) (1977) 375–389.  doi:https://doi.org/10.1093/mnras/181.3.375.
  * [2] J. J. Monaghan, Smoothed particle hydrodynamics and its diverse applications, Annu. Rev. Fluid. Mech. 44 (2012) 323–346.  doi:https://doi.org/10.1146/annurev-fluid-120710-101220.
  * [3] A. M. Zhang, P. N. Sun, F. R. Ming, A. Colagrossi, Smoothed particle hydrodynamics and its applications in fluid-structure interactions, J. Hydrodyn. 29 (2017) 187–216.  doi:https://doi.org/10.1016/S1001-6058(16)60730-8.
  * [4] T. Ye, D. Y. Pan, C. Huang, M. B. Liu, Smoothed particle hydrodynamics (SPH) for complex fluid flows: Recent developments in methodology and applications, Phys. Fluids. 31 (2019) 011301.  doi:https://doi.org/10.1063/1.5068697.
  * [5] Z. Zhang, M. Liu, A decoupled finite particle method for modeling incompressible flows with free surfaces, Appl. Math. Model. 60 (2018) 606–633.  doi:https://doi.org/10.1016/j.apm.2018.03.043.
  * [6] D. R. Vyas, S. J. Cummins, M. Rudman, P. W. Cleary, G. W. Delaney, D. V. Khakhar, Collisional SPH: A method to model frictional collisions with SPH, Appl. Math. Model. 94 (2021) 13–35.  doi:https://doi.org/10.1016/j.apm.2021.01.005.
  * [7] A. Colagrossi, M. Antuono, D. L. Touzé, Theoretical considerations on the free-surface role in the smoothed-particle-hydrodynamics model, Phys. Rev. E 79 (5) (2009) 056701.  doi:https://doi.org/10.1103/PhysRevE.79.056701.
  * [8] A. Colagrossi, M. Antuono, A. Souto-Iglesias, D. Antonio, L. Touzé, Theoretical analysis and numerical verification of the consistency of viscous smoothed-particle-hydrodynamics formulations in simulating free-surface flows, Phys. Rev. E 84 (2) (2011) 026705.  doi:https://doi.org/10.1103/PhysRevE.84.026705.
  * [9] P. N. Sun, A. Colagrossi, S. Marrone, A. M. Zhang, The $\rm{\delta}$plus-SPH model: Simple procedures for a further improvement of the SPH scheme, Comput. Methods Appl. Mech. Eng. 315 (2017) 25–49.  doi:https://doi.org/10.1016/j.cma.2016.10.028.
  * [10] C. Huang, T. Long, S. M. Li, M. B. Liu, A kernel gradient-free SPH method with iterative particle shifting technology for modeling low-reynolds flows around airfoils, Eng. Anal. Bound. Elem. 106 (2019) 571–587.  doi:https://doi.org/10.1016/j.enganabound.2019.06.010.
  * [11] M. Y. Zhang, Simulation of surface tension in 2D and 3D with smoothed particle hydrodynamics method, J. Comput. Phys. 229 (19) (2010) 7238–7259.  doi:https://doi.org/10.1016/j.jcp.2010.06.010.
  * [12] C. Sun, Z. J. Shen, M. Y. Zhang, Surface treatment technique of MPS method for free surface flows, Eng. Anal. Bound. Elem. 102 (2019) 60–72.  doi:https://doi.org/10.1016/j.enganabound.2019.02.004.
  * [13] G. A. Dilts, Moving least-squares particle hydrodynamics ii: conservation and boundaries, Int. J. Numer. Methods Eng. 48 (10) (2000) 1503–1524.  doi:https://doi.org/10.1002/1097-0207(20000810)48:10<1503::AID-NME832>3.0.CO;2-D.
  * [14] J. P. Morris, Simulating surface tension with smoothed particle hydrodynamics, Int. J. Numer. Methods Fluids. 33 (3) (2000) 333–353.  doi:https://doi.org/10.1002/1097-0363(20000615)33:3<333::AID-FLD11>3.0.CO;2-7.
  * [15] P. W. Randles, L. D. Libersky, Smoothed particle hydrodynamics: Some recent improvements and applications, Comput. Methods Appl. Mech. Eng. 139 (1) (1996) 375–408.  doi:https://doi.org/10.1016/S0045-7825(96)01090-0.
  * [16] Y. X. Lin, G. R. Liu, G. Y. Wang, A particle-based free surface detection method and its application to the surface tension effects simulation in smoothed particle hydrodynamics SPH, J. Comput. Phys. 383 (2019) 196–206.  doi:https://doi.org/10.1016/j.jcp.2018.12.036.
  * [17] A. Haque, G. A. Dilts, Three-dimensional boundary detection for particle methods, J. Comput. Phys. 226 (2) (2007) 1710–1730.  doi:https://doi.org/10.1016/j.jcp.2007.06.012.
  * [18] S. Marrone, A. Colagrossi, D. L. Touz$\rm{\acute{e}}$, G. Graziani, Fast free-surface detection and level-set function definition in SPH solvers, J. Comput. Phys. 229 (10) (2010) 3652–3663.  doi:https://doi.org/10.1016/j.jcp.2010.01.019.
  * [19] C. Huang, D. Zhang, Y. Si, Y. Shi, Y. Lin, Coupled finite particle method for simulations of wave and structure interaction, Coast. Eng. 140 (2018) 147–160.  doi:https://doi.org/10.1016/j.coastaleng.2018.07.003.
  * [20] P. P. Wang, Z. F. Meng, A. M. Zhang, F. R. Ming, P. N. Sun, Improved particle shifting technology and optimized free-surface detection method for free-surface flows in smoothed particle hydrodynamics, Comput. Methods Appl. Mech. Eng. 357 (2019) 112580.  doi:https://doi.org/10.1016/j.cma.2019.112580.
  * [21] M. Sandim, A. Paiva, L. H. de Figueiredo, Simple and reliable boundary detection for meshfree particle methods using interval analysis, J. Comput. Phys. 420 (2020) 109702.  doi:https://doi.org/10.1016/j.jcp.2020.109702.
  * [22] B. Zheng, L. Sun, P. Yu, A novel interface method for two-dimensional multiphase sph: Interface detection and surface tension formulation, J. Comput. Phys. 431 (2021) 110119.  doi:https://doi.org/10.1016/j.jcp.2021.110119.
  * [23] S. Marrone, A. Colagrossi, A. D. Mascio, D. L. Touz$\rm{\acute{e}}$, Prediction of energy losses in water impacts using incompressible and weakly compressible models, J. Fluids Struct. 54 (2015) 802–822.  doi:https://doi.org/10.1016/j.jfluidstructs.2015.01.014.
  * [24] M. Antuono, A. Colagrossi, S. Marrone, D. Molteni, Free-surface flows solved by means of SPH schemes with numerical diffusive terms, Comput. Phys. Comm. 181 (3) (2010) 532–549.  doi:https://doi.org/10.1016/j.cpc.2009.11.002.
  * [25] M. Antuono, A. Colagrossi, S. Marrone, Numerical diffusive terms in weakly-compressible SPH schemes, Comput. Phys. Comm. 183 (12) (2012) 2570–2580.  doi:https://doi.org/10.1016/j.cpc.2012.07.006.
  * [26] A. Colagrossi, M. Landrini, Numerical simulation of interfacial flows by smoothed particle hydrodynamics, J. Comput. Phys. 191 (2) (2003) 448–475.  doi:https://doi.org/10.1016/S0021-9991(03)00324-3.
  * [27] T. P. Fries, H. G. Matthies, Classification and overview of meshfree methods, Informatik-Berichte der Technischen Universität Braunschweig 2003-03 (2004).
URL https://publikationsserver.tu-braunschweig.de/receive/dbbs_mods_00001418

  * [28] R. Xu, P. Stansby, D. Laurence, Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach, J. Comput. Phys. 228 (18) (2009) 6703–6725.  doi:https://doi.org/10.1016/j.jcp.2009.05.032.
  * [29] S. J. Lind, R. Xu, P. K. Stansby, B. D. Rogers, Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves, J. Comput. Phys. 231 (4) (2012) 1499–1523.  doi:https://doi.org/10.1016/j.jcp.2011.10.027.
  * [30] H. Lyu, P. Sun, Further enhancement of the particle shifting technique: Towards better volume conservation and particle distribution in SPH simulations of violent free-surface flows, Appl. Math. Model. 101 (2022) 214–238.  doi:https://doi.org/10.1016/j.apm.2021.08.014.


11institutetext: Computer Vision Center, UAB 22institutetext: 22email:
<EMAIL_ADDRESS>
# OCR-IDL: OCR Annotations for Industry Document Library Dataset

Ali Furkan Biten 11 Rubèn Tito 11 Lluis Gomez 11 Ernest Valveny 11 Dimosthenis
Karatzas 11

###### Abstract

Pretraining has proven successful in Document Intelligence tasks where deluge
of documents are used to pretrain the models only later to be finetuned on
downstream tasks. One of the problems of the pretraining approaches is the
inconsistent usage of pretraining data with different OCR engines leading to
incomparable results between models. In other words, it is not obvious whether
the performance gain is coming from diverse usage of amount of data and
distinct OCR engines or from the proposed models. To remedy the problem, we
make public the OCR annotations for IDL documents using commercial OCR engine
given their superior performance over open source OCR models. The contributed
dataset (OCR-IDL) has an estimated monetary value over 20K US$. It is our hope
that OCR-IDL can be a starting point for future works on Document
Intelligence. All of our data and its collection process with the annotations
can be found in https://github.com/furkanbiten/idl_data.

##  1 Introduction

Analysis of masses of scanned documents is essential in intelligence, law,
knowledge management, historical scholarship, and other areas [25]. The
documents are often complex and varied in nature that can be digital or
scanned born, containing elements such as forms, figures, tables, graphics and
photos, while being produced by various printing and handwriting technologies.
Some common examples of documents comprise of purchase orders, financial
reports, business emails, sales agreements, vendor contracts, letters,
invoices, receipts, resumes, and many others [54]. Processing various document
types to user’s intent is done with manual labor that is time-consuming and
expensive, meanwhile requiring manual customization or configuration. In other
words, each type of document demands hard-coded changes when there is a slight
change in the rules or workflows of documents or even when dealing with
multiple formats.

To address these problems, Document Intelligence models and algorithms are
created to automatically structure, classify and extract information from
documents, improving automated document processing. Particularly, Document
Intelligence as a research field aims at creating models for automatically
analyzing and understanding documents, reducing the time and the cost
associated with it. From a research perspective, what makes Document
Intelligence especially challenging is the requirement of combining various
disciplines such as optical character recognition (OCR), document structure
analysis, named entity recognition, information retrieval, authorship
attribution and many more.

Recent methods on Documents Intelligence utilize deep neural networks
combining Computer Vision and Natural Language Processing. Hao et al. [17]
proposed an end-to-end training using Convolutional Neural Networks to detect
tables in documents. Several published works [46, 42, 56] exploit the advances
in object detection [40, 19] to further improve the accuracy in document
layout analysis. Even though these works have advanced the Document
Intelligence field, there are two main limitations to be recognized: (i) they
rely on a small human annotated dataset and (ii) they use pre-trained networks
that have never seen any documents, hence the interaction between text and
layout. Inspired by BERT [11],Xu et al. [54] identified these problems and
propose a pre-training strategy to unlock the potential of large-scale
unlabeled documents. More specifically, they obtain OCR annotations from an
open source OCR engine Tesseract [45] for 5 Million documents from IIT-CDIP
[25] dataset. With the introduction of pre-training strategy and advances in
modern OCR engine [1, 12, 20, 28, 34], many contemporary approaches [7, 2, 53]
have utilized even more data to advance the Document Intelligence field.

In this work, we make public the OCR annotations for 26 Millions pages using a
commercial OCR engine that has the monetary value over 20K US$. Our motivation
for releasing a massive scale documents dataset annotated with a commercial
OCR engine is two-fold. First of all, the usage of different amount of
documents and different OCR engines across the papers makes it impossible to
fairly compare their results and hence their architecture. By creating this
dataset, we hope that the works in Document Intelligence will become more
comparable and have better intuition on what the proposed architecture can
actually accomplish.

Secondly, we decide to use a commercial OCR engine, specifically Amazon
Textract111https://aws.amazon.com/textract/, over Tesseract. It is because the
performance of the OCR engines can significantly affect the model’s
performance which can be seen in fields that use OCR annotations, such as in
fine-grained classification [29, 30, 31], in scene-text visual question
answering [9, 44, 8, 13], in document visual question answering (DocVQA) [50,
33]. Apart from improving the annotation quality significantly, we want to
level the differences between research groups and companies.

(a)

(b)

(c)

(d)

(e)

Figure 1: Document images of OCR-IDL. The dataset includes a wide variety of
documents with dense text (a), tables (b), figures (c), and complex layouts
that combines different elements (d, e).

We provide the annotations for publicly available documents from Industry
Documents Library (IDL). IDL is a digital archive of documents created by
industries which influence public health, hosted by the University of
California, San Francisco Library 222https://www.industrydocuments.ucsf.edu.
IDL has already been used in the literature for building datasets: IIT-CDIP
[25], RVL-CDIP [18], DocVQA [50, 33]. Hence, our OCR annotations can be used
to further advance in these tasks.

The rest of the paper is structured as follows. First, we briefly explain all
the related works. Next, we will elaborate on our data collection and
comparison to other datasets. Finally, we will provide various statistics of
the annotations and conclude our paper.

##  2 Related Work

Document Intelligence can be considered as an umbrella term covering problems
of Key Information Extraction [10, 54], Table Detection [41, 38] and Structure
Recognition [39, 55], Document Layout Segmentation [5, 4] Document Layout
Generation [6, 36, 3, 48], Document Visual Question Answering [51, 50, 32],
Document Image Enhancement [49, 22, 47] which involves the understanding of
visually rich semantic information and structure of different layout entities
of a whole page.

Early days of Document Intelligence has relied on rule-based handcrafted
approaches mainly divided into bottom-up and top-down methods. Bottom-up
methods [15, 24, 35, 43] first detect connected components at the pixel level,
later to be fused into higher level of structure through various heuristics
and name depending on distinct structural features. While top-down methods
[16] dissects a page into smaller units such as titles, text blocks, lines,
and words.

Lately, the success of large-scale pre-training [11] in Natural Language
Processing has been integrated into Document Intelligence, resulting in
impressive performance gains. These methods follow a two step procedure where
first they pretrain the models on unlabeled documents (OCR annotations are
obtained by an off-the-shelf OCR engine), then they finetune it on specific
downstream tasks. LayoutLM [54] is one of the first works that pretrain BERT
based language model with document layout information, using masked
language/vision modeling and multi label classification. BROS [21] is built on
top of Span-BERT [23] with spatially aware graph decoder. For the pretraining
loss, they use area-masked language model. Self-Doc [27] utilizes two separate
transformer [52] encoders for visual and textual features and later to be fed
to multi-modal transformers encoder. TILT [37] tries to encode the layout
information by integrating pairwise 1D and 2D information into their models.
Uni-Doc [14] is designed to do most document understanding tasks that takes
words and visual features from a semantic region of a document image by
combining three self-supervised losses. More recent methods, Doc-Former [2]
and LayoutLMv2 [53] combine multiple pretraining losses such as image-text
alignment, learning to construct image features and multi-modal masked
language modeling together to achieve state-of-the-art results.

Yet, comparing all of these works are cumbersome since each work that performs
pretraining uses different amounts of data with diverse OCR engines. Hence,
this makes it especially hard to understand where the gain is coming from. In
other words, we can not draw clear conclusions to questions such as: “What is
the effect of the amount of pretraining data on the performance?”, “Is the
performance gain coming from a better/stronger OCR engine or from the proposed
architecture?”, “What is the effect of the pretraining loss on the downstream
tasks keeping OCR and the amount of data identical?” To help answer these
questions, we collect and annotate the largest public OCR annotated documents
dataset (OCR-IDL).

##  3 OCR-IDL Dataset

In this section, we elaborate on various details regarding OCR-IDL. Firstly,
we explain the process we follow on how we get the IDL data and use Amazon-
Textract to obtain OCR annotations. Next, we compare OCR-IDL to other datasets
that have proven useful for the document intelligence tasks. And finally, we
provide in-depth statistics on the documents we use.

Dataset # of Docs # of Pages Docs source Docs. description OCR-Text OCR-BB
Layout Doc. type IIT-CDIP [25] $6.5M$* $35.5M$* UCSF-LTD Industry documents
Unknown ✗ ✗ ✓ RVL-CDIP [18] -† $400K$ UCSF-LTD Industry documents ✗ ✗ ✗ ✓
PublayNet [56] -† $364K$ PubMedCentral Journals and articles ✗ ✗ ✓ ✗ DocBank
[26] -† $500K$ arXiv Journals and articles ✗ ✗ ✓ ✗ DocVQA [51] $6K$ $12K$
UCSF-IDL Industry documents Microsoft OCR ✓ ✗ ✓ OCR-IDL $4.6M$ $26M$ UCSF-IDL
Industry documents Amazon Textract ✓ ✗ ✓

Table 1: Summary of other Document Intelligence Datasets. *We skipped $145K$
documents that gave xml parsing errors, didn’t contain document ID or number
of pages. †No traceability between different pages of the same document.

###  3.1 Data Collection

As already mentioned, IDL is an industry documents library hosted by UCSF, its
main purpose is to “identify, collect, curate, preserve, and make freely
accessible internal documents created by industries and their partners which
have an impact on public health, for the benefit and use of researchers,
clinicians, educators, students, policymakers, media, and the general public
at UCSF and
internationally”333https://en.wikipedia.org/wiki/Industry_Documents_Library.
IDL in total contains over 70 millions documents. We use the publicly
available link444https://s3-us-
west-2.amazonaws.com/edu.ucsf.industrydocuments.artifacts/ to downlad the 4.6
Million Documents from IDL dataset which comes in the format of PDFs. Some
examples can be viewed in Figure 1. We appreciate that the documents are quite
varied in terms of layout where there are tables, figures, ads and combination
of them in a single page. Also, we see that the documents are quite text-rich
containing many OCR words. Finally, the pages contain smudges, taints and
other discoloration what is found in the in a real use-case scenario.

Figure 2: Distribution of the annotated documents in terms of document types
and dates.

We choose to annotate only 4.6M documents since annotating 13M documents not
only would be much costlier (42K dollars instead of 18K) but also in the
literature it is shown that using more data have diminishing returns on the
downstream task [7]. After obtaining the data, we preprocess the documents to
remove empty, faulty and broken pdfs. This process resulted in the elimination
of 6548 documents. Moreover, we remove also documents that have more than 2900
pages which are 71 in total. The necessity of removing such huge documents was
because Amazon-Textract OCR only accepts up to 3000 pages. After pre-
processing the documents, we feed all the documents to the OCR engine to
obtain the annotations. The annotations provided by the Textract engine
include transcription of words and lines and their corresponding bounding
boxes and polygons with text type that can be printed or handwritten.
Processing all 4.6M documents was done by a single machine with 16
parallelized cores and took about 1 month.

Figure 3: Distribution of number of pages and number of words per document
type.

###  3.2 Comparison to existing datasets

In this section, we compare the statistics of the amount of documents and
pages to other datasets that are used in Document Intelligence. To name a few,
the Illinois Institute of Technology dataset for Complex Document Information
Processing (IIT-CDIP) [25] is the biggest document dataset and it is designed
for the task of information retrieval. The Ryerson Vision Lab Complex Document
Information Processing (RVL-CDIP) [18] dataset used the IIT-CDIP metadata to
create a new dataset for document classification. PublayNet [56] and DocBank
[26] are datasets designed for layout analysis tasks and DocVQA [33, 51]
instead, is designed for Visual Question Answering task over document images.

We summarize all the key information for comparison in Table 1. First, we
stress that OCR-IDL is the second biggest dataset in amount for pre-training
and biggest dataset with annotations obtained from commercial OCR. This
provides unique opportunity for the Document Intelligence research for
utilizing the unlabeled documents in their research. Furthermore, even though
OCR-IDL uses the documents from the same source as IIT-CDIP and RVL-CDIP, it
also contains other type of industrial documents. OCR-IDL contains documents
from chemical, medical and drug industries, hence having more variety in terms
of content as well as layout information.

Figure 4: Distribution of number of words (left) and lines (right) by pages.

###  3.3 Dataset Statistics

IDL documents come with metadata that is curated by human annotators. They
include information about the date, industry, drugs, chemicals, document types
and many more. We restrict our analysis on exploring what type of documents we
have and the distribution of the dates they are created which can be found in
Figure 2. In IDL metadata, there are 35k various document types from which we
show only the most common 20 which include letters, report, email, memo, note,
etc. As can be seen in the left side of Figure 2, most of the documents’ type
is unknown. Moreover, even though we have a skewed distribution on letter,
report and memo, we also have very distinct distribution such as chart,
graphics, news articles. This is especially encouraging because it provides
diverse layouts within various contexts. Moreover, the documents are created
spanning 100 years where most of the documents are in the range of 1980-2010
as can be seen in the right side of Figure 2. The date the documents are
created contributes to the variability of the documents not only in terms of
semantics but more importantly having visual artifacts (smudges, different
resolution) with different printing technologies.

On top of the amount of documents, we give more details on the number of pages
and words for each document type. The number of pages follows the same
distribution as the amount of documents per type, as can be appreciated in
Figure 3. Also, we can see from Figure 3 that the report type is much richer
in terms of text while the rest follows more or less the same distribution.

Figure 5: Distribution of various layout blocks in OCR-IDL. Our documents
contain a lot of diversity in terms of title, figure, text block, list and
table.

We turn our attention to OCR annotation statistics. In total, we obtain OCR
annotations for 4614232 (4.6M) documents where we have 26621635 (26M) pages,
averaging 6 pages per document. Moreover, since documents are known to be a
text-rich environment, we provide extra details regarding the amount of words
and lines per page and documents. We have 166M words with 46M lines, on
average there are 62.5 words and 17.5 lines per page while 360.8 words and
101.25 lines per document. To have a better understanding of the distribution
of words and lines per page, we present Figure 4. As shown in the figure, mean
distribution is between 20 to 100 words per page while there is a significant
amount of pages that contain more than 200 words per page. The distribution
for lines follows a different distribution in which it can be observed that
most of the pages contain from 10 to 50 lines. In either case, it is clearly
observed that documents at hand are ideal for performing pretraining with
their diverse layouts and text-rich settings.

---
Figure 6: Qualitative results for segmentation of layout information in OCR-
IDL.

Finally, to quantify the diversity of the documents in terms of layout, we run
publicly available Faster-RCNN [40] trained on PubLayNet [56] to segment a
document into text block, title, figure, list, table. To obtain the
segmentation results, we randomly selected 40K pages and some segmentation
examples can be found in Figure 6. It can be appreciated from Figure 5 that
40% of the documents have at least 1 figure. We also observe that 10-20% of
the pages have at least 1 table and list, showing that documents at hand
contain very diverse layout information. Moreover, more than 45% of the pages
contain more than 1 text block and 1 title, making the documents a text rich
environment for pre-training on Document Intelligence tasks.

##  4 Conclusion

In this paper, we have presented our effort to provide OCR annotations for the
large-scale IDL document dataset called OCR-IDL. These annotations have a
monetary value over $20,000 and are made publicly available with the aim of
advancing the Document Intelligence research field. Our motivation is two-
fold, first we make use of a commercial OCR engine to obtain high quality
annotations, leading to reduce the noise provided by OCR on pretraining and
downstream tasks. Secondly, it is our hope that OCR-IDL can be a starting
point for future works on Document Intelligence to be more comparable.
Throughout this article we have detailed the process that we have followed to
obtain the annotations, we have presented a statistical analysis, and compared
them with other datasets in the state of the art. The provided analysis shows
that our contribution has a high potential to be used successfully in pre-
training strategies for document intelligence models. All the code for data
collection process and annotations can be accessed in
https://github.com/furkanbiten/idl_data.

## References

  * [1] Aberdam, A., Litman, R., Tsiper, S., Anschel, O., Slossberg, R., Mazor, S., Manmatha, R., Perona, P.: Sequence-to-sequence contrastive learning for text recognition. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 15302–15312 (2021)
  * [2] Appalaraju, S., Jasani, B., Kota, B.U., Xie, Y., Manmatha, R.: Docformer: End-to-end transformer for document understanding. arXiv preprint arXiv:2106.11539 (2021)
  * [3] Arroyo, D.M., Postels, J., Tombari, F.: Variational transformer networks for layout generation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 13642–13652 (2021)
  * [4] Biswas, S., Banerjee, A., Lladós, J., Pal, U.: Docentr: An end-to-end document image enhancement transformer. arXiv preprint arXiv:2201.11438 (2022)
  * [5] Biswas, S., Riba, P., Lladós, J., Pal, U.: Beyond document object detection: instance-level segmentation of complex layouts. International Journal on Document Analysis and Recognition (IJDAR) 24(3), 269–281 (2021)
  * [6] Biswas, S., Riba, P., Lladós, J., Pal, U.: Docsynth: a layout guided approach for controllable document image synthesis. In: International Conference on Document Analysis and Recognition. pp. 555–568. Springer (2021)
  * [7] Biten, A.F., Litman, R., Xie, Y., Appalaraju, S., Manmatha, R.: Latr: Layout-aware transformer for scene-text vqa. arXiv preprint arXiv:2112.12494 (2021)
  * [8] Biten, A.F., Tito, R., Mafla, A., Gomez, L., Rusinol, M., Mathew, M., Jawahar, C., Valveny, E., Karatzas, D.: Icdar 2019 competition on scene text visual question answering. In: 2019 International Conference on Document Analysis and Recognition (ICDAR). pp. 1563–1570. IEEE (2019)
  * [9] Biten, A.F., Tito, R., Mafla, A., Gomez, L., Rusinol, M., Valveny, E., Jawahar, C., Karatzas, D.: Scene text visual question answering. In: Proceedings of the IEEE/CVF international conference on computer vision. pp. 4291–4301 (2019)
  * [10] Carbonell, M., Riba, P., Villegas, M., Fornés, A., Lladós, J.: Named entity recognition and relation extraction with graph neural networks in semi structured documents. In: 2020 25th International Conference on Pattern Recognition (ICPR). pp. 9622–9627. IEEE (2021)
  * [11] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018)
  * [12] Fang, S., Xie, H., Wang, Y., Mao, Z., Zhang, Y.: Read like humans: Autonomous, bidirectional and iterative language modeling for scene text recognition. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 7098–7107 (2021)
  * [13] Gómez, L., Biten, A.F., Tito, R., Mafla, A., Rusiñol, M., Valveny, E., Karatzas, D.: Multimodal grid features and cell pointers for scene text visual question answering. Pattern Recognition Letters 150, 242–249 (2021)
  * [14] Gu, J., Kuen, J., Morariu, V., Zhao, H., Jain, R., Barmpalios, N., Nenkova, A., Sun, T.: Unidoc: Unified pretraining framework for document understanding. Advances in Neural Information Processing Systems 34 (2021)
  * [15] Ha, J., Haralick, R.M., Phillips, I.T.: Document page decomposition by the bounding-box project. In: Proceedings of 3rd International Conference on Document Analysis and Recognition. vol. 2, pp. 1119–1122. IEEE (1995)
  * [16] Ha, J., Haralick, R.M., Phillips, I.T.: Recursive xy cut using bounding boxes of connected components. In: Proceedings of 3rd International Conference on Document Analysis and Recognition. vol. 2, pp. 952–955. IEEE (1995)
  * [17] Hao, L., Gao, L., Yi, X., Tang, Z.: A table detection method for pdf documents based on convolutional neural networks. In: 2016 12th IAPR Workshop on Document Analysis Systems (DAS). pp. 287–292. IEEE (2016)
  * [18] Harley, A.W., Ufkes, A., Derpanis, K.G.: Evaluation of deep convolutional nets for document image classification and retrieval. In: 2015 13th International Conference on Document Analysis and Recognition (ICDAR). pp. 991–995. IEEE (2015)
  * [19] He, K., Gkioxari, G., Dollár, P., Girshick, R.: Mask r-cnn. In: Proceedings of the IEEE international conference on computer vision. pp. 2961–2969 (2017)
  * [20] He, Y., Chen, C., Zhang, J., Liu, J., He, F., Wang, C., Du, B.: Visual semantics allow for textual reasoning better in scene text recognition. arXiv preprint arXiv:2112.12916 (2021)
  * [21] Hong, T., Kim, D., Ji, M., Hwang, W., Nam, D., Park, S.: Bros: A pre-trained language model for understanding texts in document (2020)
  * [22] Jemni, S.K., Souibgui, M.A., Kessentini, Y., Fornés, A.: Enhance to read better: A multi-task adversarial network for handwritten document image enhancement. Pattern Recognition 123, 108370 (2022)
  * [23] Joshi, M., Chen, D., Liu, Y., Weld, D.S., Zettlemoyer, L., Levy, O.: Spanbert: Improving pre-training by representing and predicting spans. Transactions of the Association for Computational Linguistics 8, 64–77 (2020)
  * [24] Lebourgeois, F., Bublinski, Z., Emptoz, H.: A fast and efficient method for extracting text paragraphs and graphics from unconstrained documents. In: 11th IAPR International Conference on Pattern Recognition. Vol. II. Conference B: Pattern Recognition Methodology and Systems. vol. 1, pp. 272–273. IEEE Computer Society (1992)
  * [25] Lewis, D., Agam, G., Argamon, S., Frieder, O., Grossman, D., Heard, J.: Building a test collection for complex document information processing. In: Proceedings of the 29th annual international ACM SIGIR conference on Research and development in information retrieval. pp. 665–666 (2006)
  * [26] Li, M., Xu, Y., Cui, L., Huang, S., Wei, F., Li, Z., Zhou, M.: Docbank: A benchmark dataset for document layout analysis. In: Proceedings of the 28th International Conference on Computational Linguistics. pp. 949–960 (2020)
  * [27] Li, P., Gu, J., Kuen, J., Morariu, V.I., Zhao, H., Jain, R., Manjunatha, V., Liu, H.: Selfdoc: Self-supervised document representation learning. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 5652–5660 (2021)
  * [28] Litman, R., Anschel, O., Tsiper, S., Litman, R., Mazor, S., Manmatha, R.: Scatter: selective context attentional scene text recognizer. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 11962–11972 (2020)
  * [29] Mafla, A., Dey, S., Biten, A.F., Gomez, L., Karatzas, D.: Fine-grained image classification and retrieval by combining visual and locally pooled textual features. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 2950–2959 (2020)
  * [30] Mafla, A., Dey, S., Biten, A.F., Gomez, L., Karatzas, D.: Multi-modal reasoning graph for scene-text based fine-grained image classification and retrieval. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 4023–4033 (2021)
  * [31] Mafla, A., Rezende, R.S., Gómez, L., Larlus, D., Karatzas, D.: Stacmr: Scene-text aware cross-modal retrieval. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 2220–2230 (2021)
  * [32] Mathew, M., Bagal, V., Tito, R., Karatzas, D., Valveny, E., Jawahar, C.: Infographicvqa. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 1697–1706 (2022)
  * [33] Mathew, M., Karatzas, D., Jawahar, C.: Docvqa: A dataset for vqa on document images. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 2200–2209 (2021)
  * [34] Na, B., Kim, Y., Park, S.: Multi-modal text recognition networks: Interactive enhancements between visual and semantic features. arXiv preprint arXiv:2111.15263 (2021)
  * [35] O’Gorman, L.: The document spectrum for page layout analysis. IEEE Transactions on pattern analysis and machine intelligence 15(11), 1162–1173 (1993)
  * [36] Patil, A.G., Ben-Eliezer, O., Perel, O., Averbuch-Elor, H.: Read: Recursive autoencoders for document layout generation. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops. pp. 544–545 (2020)
  * [37] Powalski, R., Borchmann, Ł., Jurkiewicz, D., Dwojak, T., Pietruszka, M., Pałka, G.: Going full-tilt boogie on document understanding with text-image-layout transformer. In: International Conference on Document Analysis and Recognition. pp. 732–747. Springer (2021)
  * [38] Prasad, D., Gadpal, A., Kapadni, K., Visave, M., Sultanpure, K.: Cascadetabnet: An approach for end to end table detection and structure recognition from image-based documents. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops. pp. 572–573 (2020)
  * [39] Raja, S., Mondal, A., Jawahar, C.: Table structure recognition using top-down and bottom-up cues. In: European Conference on Computer Vision. pp. 70–86. Springer (2020)
  * [40] Ren, S., He, K., Girshick, R., Sun, J.: Faster r-cnn: Towards real-time object detection with region proposal networks. Advances in neural information processing systems 28, 91–99 (2015)
  * [41] Riba, P., Dutta, A., Goldmann, L., Fornés, A., Ramos, O., Lladós, J.: Table detection in invoice documents by graph neural networks. In: 2019 International Conference on Document Analysis and Recognition (ICDAR). pp. 122–127. IEEE (2019)
  * [42] Schreiber, S., Agne, S., Wolf, I., Dengel, A., Ahmed, S.: Deepdesrt: Deep learning for detection and structure recognition of tables in document images. In: 2017 14th IAPR international conference on document analysis and recognition (ICDAR). vol. 1, pp. 1162–1167. IEEE (2017)
  * [43] Simon, A., Pret, J.C., Johnson, A.P.: A fast algorithm for bottom-up document layout analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence 19(3), 273–277 (1997)
  * [44] Singh, A., Natarajan, V., Shah, M., Jiang, Y., Chen, X., Batra, D., Parikh, D., Rohrbach, M.: Towards vqa models that can read. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 8317–8326 (2019)
  * [45] Smith, R.: An overview of the tesseract ocr engine. In: Ninth international conference on document analysis and recognition (ICDAR 2007). vol. 2, pp. 629–633. IEEE (2007)
  * [46] Soto, C., Yoo, S.: Visual detection with context for document layout analysis. In: Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP). pp. 3464–3470 (2019)
  * [47] Souibgui, M.A., Biswas, S., Jemni, S.K., Kessentini, Y., Fornés, A., Lladós, J., Pal, U.: Docentr: An end-to-end document image enhancement transformer. arXiv preprint arXiv:2201.10252 (2022)
  * [48] Souibgui, M.A., Biten, A.F., Dey, S., Fornés, A., Kessentini, Y., Gomez, L., Karatzas, D., Lladós, J.: One-shot compositional data generation for low resource handwritten text recognition. In: Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. pp. 935–943 (2022)
  * [49] Souibgui, M.A., Kessentini, Y.: De-gan: A conditional generative adversarial network for document enhancement. IEEE Transactions on Pattern Analysis and Machine Intelligence (2020)
  * [50] Tito, R., Karatzas, D., Valveny, E.: Document collection visual question answering. arXiv preprint arXiv:2104.14336 (2021)
  * [51] Tito, R., Mathew, M., Jawahar, C., Valveny, E., Karatzas, D.: Icdar 2021 competition on document visual question answering. In: International Conference on Document Analysis and Recognition. pp. 635–649. Springer (2021)
  * [52] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * [53] Xu, Y., Xu, Y., Lv, T., Cui, L., Wei, F., Wang, G., Lu, Y., Florencio, D., Zhang, C., Che, W., et al.: Layoutlmv2: Multi-modal pre-training for visually-rich document understanding. arXiv preprint arXiv:2012.14740 (2020)
  * [54] Xu, Y., Li, M., Cui, L., Huang, S., Wei, F., Zhou, M.: Layoutlm: Pre-training of text and layout for document image understanding. In: Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. pp. 1192–1200 (2020)
  * [55] Zhong, X., ShafieiBavani, E., Jimeno Yepes, A.: Image-based table recognition: data, model, and evaluation. In: Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXI 16. pp. 564–580. Springer (2020)
  * [56] Zhong, X., Tang, J., Yepes, A.J.: Publaynet: largest dataset ever for document layout analysis. In: 2019 International Conference on Document Analysis and Recognition (ICDAR). pp. 1015–1022. IEEE (2019)