Rで３次元の色付き等高線図を描く

#要約
RStudioで以下のような3次元の色付き等高線図の描き方を紹介します。
![20201226-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179111/31ee72bb-29bd-20fd-0aac-cf04d5bd6b34.png)

- 等高線は、数式で与えられているものとします。上の図の場合は、点$(-1,-1,-1)$を中心とする球の式（$(x+1)^2+(y+1)^2+(z+1)^2$）で描画しています。

なお、一般的に「3次元の等高線図」というのは、下図のようなものを呼びますが、本記事では、上に示した図、強いて言い直せば「3次元空間内の等高**面**図」のことを指すこととします。
![20201226-5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179111/bdb15c80-b757-043d-9616-713295eac83e.png)
※上図の[参考Webサイト](http://cse.naro.affrc.go.jp/takezawa/r-tips/r/50.html "R-Source")

#環境
- OS: Windows10
- R: R x64 4.0.3
- RStudio: Version 1.3.1093

#説明
まずは、ソースコード（Rファイル）とその結果をお見せします。後ほど、説明をいたします。

##ソースコード

```test1.r
formula <- function(arg1, arg2, arg3){
  return((arg1+1)^2+(arg2+1)^2+(arg3+1)^2)
}

library(scatterplot3d)

n <- 301
x <- seq(-1,1,length=n)
y <- seq(-1,1,length=n)
z <- seq(-1,1,length=n)

valmin <- Inf
valmax <- -Inf

for(i in 1:n){
  for(j in 1:n){
    for(k in 1:n){
      v <- formula(x[i],y[j],z[k])
      valmin <- min(v, valmin)
      valmax <- max(v, valmax)
    }
  }
}

lvs <- pretty(c(valmin,valmax),10)
nlvs <- length(lvs)

xyz <- matrix(0,3*n^2,3)
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- min(y)
    xyz[(i-1+n)*n+j,3] <- z[i]
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}

clr <- cm.colors(nlvs-1)
clrs <- array(1:(3*n^2))
for(i in 1:length(clrs)){
  val <- formula(xyz[i,1],xyz[i,2],xyz[i,3])
  for(j in 1:nlvs-1){
    if(val < lvs[j+1]){
      clrs[i] <- clr[j]
      break
    }
  }
}

scatterplot3d(xyz,color=clrs,pch=".",
              xlab = "x",
              ylab = "y",
              zlab = "z")
```

##結果
![20201226-1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179111/31ee72bb-29bd-20fd-0aac-cf04d5bd6b34.png)

##解説
3次元の等高線図を描く関数というものは、（私の知る限り）存在しません。
今回は、3次元の**散布図**を描画する`scatterplot3d()`という関数で、無理やり3次元等高線図を描いています。

`scatterplot3d()`は、n個の空間内の点の座標を示すn行3列のデータを渡すと、3次元の散布図を描画してくれる関数です。
[こちらのWebサイト](http://www.sthda.com/english/wiki/scatterplot3d-3d-graphics-r-software-and-data-visualization#basic-3d-scatter-plots "Scatterplot3d: 3D graphics - R software and data visualization - Easy Guides - Wiki - STHDA")の例を用いると、以下のソースコード

```example1.r
scatterplot3d(iris[,1:3])
```

で、以下の3次元散布図を描画できます。
![20201226-2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179111/a4661a93-76b3-adeb-be24-987849dabf4b.png)

ここで用いた`iris`とは、Rにはじめから組み込まれているデータセットで、中身は以下のような感じです。
![20201226-3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179111/80798f2d-f752-5948-672e-11e76ae0dbbd.png)
この例では、このデータセットの最初の3列を用いて、3次元散布図を描いています。

今回は、この`scatterplot3d()`を用いて、数式の値を元に色を付けた点を、ひたすら空間内に敷き詰めた散布図を作ることで、3次元等高線図を作成します。

ところで、点を敷き詰めるのは、グラフで見えている3つの面だけで十分です。見えてないところに敷き詰めてもいいのですが、描画処理時間と、グラフを外部出力するなら出力ファイル容量の無駄です。今回のソースコードでも、見えている3つの面にのみ、色付きの点を敷き詰めています。

一方、色を付ける際には、敷き詰める点についてだけ考えるのでは不十分です。
等高線図の色をつけるためには、グラフ描画範囲における値の範囲、すなわち、最大値と最小値が必要です。値の範囲がわかることで、各値について、どの色で塗ればいいのかがわかります。
ところが、点を敷き詰める3つの面に、最大値と最小値が両方ともあるとは限りません。なので、値の範囲を調べるために、値についてはグラフ描画範囲全体について調べる必要があります。

以上を踏まえたうえで、ソースコードを見てみましょう。

```test1-1.r
formula <- function(arg1, arg2, arg3){
  return((arg1+1)^2+(arg2+1)^2+(arg3+1)^2)
}

library(scatterplot3d)

n <- 301
x <- seq(-1,1,length=n)
y <- seq(-1,1,length=n)
z <- seq(-1,1,length=n)
```

`formura`は、等高線に用いる数式です。値の範囲を調べるときと、点に色をつけるときの2回、この数式を用いるので、はじめに関数化しておきます。

`scatterplot3d()`を使うためには、`library(scatterplot3d)`でライブラリを呼び出す必要があります。なお、`scatterplot3d`がインストールされていない場合は、別途、インストールする必要があります。インストールの際には、ユーザに管理者権限がないとエラーが出ますのでご注意ください。

`n <- 301`で、1行あたりに敷き詰める点の個数を設定します。ここでは`301`で設定しています。

`x <- seq(-1,1,length=n)`、`y <- seq(-1,1,length=n)`、`z <- seq(-1,1,length=n)`で、3軸それぞれの配列を作成します。

```test1-2.r
valmin <- Inf
valmax <- -Inf

for(i in 1:n){
  for(j in 1:n){
    for(k in 1:n){
      v <- formula(x[i],y[j],z[k])
      valmin <- min(v, valmin)
      valmax <- max(v, valmax)
    }
  }
}
```
さて、ここでは、値の範囲（最大値・最小値）を計算しています。グラフ描画範囲全体（n×n×n）を走査して、値の最大値と最小値を得ます。
はじめに作成した`formula()`に、`x[i]`、`y[j]`、`z[k]`を渡して、値を取得しています。
なお、最小値の初期値に`Inf`（正の無限大）を、最大値の初期値に`-Inf`（負の無限大）を、それぞれ設定しています。Rでは無限大を変数に代入することが可能です。

```test1-3.r
lvs <- pretty(c(valmin,valmax),10)
nlvs <- length(lvs)
```

`lvs`には、等高線の水準（塗り分ける色の境界値）を代入しています。`pretty()`は、与えられた配列に対して、「いい感じの」水準配列を返してくれる関数です。このときの`pretty(c(valmin,valmax),10)`は、

```example2.r
> pretty(c(valmin,valmax),10)
 [1]  0  1  2  3  4  5  6  7  8  9 10 11 12
```

のように、`0`から`12`までの13個の整数の配列を返しています。
なお、`pretty(c(valmin,valmax),10)`の第二引数`10`は、水準数の「目安」を指定しています。あくまで目安なので、この数ぴったりの水準数になるとは限らない、という点に注意です。（今回も水準数は13個になっている）
`nlvs`には、`lvs`の長さを代入しています。

```test1-4.r
xyz <- matrix(0,3*n^2,3)
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- min(y)
    xyz[(i-1+n)*n+j,3] <- z[i]
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}
```

`xyz`には、描画する点の座標をセットします。1面あたりの点の数はn×nで、それが3面分あるので、`xyz`は`3*n^2`行`3`列の行列としています。初期値は`0`です。

二重`for`文の中では、具体的に点の座標を代入しています。
最初の3行はxy平面の、次の3行はzx平面、最後の3行はyz平面の点をx座標、y座標、z座標を代入しています。
意味合い的には以下のような3つの二重`for`文に分けた方が分かりやすいですが、処理時間短縮のために、1つの二重`for`文で済ましています。

```example3.r
#xy平面
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
  }
}

#zx平面
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- min(y)
    xyz[(i-1+n)*n+j,3] <- z[i]
  }
}

#yz平面
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}
```

さて、次は点につける色の配列を作成します。

```test1-5.r
clr <- cm.colors(nlvs-1)
clrs <- array(1:(3*n^2))
for(i in 1:length(clrs)){
  val <- formula(xyz[i,1],xyz[i,2],xyz[i,3])
  for(j in 1:nlvs-1){
    if(val < lvs[j+1]){
      clrs[i] <- clr[j]
      break
    }
  }
}

scatterplot3d(xyz,color=clrs,pch=".",
              xlab = "x",
              ylab = "y",
              zlab = "z")
```

`clrs`は、`scatterplot3d()`の引数`color`に渡す配列です。この配列は、描画する点の数と同じ長さの配列である必要があり、対応する点の色をこの配列で指定することができます。

具体例を示します。

```example4.r
scatterplot3d(rbind(c(1,1,1),c(2,2,2),c(3,3,3)),color=c("#FF0000","#00FF00","#0000FF"),pch=19)
```

このようなコードを実行すると、以下のような3次元散布図が描画されます。
![20201226-4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179111/8b34547a-25b1-b763-6fb8-cdaf5109c4ff.png)
1つ目の点`c(1,1,1)`には`#FF0000`（赤）を、2つ目の点`c(2,2,2)`には`#00FF00`（緑）を、3つ目の点`c(3,3,3)`には`#0000FF`（青）を、それぞれ対応させています。

このようにして、各点の色を指定する配列を作成します。今回は、それが`clrs`です。

`clrs`を作成するにあたり、その色サンプルとして`clr`を用意します。
`clr`には、`cm.colors(nlvs-1)`をセットします。

`for`文の中では、`formula()`関数から得た値が、各水準のどこに当てはまるのかを`if`文で判定し、`True`となったところの色を`clrs`にセットしています。

最後の`scatterplot3d()`では、点の座標データ`xyz`と、各点の色指定配列`clrs`、そして、マーカーとしてドット`.`を指定しています。また、各軸ラベルもセットしています。

マーカーとしてドット`.`を指定しているのは、最も小さな点だからです。`pch=20`が、もともと組み込まれているマーカーの中では最も小さな点ですが、これでもまだ大きいので、ドット`.`を指定しています。

#おまけ（y軸反転）
「[Rの３次元散布図（scatterplot3d()）での軸の反転方法](https://qiita.com/kaz0414/items/39b1ac17642df4042bc0 "Rの３次元散布図（scatterplot3d()）での軸の反転方法 - Qiita")」に記載の方法を使うことで、y軸を反転させて、右手系座標で描画することもできます。
##ソースコード

```test2.r
formula <- function(arg1, arg2, arg3){
  return((arg1+1)^2+(arg2-1)^2+(arg3+1)^2) #中心を(-1,1,-1)に変更
}

library(scatterplot3d)

n <- 301
x <- seq(-1,1,length=n)
y <- seq(-1,1,length=n)
z <- seq(-1,1,length=n)

valmin <- Inf
valmax <- -Inf

for(i in 1:n){
  for(j in 1:n){
    for(k in 1:n){
      v <- formula(x[i],y[j],z[k])
      valmin <- min(v, valmin)
      valmax <- max(v, valmax)
    }
  }
}

lvs <- pretty(c(valmin,valmax),10)
nlvs <- length(lvs)

xyz <- matrix(0,3*n^2,3)
for(i in 1:n){
  for(j in 1:n){
    xyz[(i-1)*n+j,1] <- x[j]
    xyz[(i-1)*n+j,2] <- y[i]
    xyz[(i-1)*n+j,3] <- max(z)
    xyz[(i-1+n)*n+j,1] <- x[j]
    xyz[(i-1+n)*n+j,2] <- max(y) #ここをmin(y)からmax(y)に変更
    xyz[(i-1+n)*n+j,3] <- z[i]
    xyz[(i-1+2*n)*n+j,1] <- max(x)
    xyz[(i-1+2*n)*n+j,2] <- y[j]
    xyz[(i-1+2*n)*n+j,3] <- z[i]
  }
}

clr <- cm.colors(nlvs-1)
clrs <- array(1:(3*n^2))
for(i in 1:length(clrs)){
  val <- formula(xyz[i,1],xyz[i,2],xyz[i,3])
  for(j in 1:nlvs-1){
    if(val < lvs[j+1]){
      clrs[i] <- clr[j]
      break
    }
  }
}

xyz[,2] <- -xyz[,2] #y成分反転

scatterplot3d(xyz[,1],xyz[,2],xyz[,3],color=clrs,pch=".",
              xlab = "x",
              ylab = "y",
              zlab = "z",
              y.ticklabs=seq(max(xyz[,2]), min(xyz[,2]), -0.5)
              )
```

##結果
![20201226-6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179111/d4b4140e-0a6f-24a7-ebab-28a042417179.png)


#姉妹記事
[Rで色付き等高線図を描く](https://qiita.com/kaz0414/items/3da341090f104fcb44c6 "Rで色付き等高線図を描く - Qiita")