#  The Simplest Viscous Flow

William Graham Hoover with Carol Griswold Hoover
Ruby Valley Research Institute
Highway Contract 60, Box 601
Ruby Valley, Nevada 89833

###### Abstract

We illustrate an atomistic periodic two-dimensional stationary shear flow,
$u_{x}=\langle\ \dot{x}\ \rangle=\dot{\epsilon}y$, using the simplest possible
example, the periodic shear of just two particles ! We use a short-ranged
“realistic” pair potential, $\phi(r<2)=(2-r)^{6}-2(2-r)^{3}$. Many body
simulations with it are capable of modelling the gas, liquid, and solid states
of matter. A useful mechanics generating steady shear follows from a special
(“Kewpie-Doll” $\sim$ “$qp$-Doll”) Hamiltonian based on the Hamiltonian
coordinates $\\{q\\}$ and momenta $\\{p\\}$ : ${\cal H}(q,p)\equiv
K(p)+\Phi(q)+\dot{\epsilon}\sum qp$. Choosing $qp\rightarrow yp_{x}$ the
resulting motion equations are consistent with steadily shearing periodic
boundaries with a strain rate $(du_{x}/dy)=\dot{\epsilon}$. The occasional $x$
coordinate jumps associated with periodic boundary crossings in the $y$
direction provide a Hamiltonian that is a piecewise-continuous function of
time. A time-periodic isothermal steady state results when the Hamiltonian
motion equations are augmented with a continuously variable thermostat
generalizing Shuichi Nosé’s revolutionary ideas from 1984. The resulting
distributions of coordinates and momenta are interesting multifractals, with
surprising irreversible consequences from strictly time-reversible motion
equations.

Shear Flow, Statistical Physics, Hamiltonian Molecular Dynamics, Periodic
Boundaries, Fractals

##  I Introduction to Isothermal Molecular Dynamics

Molecular Dynamics began in the 1950s with Fermi, Pasta, and Ulam’s
investigation of one-dimensional anharmonic chains. That work was soon
followed with Alder and Wainwright’s dynamical exploration of the
pressure/volume equations of state for two-dimensional hard disks and three-
dimensional hard spheres. While Fermi was surprised that his chains never
reached equilibrium, Alder and Wainwright found that disks and spheres
equilibrated rapidly, in just a few collision times. Alder and Wainwright also
produced evidence for first-order melting/freezing transitions for hard
particles. These earliest molecular dynamics simulations used from dozens to
hundreds of particles, solving Newton’s motion equations. By 2008 quadrillion-
atom molecular dynamics $(N=10^{12})$ was feasibleb1 . Germann and Kadau found
that a single timestep, for a short-ranged pair potential in an unstable cubic
lattice and with $N=10^{4}\times 10^{4}\times 10^{4}$ particles, took about a
minute of computer time to execute.

Periodic boundary conditions are the simplest choice, with any particle
exiting at a system boundary simultaneously introduced on the opposite side.
In the absence of phase transitions such homogeneous 100-particle simulations
are sufficient to determine the mechanical pressure/volume and thermal
energy/temperature equations of state with accuracies of order 1%.

Nonequilibrium simulations require more sophisticated boundary conditions. In
his 1974 thesis work at Livermore Ashurst sandwiched rectangular Newtonian
manybody systems between two “fluid-wall” regions.b2  The particles in these
bounding walls had specified values of the mean velocity and the mean-squared
velocity imposed at every timestep. With different boundary velocities the
resulting shear viscosity could be measured. With different boundary kinetic
temperatures, $K/N\sim T\sim\langle\ (p^{2})\ \rangle$, heat flow could be
simulated and the thermal conductivity measured. By 1980 it was well accepted
that the equilibrium and nonequilibrium properties for manybody systems
provided excellent models for the behavior of real fluids, both at, and away
from, equilibrium.

In 1984 Shuichi Nosé formulated a temperature-dependent Hamiltonian
reproducing Gibbs’ canonical-ensemble distribution of energies,
$\propto\exp[-E/kT]$. His formulationb3 ; b4  included a relaxation time
$\tau$ with rapid relaxation consistent with Ashurst’s rescaling of the mean-
squared velocity to impose temperatures within the fluid walls. The “Nosé-
Hoover” version of this workb5  is better suited to computation than was
Nosé’s original Hamiltonian-based algorithm. For $N$ particles of unit mass in
$D$ space dimensions the equations of motion have the form:

$\dot{q}=p\ ;\ \dot{p}=F(q)-\zeta p\ ;\ \dot{\zeta}=[\ (\sum p^{2}/NDkT)-1\
]/\tau^{2}\ .$

For $D=2$ $q$ and $p$ have both $x$ and $y$ components: $q=(x,y)\ ;\
p=(p_{x},p_{y})$. $k$ is Boltzmann’s constant (for convenience chosen equal to
unity) and $T$ is the specified kinetic temperature, $\langle\ p^{2}/2k\
\rangle=\langle\ p^{2}/2\ \rangle.$

The “friction coefficient” $\zeta$ imposes integral control on the kinetic
temperature in a time of order $\tau$. In the instantaneous limit,
$\tau\rightarrow 0$, integral control corresponds to Ashurst’s velocity
rescaling. In the equilibrium case, without external driving, Hoover pointed
out that the time-averaged distribution of the friction coefficient $\zeta$ is
Gaussian, with zero mean. We will see that this changes away from equilibrium
where $\langle\ \zeta\ \rangle$ corresponds to entropy production, and is
necessarily positive.

##  II Imposing Shear on Molecular Dynamics with Doll’s Tensor

With controlled boundary conditions molecular dynamics can be generalized to
nonequilibrium flows characterized by viscosity and thermal conductivity.
These material properties characterize the response of fluids or solids to
gradients in the velocity or the temperature. Viscosity is the simpler of
these possibilities. Periodic boundary conditions are natural for shear
deformation.

A straightforward viscous modification of equilibrium molecular dynamics can
be implemented by imparting a systematic motion to the neighbors of a
periodically-repeated simulation, as illustrated in Figure 1. The figure shows
the details of two particles interacting periodically. Three parallel periodic
rows, of $3N=6$ particles each are required. There is a systematic boundary
motion in the $x$ direction which increases linearly with $y$. The “Doll’s
Tensor” Hamiltonian (named after the $qp$ term reminiscent of the Kewpie Doll)
consistent with this motion was introduced by Hoover in 1980b6 . With
$(qp)\rightarrow(yp_{x})$ this Hamiltonian is

${\cal H}_{\rm Doll^{\prime}s}=\sum(p^{2}/2)+\Phi+\dot{\epsilon}\sum
yp_{x}\rightarrow\\{\ \dot{x}=+(\partial{\cal H}/\partial
p_{x})=p_{x}+\dot{\epsilon}y\ \\}\ .$

The new term proportional to the strain rate, $\dot{\epsilon}=(du_{x}/dy)$,
along with periodic boundaries in the $x$ and $y$ directions, generates a
shear flow. With just two particles, such a model provides the simplest
atomistic simulation of viscous shear flow. The Doll’s Tensor equations of
motion, applied to momenta $(p_{x},p_{y})$, with vanishing sums, produce a
systematic increase of $x$ velocity with $y$, $\dot{\epsilon}y$:

$\\{\ \dot{x}=p_{x}+\dot{\epsilon}y\ ;\ \dot{y}=p_{y}\ ;\ \dot{p}_{x}=F_{x}\
;\ \dot{p}_{y}=F_{y}-\dot{\epsilon}p_{x}\ \\}\ .$

For simplicity we consider forces $\\{\ F_{x},F_{y}\ \\}$ derived from a pair
potential $\phi(r)$: $F(r)=-\nabla\phi(r)$. In this case the total potential
is a simple sum over all pairs of particles, $\Phi=\sum\phi_{i<j}$.

With pair forces the change in the “internal” energy $E=\sum(p^{2}/2)+\Phi$ is
precisely consistent with the first-law thermodynamic relation for energy
conservation in the absence of heat flow, $\Delta{\rm
Heat}=0\longrightarrow\Delta E=\Delta{\rm Heat}-\Delta{\rm Work}=-\Delta{\rm
Work}$:

$\dot{E}=-\dot{\epsilon}\sum^{N}(p_{x}p_{y})-\dot{\epsilon}\sum^{\rm
pairs}(Fxy/r)_{i<j}=-\dot{\epsilon}P_{xy}V\ .$

Here the force-dependent pair sum includes the particle separations in the $x$
and $y$ directions:

$(xy/r)\equiv(x_{i}-x_{j})(y_{i}-y_{j})/\sqrt{[(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}]}\
.$

Notice particularly that the First Law of Thermodynamics just given is
perfectly time-reversible. In a reversed process $\dot{E}$ necessarily changes
sign, as does the strain-rate $\dot{\epsilon}$. The shear component of the
pressure tensor does not change with the direction of time, thus violating the
empirical experience that shearing a viscous fluid requires work, whether the
shearing is clockwise or counterclockwise. It is definitely paradoxical that
the time-reversible laws of microscopic mechanics provide simulations obeying
the Second Law, that work can be converted to heat, but not vice versa. We
will extend this paradox, already present in Hamiltonian mechanics, to a
nonequilibrium mechanics in which both heat and work enter, and with the same
paradoxical disagreement with the Second Law.

We adopt a particularly simple pair potential, with a range $r<2$, and
sufficiently realistic to generate gas, liquid, and solid phasesb7 . See
Figure 2 for a plot. This pair potential function is a simple polynomial in
$r$, with an attractive minimum of unity at a separation of unity,
$\phi(r=1)=-1$:

$\phi_{3}^{L}(r)=(2-r)^{6}-2(2-r)^{3}\rightarrow F=6[(2-r)^{5}-(2-r)^{2}]\ .$

Thus the potential minimum sets both the length and the energy scales. Both
the force and its first derivative vanish at the cutoff radius $r=2$,
optimizing the numerical accuracy of computer simulations.

The simplest viscous flow, illustrated in Figure 1, is a two-particle
dynamical system. With a potential range of 2 a periodic $4\times 4$ box is
the smallest for which the interaction of the two particles involves just one
of the 9 possibilities illustrated in the figure. With vanishing center-of-
mass velocity the two-particle problem is equivalent to a one-particle one,
with the boundary conditions including the steady motion of cells containing
image particles. In the two-particle version the four time-dependent variables
$(x,y,p_{x},p_{y})$ for Particle 2 exactly equal the negatives of the same
variable set for Particle 1.

Let us turn next to the integration details. We like to avoid the need for
developing special integration algorithms for the various types of problems
described by atomistic simulations. For simplicity and accuracy we prefer
fourth-order Runge-Kutta integration for all those problems that are
describable with first-order differential equationsb8 . These included models
in which constraints on the momentum or temperature or stress are imposed by
Hamiltonian or nonHamiltonian constraints. Let us review our favorite Runge-
Kutta integrator by applying it to two instructive problems. The first is the
simple harmonic oscillator. This problem is useful for code development as the
Runge-Kutta algorithm has an explicit analytic solution for the oscillator.

We introduce additional complexity and nonlinearity in our second pedagogical
example, the periodic dynamics of two two-dimensional particles using the
realistic pair potential just described, a finite-ranged very smooth
polynomial in the separation. The boundaries for that problem simplify in the
absence of the systematic velocity gradient of Figure 1. The static-boundary
problem appears to be chaotic and might even be ergodic despite its
simplicity. The problem’s description can be reduced from eight dimensions to
just three by taking the constancy of the center of mass, its motion, and the
conserved total energy into account. Three dimensions is just sufficient for
chaos. Let us consider the details of the Runge-Kutta integration algorithm
and its applications to the harmonic oscillator and the two-dimensional
periodic two-body problem next.

