Deep Eikonal Solvers

TL;DR

This paper introduces a deep learning-based method to solve the Eikonal equation more accurately and efficiently by replacing traditional local solvers with neural networks, applicable to both flat and curved geometries.

Contribution

It presents a novel neural network approach integrated into the fast marching scheme, improving accuracy and generalization for solving the Eikonal equation on various surfaces.

Findings

Smaller errors in numerical solutions

Higher orders of accuracy achieved

Effective on both flat and curved geometries

Abstract

A deep learning approach to numerically approximate the solution to the Eikonal equation is introduced. The proposed method is built on the fast marching scheme which comprises of two components: a local numerical solver and an update scheme. We replace the formulaic local numerical solver with a trained neural network to provide highly accurate estimates of local distances for a variety of different geometries and sampling conditions. Our learning approach generalizes not only to flat Euclidean domains but also to curved surfaces enabled by the incorporation of certain invariant features in the neural network architecture. We show a considerable gain in performance, validated by smaller errors and higher orders of accuracy for the numerical solutions of the Eikonal equation computed on different surfaces The proposed approach leverages the approximation power of neural networks to enhance the performance of numerical algorithms, thereby, connecting the somewhat disparate themes of numerical geometry and learning.

Figures

Equations

$$|\nabla u(x)|$$

$$u(x)$$

$$|{\nabla }_{G}u(x)|$$

$$u(x)$$

$$|\nabla u_{i,j}{|}^2\approx \max (D_{i,j}^{-x}u,-D_{i,j}^{+x}u,0{)}^2+\max (D_{i,j}^{-y}u,-D_{i,j}^{+y}u,0{)}^2.$$

$$\max (D_{i,j}^{-x}u,-D_{i,j}^{+x}u,0{)}^2+\max (D_{i,j}^{-y}u,-D_{i,j}^{+y}u,0{)}^2$$

$$q\in \text{mbox}Visited\text{if}{u}_{gt}(q)<u_{gt}(p).$$

$$\text{If }q\notin \text{mbox}Visited;u_{gt}(q)=C$$

$$L(\Theta)=\underset{i}{\sum}\Big{(}f_{\Theta}\big{(}p_{i1},u_{gt}(p_{i1}),\ldots,p_{iM},u_{gt}(p_{iM})\big{)}-u_{gt}(p_{i})\Big{)}^{2},$$

$$u_{gt}(p)=\min \{\Vert p-s\Vert ,s\in \mathfrak{s}\}.$$

$$ϵ(h)$$

$$\log (ϵ)$$

$$V_{i}=\big{(}x(p_{i})-x(p),y(p_{i})-y(p),z(p_{i})-z(p),u(p_{i})\big{)}.$$

Full text

1 Introduction

Fast and accurate computation of distances is fundamental to innumerable problems in computer science. Specifically, in the sub-fields of computer vision, geometry processing and robotics, distance computation is imperative for applications like navigating robots [14, 10, 13], video object segmentation [28], image segmentation [3] and shape matching [30]. A distance function has a gradient with unit magnitude at every point in the domain, therefore, they are computed by estimating the viscosity solutions to a Hamilton-Jacobi type non-linear PDE called the Eikonal equation. Fast marching methods are the prominent numerical schemes to estimate distance functions on discretized Euclidean domains as well as triangulated curved surfaces. Fast marching methods are known to have a $\mathcal{O}(N\log N)$ computational complexity (for a discrete domain consist in $N$ sample points). Different versions of the fast marching algorithm have been developed yielding different accuracies $\mathcal{O}(h)$ [22], $\mathcal{O}(h^{2})$ [23] and $\mathcal{O}(h^{3})$ [1] for 2D Cartesian grids and $\mathcal{O}(h)$ for triangulated surfaces [12], where $h$ is the resolution of the discretization. A fast marching method comprises of a local numerical solver and an update step. The local numerical solver approximates the gradient locally to estimate the distance of a point in the advancing unit speed wavefront. Similar to the strategy employed in the celebrated Dijkstra’s algorithm [6], the update step involves selecting the least distant point in order to propagate the wavefront. The accuracy of fast marching scheme relies heavily on the local solver which is responsible for approximations to the gradient at the point of the advancing front [26, 22].

The success of deep learning has shown that neural networks can be utilized as powerful function approximators of complex attributes governing various visual and auditory phenomena. The availability of large amounts of data and computational power, coupled with parallel streaming architectures and improved optimization techniques, have led to computational frameworks that efficiently exploit their representational power. However, the utility of neural networks to provide accurate solutions to complex PDE’s is still a very nascent topic in computational research. Notable efforts include methods like [24] which attempt to solve high-dimensional parabolic partial differential equations, whereas [18] demonstrates the learning of differential operators from data as convolution kernels.

