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.


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)

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

volumes <- vapply(trgls, vol, numeric(1L))
# 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"]]) 

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

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.