##  III Two Problems Illustrating the Runge-Kutta Algorithm

###  III.1 The Fourth-Order Runge-Kutta Integration Algorithm

For completeness we review the programming necessary to implement the fourth-
order Runge-Kutta integrator. The simplicity and accuracy of this integrator
make it our first choice whenever possible. Consider solving $N$ first-order
differential equations using the Runge-Kutta algorithm.

A main program calls the fourth-order Runge-Kutta subroutine to execute each
timestep, advancing the $N$, ${\tt N}$ in the subroutine, time-dependent
variables, stored in the vector ${\tt YNOW}$. In the subroutine, given below,
the righthandsides of the $N$ first-order differential equations are evaluated
successively for four different vectors, providing an averaged time derivative
which is accurate through the fourth power of the timestep $dt$, DT in the
program. As the timestep is specified at the beginning of each step, it can be
changed “on the fly”, if necessary. For instance, a useful diagnostic for
accuracy is the difference between the result of a single step $dt$ and that
of two successive steps half as large, $dt/2$. If the difference is too large
the step DT can be reduced to ${\tt DT/2}$. If the difference is too small
${\tt DT}$ can be doubled. The fourth-order Runge-Kutta algorithm is called at
each timestep and has the effect of converting the $N$-dimensional solution
vector YNOW(T) at time $t$ to its value at time $t+dt$: ${\tt
YNOW(T)}\rightarrow{\tt YNOW(T+DT)}$.

      CALL RHS(YNOW,YDOT)
      DO 1 J = 1,N
    1 D1(J) = YDOT(J)
      DO 2 J = 1,N
    2 YNEW(J) = YNOW(J) + DT*D1(J)/2
      CALL RHS(YNEW,YDOT)
      DO 3 J = 1,N
    3 D2(J) = YDOT(J)
      DO 4 J = 1,N
    4 YNEW(J) = YNOW(J) + DT*D2(J)/2
      CALL RHS(YNEW,YDOT)
      DO 5 J = 1,N
    5 D3(J) = YDOT(J)
      DO 6 J = 1,N
    6 YNEW(J) = YNOW(J) + DT*D3(J)
      CALL RHS(YNEW,YDOT)
      DO 7 J = 1,N
    7 D4(J) = YDOT(J)
      DO 8 J = 1,N
    8 YNOW(J) = YNOW(J) + DT*(D1(J)+2*(D2(J)+D3(J))+D4(J))/6

###  III.2 Hamiltonian Dynamics of the Harmonic Oscillator

The simplest instructive application of the Runge-Kutta integrator is the
harmonic oscillator example. Only two motion equations are involved,
$\dot{q}=+p$ and $\dot{p}=-q$. With $(q,p)$ initially $(+1,0)$ the analytic
solution has a period of $2\pi$: $q=\cos(t)\ ;\ p=-\sin(t)$. Dividing the
interval $0<t<2\pi$ into 100 steps with initial values YNOW(1) = 1, YNOW(2) =
0 gives a double-precision estimate of the algorithm’s accuracy. The error in
$q^{2}+p^{2}$ at $t=2\pi$ (twice the final total oscillator energy) is
$8.5415\times 10^{-8}$. The history of the decrease is accurately linear in
time, as is illustrated in Figure 3. Reducing the timestep twofold gives an
error nearly 32 times smaller, $2.6702\times 10^{-9}$. The error ratio, 31.99,
confirms that the energy error varies as $dt^{5}$, one order better than the
larger seperate errors in the coordinate and momentum.

Notice that the input vector YNOW is incremented by the mean value of four
estimated RHS vectors in each Runge-Kutta iteration. The derivatives D1, D2,
D3, D4 are evaluated at the initial point and three subsequent nearby points.
When the righthandsides of the differential equations are explicitly time-
dependent, as in periodic shear, the four successive evaluations of the
righthandsides in each timestep are evaluated at times T for D1, T + DT/2 for
D2 and D3, and finally, T + DT for D4.

The fourth-order Runge-Kutta solution of the oscillator equations $\dot{q}=p$
and $\dot{p}=-q$ with $[q,p]$ initially $[+1,0]$, has a relatively simple and
useful analytic form:

$[q_{n},p_{n}]=[1-(dt^{6}/72)+(dt^{8}/576)]^{n/2}[+\cos(n\lambda),-\sin(n\lambda)]\
{\rm with}$
$\lambda\equiv\arctan\frac{[dt-(dt^{3}/6)]}{[1-(dt^{2}/2)+(dt^{4}/24)]}\ .$

For sufficiently large $dt$, so that $dt^{8}/576$ is properly represented,
this result provides a useful and informative programming check for the
implementation of the Runge-Kutta integration algorithm. We can see that
$(q^{2}+p^{2})$, rather than being constant at long times $ndt$, decreases as
$ndt^{6}/72=tdt^{5}/72$, one order of $dt$ better than would be expected for a
fourth-order method. Although the oscillator amplitude is fifth-order, both
the coordinate and the momentum $[q,p]$ do exhibit errors of order $tdt^{4}$.
Further analysis shows that this difference in orders is due to a “phase
error” with the oscillator slowed to fourth order in $dt$, rather than reduced
in amplitude. Thus the dominant error at a fixed time is a negative phase
error of order $tdt^{4}$. With twenty timesteps per period a double-precision
comparison of the numerical and analytical results from Runge-Kutta up to a
time $2\pi$ gives agreement to a dozen significant figures.

With this interesting oscillator task successfully completed, confirming the
utility and accuracy of the algorithm, fully-fledged nonlinear simulations
with realistic interparticle interactions as well as both static and dynamic
boundary conditions on the coordinates provide the next steps toward our goal
of simulating “The simplest viscous flow”.

###  III.3 Equilibrium Dynamics with a Realistic Potential

The equilibrium Runge-Kutta integration of two-body dynamics with static
periodic boundaries and with interparticle forces,
$F(r<2)=-(d\phi/dr)=6(2-r)^{5}-6(2-r)^{2}$, is an excellent computational next
step. Conservation of the energy, $E=\Phi+K=\phi_{12}+(p_{1}^{2}+p_{2}^{2})/2$
is a (nearly) foolproof check on the programming. To apply fourth-order Runge-
Kutta integration the eight time-dependent variables are stored in an eight-
dimensional $yy$ vector. For the coordinates we choose

$yy(1)=x_{1}\ ;\ yy(2)=x_{2}\ ;\ yy(3)=y_{1}\ ;\ yy(4)=y_{2}$

and for the momenta

$yy(5)=p_{x_{1}}\ ;\ yy(6)=p_{x_{2}}\ ;\ yy(7)=p_{y_{1}}\ ;\ yy(8)=p_{y_{2}}\
.$

The $x$ coordinates are confined to the $4\times 4$ periodic box, both
initially and after every integration step, with four statements

    if(x1.lt.-2) yy(1) = yy(1) + 4
    if(x1.gt.+2) yy(1) = yy(1) - 4
    if(x2.lt.-2) yy(2) = yy(2) + 4
    if(x2.gt.+2) yy(2) = yy(2) - 4

and likewise for the $y$ coordinates, $y1=yy(3)$ and $y2=yy(4)$.

The eight motion equations for the two-body problem are the same for both
particles:

$\\{\ \dot{x}=p_{x}\ ;\ \dot{y}=p_{y}\ ;\ \dot{p}_{x}=F_{x}\ ;\
\dot{p}_{y}=F_{y}\ \\}\ ,$

To compute the forces the relative displacement of Particle 2 relative to
Particle 1 is first taken as the simple difference with both particles in the
central box.

    x12 = x1 - x2
    y12 = y1 - y2

Then the nearest-image distance between the two particles, smallest among the
nine possibilities visible in Figure 1, is calculated by imposing two
conditions on the $x$ separation,

    if(x12.lt.-2) x12 = x12 + 4
    if(x12.gt.+2) x12 = x12 - 4

and likewise for ${\tt y12}$. It is convenient, and perfectly permissible, to
impose the fixed center-of-mass symmetry at the conclusion of each step: yy(J)
= -yy(J-1) for J equal to 2, 4, 6, and 8.

A bit of experimentation shows that a timestep $dt=0.005$ provides good energy
conservation. At a time of 1000 the initial energy of unity has decayed to
0.9999937. This time-dependent energy decay is illustrated in Figure 3, at the
left for the oscillator and at the right for the two-body problem. The two are
also compared on a logarithmic scale in Figure 4 with similar shapes, though
the oscillator results are certainly smoother.

The jumps in the two-body energy, visible in both figures, occur at the
collisional turning points. In spite of the relative simplicity of this
problem the long-time two-body trajectories illustrated in Figure 5 suggest
that there is sufficient chaos in this problem to provide a near-ergodic
coordinate distribution. Although the apparent phase-space dimensionality here
is eight, much more than the three required for chaos, this equilibrium
problem is actually just three-dimensional ! The two-body problem, with the
center-of-mass fixed, or equivalently, expressed in relative coordinates, can
be reduced to a four-dimensional problem. It could be viewed as the constant-
energy motion of a mass-point particle within a periodic $4\times 4$ square
with four fixed particles at the corners controlling the mass point’s motion.
Conservation of energy further restricts that motion, however expressed, to a
three-dimensional energy shell, the actual minimum for chaos. In fact
isoenergetic simulations with unit energy show no tendency to approach either
of the two collisionless nonchaotic solutions:

$\\{\ x=\pm t\ ;\ y=\pm 1\ ;\ p_{x}=\pm 1\ ;\ p_{y}=0\ \\}\ ;\\{\ x=\pm 1\ ;\
y=\pm t\ ;\ p_{x}=0\ ;\ p_{y}=\pm 1\ \\}\ .$

In these hypothetical collisionless cases the motion would parallel either of
the two coordinate axes, $x$ or $y$.

The chaotic equilibrium problem, just described and pictured in Figure 5 is a
useful stepping stone toward the stationary two-body nonequilibrium shear flow
that is our goal. Let us consider next the modelling of a steadily moving
periodic boundary. Incorporating this driving mechanism will help us simulate
a “stationary”, actually time-periodic, nonequilibrium flow. The moving
boundary introduces a new variable, the phase or “strain” of the motion, which
can be reset to zero when the coordinates in all nine cells regain their
perfect square-lattice symmetry.

There is a second barrier to stationarity. Because the moving boundaries do
viscous “work” on the two particles in the central cell we will need to
“thermostat” the motion to avoid a relatively rapid divergence of the kinetic
energy. We will use a Nosé-Hoover thermostat force, $-\zeta p$ to control the
kinetic energy. With the addition of a moving boundary and a thermostat our
two-body equilibrium model is converted to a nonequilibrium steady flow. We
describe its characteristics next. For simplicity we restrict our description
to the case of unit strain rate so that the locations of the periodic-image
cells repeat at integral values of the time.

##  IV Nonequilibrium Shear Flow with Two Particles

To reach our goal of “The Simplest Viscous Flow” requires two additions to the
equations of motion and corresponding modifications of the integrator.
Applying the Runge-Kutta algorithm to shear flow requires augmenting the
motion equations to include [1] the explicit time dependence of the boundary
conditions and [2] the control of the kinetic energy. Let us consider the
modifications in more detail.

###  IV.1 Imposing a Fixed Strain Rate on the Dynamics