The main contribution of this paper is to develop a deep-learning infrastructure that enables accurate distance computation. Using the identical computational schematic of the fast marching method, the basic premise of this paper is that one can learn the local numerical solver by training a neural network and use it in conjunction with the upwind update scheme. Experiments in sections 3.2 and 4.1 demonstrate lower errors and higher order of accuracy in the solutions. Importantly, we show that our method generalizes to different geometries and sampling conditions.

2 Background

2.1 The Eikonal Equation

Given a domain $\Omega\in\mathbb{R}^{n}$ and a curve or surface $\Gamma\subset\Omega$, we wish to find the distance, denoted by $u(x)$, for each $x\in\Omega$ from the given $\Gamma$. The distance function satisfies

[TABLE]

where $\nabla$ refers to the gradient operator. Equation (1) can be generalized to arbitrary boundary conditions $u(x)=g(x),\,\,x\in\Gamma$, and also to weighted domains where the local metric at a point $x$ is defined by a positive scalar function $F:\Omega\rightarrow\mathbb{R}^{+}$. As shown in [4, 8, 9], the viscosity solution for Equation (1) is unique.

Next, we consider curved manifolds, specifically, surfaces embedded in $\mathbb{R}^{3}$. Let ${\cal M}$ be a Riemannian manifold with metric tensor $G$. The geodesic distance function $u:\mathcal{M}\rightarrow\mathbb{R}^{+}$ satisfies

[TABLE]

where $\nabla_{G}$ represents the gradient operator defined on the Riemannian manifold.

2.2 Numerical Approximation

Consider the task of numerically approximating the function $u(x)$ in Equation (1). Let $u_{i,j}\equiv u(ih,jh)$, where $h$ defines the distance between neighboring grid points along the $x$ and $y$ axis. The following numerical approximation is shown [19, 21] to pick the unique viscosity solution of Equation (1),

[TABLE]

Where $\{D^{+x}_{i,j},D^{-x}_{i,j}\},\{D^{+y}_{i,j},D^{-y}_{i,j}\}$ are the forward and backward difference operators for point $(i,j)$ along the $x$ and $y$ directions, respectively. Based on approximation (5), the solution of $u_{i,j}$, given the $u$ values at its four neighbors $\{u_{i+1,j},u_{i-1,j},u_{i,j+1},u_{i,j-1}\}$ requires a solution of the quadratic equation

[TABLE]

Typically, the approximation in Equation (5) is based on first order Taylor series expansion of $u_{i,j}$. However, more sophisticated approximations such as those using two points [23] or three points [1] away from $(i,j)$ in each direction, or using diagonal neighbors [7] yield more accurate solutions. Unlike Euclidean domains which enjoy the benefit of regular sampling, developing local approximations similar to Equation (5) for curved manifolds are evidently more challenging due to the lack of a universal regular sampling strategy.

2.3 Fast Eikonal Solvers

Fast Eikonal solvers estimate distance functions in linear or quasilinear time. They typically comprise of two computational parts: a local numerical solver and an ordering/update method. The local solver provides an approximation for the distance function at a desired point and the ordering method chooses the next point to approximate. Fast marching methods [26, 22] is the most prominent subclass of fast Eikonal solvers that use approximations similar to Equation (5) for computing the distance function at a point using the known distances of its neighbors. By employing Dijkstra’s ordering strategy, fast marching schemes display an upwind nature, where the solution can be seen to grow outward from the source, equivalent to the propagation of a unit speed wavefront (See Figure 2.2). The other prominent subclass of fast Eikonal solvers are the fast sweeping methods [32, 33, 11, 29], iterative schemes that use alternating sweeping ordering. Some variants [17] of fast sweeping establish their local solver on finite element techniques, however most of them use the same upwind difference local solver as fast marching.

Kimmel and Sethian [12] extended the fast marching scheme to approximate geodesic distances on triangulated surfaces by treating the sampled manifold as piecewise planar. In order to estimate $u(p)$ for some mesh point $p$, each triangle involving $p$ is processed independently. Given a triangle comprising of the points $p_{1},p_{2},p_{3}$, the distance $u(p_{1})$ is estimated based on $u(p_{2}),u(p_{3})$ by approximating Equation (1) on the triangle’s plane. The minimum of the solutions computed on all the triangles involving $p$ is chosen as the final distance assigned to $p$. This algorithm is analogous to the first-order fast marching method that operates on regular grids and therefore has an accuracy $\mathcal{O}(h)$.

3 Deep Eikonal Solver

