Fast Iterative Method - Numpy/Cupy

This repository implements the Fast Iterative Method on tetrahedral domains and triangulated surfaces purely in python both for CPU (numpy) and GPU (cupy). The main focus is however on the GPU implementation, since it can be better exploited for very large domains.

Details

The anisotropic eikonal equation is given by

for given boundary conditions

For a given anisotropic velocity, this can calculate the geodesic distance between a set of and all points on the domain like shown in the figure.

Note that when using multiple , they are not guaranteed to be in the final solution if they are not a valid viscosity solution. A recommended read for more details on the subject is:
Evans, Lawrence C. "Partial differential equations." Graduate studies in mathematics 19.2 (1998).

Installation

The easiest way to install the library is using pip

pip install fim-python[gpu] #GPU version

If you don't have a compatible CUDA GPU, you can install the CPU only version to test the library, but the performance won't be comparable to the GPU version (see Benchmark).

pip install fim-python #CPU version

Usage

The main interface to create a solver object to use is create_fim_solver

from fimpy.solver import create_fim_solver

#Create a FIM solver, by default the GPU solver will be called with the active list
#Set device='cpu' to run on cpu and use_active_list=False to use Jacobi method
fim = create_fim_solver(points, elems, D)

Example

The following code reproduces the above example

import numpy as np
import cupy as cp
from fimpy.solver import create_fim_solver
from scipy.spatial import Delaunay
import matplotlib.pyplot as plt

#Create triangulated points in 2D
x = np.linspace(-1, 1, num=50)
X, Y = np.meshgrid(x, x)
points = np.stack([X, Y], axis=-1).reshape([-1, 2]).astype(np.float32)
elems = Delaunay(points).simplices
elem_centers = np.mean(points[elems], axis=1)

#The domain will have a small spot where movement will be slow
velocity_f = lambda x: (1 / (1 + np.exp(3.5 - 25*np.linalg.norm(x - np.array([[0.33, 0.33]]), axis=-1)**2)))
velocity_p = velocity_f(points) #For plotting
velocity_e = velocity_f(elem_centers) #For computing
D = np.eye(2, dtype=np.float32)[np.newaxis] * velocity_e[..., np.newaxis, np.newaxis] #Isotropic propagation

x0 = np.array([np.argmin(np.linalg.norm(points, axis=-1), axis=0)])
x0_vals = np.array([0.])

#Create a FIM solver, by default the GPU solver will be called with the active list
fim = create_fim_solver(points, elems, D)
phi = fim.comp_fim(x0, x0_vals)

#Plot the data of all points to the given x0 at the center of the domain
fig, axes = plt.subplots(nrows=1, ncols=2, sharey=True)
cont_f1 = axes[0].contourf(X, Y, phi.get().reshape(X.shape))
axes[0].set_title("Distance from center")

cont_f2 = axes[1].contourf(X, Y, velocity_p.reshape(X.shape))
axes[1].set_title("Assumed isotropic velocity")
plt.show()

A general rule of thumb: If you only need to evaluate the eikonal equation once for a mesh, the Jacobi version (use_active_list=False) will probably be quicker since its initial overhead is low. Repeated evaluations with different or favor the active list method for larger meshes.
On the CPU, use_active_list=True outperforms the Jacobi approach for almost all cases.

Documentation

https://fim-python.readthedocs.io/en/latest

Citation

If you find this work useful in your research, please consider citing the paper in the Journal of Open Source Software

@article{grandits_fast_2021,
  doi = {10.21105/joss.03641},
  url = {https://doi.org/10.21105/joss.03641},
  year = {2021},
  publisher = {The Open Journal},
  volume = {6},
  number = {66},
  pages = {3641},
  author = {Thomas Grandits},
  title = {A Fast Iterative Method Python package},
  journal = {Journal of Open Source Software}
}

Benchmark

Below you can see a performance benchmark of the library for tetrahedral domains (cube in ND), triangular surfaces (plane in ND), and line networks (randomly sampled point cloud in the ND cube with successive minimum spanning tree) from left to right. In all cases, was placed in the middle of the domain. The dashed lines show the performance of the implementation using active lists, the solid lines use the Jacobi method (computing all updates in each iteration).

The library works for an arbitrary number of dimensions (manifolds in N-D), but the versions for 2 and 3D received a few optimized kernels that speed up the computations.

The steps to reproduce the benchmarks can be found in the documentation at https://fim-python.readthedocs.io/en/latest/benchmark.html

Contributing

See Contributing for more information on how to contribute.

License

This library is licensed under the GNU Affero General Public License. If you need the library issued under another license for commercial use, you can contact me via e-mail tomdev (at) gmx.net.