In introducing shear the three square cells with $-2<y<+2$ obey the usual
periodic boundary conditions in the $x$ direction. The $y$ direction is
different. The shearing motion of three square cells above and three below the
central cells is induced by giving the twelve particles in the upper and lower
cells additional velocities equal to the product of the strain rate and the
cell height: $\dot{x}=p_{x}+4$ above, $\dot{x}=p_{x}-4$ below. We choose to
restrict the cell-to-cell strain [ $\epsilon=(\Delta x/\Delta y)$ relative to
the perfect square-lattice structure ] of the nine cells of Figure 1 to lie in
the range from $-0.5<\epsilon<+0.5$. From the numerical standpoint simply
resetting the strain, $+(1/2)\rightarrow-(1/2)$, every 200 steps, prevents the
loss of significant figures in the horizontal coordinates $x$. Otherwise
simulations with billions of timesteps, corresponding to strains on the order
of millions, would erode the accuracy of the horizontal coordinates.

By limiting the strain to lie between the two extreme values $\mp(1/2)$ shown
in the figure, it is guaranteed that the nearest-neighbor image of Particle 1
lies uniquely in one of the nine cells shown. Likewise for Particle 2. Thus
the nearest-image potential and forces can both be computed for the central-
cell Particles 1 and 2 by selecting the minimum-distance particle pair out of
the nine possibilities. The unique pair $(x_{ij},y_{ij})$ which minimizes the
separation $\sqrt{(x_{ij}^{2}+y_{ij}^{2})}$ is selected from the nine
possibilities whenever the righthandsides of the differential equations are
calculated, four times within every Runge-Kutta timestep.

###  IV.2 Imposing a Fixed Time-Averaged Temperature on the Dynamics

We choose a basic time-reversible Nosé-Hoover thermostat to control the
temperature. In the time-reversed version both the momenta and the friction
coefficient $\zeta$ change signs. For simplicity we choose a thermostat
relaxation time of unity in the $\dot{\zeta}$ equation. The thermostatted
motion equations are the following, where $T$ is a specified constant target
temperature:

$\\{\ \dot{x}=p_{x}+\dot{\epsilon}y\ ;\ \dot{y}=p_{y}\ ;\
\dot{p}_{x}=F_{x}-\zeta p_{x}\ ;\ \dot{p}_{y}=F_{y}-\dot{\epsilon}p_{x}-\zeta
p_{y}\ \\}\ ;\ $ $\dot{\zeta}=K-NT\ {\rm where}\
K=\sum(p_{x}^{2}+p_{y}^{2})/2\ .$

These equations are time-reversible. To see this note first that the overdot
“$\ {}^{\bf.}\ $” signifies a time derivative so that ($\dot{x},\dot{y}$)
change signs with the direction of time while ($\dot{p}_{x},\dot{p}_{y}$) and
the forces do not. Evidently $\dot{\epsilon}$ must change sign for
reversibility, in agreement with our imagining the shear process from $-0.5$
to $+0.5$ being reversed, from $+0.5$ to $-0.5$. Likewise the sign of the
friction coefficient must be reversed just as if it were a momentum. In fact,
in Nosé’s original derivation of his thermostat equations the friction
coefficient was a momentum. Its conjugate variable $s$ is Nosé’s “time-scaling
variable”, which behaved like a coordinate. Though the motion equations just
given are patently time-reversible we will see, and come to understand, that
the reversibility is illusory, and if attempted will succeed for only a
relatively short time.

In this shear flow problem there is a competition between the driving shear
force, characterized by the strainrate $\dot{\epsilon}$, and the damping
activity of the thermostat, characterized by $\zeta$. These two forces
represent work and heat respectively. $\epsilon$ provides work and the
thermostat variable $\zeta$ on average extracts the resulting irreversible
heating, maintaining a kinetic temperature $2T=\langle\ K\ \rangle=\langle\
(p_{1}^{2}+p_{2}^{2})/2\ \rangle$. With a strainrate of unity a short
computation shows that a thermostated kinetic energy of 0.3 provides a
longtime chaotic solution.

At lower temperatures an interesting phenomenon ensues. Low temperature favors
an interparticle spacing of unity, the energy minimum. Figure 6 displays
exactly this. At lower temperatures, with kinetic energies in the range from
0.25 down to zero, the energy minimum becomes irresistable and the fractal
distributions, like those generated at 0.30 and higher, collapse to a fixed
point in the nine-dimensional phase space. Such points correspond to an
interparticle separation barely exceeding unity, and close to the potential
minimum for the interacting Particles 1 and 2. The two particles occupy fixed
points on a curve approximating a circle, with diameter slightly greater than
unity.

As a consequence, once sufficiently cooled, Particles 1 and 2 no longer can
interact with their periodic images to the sides or above or below, for which
a central-box separation of the particles greater than 2 is required. Thus
low-temperature simulations are no longer representative of boundary-driven
shear. Figure 6 shows these unphysical fixed-point solutions. With this
thermostat model the “temperature” is entirely taken up by $p_{x}^{2}$ and
exactly cancels the coordinate-dependent contribution $\dot{\epsilon}y$ to
$\dot{x}$. Meanwhile $p_{y}$ and $\dot{y}$ vanish.

$\dot{x}=p_{x}+\dot{\epsilon}\times y=0.5+(1.0\times-0.5)=0\ (\rm{For}\
2T=K=0.25\ {\rm and}\ \dot{\epsilon}=1)\ .$

Because our pair potential has its maximum value $2^{6}-2^{4}=48$ at $r=0$ it
is likewise evident that coupling with the boundary will have a much
diminished effect at sufficiently high temperatures. The high-temperature
limit is simply uniform over the $4\times 4$ square, assuming sufficient
perturbation from the thermostat to favor irrational coordinates and momenta
over the rational. Chaotic viscous flow can only exist when both the driving
and the dissipative influences on the motion are comparable, so that their
time-averaged contributions cancel.

Figure 7 pictures an $(x,y)$ projection of the fractal distribution at
$K=2T=0.30$. In the figure the six boundary-cell images of Particles 1 and 2
in the cells above and below the central cell are swept from left to right
(above) and right to left (below) accounting for the smearing of probability
at the top and bottom of the central cell. Figure 8 shows superpositions of
snapshots taken every 200 steps (corresponding to a unit increase in
$\epsilon=\dot{\epsilon}t$). Just two of the 200 different phases, at integral
and half-integral strains respectively, making up a unit time interval, are
shown in the figure. Evidently the probability density is periodic (with
period 200 steps here) in time so that a detailed microstate description would
need to include the phase of the motion as well as the the values of the
coordinates and momenta of Particles 1 and 2.

Figure 9 shows two phases of the two-body momentum distribution. The kinetic
part of the pressure tensor $P_{xy}$ is proportional to the summed values of
$p_{x}p_{y}$, identical for the two particles in our model. The figure shows
that the product is, on average, negative, as is required for agreement with
continuum mechanics where the shear stress, the negative of $P_{xy}$, is equal
to the product of the strain rate and the (positive) coefficient of shear
viscosity.

##  V Describing the Results of The Nonequilibrium Simulations

###  V.1 Paradoxical Nonequilibrium Phase-Volume Changes

Let us consider the nonequilibrium shearing motion, time-averaged over many
repetitions of the boundary time-periodicity in five-dimensional space,
$\\{x_{1},y_{1},p_{x_{1}},p_{y_{1}},\zeta\\}$. Averaging the comoving dilation
rate in this phase space for an infinitesimal five-dimensional volume element
$\otimes$ gives an empirical inevitable inequality:

$\langle\ (\dot{\otimes}/\otimes)\ \rangle=\langle\
(\partial\dot{x}_{1}/\partial x_{1})+(\partial\dot{y}_{1}/\partial
y_{1})+(\partial\dot{p}_{x_{1}}/\partial
p_{x_{1}})+(\partial\dot{p}_{y_{1}}/\partial
p_{y_{1}})+(\partial\dot{\zeta}/\partial\zeta)\ \rangle=$ $0+0-\langle\ \zeta\
\rangle-\langle\ \zeta\ \rangle+0=-2\langle\ \zeta\ \rangle<0\ .$

It is characteristic of time-reversible nonequilibrium stationary states that
their phase volume inevitably shrinks, as the time-averaged friction
coefficient $\zeta$ is invariably positive.

From the computational standpoint a negative friction coefficient would
correspond to comoving phase-space growth, numerical instability, and
numerical divergence. The inevitable positive friction coefficient,
corresponding to comoving shrinkage, and necessary for computational
stability, at first looks paradoxical in view of the formal time-reversibility
of the motion. To understand, recall that ordinary equilibrium Hamiltonian
mechanics conserves the comoving phase volume precisely. This is a consequence
of the Hamiltonian motion equations,

$\\{\ \dot{q}=+(\partial{\cal H}/\partial p)\ ;\ \dot{p}=-(\partial{\cal
H}/\partial q)\ \\}\ \longrightarrow$ $\sum[\ (\partial\dot{q}/\partial
q)+(\partial\dot{p}/\partial p)\ ]=\sum[\ (\partial/\partial
q)(+(\partial{\cal H}/\partial p))+(\partial/\partial p)(-(\partial{\cal
H}/\partial q))\ ]\equiv 0\ .$

At thermal equilibrium, imposed by Nosé-Hoover friction, the time-averaged
friction coefficient $\langle\ \zeta\ \rangle$ typically vanishes. Typically
the equilibrium $\zeta$ has a Gaussian distribution with zero mean. At
equilibrium the logarithm of the phase volume $\otimes$ corresponds to
entropy. We will see that this changes away from equilibrium where $\langle\
\zeta\ \rangle$ corresponds to entropy production, and is necessarily
positive.

###  V.2 A Microscopic Version of the Second Law of Thermodynamics

Away from equilibrium, with driving forces and a compensating nonHamiltonian
friction coefficient, the comoving dilation rate $(\dot{\otimes}/\otimes)$ is
nonzero. The internal energy $\Phi+K$ and the phase volume are no longer
constant. To prevent the divergence of both energy and phase volume the mean
friction coefficient is necessarily, though paradoxically, positive. This
positive friction, in the face of reversibility, is the microscopic version of
the Second Law of Thermodynamics, describing the conversion of work to heat.
The reversed version of this shearing process would be a violation of that
law. This consonance of reversible microscopic dynamics with macroscopic
thermodynamics was at first surprising and controversial, but ultimately
became both edifying and satisfying, a major result obtained by analyzing the
nonequilibrium simulations pioneered by Ashurst. For two early time-reversible
one-body “work-to-heat” simulations obeying the Second Law see References 9
and 10.

###  V.3 The Nature of Nonequilibrium Lyapunov Spectra

From the standpoint of dynamical systems the loss of phase volume is an
expected feature for models generating fractal distributions, but only rarely
are these models time-reversibleb11 ; b12 . An investigation of the Lyapunov
spectrum for fractal distributions invariably results in a negative sum, even
in formally time-reversible cases. The negative sum signals a collapse of
phase volume to a zero-volume fractal, limit cycle, or fixed point.

A worthy research goal is the Lyapunov analysis of the family of five-
dimensional flow equations with ranges of strain rates $\\{\ \dot{\epsilon}\
\\}$ and temperatures $\\{\ T\ \\}$:

$\\{\ \dot{x}=p_{x}+\dot{\epsilon}y\ ;\ \dot{y}=p_{y}\ ;\
\dot{p}_{x}=F_{x}-\zeta p_{x}\ ;\ \dot{p}_{y}=F_{y}-\dot{\epsilon}p_{x}-\zeta
p_{y}\ \\}\ ;\ $ $\dot{\zeta}=(K/T)-1\ {\rm where}\ K=(p_{x}^{2}+p_{y}^{2})/2\
.$

The averaged phase-space dilation rate, $\dot{\otimes}$, can be written in
terms of the five long-time-averaged Lyapunov exponents,

$\langle\ (\dot{\otimes}/\otimes)\ \rangle=-2\langle\ \zeta\ \rangle=\langle\
\lambda_{1}+\lambda_{2}+\lambda_{3}+\lambda_{4}+\lambda_{5}\ \rangle\ .$

