RGL
1. 楕円体を描く
わたしが使っているのは主に planes3d、lines3d で、単に線画を描いているにすぎません。
矩形などの基本図形を描くことも、ライティングなども知りません。
そんな中で地球儀を描こうと思って、楕円体を描くことだけは最初にやってみたことでした。
最初は、簡単に楕円体を描く関数が見つかって、何の苦もなく描けました。しばらくするとエラーが出て描画できなくなり、探すと別のパッケージに同じものがあるのを知って、パッケージ名を冠して使っていました。
しかし、いつの間にかこれも使えなくなっていました。その後、cdaパッケージにあることは分かったのですが現在何故かインストールできません。
私に必要なのは、楕円の一般的に挙げられている式の係数 a。b、c で原点を中心とした楕円が描ければ良いのです。また、回転楕円体で、赤道半径1として計算するので、a=b=1 が描ければ十分です。
止む無く、つなぎと思って以下のコードを使っています。
- ellipsoid <- function (
- x, y, z,
- a = 1, b=1, c=1,
- phi=0, theta=0, psi=0,
- subdivide = 3, smooth = TRUE, ...) {
- if(missing(x)){
- x<-0; y<-0; z<-0
- } else if(missing(y) | missing(z)){
- if(missing(y) & missing(z)){
- if(is.list(x)){ y <- x$y; z <- x$z; x <- x$x
- } else if(is.array(x)){ y <- x[,2]; z <- x[,3]; x <- x[,1] }
- else {warning("引数のy,zは省略されているがxはリストでも配列でもない")}
- } else {warning("引数が揃っていない")}
- }
- sphere <- rgl::subdivision3d(rgl::cube3d(...), subdivide)
- class(sphere) <- c("mesh3d","shape3d")
- norm <- sqrt( sphere$vb[1, ]^2 +
- sphere$vb[2, ]^2 +
- sphere$vb[3, ]^2 )
- for (i in 1:3) sphere$vb[i, ] <- sphere$vb[i, ]/norm
- sphere$vb[4, ] <- 1
- sphere$normals <- sphere$vb
- result <- rgl::scale3d(sphere, a,b,c)
- rotM <- cpp_euler_passive(phi,theta,psi)
- result <- rgl::rotate3d(result,matrix=rotM)
- result <- rgl::translate3d(result, x,y,z)
- invisible(result)
- }
|
C++で書かれている部分の代わりに
- cpp_euler_passive <- function(phi, theta, psi) {
- Rot <- matrix(rep(0,9),nrow=3,ncol=3);
- cosphi = cos(phi); cospsi = cos(psi); costheta = cos(theta);
- sinphi = sin(phi); sinpsi = sin(psi); sintheta = sin(theta);
- Rot[1,1] = cosphi*costheta*cospsi - sinphi*sinpsi;
- Rot[1,2] = sinphi*costheta*cospsi + cosphi*sinpsi;
- Rot[1,3] = -sintheta*cospsi;
- Rot[2,1] = -cosphi*costheta*sinpsi - sinphi*cospsi;
- Rot[2,2] = -sinphi*costheta*sinpsi + cosphi*cospsi;
- Rot[2,3] = sintheta*sinpsi;
- Rot[3,1] = cosphi*sintheta;
- Rot[3,2] = sinphi*sintheta;
- Rot[3,3] = costheta;
- return(Rot);
- }
|
2. 描画されているかどうかの確認
2種類の計算方法があるので、結果が同じかどうか描画して大まかな確認をしようとしました。後に指定した描画が上書きすると思い込んでいて、描画がされていないのでハマりました。
RGLの描画は、重なる部分を描画しないようになっているようです。一旦クリアして、個別に描画してみることが必要です。
3. 平面を描く
planes3d() は便利です。2Dのplot のabline() のように描画領域に合わせて表示してくれます。
一つ不便なのは、XY平面のように座標面を描くと軸が見にくくなることです。軸はアプリケーションで描いた図形の1つなので、その端も planes3d() は描画領域内と見なします。
これを回避するには矩形を描くことができれば良いのだと思います。