We present a fast Eikonal solver estimating the geodesic distance for a discrete set of points $\mathcal{S}$ from a subset of source points $\mathfrak{s}\subset\mathcal{S}$. The local solver (step 6 in Algorithm 2.2) in our scheme is represented by a neural network. The process of choosing subsequent points for evaluation is similar to Dijkstra’s algorithm [6], according to which each point $p\in\mathcal{S}$ is tagged into one of three disjoint sets,

Visited: where $u(p)$ is already computed (Black points in Figure 2.2). 2. 2.

WaveFront: where $u(p)$ calculation is in process. (Red points in Figure 2.2). 3. 3.

Unvisited: where $u(p)$ is yet to be determined (Green points in Figure 2.2).

Since we use the same update method used in the fast marching scheme, the computational complexity of our method is $\mathcal{O}(N\log N)$. In Section 3.1 we describe the general philosophy behind training the network. In Sections 3.2 and 4.1 we elaborate on the exact procedure and architecture of the network specific to the corresponding domain.

3.1 Training the Local Solver

Our basic premise is that the local solver can be learned using a neural network. The input to the local solver which is in action at a certain point $p\in{\mbox{WaveFront}}$ is the set of points in its neighborhood denoted by $\mathcal{N}(p)=\{p_{1},p_{2},\ldots,p_{M}\}$ and their corresponding distances $\{u(p_{1}),u(p_{2}),\ldots,u(p_{M})\}$. Based on the local geometry of the Visited points in $p$’s neighborhood, an estimate of the local distance is outputted by the local solver. See Figures 2.2, 1b and 1 for a visual schematic.

Following a supervised training procedure, we train the local solver with examples containing the ground truth distances. By ground truth we refer to the most accurate available solution for distance (elaborated further in sections 3.2 and 4). For a given point $p$, denote this distance as $u_{gt}(p)$. The network inputs neighborhood information: $\{p_{1},u_{gt}(p_{1}),\ldots,p_{M},u_{gt}(p_{M})\}$ and is trained to minimize the difference between its output and the corresponding ground truth distance at the point $p:u_{gt}(p)$.

To simulate the propagating unit speed wavefront in various representative scenarios, we employ the following strategy. We develop a variety of different sources and construct a dataset of local neighborhood patches at different locations relative to the sources. We allow a patch point to participate in the prediction of $u(p)$ according to the following rule.

[TABLE]

Therefore, given a point $p$ and a patch of points constituting its neighborhood: $\{p_{1},p_{2},\ldots,p_{M}\}$, the network inputs $\{p_{1},u_{gt}(p_{1}),\ldots,p_{M},u_{gt}(p_{M})\}$, where all points that are declared Visited according to rule 7 input their respective ground truth distance. For points $q\in\mathcal{N}(p)$ that do not satisfy condition 7, we simulate that the wavefront starting from the source arrives at their location later than it arrives to $p$. Hence, such points are not considered Visited and they input a constant value to the network.

[TABLE]

By employing the rules 7 and 8 stated above, we create a rich database of such local patches coupled with the desired outputs. The parameters of the network $\Theta$ are optimized to minimize the MSE loss

[TABLE]

over a dataset comprising of a variety of points $p_{i}$ and their corresponding neighborhoods.

The intuition behind our approach is twofold. First, we expect that using a larger local support leads to a more accurate estimate of the solution to the differential equation. Secondly, by training the network to follow a ground truth distance, we avoid the need to develop a sophisticated formula for locally approximating the gradient with high accuracy. Rather, we expect the network to learn the necessary relationship between the local patch and the ground truth distance. As demonstrated in our results section, our learning-based approach generalizes to different scenarios including different shapes and different datasets.

3.2 Deep Eikonal Solver for Cartesian Grids

We begin by evaluating our scheme on a Euclidean domain sampled regularly from $\Omega=[0,1]^{2}$. As the network’s input patch for point $p:\mathcal{N}(p)$, we choose all the points located within a radius of $2h$ from $p$ (yellow points in Figure 1b). We use a multilayer perceptron architecture (henceforth abbreviated as $mlp$) with four fully connected layers having number of nodes: $mlp(13,128,256,128,1)$ respectively. After each linear stage we applied ReLU function as the non-linearity. Since the grid spacing is uniform in the Cartesian case, $h$ encodes all necessary spatial information regarding neighborhood $\mathcal{N}(p)$. Accordingly, we input only the grid spacing $h$ to the network instead of the explicit coordinates $\{p_{1},p_{2},\ldots,p_{12}\}$.

We trained our network by generating 10,000 synthetic examples constructed as per the strategy enumerated in Section 3.1. We design a variety of different source configurations like points, circles, arbitrary curves etc, and construct different local patches using the ground truth distance from these sources: $\mathfrak{s}$. For the Cartesian case, we evaluate the ground truth distance as

[TABLE]