Exactly this same equation holds instantaneously though the corresponding
“local” values of the exponents depend upon the chosen coordinate system.

The instantaneous values of the five Lyapunov exponents result from the
analysis of growth rates for five orthogonal time-dependent vectors in the
space, $\\{\ \delta_{1}\ \dots\ \delta_{5}\ \\}$. Evidently the entire range
of phase-space dimensionalities from 0 (a fixed point) to 5 (ergodic coverage,
which should result for reasonable combinations of temperature with small
strain rates) could be analyzed. We heartily recommend such studies to our
readers. These simple dynamical problems are a relatively complex interesting
area ripe for study. Surely a periodic symmetric two-body problem is indeed
the simplest viscous flow.

###  V.4 The Relevance of Time-Reversible Maps

In 1996 Hoover, Kum, and Posch constructed and analyzed the results of several
time-reversible mapsb13 . Those maps were all defined in a two-dimensional
unit square with $-(1/2)<(x,y)<+(1/2)$. Shears, for instance, obeyed the maps
${\bf X(x,y)}$, similar to the shearing considered here, and ${\bf Y(x,y)}$:

${\bf X(x,y)}:x\rightarrow x^{\prime}=x+\epsilon y\ ;\ y\rightarrow
y^{\prime}=y\ ;$ ${\bf Y(x,y)}:x\rightarrow x^{\prime}=x\ ;\ y\rightarrow
y^{\prime}=y+\epsilon x\ .$

In these two-dimensional cartesian maps $x$ represents a coordinate and $y$ a
momentum so that time can be reversed by changing the sign of $y$, analogous
to reversing time at a fixed configuration. In this map modelling of time
reversibility it is convenient to introduce the time-reversal map ${\bf T}$:

${\bf T}(x,y):x\rightarrow x^{\prime}=x\ ;\ y\rightarrow y^{\prime}=-y\ .$

With this definition ${\bf T}$ and its inverse ${\bf T^{-1}}$ are one and the
same. In fact any map ${\bf M}(x,y)$ is time-reversible if its inverse is the
sequence of operations going backward in time, ${\bf TMT}$, so that ${\bf
MM^{-1}}(x,y)={\bf M^{-1}M}(x,y)=(x,y)$.

The 1996 work showed that the mapping ${\bf XYYX}$ was ergodic, with constant
density throughout the unit square. Differently, and significantly, the
mapping ${\bf XYPYX}$, where ${\bf P}$ was both time-reversible and included
regions of density expansion and contraction, produced fractal distributions,
ergodic for small deviations from the ${\bf XYYX}$ equilibrium and with an
information dimension reduced quadratically in the driving deviation away from
equilibrium. Thus the two-dimensional time-reversible maps already contain a
mechanism for irreversibility resembling that found in the present “simple”
five-dimensional viscous flows.

##  VI Acknowledgment

We thank Clint Sprott, Karl Travis, and Kris Wojciechowski for their readings
of a previous draft. In particular Karl pointed out the importance of two-
dimensional maps for the understanding of the irreversibility of the present
five-dimensional nonequilibrium flows. A particularly compelling map is the
compressible Baker Map, described in two arXiv articles, 1909.04526, “Random
Walk Equivalence to the Compressible Baker Map and the Kaplan-Yorke
Approximation to Its Information Dimension” and 1910.12642, “2020 Ian Snook
Prize Problem : Three Routes to the Information Dimensions for a One-
Dimensional Stochastic Random Walk and for an Equivalent Prototypical Two-
Dimensional Baker Map”. As no suitable entries for the latter problem were
received we again offer a Snook-Prize cash award of one thousand United States
dollars for the best work addressing the 2020 problem received during the
calendar year 2021. See our “$1000 SNOOK PRIZES FOR 2021: The Information
Dimensions of a Two-Dimensional Baker Map”, Computational Methods in Science
and Technology 21, 93-94 (2021).

##  VII Addendum of 8 July 2021

Karl suggested that we include the second-order Runge-Kutta oscillator
solution, with

QNEW = Q + P*DT - Q*(DT*DT/2) ; PNEW = P - Q*DT - P*(DT*DT/2) .

The solution is:

$[q_{n},p_{n}]=[1+(dt^{4}/4)]^{n/2}[+\cos(n\lambda),-\sin(n\lambda)]\ {\rm
with}$ $\lambda\equiv\arctan\frac{dt}{[1-(dt^{2}/2)]}\ .$

Because this oscillator amplitude now increases exponentially in time the lack
of molecular dynamics simulations using second-order Runge-Kutta algorithms is
unsurprising!

## References

  * (1) T. C. Germann and K. Kadau, “Trillion-Atom Molecular Dynamics Becomes a Reality”, Journal of Modern Physics C 19, 1315-1319 (2008).
  * (2) W. T. Ashurst, Dense Fluid Shear Viscosity and Thermal Conductivity via Nonequilibrium Molecular Dynamics (University of California, Davis, 1974) 187 pages.
  * (3) S. Nosé, “A Unified Formulation of the Constant Temperature Molecular Dynamics Method”, The Journal of Chemical Physics 81, 511-519 (1984).
  * (4) S. Nosé “Molecular Dynamics Method for Simulations in the Canonical Ensemble”, Molecular Physics 52, 255-268 (1984).
  * (5) Wm. G. Hoover, “Canonical Dynamics. Equilibrium Phase-Space Distributions”, Physical Review A 31, 1695-1697 (1985).
  * (6) W. G. Hoover, “Adiabatic Hamiltonian Deformation, Linear Response Theory, and Nonequilibrium Molecular Dynamics” in Systems Far from Equilibrium, editor L. Garrido, (Springer-Verlag, Berlin, 1980).
  * (7) K. P. Travis, W. G. Hoover, C. G. Hoover, and A. B. Hass, “What is Liquid? [in two dimensions]”, Computational Methods in Science and Technology 27, 5-23 (2021).
  * (8) Wm. G. Hoover and C. G. Hoover, “Comparison of Very Smooth Cell-Model Trajectories Using Five Symplectic and Two Runge-Kutta Integrators”, Computational Methods in Science and Technology 21, 109-116 (2015).
  * (9) B. L. Holian, W. G. Hoover, and H. A. Posch, “Resolution of Loschmidt’s Paradox: The Origin of Irreversible Behavior in Reversible Atomistic Dynamics”, Physical Review Letters 59, 10-13 (1987).
  * (10) W. G. Hoover, H. A. Posch, B. L. Holian, M. J. Gillan, M. Mareschal, and C. M. Massobrio, “Dissipative Irreversibility from Nosé’s Reversible Mechanics”, Molecular Simulation 1, 79-86 (I987).
  * (11) J. C. Sprott, “Some Simple Chaotic Flows”, Physical Review E 50, R647 (1994).
  * (12) W. G. Hoover, ‘Remark on “Some Simple Chaotic Flows” ’, Physical Review E 51, 759-760 (1995).
  * (13) W. G. Hoover, O. Kum, and H. A. Posch, “Time-Reversible Dissipative Ergodic Maps”, Physical Review E 53, 2123-2129 (1996).

Figure 1:  Steady time-dependent shearing boundary conditions are illustrated
here for unit strainrate $(du_{x}/dy)=1$ with a periodic $4\times 4$ cell
containing Particles 1 and 2. The three upper images of the central cell move
to the right at $u_{x}=+4$ while the three lower images move left with
$u_{x}=-4$. Particles 1 and 2 (red and green) occupy the central cell with a
fixed center of mass at the origin. The force on Particle 1 comes directly
from Particle 2 or its nearest image (from the eight neighboring cells) and is
nonzero whenever that separation is less than the range of the forces,
$r_{12}<2$. The force on Particle 2 is equal in magnitude and opposite in
sign, maintaining the center-of-mass position at the origin, the center of the
central cell. Whenever a horizontal cell boundary is crossed the systematic
horizontal velocity component $\dot{x}$ changes by $\pm 4$ and the $x$
coordinate undergoes a change equal to the fractional part of the strain,
$\dot{\epsilon}t-[\dot{\epsilon}t]_{\rm integer}$, where $t$ is the time. The
crossing of a vertical cell boundary changes only the value of the $x$
coordinate, by $\pm 4$, just as in the usual application of periodic boundary
conditions. To avoid the unnecessary loss of significant figures whenever the
strain reaches $+0.5$, as shown at the right, the strain is decreased by unity
to $-0.5$, the strain shown at the left.

Figure 2:  Pair potential for molecular dynamics simulation of shear,
$\phi(r)=(2-r)^{6}-2(2-r)^{3}$. The minimum at $r=1$ sets the energy scale and
the distance scale. The vanishing of both the force and its first derivative
at $r=2$ provides highly-accurate Runge-Kutta trajectories.

Figure 3:  At the left we see the steady growth of energy error for the
harmonic oscillator, with initial conditions $(q,p)=(1,0)$, shown for forty
time steps of $\pi/10$. The decay of the energy is accurately linear in time
with the error varying as $dt^{5}$ at a fixed time. The analytic solution to
this problem is given in Wm. G. Hoover’s 1991 Elsevier book, Computational
Statistical Mechanics. At the right we see the analogous error growth for the
Hamiltonian dynamics of two two-dimensional particles with periodic boundary
conditions for 200,000 steps of 0.005 each. Though there are collisional
fluctuations the growth rate is still roughly linear and varies likewise as
$dt^{5}$.

Figure 4:  Comparison of the time-dependent error in $q^{2}+p^{2}$ for the
harmonic oscillator (below) and the periodic two-body equilibrium problem
(above), with an initial energy equal to unity in the latter case. The
analytic oscillator solution shows that the energy error at a fixed time
varies as $dt^{5}$, which the numerical work for two bodies nicely confirms.
Although the two-body error increases less regularly, being dominated by the
collisions, it is clear that the same error analysis applies in that case.
Replicating these results is a useful debugging exercise for both these
equilibrium Hamiltonian problems.  Figure 5:  200,000 points along the
equilibrium Hamiltonian $(x,y)$ trajectories for Particle 1 (at the left) and
Particle 2 (at the right). Notice the inversion symmetry resulting from the
center-of-mass constraints described in the text. The total energy of unity is
conserved here. The eight neighboring cells in this case are arranged in a
fixed square lattice corresponding to zero strain. Initially the two particles
are separated as in Figure 1 with $r=2$, $p_{x}=\mp 0.8$, and $p_{y}=\mp 0.6$,
giving $E=1$. The Runge-Kutta timestep is 0.005.  Figure 6:  Particle
locations for the fixed points at $p^{2}$-based kinetic energies from 0.00 to
0.25 in steps of 0.05 with strainrate unity. At a fixed point the horizontal
motion vanishes so that $\dot{x}=p_{x}+y=0$. For Particle 1 the $y$ coordinate
is $-\sqrt{K}$, precisely cancelling the horizontal momentum. As the vertical
momentum vanishes too $p_{x}^{2}=K=2T$. Convergence at $T=0$ is poor while at
kinetic energies from 0.05 to 0.25 convergence to machine accuracy is rapid
and accurate. The timestep is 0.005 and the friction coefficient obeys the
motion equation $\dot{\zeta}=K-2T$. Figure 7:  $(x_{1},y_{1})=(-x_{2},-y_{2})$
coordinates at unit shear strain rate with $\langle\ K\ \rangle=0.3$ and
$dt=0.005$. The motion shown here covers a simulation from time 50,000 to time
100,000 as represented by one million equally-spaced points.

Figure 8:  $(x_{1},y_{1})=(-x_{2},-y_{2})$ coordinates at integral (right) and
half-integral (left) shear strains from the last half of a 200,000,000
timestep simulation at unit strainrate and averaged kinetic energy 0.3. One
million points are shown.

