Delaunay Tesselation Field Estimator (DTFE) is a mathematical tool for the reconstruction of the density and velocity field of a discrete point set developed by Willem Schaap and Rien van de Weijgaert (Schaap and Weijgaert 2004, Willem Schaap 2007). The DTFE method dynamically adapts to the variation of the density and geometry of the point set. The method is used in cosmological datasets (see for example the astronomy picture of the day, 7 November 2007) as it can simultaneously capture the geometry of the voids, walls, filaments and clusters and preserve the total mass of the pointset. The total mass is not preserved in the closely related natural neighbor interpolation based on the Voronoi tessellation.
I will here describe the method and provide a Python implementation using the scipy library in two dimensions and three dimensions. The complete code including an example can be found in the repository. For an efficient implementation of the DTFE method in c++ see the code by Marius Cautun.
Consider the pointset \(\mathcal{P}\) consisting of \(N\) labaled points \(\boldsymbol{p}_i \in \mathbb{R}^d\), the velocity set \(\mathcal{V}\) consisting of the \(N\) velocities \(\boldsymbol{v}_i \in \mathbb{R}^d\), and the mass set \(\mathcal{M}\) consisting of the masses \(m_i \in \mathbb{R}\) corresponding to the points in \(\mathcal{P}\), with \(i=1,\dots,N\). Using the points, we reconstruct the density field \(\rho:\mathbb{R}^d \to \mathbb{R}\). Using both the points and the velocities, we construct the velocity field \(\boldsymbol{v}:\mathbb{R}^d\to \mathbb{R}\).
Given the pointset \(\mathcal{P}\), consider the Delaunay tesselation \(\mathcal{D}\) consisting of \(N_T\) labelled simplices \(D_i\), with \(i=1,\dots, N_T\). In two dimensions, a simplex is a triangle spanning three points in \(\mathcal{P}\). In three dimensions, a simplex is a tetrahedron spanning four points in \(\mathcal{P}\).
Let’s assume we can associate a (to be determined) density estimate \(\rho_i\) to each point in \(\mathcal{P}\). Given a symplex \(D \in \mathcal{D}\) spanned by the vertices \(p_{l_0},\dots, p_{l_d} \in \mathcal{P}\), with the corresponding densities \(\rho_{l_0}, \dots, \rho_{l_d}\), we construct a linear interpolation of the density field in the simplex
\[\rho(\boldsymbol{x}) = \rho_{l_0} + [\nabla \rho] (\boldsymbol{x}-\boldsymbol{p}_{l_0}),\]
with \(\boldsymbol{x} \in D\) and the gradient vector \(\nabla \rho\in \mathbb{R}^d\) associated to the simplex determined by the \(d\) linear consistency relations \(\rho(\boldsymbol{p}_{l_i}) = \rho_{l_i}\) for \(i=1,\dots,d\). In matrix notation,
\[ \nabla \rho = \begin{pmatrix} \boldsymbol{p}_{l_1}-\boldsymbol{p}_{l_0}\\ \vdots\\ \boldsymbol{p}_{l_d}-\boldsymbol{p}_{l_0}\\ \end{pmatrix}^{-1} \begin{pmatrix} \rho_1-\rho_0\\ \vdots\\ \rho_d-\rho_0\\ \end{pmatrix}. \]
This system of equations is solvable when non of the points \(\boldsymbol{p}_{l_0}, \dots, \boldsymbol{p}_{l_d}\) are collinear. The integral over the linear interpolation yields
\[\int_D \rho(\boldsymbol{x})\mathrm{d}\boldsymbol{x} = \frac{V(D)}{d+1} \sum_{i \in D} \rho_i ,\]
with \(V(D)\) the volume of the simplex \(D\) and \(i\) the labels of the vertices of \(D\) in \(\mathcal{P}\). Note that the volume of a simplex \(D\) spanned by \(\boldsymbol{p}_{l_0}, \dots, \boldsymbol{p}_{l_d}\) can be expressed as the determinant
\[ V(D) = \frac{1}{d!} \begin{vmatrix} \boldsymbol{p}_{l_1}-\boldsymbol{p}_{l_0}\\ \vdots\\ \boldsymbol{p}_{l_d}-\boldsymbol{p}_{l_0}\\ \end{vmatrix}. \]
The integral over density – assuming the density vanishes outside of the convex hull of the Delaunay tesselation – takes the form
\[ \begin{align} \int \rho(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} &= \sum_{i=1}^{N_T} \int_{D_i} \rho(\boldsymbol{x})\mathrm{d}\boldsymbol{x}\\ &= \frac{1}{d+1} \sum_{i=1}^{N_T} V(D_i) \sum_{j \in D_i} \rho_j, \end{align} \]
where the first sum runs over the simplices of the tessellation and the second sum runs over the vertices of a given simplex. Note that \(\rho_i\) enters the sum for each simplex for which it is a vertex. These simplices form the star \(W_i\) of the point \(\boldsymbol{p}_i\)
Using this observation, we reorder the double sum, by collecting the terms involving \(\rho_i\) leading to the terms \(\rho_i(V(D_{l_0}) + \dots + V(D_{l_n})) = \rho_i V(W_i)\), with the \(D_{l_i}\)’s forming the star of \(\boldsymbol{p}_i\). The integral over the density is now a single sum over the points in \(\mathcal{P}\), i.e.,
\[ \int \rho(\boldsymbol{x}) \mathrm{d}\boldsymbol{x} =\frac{1}{d+1} \sum_{i=1}^N \rho_i V(W_i). \]
The key observation in DTFE is that when we chose the natural estimate of the density,
\[ \rho_i = \frac{(d+1) m_i}{V(W_i)}, \]
by which the density at \(\boldsymbol{p}_i\) only depends on the mass at the point and the local geometry, the integral over the density reduces to the total mass of the point set
\[ \int \rho(\boldsymbol{x})\mathrm{d}\boldsymbol{x} = \sum_{i=1}^N m_i \]
as one might expect of a density reconstruction method.
Now, given a point \(\boldsymbol{x} \in \mathbb{R}^d\), we can reconstruct the density field at this point by finding the simplex in which the point lays and evaluating the linear interpolation of the density in the simplex.
The velocity field of the dataset can be analogously reconstructed reconstructed with an linear interpolation in the Delaunay cells. Given a simplex \(D\) spanned by the vertices \(\boldsymbol{p}_{l_0}, \dots, \boldsymbol{p}_{l_d}\) and the associated velocities \(\boldsymbol{v}_{l_0},\dots,\boldsymbol{v}_{l_d}\), we write the velocity field in the simplex \(D\) as
\[ \boldsymbol{v}(\boldsymbol{x}) = \boldsymbol{v}_{l_0} + [\nabla \boldsymbol{v}] (\boldsymbol{x} - \boldsymbol{p}_{l_0}) \]
with the velocity gradient \(\nabla \boldsymbol{v}\) associated to the simplex determined by the linear consistency relations \(\boldsymbol{v}(\boldsymbol{p}_{l_i}) = \boldsymbol{v}_{l_i}\) for \(i=1,\dots,d\). In matrix notation,
\[ \nabla \boldsymbol{v} = \begin{pmatrix} \boldsymbol{p}_{l_1}-\boldsymbol{p}_{l_0}\\ \vdots\\ \boldsymbol{p}_{l_d}-\boldsymbol{p}_{l_0}\\ \end{pmatrix}^{-1} \begin{pmatrix} \boldsymbol{v}_1-\boldsymbol{v}_0\\ \vdots\\ \boldsymbol{v}_d-\boldsymbol{v}_0\\ \end{pmatrix}. \]
Now, given a point \(\boldsymbol{x} \in \mathbb{R}^d\), we can reconstruct the velocity field at this point by finding the simplex in which the point lays and evaluating the corresponding linear interpolation in the simplex.
Note that the gradient \(\nabla \boldsymbol{v}\) is a piecewise constant function. Given the gradient \(\nabla \boldsymbol{v}\) for each simplex, it is natural to evaluate velocity deformation modes. In two dimensions, we evaluate the divergence \(\theta\) and the curl \(\omega\) defined by
\[ \begin{align} \theta &= \nabla \cdot \boldsymbol{v} = \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y},\\ \omega &= \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y}. \end{align} \]
In three dimensions, we evaluate the divergence \(\theta\), the shear \(\sigma_{ij}\) and the vorticity \(\boldsymbol{\omega} = \epsilon^{ijk} \omega_{ij}\) defined by
\[ \begin{align} \theta &= \nabla \cdot \boldsymbol{v} = \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z},\\ \sigma_{ij} &= \frac{1}{2} \left[\frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i}\right] - \frac{1}{3} (\nabla \cdot \boldsymbol{v}) \delta_{ij},\\ \omega_{ij} &=\frac{1}{2} \left[\frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i}\right]. \end{align} \]
with the Kronecker delta \(\delta_{ij}\) and the Levi-Civita symbol \(\epsilon^{ijk}\).
For a point \(\boldsymbol{x} \in \mathbb{R}^d\), the reconstruction amounts to looking up the corresponding simplex and velocity deformation modes.
For simplicity, we will here only discuss the two-dimensional case. The code generalizes straightforwardly to higher dimensions.
#Load the numpy and scipy libraries
import numpy as np
from scipy.spatial import Delaunay
#The Delaunay tesselation
delaunay = Delaunay(points)#Area of a triangle
def volume(sim, points):
return abs(np.linalg.det(np.array([points[sim[1]] - points[sim[0]],
points[sim[2]] - points[sim[0]]]))) / 2
#The density estimate
rho = np.zeros(len(points))
for sim in delaunay.simplices:
vol = volume(sim, points)
for index in sim:
rho[index] += vol
rho = (2 + 1) * m / rho
#The gradients
Drho = np.zeros([len(delaunay.simplices), 2])
Dv = np.zeros([len(delaunay.simplices), 2, 2])
for i in np.arange(len(delaunay.simplices)):
[p0, p1, p2] = points[delaunay.simplices[i]]
[r0, r1, r2] = rho[delaunay.simplices[i]]
[v0, v1, v2] = velocities[delaunay.simplices[i]]
A = np.array([p1 - p0, p2 - p0])
Drho[i] = np.linalg.inv(A) @ np.array([r1 - r0, r2 - r0])
Dv[i] = np.linalg.inv(A) @ np.array([v1 - v0, v2 - v0])def density(x, y, rho, Drho, delaunay):
simplexIndex = delaunay.find_simplex([x, y])
pointIndex = delaunay.simplices[simplexIndex][0]
return rho[pointIndex] + Drho[simplexIndex] @ ([x,y] - delaunay.points[pointIndex])
def v(x, y, velocities, Dv, delaunay):
simplexIndex = delaunay.find_simplex([x, y])
pointIndex = delaunay.simplices[simplexIndex][0]
return velocities[pointIndex] + Dv[simplexIndex] @ ([x, y] - delaunay.points[pointIndex])The velocity deformation modes follow similarly:
def theta(x, y, velocities, Dv, delaunay):
simplexIndex = delaunay.find_simplex([x, y])
return Dv[simplexIndex][0,0] + Dv[simplexIndex][1,1]
def omega(x, y, velocities, Dv, delaunay):
simplexIndex = delaunay.find_simplex([x, y])
return Dv[simplexIndex][1,0] - Dv[simplexIndex][0,1]To illustrate the method, we apply DTFE to a two-dimensional point set consisting of \(256^2\) particles and their velocities (see code). Consider the dataset
leading to the Delaunay tesselation
Using the DTFE method, we reconstruct the density field
