如何计算由离散数据定义的表面下的体积? [英] How to calculate the volume under a surface defined by discrete data?

问题描述

我需要确定由离散数据点表示的一系列表面下的体积.在我的数据中，每个样本都作为单独的数据帧存储在数据帧列表中.这是一些(小的)示例数据:

I need to determine the volume beneath a series of surfaces represented by discrete data points. In my data, each sample is stored as a separate data frame within a list of data frames. Here is some (small) example data:

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)，或者答案不包含用于解决该问题的有用脚本(

Answers to similar questions are either in other languages (matlab, python), or the answers do not contain useable script to address the problem (as here). I can think of two acceptable ways to estimate the volume beneath each surface: 1) write out a discretized version of simpson's rule as a function in R that is applied across the list of data frames (DF); 2) calculate an arbitrary relationship between x, y, and z and use multivariate numerical integration to find the volume under the surface (with functions like simpson2d / quad2d in the package pracma or adaptIntegrate in cubature).

关于第一种方法，复合辛普森规则(我想使用)的公式为

Regarding the first approach, the formula for the composite simpson's rule (that I would like to use) is here, but due to its complexity, I have been unsuccessful in writing a working double summation function. In this expression, I(lambda(em) lambda(ex)) is equal to z in the above datasets at each x,y grid point, and Delta(em) and Delta(ex) represent the interval between x and y points.

第二种方法从本质上扩展了方法在此处找到以进行多元样条拟合，然后将预测的z值作为积分函数传递.到目前为止，这是我尝试过的方法:

The second approach would essentially extend the approach found here to multivariate spline fits and then pass the predicted z values as a function for integration. Here's what I have tried so far for this approach:

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

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

But this does not yield useful results.

实际上，我为单个样本提供了将近100个数据帧的列表，因此我真的需要能够将解决方案表示为一系列lapply函数，这些函数可以自动对列表中的所有数据帧进行这些计算.一个例子看起来像这样:

In reality, I have a list of almost 100 data frames for individual samples, so I really need to be able to express the solution as a series of lapply functions that automate these calculations across all data frames in the list. An example looks something like this:

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)能否产生更好的拟合度，并可以传递给数值积分?

Unfortunately, this produces an exception for missing values and Infs. I'm extremely open to any suggestions for how to successfully implement one of these strategies, or to utilize a different (simpler?) approach. Could a kriging function (like km in the DiceKriging package) produce a better fit that could be passed on for numerical integration?

推荐答案

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

I am assuming that the volume surface mesh is defined by connecting points via straight lines. Then you can find the volume beneath that surface via

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

1. triangular tessellation of the (x,y) grid into triangles T_i with area A_i
2. finding the corresponding z values Z_i for each of the triangles T_i
3. calculating the volume V_i of the truncated prisms (defined by T_i and Z_i) via V_i=A_i*sum(Z_i)/3 (see https://en.wikipedia.org/wiki/Prism_(geometry) and https://math.stackexchange.com/questions/2371139/volume-of-truncated-prism)
4. summing up all truncated prism volumes V_i

但是请记住，音量确实取决于镶嵌，并且镶嵌不是唯一的.但是您的问题并未完全定义，因为它没有描述应如何在点之间进行插值.因此，任何计算体积的方法都必须做出其他假设.

Keep in mind, however, that the volume does depend on your tessellation and that the tessellation is not unique. But your problem is not fully defined in the sense that it does not describe how one should interpolate between points. So any approach to calculate a volume will have to make additional assumptions.

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

Going back to my solution approach, points 1 and 2 can be achieved via the geometry package. Here some code

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.

Since it's hard to check whether the results are consistent, here a simple example where we know the right answer. It's a cube with side length 2 where we have cut out a wedge along the x axis. The cut-out region is 1/4 of the total volume.

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 值均设置为max z 的最大值(在您的情况下， max(z)= 9 ).对于两种情况，简单的长方体体积均为72.不让我们定义补体曲面并总结体积和补体体积

Another sanity check can be performed by calculating the "complement" of the volume w.r.t. to a simple cuboid where all z values are set to the max z value (in your case max(z)=9 in both cases). The simple cuboid volumes are 72 for both of your cases. Not let's define the complement surfaces and sum up volume and complement volume

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

因此，在两种情况下，体积和补体体积都能提供正确的简单长方体体积.

So volume and complement volume give the right simple cuboid volume in both cases.

这篇关于如何计算由离散数据定义的表面下的体积?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！