Figure 9:  $(p_{x_{1}},p_{y_{1}})=(-p_{x_{2}},-p_{y_{2}})$ momenta at integral
(right) and half-integral (left) strains from shear strains beginning at
$10^{5}$ and finishing at $10^{6}$. The timestep is 0.005, the averaged
kinetic energy $\langle\ K\ \rangle=0.3=2T$, and the shear strain rate
$\dot{\epsilon}=(du_{x}/dt)$ is unity. $\dot{\zeta}=K-0.3$. Evidently the
velocity distribution is time-periodic and sensitive to the phase of the
boundary cells. Here 900,000 points are shown in each panel, selected from the
time interval from 100,000 to 1,000,000.


# Robust Estimation of Surface Curvature Information from Point Cloud Data

Jared Spang Department of Applied Physics and Applied Mathematics, Columbia
University, New York, NY 10027<EMAIL_ADDRESS>
(May 10, 2023)

###### Abstract

This paper surveys and evaluates some popular state of the art methods for
algorithmic curvature and normal estimation. In addition to surveying existing
methods we also propose a new method for robust curvature estimation and
evaluate it against existing methods thus demonstrating its superiority to
existing methods in the case of significant data noise. Throughout this paper
we are concerned with computation in low dimensional spaces ($N<10$) and
primarily focus on the computation of the Weingarten map and quantities that
may be derived from this; however, the algorithms discussed are theoretically
applicable in any dimension. One thing that is common to all these methods is
their basis in an estimated graph structure. For any of these methods to work
the local geometry of the manifold must be exploited; however, in the case of
point cloud data it is often difficult to discover a robust manifold structure
underlying the data, even in simple cases, which can greatly influence the
results of these algorithms. We hope that in pushing these algorithms to their
limits we are able to discover, and perhaps resolve, many major pitfalls that
may affect potential users and future researchers hoping to improve these
methods.

##  1 Introduction

The robust and accurate estimation of surface curvature information from point
cloud data is a crucial task in many fields such as computer graphics,
robotics, and manufacturing. For example, in the realm of 3D surface analysis
we would ideally work with mesh data; however, oftentimes the only information
available is that of point cloud data from which we must then estimate the
underlying surface / mesh [HKL20]. A particularly salient example
demonstrating the need for estimates associated with point clouds is that of
LiDAR technologly. LiDAR is a popular remote sensing technology that generates
data in the form of point clouds by utilizing light pulses to estimate the
distance to a target surface which is currently being used in a variety of
applications from autonomous vehicle navigation to atmospheric science [LT19,
KBP+04].

Although there are numerous examples of instances in which we lack explicit
surface information, the precise location of highly curved regions can be
difficult to determine from just point cloud data; however, there have been
great strides made in this field of curvature estimation in the last few
decades. A particularly interesting challenge associated with point cloud data
is the inherent uncertainty associated with the data. We might imagine a
scenario in which the point cloud data at hand was scanned from an inherently
noisy sensor and, as a result, there is great uncertainty in the location of
the points themselves. In addition to this common complication there are often
numerous parameters to tune in manifold learning models which effectively
introduce additional noise to the algorithm as there is no inherent ”best”
parameter choice.

As a result, designing an algorithm capable of reliably and robustly
estimating surface curvature from nothing but point cloud data is an
incredibly difficult task. There have been several papers which have sought to
solve this problem; however, few provide strict guarantees on notions such as
error or convergence. It is worth noting that there have been several papers
that take a geometric measure theory approach to solving similar problems
which are typically able to prove convergence; however, many of these results
do not hold in the case of highly irregular surfaces and are still highly
sensitive to the issues discussed [CCSLT08, CSM06, BLM18].

In the field of geometry it is well known that there are a slew of different
notions of curvature; however, here we will be primarily concerning ourselves
with the the computation of the Weingarten map and principle curvatures.
Specifically, we will be implementing, building on, and analyzing the results
of [CLS+21] where we directly estimating the Weingarten map and of [MOG11]
which performs the estimation of curvature related quantities through the
calculation of Voronoi Covariance Measure. The benefit of these methods are
that they claim robustness and are able to perform estimation in arbitrary
dimension. There are of course numerous other methods for computing curvature
(see [MSR07]); however, most of these methods are specific to
$\mathbb{R}^{3}$. It is our hope that the methods discussed are diverse
enough, both methodologically and temporally, to provide an interesting avenue
for comparison. It is the purpose of this paper to evaluate the hardiness of
these methods and work to evaluate some potential improvements on them.

##  2 Background

The underlying assumption we will be working with is one common to the field
of manifold learning; although accumulated data may be presented in a high
dimensional Euclidean space $\mathbb{R}^{N}$, it actually lies within a low
dimensional latent manifold $\mathcal{M}\subseteq\mathbb{R}^{N}$. This
manifold $\mathcal{M}$ is the surface from which we hope to determine our
curvature information. Here we will review the necessary mathematical
background for the understanding of the algorithms to be presented. We start
by defining the Weingarten map and its relationship to the second fundamental
form as well as its connection to the mean and Gaussian curvature. As the
algorithms discussed in this paper are theoretically viable for arbitrary
dimension we present these concepts in their full generality using notions
from differential geometry. Following this we develop the theory surrounding
Voronoi cells and both establish the Voronoi Covariance measure and its
relevant variants. Before we proceed, however, we must establish some
ubiquitous notation. Throughout this paper we let
$\mathcal{M}\subseteq\mathbb{R}^{N}$ be our $m$-dimensional manifold which is
embedded in $N$-dimensional space. Hence, for any point $x_{i}\in\mathcal{M}$
our tangent space $T_{x_{i}}$ is $m$-dimensional. Here we are following the
manifold hypothesis where $\mathcal{M}$ is the ”low” dimensional manifold in
which the data rests.

###  2.1 Curvature Estimation and PCA

When solving the problem of curvature estimation one often first simplifies
the problem to that of estimating normal and tangent spaces. This makes
intuitive sense since, in $\mathbb{R}^{3}$, for any particular point
$x_{i}\in\mathcal{M}$ we may associate a normal vector $n$ to that point.
Since curvature is a metric of how much a surface curves in a different
directions we can estimate the curvature by observing how the surface normal
changes across the surface. For example, if a small perturbation $\varepsilon$
in $x_{i}$ s.t. $x_{i}+\varepsilon\in\mathcal{M}$ yields a large change in the
surface normal then, intuitively, we would hope the surface to be highly
curved at that point. So, one can imagine a valid method of determining
curvature is to compute the surface normals and observe how they change. One
common method of determining the surface normals and tangent space is through
that of principal component analysis (PCA). The basic algorithm can be defined
as follows:

Algorithm 1 PCA Normal and Tangent Estimation

1:function PCA_NT_Est($\varepsilon$, $x_{i}$)

2: $N(x_{i})=\\{x:\|x-x_{i}\|<\varepsilon\\}$

3: $C_{x_{i}}=\sum_{x_{j}\in
N(x_{i})}(x_{j}-\bar{x}_{i})\otimes(x_{j}-\bar{x_{i}})$ $\triangleright$
$\bar{x_{i}}$ is the mean of $x_{j}\in N(x_{i})$

4: $\vec{\lambda},E=\text{Eig}(C_{x_{i}})$ $\triangleright$ $E$ is the matrix
of eigenvectors

5: $n=\\{E_{1},E_{2},\ldots E_{m}\\}$ $\triangleright$ First $m$ eigenvectors

6: $t=\\{E_{m},E_{m+1},\ldots E_{N}\\}$ $\triangleright$ Last $N-m$
eigenvectors

7: return $n$, $t$

8:end function

Here it is worth noting that the eigenvalues $\vec{\lambda}$ are sorted from
smallest to largest so the normal space is eigenvectors corresponding to the
first $m$ smallest eigenvalues and the tangent space consists of the
eigenvectors corresponding the $N-m$ largest eigenvalues. So with this
algorithm we have effectively associated an estimated tangent and normal space
to each point $x_{i}$ which can then be used in the curvature computation.

###  2.2 The Weingarten Map and Curvature

The Weingarten map, also known as the shape operator or second fundamental
form, is a mathematical quantity of focus in this paper. Before we define the
Weingarten map we first endow our manifold $\mathcal{M}$ with some additional
structure; let our manifold $\mathcal{M}$ be a Riemannian manifold
$(\mathcal{M},d)$ endowed with the standard Levi-Civita connection. We then
define $\mathcal{S}$ to be the Weingarten map which serves to associate a
linear map with each point $x_{i}\in\mathcal{M}$ on a surface that details how
the surface curves in the ambient space. Specifically, the Weingarten map
tells us how a normal vector transforms as we move along a chosen tangent
vector. Let $x_{i}$ be a point on $\mathcal{M}$ and $v\in
T_{x_{i}}\mathcal{M}$ then we may more precisely define $\mathcal{S}$ at point
$x_{i}$ as

$\mathcal{S}_{x_{i}}:T_{x_{i}}\mathcal{M}\mapsto T_{x_{i}}\mathcal{M}$ (1)
$\mathcal{S}_{x_{i}}(v)=-\nabla_{v}\zeta(x_{i})$ (2)

where $\zeta$ is a normal vector field to $\mathcal{M}$ and $\nabla_{v}$ is
the covariant derivative. It is also useful to write the Weingarten map
strictly in terms of the connection. For an arbitrary immersion
$\imath:\mathcal{M}\mapsto A^{N}$ which maps an $n$-dimensional manifold
$\mathcal{M}$ into an ambient space $A$ we let $\bar{\nabla}$ be the
connection on $A$ and $\nabla$ be the induced Riemannian connection on $M$.
Then for any tangent fields $X,Y$ to $\mathcal{M}$ we may decompose the
connection in terms of tangential ($\top$) and normal ($\bot$) projections
[Spi05]

$\begin{gathered}\bar{\nabla}_{X}Y=\top\bar{\nabla}_{X}Y+\bot\bar{\nabla}_{X}Y\\\
\bar{\nabla}_{X}Y=\nabla_{X}Y+\mathrm{I\\!I}(X,Y)\end{gathered}$ (3)

where $\mathrm{I\\!I}$ is the second fundamental form which serves to define
the Weingarten map as

$\mathrm{I\\!I}(v,w)=\langle\mathcal{S}(v),w\rangle\zeta=\langle-\nabla_{v}\zeta,w\rangle\zeta$
(4)

So we see the mathematics supports our initial intuition that the change in
the normal vector is essential in defining the Weingarten map and curvature.
In fact, the relation between the Weingarten map and curvature can be made
clear by consider both mean curvature and Gaussian curvature.

In the case of mean curvature we are interested in the average curvature of
the surface at a given point which is given by the average of the principal
curvatures at the desired point. Similar to mean curvature, we can compute the
Gaussian curvature as the product of the principal curvatures at a given
point. Luckily we can easily express both of these curvatures in terms of the
Weingarten map

$\begin{gathered}H_{x_{i}}=\frac{1}{n}\text{tr}(S_{x_{i}})\\\
K_{x_{i}}=\det(S_{x_{i}})\end{gathered}$ (5)

where $H$ is the mean curvature and $K$ is the Gaussian curvature.
Additionally, the directions of principal curvature can be computed as the
eigenvectors of $S_{x_{i}}$. As a result, computation of the Weingarten map
$S_{x_{i}}$ allows us to easily compute both the mean curvature and the
Gaussian curvature.

