Rを用いた多次元立方体の描画

#概要
　以下のツイートがややバズったのでソースコードを含めて解説します。
<blockquote class="twitter-tweet" data-lang="ja"><p lang="ja" dir="ltr">VR界で4次元立方体がバズってるから私が作った5次元立方体ももうちょい流行ってもよくない？ <a href="https://t.co/unagNhuoqg">pic.twitter.com/unagNhuoqg</a></p>&mdash; Hibit (@hibit_at) <a href="https://twitter.com/hibit_at/status/1037664390214897666?ref_src=twsrc%5Etfw">2018年9月6日</a></blockquote>

<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>
　使用する言語はR。
#ソースコード
　とりあえずソース見せろという方へ。

```R:
library(rgl)

col_ind = c("red","green","blue")
snap = 1

for (count in 1:10){

t = 0

while (t < 4*pi){ 

x = numeric(4)

#rot行列を定義
n = 5
rot = array(1,dim=c(n,n,n*(n-1)/2))

for (k in 1:(n*(n-1)/2)){
rot[,,k]=diag(n)
}

k = 1

for (i in 1:(n-1)){
for (j in (i+1):n){
rot[i,i,k] = cos(t)
rot[i,j,k] = sin(t)
rot[j,i,k] = -sin(t)
rot[j,j,k] = cos(t)
k = k + 1
}
}

#ピッチ行列を決定
pitch = 0.01
delta = matrix(1, nrow=5, ncol=5)
for (i in 1:5){
delta[i,i] = pitch
}

#こっから描画
for (i in 1:5){
x1 = -1/2
x2 = -1/2
x3 = -1/2
x4 = -1/2
x5 = -1/2
while(x1 <= 1/2){
while(x2 <= 1/2){
while(x3 <= 1/2){
while(x4 <= 1/2){
while(x5 <= 1/2){
z = rot[,,count] %*% c(x1,x2,x3,x4,x5)
y = z[1:3]*2/(3/2-z[4])+z[5]*1
y[3] = z[3]*2/(3/2-z[4])+z[5]*0.7

#グラフの色付け
if (x5 == -1/2){
y = c(y,1)
}else{
if (i == 1){
y = c(y,3)
}else{
y = c(y,2)
}
}

x = rbind(x,y)
x5 = x5 + delta[i,1]
} #while-x5
x5 = -1/2
x4 = x4 + delta[i,2]
} #while-x4
x4 = -1/2
x3 = x3 + delta[i,3]
} #while-x3
x3 = -1/2
x2 = x2 + delta[i,4]
} #while-x2
x2 = -1/2
x1 = x1 + delta[i,5]
} #while-x1
} #for-i

x = x[-1,]
clear3d()
part = c(-2,2)
plot3d(x,type="s",col=col_ind[x[,4]],xlim=part,ylim=part,zlim=part,size=1.5)

rgl.snapshot(paste("C:/Users/",count,"軸",snap,".png",sep=""))
snap = snap + 1
t = t + pi/40

} #while-t

} #for-count
```

##注意点１
　このプログラムはrglライブラリを用いているため、事前にインストールしておく必要があります。
　インストールされていない方は、

```R:
install.packages("rgl")
```
　でインストールされたし。
##注意点2
```R:
rgl.snapshot(paste("C:/Users/",count,"軸",snap,".png",sep=""))
```
　を見ればわかる通り、このソースコードをコピペして実行した場合、貴方のユーザーディレクトリ直下に大量（16,000枚）のpngが投下されます。作業用のフォルダを指定してそこに生成するようにしましょう。

#解説

##rot行列
```R:
n = 5
rot = array(1,dim=c(n,n,n*(n-1)/2))

for (k in 1:(n*(n-1)/2)){
rot[,,k]=diag(n)
}

k = 1

for (i in 1:(n-1)){
for (j in (i+1):n){
rot[i,i,k] = cos(t)
rot[i,j,k] = sin(t)
rot[j,i,k] = -sin(t)
rot[j,j,k] = cos(t)
k = k + 1
}
}
```
　ここでは回転行列（正確には3次元行列）を生成しています。
　2次元だと

```math
\begin{bmatrix}
cosθ & -sinθ \\
sinθ & cosθ 
\end{bmatrix}
```

