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)