A final fact of note is the relationship between the Weingarten map and the
Gauss map which, for any manifold of dimension $N-1$, maps
$x_{i}\in\mathcal{M}$ to the unit hypersphere
$\mathbb{S}^{N-1}\subseteq\mathbb{R}^{N}$. Intuitively this map allows us to
associate a normal $\zeta(x_{i})$ with each point $x_{i}\in\mathcal{M}$. Since
we have shown the Weingarten map to be the covariant derivative of the normal
field we immediately see the connection between the Gauss map $g$ and $S$,

$\begin{gathered}g_{x_{i}}:T_{x_{i}}\mathcal{M}\mapsto\mathbb{S}^{N-1}\\\
S_{x_{i}}(X)=-Dg_{x_{i}}(X)=\top(-Dg_{x_{i}}(X))=-\nabla_{X}\zeta(x_{i})\end{gathered}$
(6)

Where $X\in T_{x_{i}}\mathcal{M}$ and so we have explicitly recovered our
equation (2) for the Weingarten map as the differential of the Gauss map.

###  2.3 Voronoi Cells

Let $K\subset\mathbb{R}^{N}$ be a general compact subset of $R^{N}$ then we
define the distance metric $d_{K}$ for any $x\in\mathbb{R}^{K}$ as

$d_{K}(x)=\min_{y\in K}\|x-y\|$ (7)

which just measures the distance to the closest point in $K$ from $x$. Now let
$K\subset\mathbb{R}^{N}$ be a finite set of points (a point cloud) rather than
a general compact set; then we define the Voronoi cell, induced by a
particular norm $\|\cdot\|$, of point $x_{i}$ to be

$\text{Vor}(x_{i})=\\{y:d_{K}(y)=\|x_{i}-y\|\\}$ (8)

which is the set of points in the ambient space $\mathbb{R}^{N}$ that are
closer to $x_{i}\in K$ than to any other $x_{j}\in K$. It is worth noting that
$\\{\text{Vor}(x_{i})\\}$ for all $x_{i}\in K$ forms a partition of the space
$R^{N}$ which is identical to that of the K-NN algorithm as Voronoi cells are
fundamentally solving the same problem. It is also useful to define the medial
axis $\mathbb{M}(K)$ of $K$. The medial axis is the set of points that have
more than one nearest Voronoi point and forms the ”skeleton” of the Voronoi
cell diagram. The key difference between the medial axis and the general
boundary of the Voronoi cell is that the medial axis must be equidistant to
two or more of the Voronoi points.

From this definition of medial axis we also define the projection function
$p_{K}:\mathbb{R}^{N}\setminus\mathbb{M}(K)\mapsto K$ which projects any
$x\in\text{Vor}(x_{i})\subseteq\mathbb{R}^{N}$ to the closest point in $K$ (in
our case $K$ will be our point cloud). The projection function will be used as
in [MOG11] for which the authors note that the function is well behaved almost
everywhere which is sufficient for the use case of (2) as we are interested
purely in the integral of $p_{K}$.

###  2.4 Voronoi Covariance Measure

Here we define the Voronoi Covariance Measure (VCM) which is a tensor-valued
measure defined for any compact set $K\subseteq R^{N}$.

$\mathcal{V}_{K,R}(B)=\int_{K^{R}\cap p_{k}^{-1}(B\cap
K)}(x-p_{K}(x))\otimes(x-p_{K}(x))dx$ (9)

The VCM can be thought of as a generalized covariance matrix since replacing
$p_{K}(x)$ with a point $p$ and adjusting the domain yields the definition of
covariance. Instead of a single point $p$ we are now integrating over the
Voronoi cells of a ”curve” of points. This notion of VCM is good because it
encodes normal information; consider the set
$\\{x_{i}-p:p_{K}(x_{i})=p,x_{i}\in K^{R}\\}$ which takes the set of points
that can be projected to $x_{i}$ and does so, then this set is the normal cone
for the point $x_{i}$. So for a small $\varepsilon$-ball around $x_{i}$ the
VCM is the sum of the covariance matrices (across an infinitesimally small
volume) of the normal cones about that point. Thus, the VCM has a clear
relation with the variation in the normal of a surface which helps us
intuitively establish a connection to the Weingarten map and the curvature of
the surface.

In addition to the standard VCM we can also convolve the measure by a
convolution kernel function $\chi:\mathbb{R}^{N}\mapsto\mathbb{R}^{+}$ which
allows us to locally smooth the VCM in the case of highly noisy observations.
The convolved VCM can be written as,

$\mathcal{V}_{K,R}*\chi(p)=\int_{K^{R}}(x-p_{K}(x))\otimes(x-p_{K}(x))\chi(x-p_{K}(x))dx$
(10)

Typically we are interested in performing a convolution over a ball of radius
$r$ local to the desired point $x_{i}$ and as such we set $\chi$ to be the
indicator function of the $r$-ball which yields
$\mathcal{V}_{K,R}*\chi(p)=\mathcal{V}_{K,R}(B(p,r))$. This provides a
convenient way to add robustness to noise via local smoothing.

As we are not working with continuous sets but rather a discrete point cloud
what we are primarily interested in is the discrete analog of these functions.
Luckily this is quite straightforward and we get the following for the VCM of
a point $x_{i}$ in a point cloud $K$

$\mathcal{V}_{K,R}(\\{x_{i}\\})=\int_{K^{R}\cap
p_{k}^{-1}(x_{i})}(y-x_{i})\otimes(y-x_{i})dy$ (11)

which amounts computing the covariance of $x_{i}$ across its Voronoi cell (and
within the radius $R$). Additionally, we can extend this definition to any set
of discrete points by simply summing them together which leads us to the
Convolved VCM for point clouds,

$\begin{gathered}\mathcal{V}_{K,R}*\chi(\\{x_{i}\\})=\int_{K^{R}}(y-x_{i})\otimes(y-x_{i})\chi(y-x_{i})dy\\\
\mathcal{V}_{K,R}*\chi(\\{x_{i}\\})=\sum_{x_{j}\in B(x_{i},r)\cap
K}\text{cov}(\text{Vor}(x_{j})\cap B(x_{j},R),x_{j})\end{gathered}$ (12)

This provides us with a syntactically simple way to express the ”smoothed”
VCM; however, it can still be difficult to compute the covariance across an
arbitrary Voronoi cell which is why we resort to Monte-Carlo integration in
(2).

###  2.5 Theoretical Guarantees of VCM

A particularly useful facet of VCM based algorithms is that VCM and convolved
VCM are provably robust. Essentially, if we consider two sets $K$ and
$K^{\varepsilon}$ where $K^{\varepsilon}$ is some ”noised” version of $K$ then
if $K$ and $K^{\varepsilon}$ are close in Hausdorff distance then they are
close in convolved VCM. We can make this formal using the results of [MOG11].
Let $d_{H}(X,Y)$ denote the Hausdorff distance between two sets $X$ and $Y$.
If for every $K\subseteq\mathbb{R}^{N}$ and $R>0$,
$\chi:\mathbb{R}^{N}\mapsto\mathbb{R}$ is bounded and $k$-Lipschitz then then
there exists a constant $C$ s.t. for any other
$K^{\varepsilon}\subseteq\mathbb{R}^{N}$

$\|\mathcal{V}_{K,R}*\chi-\mathcal{V}_{K^{\varepsilon},R}*\chi\|_{\infty}\leq
Cd_{H}(K,K^{\varepsilon})^{1/2}$ (13)

Where $C$ is a constant that depends on $N,K$ and $R$. So if the Hausdorff
distance between $K$ and $K^{\varepsilon}$ is $\varepsilon$ we get

$\|\mathcal{V}_{K,R}*\chi-\mathcal{V}_{K^{\varepsilon},R}*\chi\|_{\infty}\leq
C\varepsilon^{1/2}$ (14)

Interpreting this result we see that for a fixed compact set (or point cloud)
$K$, ambient dimension $N$, and $R$ we have that the maximum difference in the
convolved VCM between $K$ and $K^{\varepsilon}$ is bounded by a constant times
the square root of the Hausdorff distance. This essentially let us know that
as the Hausdorff distance between two sets shrinks (e.g. because there is less
noise) so does the difference in their convolved VCM. Although this result
proves that the convolved VCM converges in Hausdorff distance, a common
complaint lodged against this method is that we are not often faced with
”Hausdorff noise” in practice. As a result, some work has been done to
generalize this method in a way that provides more robust convergence
properties in the case of more realistic noise; however, we will not involve
ourselves with these variants in this paper [CLMT15, CCSM11].

##  3 Numerical Estimators

###  3.1 Monte-Carlo VCM Estimator

The Monte-Carlo Voronoi Covariance Measure Estimator (MCVCM) is a randomized
algorithm for estimating the Voronoi Covariance Measure. This algorithm has
been shown to, with high probability, converge to an
$\varepsilon$-approximation of the true VCM when run for
$N=\mathcal{O}(|K|\ln(1/\varepsilon)/\varepsilon^{2})$ iterations. It’s worth
highlighting here that the VCM is a tensor field which, in $\mathbb{R}^{N}$,
associates an $N\times N$ matrix with any point $x\in R^{N}$. So if we wish to
compute the VCM for a point cloud $K$ with $P$ points at every point $p\in K$
we will have a tensorial array of size $P\times N\times N$. We also note that
throughout this paper we perform various calculations using the MCVCM
algorithm (2) and in every experiment we choose $N$ large enough to establish
an estimation within a 5% $\varepsilon$-perturbation of the true solution with
high probability.

Algorithm 2 MCVCM

1:function VCM($K$, $R$, $N$)

2: $V(x)\leftarrow\underline{0}$ $\forall x\in K$ $\triangleright$ Initialize
VCM to zero matrix

3: $M\leftarrow 0$

4: for $1$ to $N$ do

5: $x\leftarrow$ random point in $K$

6: $s\leftarrow$ uniformly sampled random point in $B(x,R)$

7: $k\leftarrow\\#(B(x,R)\cap C)$ $\triangleright$ Number of points within
radius R of x

8: $p^{\prime}\leftarrow$ NearestNeighbor(x,C)

9: $V(p^{\prime})\leftarrow
V(p^{\prime})+\frac{1}{k}(x-p^{\prime})\otimes(x-p^{\prime})$

10: $M\leftarrow M+\frac{1}{k}$

11: end for

12: $V\leftarrow V/M$

13: return $V$

14:end function

###  3.2 WME Estimator

Once we are able to accurately estimate the normals and tangent bundle of our
point cloud via a method such as (1) we can estimate the Weingarten map via
the WME estimator (3). The WME estimator can be derived as follows, first
consider our manifold $\mathcal{M}\subseteq\mathbb{R}^{N}$. Let
$p\in\mathcal{M}$ and $q\in\mathcal{M}$ is a point s.t.
$d_{G}(p,q)<\varepsilon$ where $d_{G}$ is the geodesic distance on
$\mathcal{M}$. As shown in [CLS+21] we have that

$\top(\zeta_{p}-\zeta_{q})=-\mathcal{S}_{p}(\top(p-q))+\mathcal{O}(\|p-q\|^{2})$
(15)

where $\top$ is again the projection to the tangent space. This result is
quite interesting as it allows us to directly model the local change in the
Weingarten map as a change in normals. This relation can be derived by
considering the exponential map $\bf{r}:\mathbb{R}^{N}\mapsto\mathcal{M}$ s.t.
$\bf{r(0)}=p$ (i.e. the geodesic spawning from $p$). Additionally, let
$\bf{u}\in\mathbb{R}^{N}$ then the Taylor expansion of
$\zeta(\bf{r}(u))-\zeta(\bf{r}(0))=\zeta(\bf{r}(u))-\zeta(p)$ results in the
following,