We pipeline each example by subtracting $\min\{u_{gt}(q),q\in\mathcal{N}(p)\}$ from each distance and by scaling the example to mean magnitude 0.5, thereby simplifying the diversity of inputs prior to its feeding into the network. Note, that the subtraction is equivalent to solving Equation (1) with corresponding initial condition $g^{\prime}(x)=g(x)-\min\{u_{gt}(q),q\in\mathcal{N}(p)\}$ and the scaling is equivalent to solving Equation (1) on correspondingly scaled coordinates. Before proceeding with the update step, the network’s output is re-scaled and the bias is added to achieve the distance estimate.

We test our method using the benchmark shown in [23] comprising of two point sources. The learned scheme shows a superior and more accurate approximation of the distance function compared to fast marching local solvers (See Figure 1b). We measure the order of accuracy $r$ of the proposed scheme by evaluating the error $\epsilon$ as a function of the grid spacing $h$ [15]. Let $u_{gt}$ be the ground truth and $u_{h}$ represent the solution of a numerical scheme on a grid spacing $h$.

[TABLE]

where $C$ is a constant. We evaluated our scheme for a range of grid sizes $h$ and plotted $\log\epsilon$ as a function of $\log h$. The slope of this line gives the order of accuracy $r$. From Figure 3.2 we can see that our local solver has a higher order of accuracy as compared to the classical fast marching local solvers in addition to lower errors.

4 Deep Eikonal Solver for Triangulated Meshes

For triangulated meshes, we chose the second-ring neighbors of a point $p$ to constitute the local patch $\mathcal{N}(p)$ (See Figure 1). In Euclidean domains, the regular sampling allows us to encode the spatial information of a patch using only the grid spacing $h$. However, the unordered nature of the sampled points within a mesh patch demands an explicit description of the spatial position associated with each point in a given patch. For each point $p_{i}$, we construct a vector comprising of four values: The 3D coordinates along with the distance $u(p_{i})$, namely,

[TABLE]

The architecture we employ (see Figure 1) has the following structure. It comprises of $M$ input vectors, each encoding the information about a single point in the neighborhood point of some patch as shown in Figure 1. Each vector $V_{i}$ is processed independently in linear layers $mlp(4,64,128,512,1024)$ which produces a 1024-dimensional feature vector for each point. These 1024-dimensional feature vectors are max-pooled to generate a single vector of the same dimension and finally passed through a $mlp(512,256,1)$ which generates the desired output distance $u(p)$. Our architecture is motivated from previous deep learning frameworks for point clouds like [20]. This scheme of learning functions over non-Euclidean domains was analytically validated in [31].

4.1 Experimental results

We create $100,000$ training examples from 8 TOSCA shapes [2] using the methodology described in Section 3.1. In order to train for rotation invariance, we augment the data by multiplying each example by a random rotation matrix. As a pre-processing stage, we practice the same pipeline described in 3.2.

We compare our Deep Eikonal solver to two different axiomatic approaches: our direct competitor, the fast marching method and recent geodesic approximation, the heat method [5]. The heat method uses the heat equation solutions to estimate geodesics on surfaces, motivated by Varadhan’s formula [27], that connects the two. As shown in Figures 4.1 and 4.1, the Deep Eikonal solver approximates the geodesic distance more accurately than the axiomatic approaches. Fast marching and heat method yield smooth but inaccurate approximations in comparison to the polyhedral scheme [25] which is known to be the most accurate scheme for computing distances on surfaces. Tables 4.1 and 4.1 show quantitative results on TOSCA [2] and SHREC [16] databases respectively. The error at point $p$ was computed as $|\frac{u(p)-u_{gt}(p)}{u_{gt}(p)}|$ where the ground truth is taken as the polyhedral scheme solution. To evaluate our method’s order of accuracy as described in Section 3.2, we use different resolutions of unit sphere, since the geodesic distance on a sphere can be computed analytically. Figure 2 demonstrates that the order of accuracy of the Deep Eikonal solver is comparable with the polyhedral scheme accuracy. Finally, we test the robustness of our method to noise in the sampled manifold. Figure 3 demonstrates that the isolines for our method remain robust to noise in comparison to the fast marching.

5 Conclusion

In this paper, we introduce a new methodology by which we train a neural network to numerically solve the Eikonal equation on regularly and non-regularly sampled domains. We develop a numerical scheme which uses a data-centric learning-based computation in conjunction with a well-established axiomatic update method. In comparison to the axiomatic counterparts, we demonstrate that our hybrid scheme results in a considerable gain in performance measured by lower errors and larger orders of accuracy tested on a variety of different settings. Our approach advocates combining the best of both worlds: approximation power of neural networks combined by meaningful axiomatic numerical computations in order to achieve performance as well as better understanding of computational algorithms.