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

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.

• You can explore the Point Cloud Library (PCL). For example, pcl.readthedocs.io/projects/tutorials/en/latest/… . Dec 15 '20 at 0:43
• The TRIGRID function in IDL software (currently owned and maintained by Harris Geospatial) seems like what you could use here. The algorithm is based on Delanay triangulation. Here is the Harris webpage, l3harrisgeospatial.com/docs/TRIGRID.html and here is a webpage with good explanations and example idlcoyote.com/tips/grid_surface.html Dec 15 '20 at 5:52
• Thanks @AloneProgrammer . It is an open surface. I goint to study the algorithms that you suggest me. Dec 18 '20 at 16:44

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

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)