rgl パッケージ

R言語の3Dには等高線のようなものからインタラクティブに操作できる高水準のものまであります。おそらくrgl パッケージは最も高レベルで、描画しながら考えるためのツールです。マウス操作でズーミングや回転がストレスなしに可能です。また、スクリプトでインタラクティブに書き加えることも出来ます。

library(rgl)
library(mapdata)
# mapdataの図形を読み込んでXYZ座標系に変換
d <- map("world",plot=F)
xyz <- mapply(XYZ,d$x,d$y)
# 描画
open3d()
plot3d(xyz[1,],xyz[2,],xyz[3,],t="l",
col="blue",axes=F,add=T)

左図の地球儀は mapdataパッケージの "world" データセットを描画したものです。緯度、経度の値から、XYZ直交座標系に変換をして plot3d() を呼び出せば描画できます。XYZ() は関数で「地図、地球儀」にあります。地表を黄色にしたり、赤道の描画もあります。
rlg で描画した地球儀は、マウスホイールでズーミングや回転ができます。
この地球儀は、ほとんど rlg を知らなくても描けますが、山脈を描こうと考えて調べました。

3D表示について

3Dは計算機でも容易に扱えないことが多かったので便法が今でも残っている分野だと思います。逆に計算機の能力不足から実用にされていないが今では利用可能なアルゴリズムもあるのだと思います。
等高線や鳥瞰図も3Dに違いありません。例えば、contour(x,y,z) で等高線を描く場合、z（高さ）の個数は、x の個数と y の個数の積になります。(x,y,z) の組でデータを表しているのではなく、XY平面状の（x,y）で表される格子状に z があります。
rgl の polygon3d() は fill=TRUE （塗り潰し）が指定されると、（x,y）から三角形のメッシュを生成するようです。したがって、XY平面に垂直な面は描画不能なようです。

地球儀を考えると、データは本来、[経度、緯度、半径]からなります。「半径」は良い呼び名が見つからないからで、地心からの距離と言うことです。正確には一般に使用されている地理緯度は地心を基準としていないので地心緯度に変換しての話になります。
しかし、地球の半径に比べで地表の起伏は一般に極端に小さく「半径」として「標高データ」が示されることは実際にはないようです。

3D表示にはプロジェクション（投影）があります。また、サーフェイス（図形の表面）の概念があります。
物の見え方には、遠近法、陰影、光沢と言った要素があります。これらは、１）照明の色や強さ、種類（点光源、面光源）、位置、２）視点の位置、３）サーフェイスの材質として表されます。
light3d() や par3d()によって照明や視点は変えられるようですが、rgl は描画される図形によって自動的に最適な viewport を与えています。図が行方不明にならないのは大変助かります。

3Dでは色は直接的に指定できないことになります。照明の色と向き、表面の材質によって動的に演算されます。たとえば球は何色でもスペキュラと呼ばれる白い部分が描かれます。
たとえば色を指定した線を引いて色が黒なので不思議に思ったりしますが回転してみるとちゃんと色が付いていると言ったことが置きます。影は黒いと言うことです。
面には表裏があります。裏は光が当たらなければ黒ですが、図形を回転して光源の方に向けられます。表と裏の色は別々に存在します。

点、線、面、文字

他を知らなければ rgl は当たり前にしか見えないと思います。しかし過去の作図を考えると rgl が点、線、面、文字を区別していることは画期的です。多くがビットマップを生成するものだったので拡大すると線幅も拡大されました。
逆に言うと着色する箇所は面として描画しなければならないと言うことです。太い線を引いたり、点を打って一見塗り潰したように見せてもだめです。
左図は、同じ図をマウスで回転拡大縮小して見たものです。
文字は拡大縮小の影響を受けません。図形の回転によって位置は移動しますが常に正面を向いています。試していませんが照明の影響も受けないのだと思います。
点は拡大縮小によって大きさが変わりません。線は幅が変わりません。面は拡大縮小に依らず常に塗り潰されます。面の色は照明と視点の位置関係から調整されています。

open3d()
x <- c(-1,4)
y <- c(-1,4)
z <- c(-1,4)
points3d(x,y,z,col="red")
texts3d(x,y,z,texts=c("P1","P2"),adj=1.2)
x <- c(0,2,0, 4,0,4)
y <- c(0,2,0, 2,2,0)
z <- c(0,0,3, 3,3,0)
triangles3d(x,y,z,col="blue")
texts3d(x[2],y[2],z[2],texts=c("blue-tri"))
x <- c(0,1,1,0,0,1)
y <- c(0,0,0.5,0.5,0,0.5)
z <- c(0,0,-0.5,-0.5,0,-0.5)
lines3d(x,y,z,col="red",lwd=10)