　で済むアレですが、アレを多次元に拡張する必要があります。と言っても話はほとんど変わらなくて、要は「2つの成分だけ上の行列を作用させて後はそのまま」なので、

```math
\begin{bmatrix}
cosθ & 0 & 0 & 0 & -sinθ \\
0 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 \\
sinθ & 0 & 0 & 0 & cosθ
\end{bmatrix}
```
　みたいな行列を用意すればいいだけです。上の行列は1次元目と5次元目の回転行列ですが、実際は$_nC_2$通りつまり$n=5$では10通り作る必要があるので、上のコードではそれを生成しています。<font color="Red">（2018/09/08　指摘を基に、上の行列のミスを直しました。あーあーあ～さんありがとうございます。）</font>

　後で、

```R:
z = rot[,,count] %*% c(x1,x2,x3,x4,x5)
```
　という形で使います。

##ピッチ行列
```R:
pitch = 0.01
delta = matrix(1, nrow=5, ncol=5)
for (i in 1:5){
delta[i,i] = pitch
}
```
　ここではピッチ行列を生成しています。このプログラムでは、大量の球体を描画することによって立方体を表現していますが、ある次元は0.01刻みで描画し、別の次元は1刻みで描画することにより辺を表現しています。スマートではないですね。スマートなやり方があったら教えてください。
　この刻みも行列で制御します。

```math
\begin{bmatrix}
0.01 & 1 & 1 & 1 & 1 \\
1 & 0.01 & 1 & 1 & 1 \\
1 & 1 & 0.01 & 1 & 1 \\
1 & 1 & 1 & 0.01 & 1 \\
1 & 1 & 1 & 1 & 0.01
\end{bmatrix}
```

　というような行列を作ります。<font color="Red">（2018/09/08　指摘を基に、上の行列のミスを直しました。あーあーあ～さんありがとうございます。）</font>
　後で、

```R:
x5 = x5 + delta[i,1]
} #while-x5
x5 = -1/2
x4 = x4 + delta[i,2]
} #while-x4
x4 = -1/2
x3 = x3 + delta[i,3]
} #while-x3
x3 = -1/2
x2 = x2 + delta[i,4]
} #while-x2
x2 = -1/2
x1 = x1 + delta[i,5]
} #while-x1
```

　という形で使います。
##投影
　5次元のままでは3次元で描画できないので3次元に投影する必要があります。投影の方法としては[ステレオ投影](https://ja.wikipedia.org/wiki/%E3%82%B9%E3%83%86%E3%83%AC%E3%82%AA%E6%8A%95%E5%BD%B1)が有名ですが、これだけだと1次元しか集約できないので、その他に[斜投影](https://ja.wikipedia.org/wiki/%E6%8A%95%E5%BD%B1%E5%9B%B3#%E6%96%9C%E6%8A%95%E5%BD%B1)を使います。
　ソースコードでは、

```R:
z = rot[,,count] %*% c(x1,x2,x3,x4,x5)
y = z[1:3]*2/(3/2-z[4])+z[5]*1
y[3] = z[3]*2/(3/2-z[4])+z[5]*0.7
```
　の部分になります。ステレオ投影の中心軸を（3/2）にしたり、斜投影の角度をずらしていたり（そうしないと他の立方体の辺と重複してわかりづらくなってしまうので）するので、無駄に複雑に見えますが、要は$(a,b,c,d,e)\in\mathbb{R^5}\mapsto(a',b',c') \in\mathbb{R^3}$を考えた時に、

```math
a'=\frac{a}{const(a')-d}+const'(a')e\\
b'=\frac{b}{const(b')-d}+const'(b')e\\
c'=\frac{c}{const(c')-d}+const'(c')e\\
```
　なる変換を行っているだけです。※$const,const'$は$a',b',c'$に応じた定数
　後は、各辺を色付けして、rglのplot3dコマンドを用いてグラフに描画するだけです。

#補足

　動画への書き出しは[AVIUtl](http://spring-fragrance.mints.ne.jp/aviutl/)を用いています（pngのファイル名を連番にしているのはそのためです）。Rでも出来ると思いますが、詳しくは調べていません。

　それでは皆様も良い多次元ライフを。

（2018/09/15追記）[ブログ](https://deux-hibi.hatenablog.com/entry/2018/09/14/203205)に多次元空間や投影についての補足記事を書きました。