如何计算由离散数据定义的曲面下的体积？

Question

我需要确定由离散数据点表示的一系列表面下方的体积。在我的数据中，每个样本都存储为数据框列表中的单独数据框。这是一些（小）示例数据：

df1 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6),
                  y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3),
                  z=c(0,2,0,4,6,7,3,2,1,2,7,8,9,4,2))

df2 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6),
                  y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3),
                  z=c(1,1,2,3,5,6,2,1,3,3,8,9,8,3,1))

DF <- list(df1,df2)

类似问题的答案要么使用其他语言（matlab、python），要么答案不包含解决问题的可用脚本（as here）。我可以想到两种可接受的方法来估计每个表面下的体积：1）写出辛普森规则的离散版本，作为 R 中的一个函数，应用于数据帧列表（DF）； 2) 计算 x、y 和 z 之间的任意关系，并使用多元数值积分来找到表面下的体积（使用 pracma 包中的 simpson2d / quad2d 或 cubature 中的 adaptIntegrate 等函数）。

关于第一种方法，复合辛普森规则的公式（我想使用）是here，但由于其复杂性，我未能成功编写一个有效的双重求和函数。在这个表达式中，I(lambda(em) lambda(ex)) 等于上述数据集中每个 x,y 网格点的 z，Delta(em) 和 Delta(ex) 表示 x 和 y 点之间的间隔。

第二种方法本质上将found here 方法扩展为多元样条拟合，然后将预测的 z 值作为积分函数传递。以下是我迄今为止尝试过的这种方法：

require(pracma)

df1.loess <- loess(z ~ x + y, data=DF[[1]])
mod.fun <- function(x,y) predict(df1.loess, newdata=x,y)

simpson2d(mod.fun, x=c(2,6), y=c(1,3))

但这不会产生有用的结果。

实际上，我有一个包含近 100 个单个样本的数据帧的列表，因此我确实需要能够将解决方案表示为一系列 lapply 函数，这些函数可以跨列表中的所有数据帧自动执行这些计算。一个示例如下所示：

require(akima)
DF.splines <- lapply(DF, function(x,y,z) interp(x = "x", y = "y", z = "z",
                                                linear=F, nx=4, ny=2))

不幸的是，这会产生缺失值和 Inf 的异常。对于如何成功实施这些策略之一或使用不同（更简单？）方法的任何建议，我都非常愿意接受。克里金函数（如 DiceKriging 包中的 km）能否产生更好的拟合，可以传递给数值积分？

Answer 1

我假设体积表面网格是通过直线连接点来定义的。然后您可以通过

找到该表面下方的体积

1. 将(x,y) 网格三角形镶嵌成三角形T_i，面积为A_i
2. 为每个三角形T_i找到对应的z值Z_i
3. 通过V_i=A_i*sum(Z_i)/3计算截断棱镜的体积V_i（由T_i和Z_i定义）（见https://en.wikipedia.org/wiki/Prism_(geometry)和https://math.stackexchange.com/questions/2371139/volume-of-truncated-prism）
4. 总结所有截断棱镜体积V_i

但请记住，体积确实取决于您的细分，并且细分不是唯一的。但是从某种意义上说，您的问题并没有完全定义，因为它没有描述应该如何在点之间进行插值。因此，任何计算体积的方法都必须做出额外的假设。

回到我的解决方法，第 1 点和第 2 点可以通过 geometry 包实现。 这里有一些代码

library(geometry)

getVolume=function(df) {
  #find triangular tesselation of (x,y) grid
  res=delaunayn(as.matrix(df[,-3]),full=TRUE,options="Qz")
  #calulates sum of truncated prism volumes
  sum(mapply(function(triPoints,A) A/3*sum(df[triPoints,"z"]),
             split.data.frame(res$tri,seq_along(res$areas)),
             res$areas))
}

sapply(DF,getVolume)
#[1] 32.50000 30.33333

由于很难检查结果是否一致，这里举一个我们知道正确答案的简单示例。这是一个边长为 2 的立方体，我们在其中沿 x 轴切出一个楔形。裁剪区域为总体积的 1/4。

cutOutCube=expand.grid(c(0,1,2),c(0,1,2))
colnames(cutOutCube)=c("x","y")
cutOutCube$z=ifelse(cutOutCube$x==1,1,2)

sapply(list(cutOutCube),getVolume)
#[1] 6

这是正确的，因为2^3*(1-1/4)=6。

可以通过计算卷 w.r.t. 的“补码”来执行另一个完整性检查。到一个简单的长方体，其中所有z 值都设置为最大z 值（在这两种情况下都是max(z)=9）。两种情况下的简单长方体体积均为 72。不让我们定义补体曲面并总结体积和补体体积

df1c=df1
df1c$z=max(df1c$z)-df1c$z
df2c=df2
df2c$z=max(df2c$z)-df2c$z
DFc=list(df1c,df2c)

sapply(DFc,getVolume)+sapply(DF,getVolume)
#[1] 72 72

所以体积和补体积在这两种情况下都给出了正确的简单长方体体积。

Answer 2

您可以通过 pracma 包中的函数 barylag2d 中实现的“重心拉格朗日”方法来近似曲面。然后，为避免任何矢量化问题，请明确应用高斯正交规则。

library(pracma)

df1 <- data.frame(x=c(2,2,2,3,3,3,4,4,4,5,5,5,6,6,6),
                  y=c(1,2,3,1,2,3,1,2,3,1,2,3,1,2,3),
                  z=c(0,2,0,4,6,7,3,2,1,2,7,8,9,4,2))

# Define the nodes in x- and y-direction
xn <- df1$x[c(1,4,7,10,13)]
yn <- df1$y[1:3]

# Define the matrix representing the function
m1 <- matrix(df1$z, nrow=5, byrow=TRUE)
f <- function(x, y) 
        c(pracma::barylag2d(m1, xn, yn, x, y))

# 32 nodes in integration intervals
n <- 32
xa <- 2; xb <- 6; ya <- 1; yb <- 3

# Apply quadrature rules explicitely
cx <- gaussLegendre(n, xa, xb)
x <- cx$x; wx <- cx$w
cy <- gaussLegendre(n, ya, yb)
y <- cy$x; wy <- cy$w

# Sum weights * values over all nodes
I <- 0
for (i in 1:n) {
  for (j in 1:n) {
    I <- I + wx[i] * wy[j] * f(x[i], y[j])
  }
}
I  # 40.37037

考虑到数据，40 的整数值似乎是合理的。 simpson2d 或 quad2d 在此设置中不起作用。

您可以尝试adaptIntegrate 是否可以与如此定义的函数f 一起使用。