描画関数

使い方の推測がつくものを試してみました。
plot3d() は、描画のスーパーセットを意図したものだと思いますが、t(ype) = "l" を指定すると線分を描きます。
lines3d() と同じだと思います。各ポイント間が直線で結ばれます。ポイントはNAを使用して分離して複数の図形を指定することが出来ます。
segments3d() も同様の描画が出来ますが、2ポイントで出来た線分の繰り返しで指定します。

library(rgl)
open3d()
#----------------------------------------------------------------------
# 1x2x3の箱の座標データ
x <- c(0,1,1,0,0,0,1,1,0,0,0,NA,1,1,NA,1,1,NA,0,0)
y <- c(0,0,2,2,0,0,0,2,2,0,0,NA,0,0,NA,2,2,NA,2,2)
z <- c(0,0,0,0,0,3,3,3,3,3,0,NA,0,3,NA,0,3,NA,0,3)
# 赤い箱
plot3d(x,y,z,col="red",t="l",axes=F,
xlim=c(-1,4),ylim=c(-1,4),zlim=c(-1,4))
# 青い箱
lines3d(x+0.5,y+0.5,z+0.5,col="blue")
# 緑の箱
x <- c(0,1, 1,1, 1,0, 0,0, 0,0, 0,1, 1,1, 1,0, 0,0, 0,0, 1,1, 1,1, 0,0)
y <- c(0,0, 0,2, 2,2, 2,0, 0,0, 0,0, 0,2, 2,2, 2,0, 0,0, 0,0, 2,2, 2,2)
z <- c(0,0, 0,0, 0,0, 0,0, 0,3, 3,3, 3,3, 3,3, 3,3, 0,0, 0,3, 0,3, 0,3)
segments3d(x-0.5,y-0.5,z-0.5,col="green")
#----------------------------------------------------------------------
# 赤い面（多角形）
x <- c(0,2,1,0)
y <- c(0,0,1,3)
z <- c(0,4,2,0)
polygon3d(x,y,z,col="red")
# 青い面（三角形）2つ
x <- c(0,2,0, 4,0,4)
y <- c(0,2,0, 2,2,0)
z <- c(0,0,3, 3,3,0)
triangles3d(x,y,z,col="blue")
# 緑の面（矩形）
x <- c(0,2,2,0, 4,0,0,4)
y <- c(0,2,2,0, 2,2,2,0)
z <- c(0,0,3,0, 3,3,0,0)
quads3d(x,y,z,col="green")
#----------------------------------------------------------------------
# 赤い点
x <- c(-1,4)
y <- c(-1,4)
z <- c(-1,4)
points3d(x,y,z,col="red")

points3d() は、矩形のマークが描かれ大きさは size で指定できます。形状については分かりません。

polygon3d() は多角形を描いて塗り潰しを行います。この関数は lines3d() と shade3d() のスーパーセットです。
各ポリゴンは始端と終端を同じになっていない前提で出来ています。図形を閉じる処理は関数内部で行っていて塗り潰しを指定した場合には始端と終端を同じに指定するとエラーになります。
fill = FALSEの時には、1図形の最後に開始点の座標を補って図形を閉じて lines3d() に渡しています。
fill = TRUEの時には、NA で区切って複数図形を指定することは出来ないようです。この処理は x、y から図形を三角形に分割してtmesh3dを作ってshade3d() で描画します。したがってXY平面に垂直な図形はエラーになります。ただし、coords に依存します。coordsは省略すると、1,2となり x、y が使われます。例えば、1、3 を指定すれば x、z で同じことをします。
三角形に分割するのは triangulate() ですが、上手く使えません。左図は地球儀にヒマラヤを描こうとしたときにエラーになったポリゴンです。ヒマラヤは3700ほどのポリゴンになりましたが多数のポリゴンがエラーになって塗り潰すことが出来ませんでした。
triangulate() は tmesh3d() の入力となるインデックスを返します。これは既存のポイントを使用することを意味し、常に成り立つのかどうかは良く分かりません。また、局面を近似する目的なので、ポイントの間隔はあまり長くないことが前提になるのだと思います。

三角形メッシュ、四角形メッシュ

3D 描画では面を三角形のメッシュに分割して塗るのは基本的なことのようです。三点は必ず平面を一義に指定します。
サンプルを見ると矩形のメッシュ qmesh3d() も多く使用されています。
これらは局面が自然に見える程度に分割してくれる関数がありそうに思いましたがそうではないようです。
三角形メッシュにも四角形メッシュにも小さなメッシュを大量に描く仕組みがあるわけではないようです。単に三角形や四角い面を描く関数と考えて良さそうです。