$\top(\zeta(\bf{r}(u))-\zeta(p))=\top(\sum_{i=1}^{N}u^{i}\frac{\partial\zeta}{\partial
u^{i}})+\mathcal{O}(\|p-q\|)^{2}=-\mathcal{S}_{p}(\top(p-q))+\mathcal{O}(\|p-q\|^{2})$
(16)

We would initially expect the higher order term to involve some square of the
exponential map; however, by evaluating the exponential map Taylor expansion
it can be shown that this higher order term is on the same order as the
euclidean distance [CLS+21, MMASC14]. Consider the tangent basis matrix
$E_{x_{i}}=[e_{x_{i}}^{1},\ldots,e_{x_{i}}^{m}]$ then from (16) we have

$(\zeta_{x_{j}}-\zeta_{x_{i}})^{\intercal}E_{x_{i}}=-(x_{j}-x_{i})^{\intercal}E_{x_{i}}S_{x_{i}}$
(17)

Then we hope to find a matrix $S_{x_{i}}$ which minimizes the residual of this
equation,

$S_{x_{i}}\approx\min_{S}\sum_{j=1}^{N}\chi(|x_{j}-x_{i}|<\varepsilon)\|(\zeta_{x_{j}}+\zeta_{x_{i}})^{\intercal}E_{x_{i}}-(x_{j}-x_{i})^{\intercal}E_{x_{i}}S\|^{2}$
(18)

where $\chi(|x_{j}-x_{i}|<\varepsilon)$ is the indicator function for if
$x_{j}$ is in the $\varepsilon$ ball of $x_{i}$. This indicator is necessary
as our approximation only holds for points near $x_{i}$ on the manifold;
however, we could replace it with some other measure of closeness such as the
K nearest neighbors of $x_{i}$. Thus, we have derived an approximation of the
Weingarten map which is solvable as a fairly simple least-square problem for
matrices.

Algorithm 3 WME

1:function WME($\varepsilon$, $K$, $i$)

2: $\zeta,E\leftarrow\text{PCA\\_NT\\_EST}(\varepsilon,K)$ $\triangleright$
Get normal and tangent basis

3: $\Delta\leftarrow E^{\intercal}[x_{1}-x_{i},\ldots,x_{n}-x_{i}]$

4: $\Xi\leftarrow
E^{\intercal}[\zeta_{x_{1}}-\zeta_{x_{i}},\ldots,\zeta_{x_{n}}-\zeta_{x_{i}}]$

5: $W\leftarrow\text{diag}\\{\chi(x_{1}-x_{i}),\ldots,\chi(x_{n}-x_{i})\\}$
$\triangleright$ KNN mask matrix

6: $S=-\Xi W\Delta^{\intercal}(\Delta W\Delta^{\intercal})^{-1}$
$\triangleright$ Closed form to LLS

7: return $S$

8:end function

This algorithm has been theoretically shown to be consistent with bias
$\mathcal{O}(h\varepsilon^{4})$ and variance
$\mathcal{O}(\frac{1}{n\varepsilon^{m}})$. Here we define $\varepsilon$ as the
”bandwidth” of the model, which essentially serves as the size of the
neighborhood surrounding $x_{i}$ we take points from (this is $\varepsilon$ in
(18)). Additionally, we have $m$ as the dimension of the manifold and $n$ as
the number of points. Thus it has been shown that

$MSE=\mathcal{O}(\varepsilon^{4}+\frac{1}{n\varepsilon^{m}})$

By considering the formulas for curvature in (5) we can then easily extend
computation of $S_{x_{i}}$ via (3) to computation of mean and Gaussian
curvature.

###  3.3 VWME

Based on the robustness and superior orientability (see Discussion) associated
with VCM normal estimation compared to that of PCA-based normal estimation we
propose to modify the WME algorithm by augmenting it with the VCM.
Specifically, we replace the PCA-based normal and tangent estimation with the
normal and tangent estimation found by computing the eigenbasis of the VCM (or
locally convolved VCM in this case). We title this variation on the WME
algorithm the Voronoi Weingarten Map Estimator (VWME) and refer to it as such
throughout the remainder of the paper. The exact details are shown in (4)
where we see VWME amounts to this simple mix of VCM and WME algorithm.

Algorithm 4 VWME

1:function VWME($\varepsilon$,$r$,$VCM$,$K$,$i$)

2: $N(x_{i})\leftarrow\\{x:\|x-x_{i}\|<\varepsilon\\}$ $\triangleright$
$\varepsilon$-ball or KNN

3: $CVCM\leftarrow\sum_{x_{j}\in N(x_{i})}VCM(x_{i})$ $\triangleright$ Compute
Convolved VCM

4: $\vec{\lambda},E\leftarrow\text{Eig}(CVCM)$ $\triangleright$ $E$ is the
matrix of eigenvectors

5: $\zeta\leftarrow\\{E_{1},E_{2},\ldots E_{m}\\}$ $\triangleright$ First $m$
eigenvectors

6: $E\leftarrow\\{E_{m},E_{m+1},\ldots E_{N}\\}$ $\triangleright$ Last $N-m$
eigenvectors

7: $\Delta\leftarrow E^{\intercal}[x_{1}-x_{i},\ldots,x_{n}-x_{i}]$

8: $\Xi\leftarrow
E^{\intercal}[\zeta_{x_{1}}-\zeta_{x_{i}},\ldots,\zeta_{x_{n}}-\zeta_{x_{i}}]$

9: $W\leftarrow\text{diag}\\{\chi(x_{1}-x_{i}),\ldots,\chi(x_{n}-x_{i})\\}$
$\triangleright$ KNN mask matrix

10: $S=-\Xi W\Delta^{\intercal}(\Delta W\Delta^{\intercal})^{-1}$
$\triangleright$ Closed form to LLS

11: return $S$

12:end function

##  4 Error Analysis

(a) Torus Sectional

(b) Torus

Figure 1: Applying the WME algorithm (3) to estimate the mean curvature of the
torus (right) and torus sectional (left). This plot demonstrates the baseline
capability of the WME algorithm to estimate curvature on some simple 3D shapes
with no noise and fairly dense point clouds (2000 points in each example).

(a) Torus Sectional

(b) Torus

Figure 2: Mean squared error in the WME algorithm (3) mean curvature
computation of figure 1. Again we have 2000 points and here we have rotated
the view of the to get a better view of high error areas. We see that in the
case of the torus the errors are fairly randomly distributed across the
surfaces; however, in the torus sectional we get localized errors around the
edges of the surface as we may expect.

###  4.1 Normal Estimators

As many algorithms for curvature estimation utilize the change in the normals
for that estimation it is nature to evaluate the estimation of the surface
normals prior to that of the entire algorithm. Luckily for many simples
surfaces in $\mathbb{R}^{3}$ we can construct point clouds with known normal
information which allows us to compare our estimates of the normal basis to
that of our known surfaces. The most common method for estimating normals is
that of the PCA based method detailed in (1) which we also evaluate. In
addition to the PCA estimator one of the primary use cases propounded for the
Voronoi Covariance Measure is as an alternative method for estimating surface
normals. As a result, we thoroughly compare the robustness of these two
algorithms and how they compare to the true normals in known surfaces.

In the case of computing normal vectors we have little interest in the
magnitude of the vectors but rather of its angle as we are concerned with the
proper decomposition of the local space around a point $x$ into normal and
tangent spaces. In order to accommodate this requirement we compute a proxy
for the error (or lack thereof) as the cosine similarity in the computed
normals to that of the true normals. In the case of increasing Gaussian noise
applied to a sectional torus point cloud, we plot both the mean cosine
similarity and the mean absolute cosine similarity in figure 3. Specifically,
to determine the cosine similarity ($S_{C}(A,B)$) between two vectors $A$ and
$B$ we compute $S_{C}(A,B)=\frac{\langle A,B\rangle}{\|A\|\|B\|}$ and the
absolute cosine similarity as $|S_{C}(A,B)|$.

The reason for the inclusion of both the absolute and regular cosine
similarity is that when viewing the distribution of errors for both these
methods, e.g. in figure 5, we see that a large number of errors crop up from
simply orienting the normal incorrectly. Merely orienting the normal in the
wrong direction is not usually a huge issue in practice as the vector will
still be orthogonal to the tangent space. In addition, there exist efficient
algorithms for establishing a consistent normal orientation across a surface
using minimal spanning trees [HDD+92].

(a) Cosine Similarity

(b) Absolute Cosine Similarity

Figure 3: Cosine and Absolute Cosine Similarity between ground truth and the
PCA (1) and VCM (2) normal estimators. Here we evaluate the models on the
torus sectional point cloud with 1000 points with increasing additive Gaussian
noise generated with variance $\sigma^{2}$ (x-axis). The parameters for the
convolved VCM estimator are $r=0.2$ and $R=0.5$.

(a) Cosine Similarity

(b) Absolute Cosine Similarity

Figure 4: Cosine and Absolute Cosine Similarity between ground truth and the
PCA (1) and VCM (2) normal estimators. Here we evaluate the models on the
torus sectional point cloud with 1000 points with increasing additive uniform
noise generated in a range $[-\alpha,\alpha]$ for a parameter $\alpha$
(x-axis). The parameters for the convolved VCM estimator are $r=0.2$ and
$R=0.5$.  Figure 5: Here we plot a histogram of the Cosine Similarity between
the PCA and VCM normal estimators under moderate additive Gaussian noise
($\sigma=0.4$). Here we use 50-NN for both PCA and VCM on the torus sectional
point cloud with 1000 points where the convolved VCM uses the parameters
$r=0.2$ and $R=0.5$.

###  4.2 WME Estimator

Although [CLS+21] does perform a series of error analyses through a variation
in the number of sample points on the surface not much is done in the way of
varying the construction of the approximation to the manifold itself. As a
result, we benchmark the performance of this algorithm under various changes
in the type of neighborhood embedding used (K-NN, $\varepsilon$-ball, and
Gaussian kernel) as well as the parameters for each of these embeddings. In
addition, the performance of the algorithm under various different types of
noise is examined for a variety of surfaces.

One thing that is particularly nice about the WME estimator is that, unlike
many estimators used in computer graphics, it is directly applicable to
surfaces of arbitrary dimension. As a result, although we cannot visualize it,
we can evaluate how the algorithm performs on high dimensional shapes such as
the hypersphere. A particular example evaluating the performance of the
algorithm on the hypersphere can be seen in figure 10.

Figure 6: Mean Squared Error (MSE) of Mean Curvature approximation via the WME
algorithm (3) on the torus when varying both point cloud size and variance of
the additive Gaussian noise. For consistency, in every MSE Mean Curvature
approximation the PCA (1) kNN was computed with 50-NN and the WME-$\chi$ was
computed with 30-NN. Figure 7: Mean Squared Error (MSE) of Mean Curvature
approximation via the WME algorithm (3) on the torus when varying both point
cloud size and range of additive uniform noise. For consistency, in every MSE
Mean Curvature approximation the PCA (1) kNN was computed with 50-NN and the
WME-$\chi$ was computed with 30-NN.

(a) Noiseless

(b) Heavy Gaussian Noise ($\sigma=0.7$)

Figure 8: This plot explores the effect on Mean Curvature MSE of modifying the
neighborhood size (chosen by KNN) for a torus point cloud with 2000 points in
both the of significant ($\sigma=0.7$) Gaussian noise and the case of no
noise. The PCA line depicts the MSE as we vary the kNN size of the PCA
algorithm (1) whereas the WME line depicts the MSE as we vary the kNN size of
the WME-$\chi$ parameter. In both cases when varying the neighborhood size the
$k$ parameter the other $k$ value was kept constant at the proposed ”optimal”
value which is $k=\mathcal{O}(\log n)$ for a set of $n$ points [CLS+21].

(a) Noiseless

