R中表面下体积的计算方法

问题描述

我正在尝试计算R中3d表面下的体积。

我的数据dat如下：

          0.003     0.019     0.083      0.25       0.5         1
0     1.0000000 0.8884265 0.8603268 0.7719994 0.7443621 0.6571405
0.111 0.6909722 0.6775000 0.6443750 0.6243750 0.5914730 0.5698242
0.25  0.5847205 0.6022367 0.5572917 0.5432991 0.5170673 0.4835819
0.429 0.5210938 0.5139063 0.4995312 0.4864062 0.4648636 0.4163698
0.667 0.4363103 0.4526562 0.4321859 0.4027519 0.4046011 0.3661616
1     0.3958333 0.4167468 0.3964428 0.3810459 0.3486328 0.3487930

其中x = rownames(dat)，y = colnames(dat)和z = dat。

我看过here，here和here，但是似乎无法弄清楚如何将这些答案应用于我的用例。

这是我数据的可复制版本：

dat = structure(c(1,0.690972222222222,0.584720477386935,0.52109375,0.436310279187817,0.395833333333333,0.888426507537688,0.6775,0.602236675126904,0.51390625,0.45265625,0.416746794871795,0.860326776649746,0.644375,0.557291666666667,0.49953125,0.432185913705584,0.396442819148936,0.771999378109453,0.624375,0.543299129353234,0.48640625,0.402751865671642,0.381045854271357,0.744362113402062,0.591472989949749,0.517067307692308,0.464863578680203,0.404601130653266,0.3486328125,0.657140544041451,0.56982421875,0.483581852791878,0.41636981865285,0.366161616161616,0.348792989417989),.Dim = c(6L,6L),.Dimnames = list(c("0","0.111","0.25","0.429","0.667","1"),c("0.003","0.019","0.083","0.5","1")))

解决方法

如果矩阵是必需的数据框格式，则可以使用链接here的答案中提供的getVolume()函数。

以下是制作该数据框的代码：

df <- expand.grid(x = as.numeric(rownames(dat)),y = as.numeric(colnames(dat)))
df$z = as.vector(dat)

然后定义函数并应用：

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

getVolume(df)

[1] 0.4714882

auc r r surface