Software to build a mesh of a surface from points on the surface

Question

I have a set of points $(x_i,y_i,u(x_i,y_i))\in\mathbb{R}^3$, $i=1,\dots N$, over a surface $S$ (from experimental data). I need to calculate the integral of a function $F$ over that surface.

If the points were points over a volume I could use some mesh software (tetgen, for example) and build a mesh, and after that calculate everything. My problem is that it is a surface only, so if I try to use some mesh software I am going to calculate a volume...

How can I build, starting with the given points, a 2D (face) mesh, then iterate over each face in order to calculate the integral?

At the moment, I just need some suggestions about how to compute that "2D mesh". If there exists software that generates the mesh that would be perfect.

Answer 1

Is it supposed to be a closed surface or not? If yes Poisson surface reconstruction from VTK library is your best bet and easiest way to construct such surface, see this example here: https://lorensen.github.io/VTKExamples/site/Cxx/Points/PoissonExtractSurface/

But, if it's not a closed surface, your problem is much more difficult to solve, and you need this advancing front surface reconstruction algorithm from CGAL library: https://doc.cgal.org/latest/Advancing_front_surface_reconstruction/index.html

Answer 2

Here is a way using a Delaunay triangulation. It is performed in R with the help of the deldir package.

f <- function(x, y){
  exp(-(x^2+y^2)) # integrate to pi
}

x <- seq(-5, 5, length.out = 100)
y <- seq(-5, 5, length.out = 100)
grd <- transform(expand.grid(x=x, y=y), z = f(x,y)) # data (x_i, y_i, z_i)

library(deldir)
dd <- deldir(grd[["x"]], grd[["y"]], z = grd[["z"]]) # Delaunay
trgls <- triang.list(dd) # extracts all triangles

vol <- function(trgl){ # calculates volume under a triangle
  with(
    trgl, 
    sum(z)*(x[1]*y[2]-x[2]*y[1]+x[2]*y[3]-x[3]*y[2]+x[3]*y[1]-x[1]*y[3])/6
  )
}

volumes <- vapply(trgls, vol, numeric(1L))
sum(volumes)
# result: 3.141593, approx pi!

And you can plot the triangulated surface:

x <- seq(-3, 3, length.out = 20)
y <- seq(-3, 3, length.out = 20)
grd <- transform(expand.grid(x=x, y=y), z = f(x,y))
dd <- deldir(grd[["x"]], grd[["y"]], z = grd[["z"]]) 

library(rgl)
persp3d(dd, front = "lines", back = "lines", col = "blue")
aspect3d(2, 2, 1)