(b) Heavy Gaussian Noise ($\sigma=0.5$)

Figure 9: Here we plot just the noiseless case explores the effect on Mean
Curvature MSE of modifying the neighborhood size (chosen by KNN) for a torus
point cloud with 2000 points in both the of significant ($\sigma=0.7$)
Gaussian noise and the case of no noise. The PCA line depicts the MSE as we
vary the kNN size of the PCA algorithm (1) whereas the WME line depicts the
MSE as we vary the kNN size of the WME-$\chi$ parameter. In both cases when
varying the neighborhood size the $k$ parameter the other $k$ value was kept
constant at the proposed ”optimal” value which is $k=\mathcal{O}(\log n)$ for
a set of $n$ points [CLS+21].

(a) Noiseless

(b) Heavy Gaussian Noise ($\sigma=0.7$)

Figure 10: Here we plot the MSE of the computation of the Mean Curvature
varying both the PCA neighborhood size and the WME-$\chi$ neighborhood size
both approximated via kNN. Here we plot both the noiseless case and the case
with $\sigma=0.7$ additive Gaussian noise. These are both evaluated across the
$\mathbb{S}^{4}$ unit hypersphere sampled with 3000 points.

###  4.3 VWME Estimator

In addition to performing an error analysis on the WME estimator itself, we
also perform an error analysis on our proposed VWME estimator and compare the
results to that of the base WME algorithm. Specifically, we evaluate the mean
square error in curvature reconstruction across a variety of point cloud sizes
with the introduction of increasing additive noise. We consider the case of
both Gaussian additive noise (fig 12) and uniform additive noise (fig 13) just
as we have done for the standard WME case. This analysis allows us to compare
how robust our proposed VWME algorithm is to noise when compared to the
standard WME algorithm. Not only have we evaluated the VWME algorithm across a
variety of noise levels and point cloud sizes but we also directly compare it
to the WME algorithm in figure 11. Here we consider two cases; a noiseless
case and a noised case; in both cases we directly evaluate the VWME algorithm
against the WME algorithm when varying the $kNN$ neighborhood size for the
normal estimator used. This allows us to directly see how the VWME algorithm
holds up against the WME algorithm on a simple test surface.

(a) Noiseless

(b) Heavy Gaussian Noise ($\sigma=0.5$)

Figure 11:  Here we plot the mean squared error in the calculation of mean
curvature using the VWME algorithm (4). Specifically, we compare the
performance of the VWME algorithm to the standard WME algorithm (3) as the
size of the KNN neighborhood size changes. This comparison is made across a
torus point cloud with 1000 points with a VCM (2) R value of 0.5 and a
constant VWME-$\chi$ KNN mask with 30-NN. Additionally, this progression was
plotted both both the noiseless case and the case of heavy additive Gaussian
noise ($\sigma=0.7$).  Figure 12: Mean Squared Error (MSE) of Mean Curvature
approximation via the VWME algorithm (4) on the torus when varying both point
cloud size and range of additive Gaussian noise. For consistency, in every MSE
Mean Curvature approximation the VCM (2) kNN was computed with 50-NN and the
VWME-$\chi$ was computed with 30-NN. Figure 13: Mean Squared Error (MSE) of
Mean Curvature approximation via the VWME algorithm (4) on the torus when
varying both point cloud size and range of additive uniform noise. For
consistency, in every MSE Mean Curvature approximation the VCM (2) kNN was
computed with 50-NN and the VWME-$\chi$ was computed with 30-NN.

##  5 Discussion

As normal / tangent space estimators are the base for our curvature algorithms
we began by evaluating these algorithms as in figure 3. From this graph we see
that the introduction of noise in the data certainly has a significant effect
on the accuracy (cosine similarity) of the normal estimator. In particular,
here we compare the PCA normal estimator to that of the VCM estimator across
both cosine similarity and absolute cosine similarity. When viewing these
graphs we see that the difference between the two algorithms is much more
stark when just comparing cosine similarity. The reason for this can be seen
more clearly in figure 5 where we plot the distribution of the cosine
similarity for both methods. From this we see that the PCA based estimator
performs comparably to the VCM algorithm; however, it has a much greater
tendency to flip the direction of the normal thus causing it to yield low
overall cosine similarity but decent absolute cosine similarity to the true
solution.

We also tested running these algorithms on graphs constructed both by KNN and
via $\epsilon$-balls; however, when evaluating noised data on the
$\epsilon$-ball graphs a common issue that was encountered was the existence
of ”singleton” points which were not connected to the rest of the graph. If
there is a single singleton point then the algorithms can not be run on that
point because they require local information but when the point is
disconnected from the main graph there is no local information. One possible
way to resolve this would be to generate a more dense sampling of points;
however, this is both not always realistic when using real world data and
becomes computationally intractable rather quickly. As a result, graphs
constructed via KNN over $\epsilon$-balls tend to be superior when one is
concerned with the practical application of these algorithms; however, KNN
based graphs do not have the same theoretical guarantees as $\epsilon$-balls
or kernel graphs (which are even more computationally intractable). One
example of this can be seen in figure 9 where smaller $\varepsilon$ values had
to be cut from the noised tests. This example is particularly interesting
because shows how the addition of noise can directly affect our neighborhood
parameter choice. In graph 9(a) we see a local error minima for the PCA
$\varepsilon$ choice around 1.5; whereas, in 9(b) the error doesn’t flatten
out until around $\varepsilon=2.2$ and remains low afterwards. It is possible
that further increasing $\varepsilon$ could lead to another rise in error;
however, in this particular example the minor radius of the torus is reached
at $\varepsilon=3$ at which point we obtain meaningless results.

In addition to the tests run the WME algorithm we also ran several error
evaluations on our proposed VWME algorithm (4). As we are mostly concerned
with robustness we largely compared the two algorithms in the case of noisy
point clouds. After viewing these results we see that, ostensibly, the VWME
algorithm outperforms the WME in terms of robustness for sufficiently sampled
surfaces ($n\geq 1000$ in the case of figure 12 and figure 13). In fact, as we
increase the point count we may compare our VWME results to that of figure 6
and figure 7 and see that the VWME significantly outperforms the WME algorithm
in terms of MSE on this surface. Of course, all this demonstrates is that the
VWME algorithm has the capability of outperforming the WME algorithm in the
case of significant noise. We should keep in mind that the VWME algorithm has
an additional parameter choice of $R$ in its normal estimation that the WME
algorithm does not have which has the potential to significantly effect the
final result. However, in figure 11 we see that for a fixed $R$ there is a
tendency to outperform the WME algorithm when varying the neighborhood size
for the normal estimation.

##  6 Conclusion

Throughout this paper we have developed the theory for several modern normal
and curvature estimators that are applicable to point clouds of arbitrary
dimension. Following the introduction of relevant theory we presented these
algorithms and evaluated them. In particular, we have subjected both PCA-based
(1) and Voronoi-based (2) curvature estimators to trials of robustness through
varying point cloud sizes, parameter sweeping, and a variety of additive
noise. From this we unveiled the behavior of these algorithms and discovered
several interesting facets of these algorithms including both the propensity
of the PCA estimator to incorrectly orient normals and the finicky nature of
$\varepsilon$-ball neighborhood approximations. In addition to the evaluation
of existing algorithms we also proposed a new variant of the WME algorithm we
call VWME (4) which replaces the PCA-based normal estimation of the WME
algorithm with a more robust Voronoi-based normal estimation. The intuition
behind this new algorithm is that, since the WME algorithm is proven to have
low error using true normals, a more robust estimation of the normal vector
should yield a more robust overall estimator. We have evaluated this VWME
algorithm across a variety of point cloud densities, parameter choices, and
noise values and have shown that it is, in general, more robust to noise than
the standard PCA-based method for curvature estimation.

Although we show this algorithm to be more robust some future work still needs
to be done. In particular, the algorithm should be tested across a larger and
more complex variety of surfaces sampled with a greater number of points. We
currently use a Monte-Carlo variation of the VCM algorithm (2); however, there
exists a much more efficient method for computing the VCM involving 3D
Delaunay triangulation which could be implemented and used to run these tests
[MOG11]. This would be quite beneficial as the VWME algorithm takes
significantly longer to run than the WME algorithm and does not scale as well
due to its reliance on the VCM. We note that although an efficient variation
of VCM exists we stuck with the Monte-Carlo version as, in the spirit of this
paper, it is applicable to arbitrary dimension whereas the efficient version
only applies to surfaces embedded in $\mathbb{R}^{3}$.

## References

  * [BLM18] Blanche Buet, Gian Leonardi, and Simon Masnou.  Discretization and approximation of surfaces using varifolds.  Geometric Flows, 3:28–56, 03 2018.
  * [CCSLT08] Frédéric Chazal, David Cohen-Steiner, André Lieutier, and Boris Thibert.  Stability of curvature measures, 2008.
  * [CCSM11] Frédéric Chazal, David Cohen-Steiner, and Quentin Mérigot.  Geometric inference for probability measures.  Foundations of Computational Mathematics, 11:733–751, 12 2011.
  * [CLMT15] Louis Cuel, Jacques-Olivier Lachaud, Quentin Mé rigot, and Boris Thibert.  Robust geometry estimation using the generalized voronoi covariance measure.  SIAM Journal on Imaging Sciences, 8(2):1293–1314, jan 2015.
  * [CLS+21] Yueqi Cao, Didong Li, Huafei Sun, Amir H. Assadi, and Shiqiang Zhang.  Efficient weingarten map and curvature estimation on manifolds.  Machine Learning, 110(6):1319–1344, may 2021.
  * [CSM06] David Cohen-Steiner and J. M. Morvan.  Second fundamental measure of geometric sets and local approximation of curvatures.  Journal of Differential Geometry, 74, 11 2006.
  * [HDD+92] Hugues Hoppe, Tony DeRose, Tom Duchamp, John McDonald, and Werner Stuetzle.  Surface reconstruction from unorganized points.  In Proceedings of the 19th annual conference on computer graphics and interactive techniques, pages 71–78, 1992.
  * [HKL20] Yuchen He, Sung Ha Kang, and Hao Liu.  Curvature regularized surface reconstruction from point cloud.  CoRR, abs/2001.07884, 2020.
  * [KBP+04] Grady Koch, Bruce Barnes, Mulugeta Petros, Jeffrey Beyon, Farzin Amzajerdian, Jirong yu, Richard Davis, Syed Ismail, Stephanie Vay, Michael Kavaya, and Upendra Singh.  Coherent differential absorption lidar measurements of co2.  Applied optics, 43:5092–9, 10 2004.
  * [LT19] Hazel Si Min Lim and Araz Taeihagh.  Algorithmic decision-making in avs: Understanding ethical and technical concerns for smart cities.  CoRR, abs/1910.13122, 2019.
  * [MMASC14] Maria G Monera, A Montesinos-Amilibia, and Esther Sanabria-Codesal.  The taylor expansion of the exponential map and geometric applications.  Revista de la Real Academia de Ciencias Exactas, Fisicas y Naturales. Serie A. Matematicas, 108:881–906, 2014.
  * [MOG11] Quentin Mérigot, Maks Ovsjanikov, and Leonidas J. Guibas.  Voronoi-based curvature and feature estimation from point clouds.  IEEE Transactions on Visualization and Computer Graphics, 17(6):743–756, 2011.
  * [MSR07] Evgeni Magid, Octavian Soldea, and Ehud Rivlin.  A comparison of gaussian and mean curvature estimation methods on triangular meshes of range image data.  Computer Vision and Image Understanding, 107(3):139–159, 2007.
  * [Spi05] Michael Spivak.  A comprehensive introduction to differential geometry.  Publish or Perish, Inc., 2005.