10

いくつかの平面をプロットしたいと思いますが、それぞれが不等式です。すべての平面をプロットしたら、それらを結合して、これらの線の内側の領域に色を付けたいと思います。多くの 3D 線を描画し、内側の領域を着色するイメージ - これが私がやろうとしていることです。

私のデータは次のようになります。

df <- structure(list(z = c(0, 0.06518, 0.08429, -0.01659, 0, 0.06808, 
0.12383, -1, -0.01662, 0.28782, 0, -0.09539, 0.04255, 0.09539, 
-0.13361, -0.28782, -0.14468, -0.19239, 0.10642), x = c(1, 0.02197, 
0.03503, -0.02494, 0, 0.04138, 0.17992, 0, -0.02482, 0.1122, 
0, 0.01511, 0.0011, -0.01511, -0.06699, -0.1122, -0.06876, 0.12078, 
0.10201), y = c(0, 0.08735, 0.09927, 0.03876, -1, 0.22114, -0.00152, 
0, 0.03811, -0.07335, 0, -0.03025, 0.07681, 0.03025, -0.23922, 
0.07335, -0.25362, -0.09879, 0.05804), value = c(5801L, 135L, 
162L, 109L, 4250L, 655L, 983L, 4500L, 108L, 1594L, 4400L, 540L, 
147L, 323L, 899L, 1023L, 938L, 1627L, 327L)), .Names = c("z", 
"x", "y", "value"), class = "data.frame", row.names = c(NA, -19L
))

各行は、次の形式の方程式を表しますz + x + y < value。xは水平方向の値、yは垂直方向、zは深さです。z は次のように解決できます: -x - y + 値 > z。

座標系の制限は次のとおりです。

x <- z <- seq(-6000, 6000, by = 1)
y <- seq(-4000, 4000, by = 1)

ですから、各行から平面を描きたいと思います。次に、これらすべての平面を組み合わせて、線の内側に値を入力したいと思います。結果は、複数面の不等サイコロのように見えるはずです。または醜いカットのダイヤモンド。

rgl関数と関数の両方で遊んでいますがpersp、どこから始めればよいかわかりません。私は他のソフトウェアの推奨事項を受け入れます。

からの例の1つに頼るpersp3d：

x <- seq(-6000, 6000, by = 1)
z <- x 
y <- seq(-4000, 4000, by = 1)

f <- function(x, y) <- { r <- -x - y + value > z } # stuck here, can you handle an inequality here?
z <- outer(x, y, f)
open3d()
bg3d("white")
material3d(col = "black")
persp3d(x, y, z, col = "lightblue",
        xlab = "X", ylab = "Y", zlab = "z")

私はそれがかなり大きな限界であることを認識しています。それらを減らすのに役立つ場合は、遠慮なく、またはsequence(..., by = ).

4

2 に答える 2

4

行列の乗算を利用することで、計算時間を大幅に節約できます。

library(dplyr)
library(geometry)
library(rgl)

# define point grid
r <- 50   # resolution
grid <- expand.grid(
  x = seq(-6000, 6000, by = r),
  y = seq(-4000, 4000, by = r),
  z = seq(-6000, 6000, by = r))  # data.table::CJ(x,y,z) if speed is a factor

# get points satisfying every inequality
toPlot <- df %>% 
  select(x, y, z) %>% 
  data.matrix %>% 
  `%*%`(t(grid)) %>%
  `<`(df$value) %>% 
  apply(2, all)

## Alternative way to get points (saves time avoiding apply)
toPlot2 <-
  colSums(data.matrix(df[, c('x', 'y', 'z')]) %*% t(grid) < df$value) == nrow(df)

内部は必要ないので、ポイントを凸包に減らし、表面をプロットするだけです。

# get convex hull, print volume
gridPoints <- grid[toPlot, ]
hull <- convhulln(gridPoints, "FA")
hull$vol
#> 285767854167

# plot (option 1: colors as in picture)
apply(hull$hull, 1, function(i) gridPoints[i, ]) %>% 
  lapply(rgl.triangles, alpha = .8, color = gray.colors(5))

## plot (option 2: extract triangles first - much faster avoiding apply)
triangles <- gridPoints[c(t(hull$hull)), ]
rgl.triangles(triangles, alpha=0.8, color=gray.colors(3))

この奇妙な溶けた角氷のようなものを与える：

奇妙なとろけるアイス キューブの事

于 2015-12-18T04:09:59.983 に答える
2

df不等式 ( )の data.frame を受け取り、次のような関数を返す関数ファクトリを作成できます。

  1. 3 つの値 (x、y、z) を取ります
  2. ロードされた不等式のいずれかが成り立たない場合は false を返します

工場：

eval_gen <- function(a,b,c,d){
  force(a); force(b); force(c); force(d)
  check <- function(x,y,z){
    bool <- T
    for (i in 1:length(a)){
      bool <- bool && (a[i] * x + b[i] * y + c[i] * z < d[i])
    }
    return(bool)
  }
}

次に、不等式データフレームを読み込み、関数を作成します。

ueq_test <- eval_gen(df$x,df$y,df$z,df$value) #load the inequalities

あとは、グリッドを作成し、すべての不等式が成り立つ場合に色を付けるだけです。

library(data.table)
library(rgl)

#Note, you can change the resolution by changing the `by` argument here, I've set to 100 to keep computation time and object size manageable

lx <- lz <- seq(-6000, 6000, by = 100)
ly <- seq(-4000, 4000, by = 100)
df_pixels <- data.table(setNames(expand.grid(lx, ly, lz), c("x", "y", "z")))
df_pixels[, Ind := 1:.N]
df_pixels[, Equal := ueq_test(x,y,z), by = Ind]
df_pixels[Equal == T, colour := "red"]

rgl にプロット:

with(df_pixels[Equal == T, ], plot3d(x=x, y=y, z=z, col= colour, type="p", size=5, 
                                     xlim = c(-6000,6000),
                                     ylim = c(-4000,4000),
                                     zlim = c(-6000,6000)
                                     ))

これにより、次のことが得られます。

ここに画像の説明を入力

于 2015-12-17T19:36:17.907 に答える