ただし、メッシュを描く関数の位置指定には特徴があります。（x、y、z）で示されるポイントと、その接続関係が分離されていることです。lines3d() などではポイントの並びの順番が接続関係を示しています。メッシュデータではポイントの並びは接続順序とは無関係です。
メッシュ関連の関数の引数ではバーテックス（vertices）と記されていて、頂点の並びの順序は接続順では無いことを示しているようです。
また、バーテックスの各ポイントは（x、y、z、S）と4つの値の組になっています。最後のSが何なのかは良く分かりません。

(0,0),(1,0),(1,1),(0,1)で表されるXY平面状の四角形を triangulate() に入力すると、2組のインデクスが得られます。
[2,3,4] と[4,1,2] で、四角形を2つの3角形として表しています。
ただし triangulate() は2D用の関数で、XY平面と高さのデータを扱う場合には便法として使用できるということだと思います。XY平面に垂直な面が処理できません。それでも polygon3d() で使用されていることから3Dの triangulate() は無いのだろうと思います。

バーテックスの4番目のSの影響はスケーリングのようです。

x <- c(0,1,1,0)
y <- c(0,0,1,1)
z <- c(0,0,0,0)
S <- c(1,1,1,1)
vertices <- rbind(x,y,z,S)
indices <- c(2,3,4, 4,1,2)
mesh <- tmesh3d(vertices, indices=indices)
open3d()
shade3d(mesh,color="green")
mesh <- tmesh3d(points+0.5, indices=tri)
wire3d(mesh,color="red")
mesh <- tmesh3d(points+1, indices=tri)
wire3d(mesh,color="blue")

S=1 として立方体の頂点のバーテックスを作って描画して見ます。

go <- 100; gw <- 100
op <- function(){open3d(windowRect=c(go,100,go+gw,100+gw));go<<-go+gw}
#
vertices <- c( 0,0,0,1, 1,0,0,1, 1,1,0,1, 0,1,0,1,
0,0,1,1, 1,0,1,1, 1,1,1,1, 0,1,1,1)
indices <- c(2,3,4, 4,1,2, 6,7,8, 8,5,6)
mesh <- tmesh3d(vertices,indices)
op(); wire3d(mesh,color="red"); title3d("wire3d")
op(); shapelist3d(mesh,color="blue"); title3d("shapelist3d")
op(); shade3d(mesh,color="green"); title3d("shade3d")
op(); dot3d(mesh,color="blue"); title3d("dot3d")
op(); plot3d(mesh,t="wire",col="red",axes=F,xlab="",ylab="",zlab="")
title3d("plot3d - wire")
op(); plot3d(mesh,t="shade",col="red",axes=F,xlab="",ylab="",zlab="")
title3d("plot3d - shade")
op(); plot3d(mesh,t="dots",col="red",axes=F,xlab="",ylab="",zlab="")
title3d("plot3d - dots")

バーテックスとインデクスから tmesh3d() で作った mesh（三角形メッシュデータ）は、いろいろな描画関数の入力になります。
おそらく全て plot3d() で全て描画できると言うことなのだと思います。plot3d() は自ら open3d() や rgl.clear() も行います。
窓がないと新たな窓を開きます。開いた窓があると、クリアーして描画します。add=T が指定されているとクリアせずに描画します。

メッシュデータは、バーテックスとインデクスが基本的な内容で、material などのプロパティを設定できるようになっているようです。
メッシュデータは特殊なものではないようです。これを生成すれば「曲面」を描画できると考えて良さそうです。

rgl.ellipsoid() は楕円体のメッシュデータを生成します。この関数は cda パッケージにあります。
rgl.ellipsoid(a=1,b=1,c=1) は楕円体の式の係数を全て1に設定することを指定するもので半径1の球と言うことです。
この結果は386のバーテックスと384のインデクスを生じます。立方体の各面を正方形の升目に分割して膨らませたような四角形のメッシュです。経線のように1点に集まるような分割の仕方ではないようです。

四角形が384あって、その頂点は、その4倍ですが、1頂点は4つの矩形で重複すると言うことだと思います。

