--- title: "Volume under surface from points" author: "Stéphane Laurent" date: '2022-02-26' tags: R, maths, geometry, rgl rbloggers: yes output: md_document: variant: markdown preserve_yaml: true html_document: highlight: kate keep_md: no highlighter: pandoc-solarized --- {r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, collapse = TRUE)  ## The problem Suppose you want to get the volume under a surface but you only have some points on this surface. For the illustration, I will take the surface defined by$z = \exp\bigl(-(x^2 + y^2)\bigr)$on the square$[-5, 5] \times [-5, 5]$. Then the volume we're looking for is close to$\pi$(the integral on$[-\infty, +\infty] \times [-\infty, +\infty]$is exactly$\pi$). {r fxy} f <- function(x, y){ exp(-(x*x + y*y)) }  Now let's define a grid on$[-5, 5] \times [-5, 5]$and the value of$z\$ for each point on this grid: {r grid} x <- seq(-5, 5, length.out = 100) y <- seq(-5, 5, length.out = 100) grd <- transform( # data (x_i, y_i, z_i) expand.grid(x = x, y = y), z = f(x, y) )  ## Elevated Delaunay tessellation - using 'deldir' A solution consists in constructing a Delaunay tessellation of the surface and then to sum the volumes under the Delaunay triangles. The **deldir** package allows to construct such a Delaunay tessellation (which I call an *elevated Delaunay tessellation*). {r deldir, message=FALSE} library(deldir) del <- deldir( # Delaunay x = grd[["x"]], y = grd[["y"]], z = grd[["z"]], rw = c(-5, 5, -5, 5), round = FALSE ) trgls <- triang.list(del) # extracts all triangles  The function below calculates the volume under a triangle: {r fun_volume} volume_under_triangle <- function(trgl){ with( trgl, sum(z) * (x[1L]*y[2L] - x[2L]*y[1L] + x[2L]*y[3L] - x[3L]*y[2L] + x[3L]*y[1L] - x[1L]*y[3L]) / 6 ) }  So here is our approximation of the volume: {r volume_approx} volumes <- vapply(trgls, volume_under_triangle, numeric(1L)) sum(volumes)  ## Using 'RCGAL' If you ran the above code, you noticed that the deldir function as well as the triang.list function are a bit slow. My package [RCGAL](https://laustep.github.io/stlahblog/posts/SurfaceReconstruction.html) (not on CRAN) can construct an elevated Delaunay tessellation, and it is faster. {r RCGAL} library(RCGAL) points <- as.matrix(grd) del <- delaunay(points, elevation = TRUE)  You can directly get the volume: {r RCGAL_volume} del[["volume"]]  And you can easily plot the elevated Delaunay tessellation with the help of the **rgl** package: {r rgl_RCGAL, eval=FALSE} mesh <- del[["mesh"]] library(rgl) open3d(windowRect = c(50, 50, 562, 306), zoom = 0.5) aspect3d(1, 1, 3) shade3d(mesh, color = "limegreen") wire3d(mesh)  ![](figures/rgl_elevated_delaunay.png) ## Update 2022-03-02: using 'tessellation' The elevated Delaunay tessellation is now available in my package [tessellation](https://github.com/stla/tessellation). The command to get it is the same as the 'RCGAL' command and the output is similar. {r tessellation} del <- tessellation::delaunay(points, elevation = TRUE) del[["volume"]]  ## Update 2022-03-10: using 'RCDT' The elevated Delaunay triangulation is now available in my package [RCDT](https://github.com/stla/RCDT). {r RCDT} del <- RCDT::delaunay(points, elevation = TRUE) del[["volume"]]  ## Interactive plot with 'deldir' The **deldir** also allows to get an interactive graphic from the elevated Delaunay tessellation. This requires the **rgl** package. I do it below with a less fine grid, otherwise the visualization is not nice (too dense): {r deldir_for_plot} 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) ) del <- deldir( x = grd[["x"]], y = grd[["y"]], z = grd[["z"]], rw = c(-3, 3, -3, 3), round = FALSE )  {r rgl, eval=FALSE} library(rgl) persp3d(del, front = "lines", back = "lines", col = "blue") aspect3d(2, 2, 1)  ![](figures/deldir_elevated_delaunay.gif)