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

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

您可以通过在 pracma 包中的函数 barylag2d 中实现的 "barycentric Lagrangian" 方法来近似曲面。然后，为避免任何矢量化问题，请明确应用高斯求积规则。

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。

Answer 2

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

找到那个表面下面的体积

将 (x,y) 网格三角形细分为面积 A_i

T_i

为每个三角形找到相应的 z 值 Z_i T_i
通过V_i=A_i*sum(Z_i)/3（见https://en.wikipedia.org/wiki/Prism_(geometry) and https://math.stackexchange.com/questions/2371139/volume-of-truncated-prism）[=49=计算截棱柱（由T_i和Z_i定义）的体积V_i ]
汇总所有截断的棱镜体积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 的 "complement" 来执行另一项健全性检查。到一个简单的长方体，其中所有 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

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

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

How to calculate the volume under a surface defined by discrete data?

interpolation

r

spline

numerical-integration