6

rgl パッケージを使用して、データの 3D プロットを作成しています。いくつかの理由 (3D PCA バイプロット) で、ベクトル (矢印付きの線分) が必要です。3Dコーンを矢印の頭にしたいので、私は立ち往生しています。

どういうわけか、私は老人の心を問題の幾何学に巻き込むことができません。たとえば、次のようにベクトルを描画します

segments3d( rbind( c( 0, 0, 0 ), c( 3, 3, 3 ) ) )

つまり、ユーザー座標系の原点から [3,3,3] までのベクトルです。

先端が [3,3,3] の円錐を作成したいと思います。円錐の底面は円で形成できます。xz 平面 (y 平面に垂直) に半径 r の円を描くのは簡単です。

n <- 10
sin.t <- sin( seq( 0, 2 * pi, len= n ) )
cos.t <- cos( seq( 0, 2 * pi, len= n ) )
r <- 0.1 
xv <- x + r * sin.t
yv <- rep( y, n )
zv <- z + r * cos.t

しかし、円がベクトルに対して垂直になるようにこれらの点を変換するにはどうすればよいでしょうか? そして、その中心は先端からベクトル方向に沿って 0.2 ですか? この変換が完了したら、triangles3d関数を使用して三角形を描画します。各三角形は、頂点に 1 つの角があり、円の点内に 2 つの頂点があります。

これは基本的な数学であり、18 歳の私 (または 28 歳の私でさえ) には問題がないことを私は知っています。（魚とは対照的に）任意のフックをいただければ幸いです。

4

1 に答える 1

8

rgl のデモには、cone3d 関数があります。ベースとチップを別々に取ります。いずれにしても、次のようなことができます。

vec=rbind( c( 0, 0, 0 ), c( 3, 3, 3 ) )
segments3d( vec )


cone3d(base=vec[2,]-(vec[1,]+vec[2,]/6), 
     #this makes the head go 1/6th the length of the arrow
       rad=0.5,
       tip=vec[2,],
       col="blue",
       front="lines",
       back="lines")

ここに、cone3d 関数があります。

   cone3d <- function(base=c(0,0,0),tip=c(0,0,1),rad=1,n=30,draw.base=TRUE,qmesh=FALSE,
                    trans = par3d("userMatrix"), ...) {
   ax <- tip-base
   if (missing(trans) && !rgl.cur()) trans <- diag(4)
   ### is there a better way?
   if (ax[1]!=0) {
     p1 <- c(-ax[2]/ax[1],1,0)
     p1 <- p1/sqrt(sum(p1^2))
     if (p1[1]!=0) {
       p2 <- c(-p1[2]/p1[1],1,0)
       p2[3] <- -sum(p2*ax)
       p2 <- p2/sqrt(sum(p2^2))
     } else {
       p2 <- c(0,0,1)
     }
   } else if (ax[2]!=0) {
     p1 <- c(0,-ax[3]/ax[2],1)
     p1 <- p1/sqrt(sum(p1^2))
     if (p1[1]!=0) {
       p2 <- c(0,-p1[3]/p1[2],1)
       p2[3] <- -sum(p2*ax)
       p2 <- p2/sqrt(sum(p2^2))
     } else {
       p2 <- c(1,0,0)
     }
   } else {
     p1 <- c(0,1,0); p2 <- c(1,0,0)
   }
   degvec <- seq(0,2*pi,length=n+1)[-1]
   ecoord2 <- function(theta) {
     base+rad*(cos(theta)*p1+sin(theta)*p2)
   }
   i <- rbind(1:n,c(2:n,1),rep(n+1,n))
   v <- cbind(sapply(degvec,ecoord2),tip)
   if (qmesh) 
     ## minor kluge for quads -- draw tip twice
     i <- rbind(i,rep(n+1,n))
   if (draw.base) {
     v <- cbind(v,base)
     i.x <- rbind(c(2:n,1),1:n,rep(n+2,n))
     if (qmesh)  ## add base twice
       i.x <-  rbind(i.x,rep(n+2,n))
     i <- cbind(i,i.x)
   }
   if (qmesh) v <- rbind(v,rep(1,ncol(v))) ## homogeneous
   if (!qmesh)
     triangles3d(v[1,i],v[2,i],v[3,i],...)
   else
     return(rotate3d(qmesh3d(v,i,material=...), matrix=trans))
 }
于 2013-03-01T20:24:43.873 に答える