> ell <- rgl.ellipsoid(a=1,b=1,c=1)
> wire3d(ell)
> str(ell)
List of 5
$ material : list()
$ ib : num [1:4, 1:384] 1 195 99 196 51 197 99 195 ...
$ vb : num [1:4, 1:386] -0.577 -0.577 ...
$ primitivetype: chr "quad"
$ normals : num [1:4, 1:386] -0.408 -0.408 -0.408 ...
- attr(*, "class")= chr [1:2] "mesh3d" "shape3d"

現在の3Dでメッシュは面を表現するアトミックな単位だと思います。大きく拡大するなら最初から細かく用意する必要があります。自動的に細分化することは可能に思えます。

open3d()、rgl.clear()、plot3d()

R言語ではウィンドウ（窓）をデバイスと読んでいます。plot() は windows() で作ったデバイスに図を描きますが、rgl のプロットは windows() で開いたデバイスを使用しません。 rgl は、open3d()を呼び出すか、最初の plot3d() の呼び出しで専用のデバイスを開きます。
rgl のデバイスが開いていると、plot3d() は既存のデバイスを消去（rgl.clear()）して描画します。
plot3d() には add=TRUE と言うオプションが指定でき、この場合は消去を行わず上書きされます。
wire3d() や shade3d() もデバイスがないとデバイスを生成します。既にデバイスがある場合は重ね書きをします。消去して書くことはないようです。
複数のデバイスがある場合、使用されるデバイスはFocus状態のデバイスです。常にデバイスのどれか1つがFocusになっています。
デバイスは番号（シーン）で識別されます。rgl.set(シーン) で、デバイスを Focus に設定できます。
open3d(windowRect=c(x1,y1,x2,y2)) のように指定して、デバイスの位置とサイズを指定できます。R言語のGUIの表示領域の左上隅が(0,0)で、デバイスのウインドウバーを含まない左上隅が合わせられます。したがって（0,0,x2,y2）と指定するとウインドウバーが表示されず閉じることが出来なってしまいます。rgl.close() とすれば閉じることが出来ます。

surface3d

open3d()
#(xの長さ) x (ｙの長さ)
z <- matrix(rnorm(100),ncol=10)
surface3d(1:10, 21:30, z, col="white")

データは（x、y、z）ではなく、x、y は矩形メッシュの刻みを表しています。
z は、x、y の長さの積の個数あります。
この方法は基本的にはXY平面状に高さを表すもので、球体の表面の凹凸と言ったものは描画できないのだと思います。

定番の図形

open3d()
rgl.material(color="yellow")
spheres3d(0,0,0)
# 以下ではrgl.material()が作用しない
open3d()
rgl.spheres(0,0,0)
# 以下の色指定は有効
open3d()
rgl.spheres(0,0,0,color="yellow")

楕円体は、メッシュデータを生成して、それを描画するようです。

open3d()
shade3d(rgl.ellipsoid(a=1,b=1,c=1),col="khaki")
open3d()
x <- matrix(c(1,0,0,0,1,0,0,0,1),3,3)
shade3d(ellipse3d(x),color="pink")

O（オー）オブジェクトは図のような図形を指すようです。これもメッシュデータを作って描画します。

open3d()
shade3d(oh3d(), color="green")

直方体は少し考え方が違うようです。
cube3d()は一辺が2の立方体の基本図形のバーテックスとインデクスを作る関数と考えると良いようです。
このデータから qmesh3d() で四角形メッシュを作って使用します。座標変換は、元のバーテックスに対して行うことも出来ると思いますがメッシュデータに対して演算する関数があります。

> str(cube3d())
List of 6
$ vb : num [1:4, 1:8] -1 -1 -1 1 1 -1 -1 1 -1 1 ...
$ ib : num [1:4, 1:6] 1 3 4 2 3 7 8 4 2 4 ...
$ primitivetype: chr "quad"
$ material : list()
$ normals : NULL
$ texcoords : NULL
- attr(*, "class")= chr [1:2] "mesh3d" "shape3d"
> cube3d()$vb
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] -1 1 -1 1 -1 1 -1 1
[2,] -1 -1 1 1 -1 -1 1 1
[3,] -1 -1 -1 -1 1 1 1 1
[4,] 1 1 1 1 1 1 1 1

回転、配置、大きさを変えるのは以下のようにすれば良いようです。

open3d()
origin <- cube3d() # 基本の立方体
cube <- qmesh3d(origin$vb,origin$ib)
plot3d(rotate3d(cube,pi/4,0,0,1),color="green")
plot3d(translate3d(cube,1,2,3),color="blue",add=T)
plot3d(scale3d(cube,0.5,0.5,3),color="red",add=T)

