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

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))

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

Answer 1

我假设体积曲面网格是通过直线连接点定义的。然后你可以通过以下方式找到该表面下方的体积

1. (x,y) 网格的三角形分割为面积为 A_i 的三角形 T_i
2. 为每个三角形 T_i 找到相应的 z 值 Z_i
3. 通过 V_i=A_i*sum(Z_i )/3(参见 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

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