3D point density calculation: a C++ program

Source code and sample result: https://github.com/mauroalberti/InterpDensity3D

Report bugs and suggestions to: alberti.m65@gmail.com

In this post, I present a C++ program that computes the density of points in 3D. This program generates a 3D grid representing the smoothed density of points. While most GIS softwares calculate the 2D point density, the same is not true in the 3D case. Below we compare the two calculations: in the 2D case, the z (vertical component) is not taken into account and vertical variations are not resolved. 3D densities evidence vertical changes.

Visualization of point density for seismic events during the Umbria-Marche 1997 sequence: left, 2D calculation, right, 3D calculation with the C++ program (quartic distribution function), visualised with ParaView. Data from INGV.


Two different kernel functions for density calculation are available in this program: the normal distribution and the quartic functions. The CrimeStat software manual describes the basic theory and these functions (chapter 8, CrimeStat III manual). In the normal distribution function, the contribution of a point to the density is equal to (eq. 8.1 in chapter 8 of CrimeStat III manual):

where h is the standard deviation (kernel bandwidth) and d_ij is the distance between the cell centre and the considered point. For the quartic function, we distinguish two cases. For points outside the search radius (i.e., the kernel bandwidth), the contribution is:

while, for points inside the search radius, the formula is:

where h is the search radius (kernel bandwidth) and d_ij is the distance between the cell centre and the considered point.

Higher bandwidth values give rise to smoother results.

The animations below represent the densities starting from the deepest events. Chosen cell sizes are 500 meters and a kernel bandwidth of 5 cells. This bandwidth is used in both the normal distribution and the quartic functions. We can observe that the normal distribution function produces results that are much smoother than in the quartic function. From a seismological point of view, events are more connected at depth, while at shallow depths there are local event concentrations with oblique alignment with respect to the general fault trend (that is NNW-SSE).

3D densities calculated with normal distribution function.

3D densities calculated with quartic function.

Technical details

The input consists in a text file storing x-y-z coordinates of points (with an initial header line), where the x axis points to the East, the y axis to the North, and the z axis points upwards, i.e. following the right-hand rule. Commas separate coordinate values, e.g.:

"x","y","z"
327572,4765486,-5466
328216,4765470,-4870
325835,4766420,-4862
321250,4775017,-2100
332530,4765174,-3059
326421,4769614,-1813
322404,4766816,-5469
325752,4767419,-5453
331115,4759864,-1276
330473,4767684,-913
325913,4769130,-2164

The program requires a file with analysis parameters, i.e., input and output file names, cell size, kernel functions, kernel bandwidth. An example is:

colfiorito_ipo.txt # input 3D point file
densita_ipo_11.vtk # output file name
500 # cell size
qn # kernel functions: n stands for normal distr., q for quartic distr.
5 # kernel bandwidth

The critical analysis parameters are the cell size of the resulting 3D grid, and the kernel bandwidth (expressed as multiples of the cell size) that regulates the smoothness of the resulting density.

Calculation times can be of the order of tens of seconds- minutes. Since this is a first version, it is not optimised for speed or large data numbers.

The result is a VTK file (name provided in the parameter file) that can be visualized with ParaView and VisIt.
Density values are stored within this file.

Source code compiled with g++ under Windows.