直方体の表面の色を一色で塗るのは上手くいきます。各面を塗り分けるには工夫がいるようです。
6面あるので6色して見ます。
各面は同じようなグラデーションになりました。どうやら各面とも向かい合う2つの頂点からのグラデーションで4色で1組のようです。各面を1色で塗るには同じ色を4つ指定する必要があるようです。
3面分を赤、青、緑で指定してみます。他の面は黄色をセットします。

open3d()
origin <- cube3d() # 基本の立方体
cube <- qmesh3d(origin$vb,origin$ib)
plot3d(cube,color=c("red","red","red","red",
"blue","blue","blue","blue",
"green","green","green","green",
rep("yellow",4*3))

図の手前の頂点は (1,1,0) です。底に当たるXY平面の部分が赤になりました。
色指定の2番目と4番目の4つの組は、XY平面に垂直な4つ面のうちの2つでした。

以下のようにして色を指定することも出来ます。結果から最後の4つの緑が上の面の色でした。

open3d()
origin <- cube3d(color=c(rep("yellow",4*3),
"red","red","red","red",
"blue","blue","blue","blue",
"green","green","green","green"))
plot3d(origin)

この例は qmesh3d() を使用していません。qmesh3d() を呼び出すことは必要ではないようです。
各段階で座標変換マトリックスが指定できますが、cube3d() で設定した値は、qmesh3d()の出力に反映させる必要があります。

open3d()
origin <- cube3d(color=c(rep("yellow",4*3),
"red","red","red","red",
"blue","blue","blue","blue",
"green","green","green","green"))
cube <- qmesh3d(origin$vb,origin$ib)
cube$material$color <- origin$material$color
plot3d(cube)

座標変換マトリックス

前述のように cube3d() の出力やメッシュデータの回転、配置、大きさを変換する関数があります。
マトリクスを生成する基本的な関数を使うことも出来ます。

m <- scaleMatrix(2, 3, 4) %*%
rotationMatrix(pi/4, 0, 0, 1) %*%
translationMatrix(2, 3, 4)
# -1,1の立方体をZ軸に45度回転し辺を2,3,4に拡大し、
# 2,3,4移動すると
cube <- cube3d(trans=m,color="yellow")
plot3d(cube)

ユークリッド座標と同次座標

バーテックスで使用されている（x,y,z,S）は、homogeneous coordinate（同次座標、斉次座標、ホモグラフィ）と言うようです。
両者は左のリストのように変換されます。
Sは、x、y、z を割る効果があります。Sが2なら x、y、z は半分になります。
なぜ3Dの座標が同次座標で扱われるのかは良く分かりません。
（1,1,1）は同次座標では（1,1,1,1）とも（2,2,2,2）、（0.5,0.5,0.5,0.5）とも表せます。

行列演算で座標変換を行うことをアフィン変換と言うようです。affineはラテン語由来で neighboring, adjacent, next, bordering, related,marriage, akin, connected などが語義のようです。主に同次座標を対象にすることを指すものと思います。

前述の scaleMatrix() 、rotationMatrix(pi/4, 0, 0, 1)、translationMatrix() を組み合わせた cube の変換の例を考えます。

library(heplots)
arrow3d(c(4,4,8),cube$vb[1:3,8],barblen=.2,
lwd=3, col="red")
text3d(4,4,8.1,paste("(1,1,1) => ",
paste(round(cube$vb[1,8],3),
round(cube$vb[2,8],3),
round(cube$vb[3,8],3),sep=",")),
adj=c(1,1.5),color="#007700",cex=0.8)

cube3d() が生成する基本の立方体の頂点（1,1,1）がどこに移動したかを書き加えると図のように、（5.536,3.707,8）でした。

%*% は、R言語の行列の積の演算子です。* は、要素ごとの積になります。

> N <- matrix(c(1.414,-1.414,0,2.121,2.121,0,0,0,4),nrow=3,byrow=T)
> A <- scaleMatrix(2, 3, 4)[1:3,1:3]
> B <- rotationMatrix(pi/4, 0, 0, 1)[1:3,1:3]
> C <- translationMatrix(2, 3, 4)[4,1:3]
> A
[,1] [,2] [,3]
[1,] 2 0 0
[2,] 0 3 0
[3,] 0 0 4
> B
[,1] [,2] [,3]
[1,] 0.7071 -0.7071 0
[2,] 0.7071 0.7071 0
[3,] 0.0000 0.0000 1
> C
[1] 2 3 4
> N <- A %*% B
> t(N) %*% matrix(c(1,1,1),nrow=3) + C
[,1]
[1,] 5.536
[2,] 3.707
[3,] 8.000

mikeo_410