Perlinノイズっぽいものをpythonで、ついでにGPU化

# 概要

https://postd.cc/understanding-perlin-noise/

こちらを読んで理解した範囲でまとめ、pythonで実装してGPUで走らせます。
でも途中で僕の都合で違う式を使ったりもしています。

Perlinノイズというより、numpy配列の練習とcupyによるGPU化の実験をしたかったのかもしれません。

# 理屈

図は一応用意しましたが、参考URLにあります図を踏まえて書いていきます。

Perlinノイズのポイントは2点。

## １.グリッドによるノイズ関数

格子点に勾配ベクトルを仮定し、これとノイズを計算したい点の位置ベクトルとの内積を求めることで、実数座標での関数を実現しています。

ここでは座標を与えるとその座標に応じたノイズを返す関数を$f(\vec{x}) = v$と書きます。 

### 格子

$N_d$次元空間において実数での座標ベクトル$\vec{x}$が与えられたとき、それを含む格子の格子点$i\ (i=1, 2, \cdots, 2^{N_d})$について、$\vec{x}$と位置$\vec{p_i}$のベクトル

```math
\vec{d_i} = \vec{x} - \vec{p_i}
```
を得ます。ここで格子点$i$の座標$p_i$は床関数を使って$\lfloor x\rfloor + \vec{e_i}$ （$\vec{e_i}$は各次元が0または1のベクトル） と求められます。
以降、$\vec{x} \leftarrow \vec{x} - \lfloor \vec{x} \rfloor$として１つの格子内に注目し、式を簡単にします。

### 勾配

![2dgrid.png](https://qiita-image-store.s3.amazonaws.com/0/258019/fbf29bde-5ce8-2e4b-922a-19e305cd601f.png)


それぞれの格子点にランダムな勾配$\vec{g_i}$が割り当てられます。これらをもとに、各格子点で内積

```math
\vec{d_i}\cdot\vec{g_i}
```

を計算します。$\vec{g_i}$はランダムですが、格子点の座標が決まると一意に決まるべき値ですので、「ハッシュ」と呼ばれる値$h_i$に基づいて決定されます。

あとは格子中の$\vec{x}$の値に基づいて各格子点上の内積$\vec{d_i}\cdot\vec{g_i}$の値を補間して$f(\vec{x})$を得ます。

### 補間

補間ですが、参考URLでは線形補間`lerp`を各次元で順番に行ってますが、今回はBilinear補間を拡張したものを使ってみます。
すなわち、格子点$i$の重みは、その$i$の対角にある格子点$i_c$と$\vec{x}$で作られる直方体の体積となります。

![3dgrid.png](https://qiita-image-store.s3.amazonaws.com/0/258019/ad89490b-b6bb-220f-6e56-d394871486cc.png)


この定義だと、格子点$i_c$の座標は$p_{i_c} $で、$i_c$と$\vec{x}$で張る直方体の体積は対角線のベクトルの全成分の絶対値の積となります。numpyでは`prod()`で実装されているので、$\mathrm{prod}(\cdot)$と書きます。

```math
w_i = \left|\mathrm{prod}(\vec{e_{i_c}} - \vec{x})\right|
```

この時、補間するのに用いる座標をそのまま使うと格子をまたがるときに導関数が不連続になります。$y=|x|$のグラフを考えてみてください。
なので、これにfade関数$f_{fade}(x) = 6x^5 - 15x^4 + 10x^3$を$\vec{x}$の各座標で適用してできる新たな座標$\vec{u} = f_{fade}(\vec{x})$を使って補間します。
fade関数はなんでもいいのですが、

* $f_{fade}(0) = 0$, $f_{fade}(1) = 1$,
* 0, 1で導関数を0に：$f_{fade}^{(1)}(0) = f_{fade}^{(1)}(1) = 0$

という性質があることが重要です。

### まとめ

これらをまとめると、

```math
f(\vec{x}) = \sum_iw_i(\vec{d_i}\cdot\vec{g_i}) =\\
 \sum_i\left|\mathrm{prod}(\vec{e_{ic}} - f_{fade}(\vec{x} - \lfloor \vec{x} \rfloor))\right|\left[\left(\vec{x} - \lfloor \vec{x} \rfloor + \vec{e_i}\right)\cdot\vec{g_i}\right]
```

## 2. ノイズ関数の合成

大局的な構造と小さなノイズを両立させます。

```math
F(\vec{x}) = \sum_i a_if(\omega_i\vec{x})
```

たとえば、

```math
F(\vec{x}) = 64*f(1*\vec{x}) + 16*f(4*\vec{x}) + 16*f(4*\vec{x}) + 1*f(64\vec{x})
```

という感じです。$\vec{x}$に大きい数字がかかっている項は座標軸方向に圧縮する効果がありますので、局所的なギザギザを反映します。

### 注意

ホワイトノイズをぼかしたものだとバリューノイズというそうですが、それによる$f(x)$を使ってこの合成を以てperlinノイズとしている記事がちらほらありました。

# 実装

まずグリッドを用いたノイズ関数 $f$ を構築します。
以降、以下のようにnumpyをimportしています。

```python
import numpy as np
```

## グリッドの準備

ここでは$\vec{x}$によらない$f$の性質を決定します。クラスで実装するならコンストラクタに書きます。

### 格子点

実数で座標 $\vec{x}$が与えられたとき、床関数を適用すれば、$\vec{x}$を含む格子の原点側の格子点の座標が得られます。
これに各次元について１を足したり足さなかったりすれば別の格子点の座標も得られますが、これは各次元について足す・足さないの2通りで、$N_d$次元においては計$2^{N_d}$個の格子点となります。

この足す・足さないのベクトルをnumpyでまとめて得るのが以下の式です。

```python
unit = np.array(np.meshgrid(*([[0, 1]]*ndim))).T.reshape(-1, ndim)
```
`meshgrid`の中身は`meshgrid([0, 1], [0, 1], ・・・ )`と展開されます。

これの結果は、3次元(`ndim = 3`)だと

```python
array([[0, 0, 0],
       [0, 1, 0],
       [1, 0, 0],
       [1, 1, 0],
       [0, 0, 1],
       [0, 1, 1],
       [1, 0, 1],
       [1, 1, 1]])
```

となります。さきの理屈では$\vec{e_i}$に当たるベクトル群です。

### 対面のインデックス

`unit[i]`の対面の番号は、`unit`から`1-unit[i]`に一致する成分を探すことに相当します。

まず`i`と`i`の対面の`unit`ベクトルは足すと全成分が１になります。これを探します。

```python
np.all((unit + unit[2]) == 1, axis=1)
```

```python
array([False, False, False, False, False,  True, False, False])
```

これは各頂点について`unit[2]`の対面としての条件を満たすかどうかを判定していったものになります。今欲しいのはこのベクトルから何番目が`True`かという番号ですが、これは`np.where`を使って取得できます。

```python
np.where(np.all((unit + unit[2]) == 1, axis=1))
```

```python
(array([5]),)
```

これを全頂点についてスキャンします。ついでにnumpy配列にして、使いやすいように形を整えます。

```python
opposite = np.array([np.where(np.all((unit + ui) == 1, axis=1)) for ui in unit]).reshape(-1)
```

```python:oppositeの中身
array([7, 6, 5, 4, 3, 2, 1, 0])
```

### ハッシュ配列

0から定数まで、ここでは255までとして256までの数字をランダムに並べ替えた配列を用意します。
この使い方は、

* 格子点の座標を0 ~ 255の整数に正規化
* `h = hash[x]`
* `h = hash[h + y]`
* `h = hash[h + z]`

というようにして、全次元の座標を乱数決定に関わらせます。ここで0～255の値を取る`h`や`x, y, z`が引数に足し算で入ってきますので、長さを511までにしておく必要があります。
ここでは0～255までの配列を2つ合わせてシャッフルします。

```python
hash = np.arange(256 * 2, dtype=int) % 256
np.random.shuffle(hash)
```

### 勾配ベクトル

格子点にハッシュ値が割り当てられたら、それをもとに勾配を決定します。
式を見るとどんなのでもいいのですが、普通は各成分が`-1, 0, 1`のどれかになるようにするようです。
すると内積の計算が$\vec{d_i}$の成分の足し算引き算で大丈夫になります。

```python
sign = np.array(np.meshgrid(*([[-1, 0, 1]]*ndim))).T.reshape(-1, ndim)
```
`ndim = 3`だと

```python
array([[-1, -1, -1],
       [-1,  0, -1],
       [-1,  1, -1],
       [ 0, -1, -1],
       [ 0,  0, -1],
       [ 0,  1, -1],
       [ 1, -1, -1],
以下略
```
となります。参考URLではさらに、このうち2成分が`1`または`-1`で残りひとつが`0`のベクトルのみとしています。

```python
sign = sign[np.sum(np.abs(sign), axis=1) == 2]
```

```python:sign
array([[-1,  0, -1],
       [ 0, -1, -1],
       [ 0,  1, -1],
       [ 1,  0, -1],
       [-1, -1,  0],
       [-1,  1,  0],
       [ 1, -1,  0],
       [ 1,  1,  0],
       [-1,  0,  1],
       [ 0, -1,  1],
       [ 0,  1,  1],
       [ 1,  0,  1]])
```

numpyが1.12以降だと

```python
sign = sign[np.count_nonzero(sign, axis=1) == 2]
```

と、`count_nonzero`が`axis`引数を受け付けられるようになります。

この2成分というのが「次元数-1」なのか、「次元数に関わらず２」なのかが分かりませんでしたが、後者だと考えておきます。

あと、ハッシュ値からの勾配ベクトルの決定ですが、参考URLでは下位4ビットを元に12通りから決定しす。
今回は`sign`配列の長さでハッシュ値を割った余りをもとに決定することにします。

## ノイズ関数の実装

以降、$\vec{x}$から計算される部分を実装していきます。これはクラスで実装する場合は`__call__`などのメソッドに書くと良いと思います。

### パラメータ

* 次元数

`ndim`

* 格子数

`ntick`
各次元に格子は`ntick`個並んでいます。与えられた座標が0以下もしくは`ntick`以上になったら剰余を取ってこの範囲に収めます。

### 格子の取得

* $\lfloor \vec{x} \rfloor$

```python
ix = np.array(np.floor(x), dtype=int)
```

入力座標である`x.shape`は`(ndim)`、つまりndim次元のベクトルとして用意します。
`ix`は`x`と同じshape、`(ndim)`となります。

*  $p_i = \lfloor \vec{x} \rfloor + \vec{e_i}$ 

```python
grid_coord = ix + unit
```

`unit.shape`は`(2**ndim, ndim)`ですので、足し算で`ix`がブロードキャスト（次元追加とtiling）されます。
`grid_coord.shape`は`(2**ndim, ndim)`となります。


* $\vec{d_i} = \vec{x} - \vec{p_i}$

```python
dx = x - grid_coord
```

同様に`x`がブロードキャストされ、`dx.shape`は`(2**ndim, ndim)`です。

### 勾配とその内積

* ハッシュ値

まずは$\vec{x}$を取り囲んでいる格子点座標（整数値）を得ます。

```python
surround = grid_coord % ntick
```

`surround.shape`は`(2**ndim, ndim)`です。
この整数座標をもとに周りの格子点のハッシュ値を求めます。

```python
surround_hash = 0
for idim in range(ndim):
    surround_hash = hash[surround_hash + surround[:, idim]]
```
`surround[:, idim]`は`(2**ndim)`の配列です。
このshapeでそのまま1次元配列の`hash`のindexingに利用されるため、`surround_hash.shape`は`(2**ndim)`と、1次元配列です。

* $\vec{g_i}$

```python
grad_vec = sign[surround_hash % self.sign.shape[0]]
```

`sign`はベクトルが複数並んでいます。`sign`のindexに整数が渡されるとベクトルが得られるのですが、渡されたのが配列のため、そのまま渡された配列の数値に応じたベクトルを並べた配列が得られます。詳しくは「numpy indexing」でお調べ下さい。
とりあえず、`grad_vec.shape`は`(2**ndim, ndim)`です。`grad_vec[i]`は格子点`i`に割り当てられた勾配ベクトル$\vec{g_i}$です。

* $\vec{d_i}\cdot\vec{g_i}$

上で求めたベクトルのドット積を計算すればOKです。

```python
grad = np.sum(grad_vec * dx, axis=1)
```

`grad.shape`は先頭の次元のみ残すため、`(2**ndim)`です。
参考URLでは、ここで`grad_vec`の成分が`1, -1, 0`のみであるということを使って、`dx`の成分の足し算引き算のみで書いています。

### 補間

* $\vec{u}$と$w_i$

各頂点からの相対位置ベクトル`dx`を使います。これに`fade`関数を適用するにあたり、体積の計算をするために絶対値をとっておきます。

```python
u = fade(np.abs(dx))
```
これから、重みとして使うのは反対側の格子点での直方体ですので、添え字も並べ替えておきます。

```python
u = u[opposite, :]
#あるいは一行にまとめて　u = fade(np.abs(dx))[opposite, :]
```

ここで`u.shape`は`(2**ndim, ndim)`です。
このまま体積を計算します。

```python
weight = np.prod(u, axis=1)
```

`weight.shape`は`(2**ndim)`になります。

### まとめ

* $f(\vec{x}) = \sum_iw_i(\vec{d_i}\cdot\vec{g_i})$

いよいよ結論です。

```python
val = np.sum(grad*weight)
val = (val + 1.0)/2.0
```
スカラー値の`val`が得られました。これが$f(\vec{x})$の値になります。
最後の１足して２で割る行はあってもなくてもいいです。

お疲れ様でした。

# 複数点を一度に

ところが、この実装では使い物になりません。
どういうことかというと、これまでの実装では、`N`個の点があった時に各点での値がすべてほしいとなった時、

```python
xs = np.array(["N個の３次元ベクトル"])
v = [f(xi) for xi in xs]
```

とするしか無いですが、pythonのforは非常に重いです。実際、上の実装では10kオーダーの点数になると現実的でないです。
実装を`xs`のまま処理できるようにしていきます。

## ポイント

ニューラルネットワークで、numpy配列の最初の次元がミニバッチのサンプル番号に相当する配列となっているかと思いますが、同様の発想で、単純に先頭に１次元足してやります。

### 入力ベクトル

とりあえず、`xs.shape`は`(npoint, 1, ndim)`と真ん中に１を入れてやります。これは、

```python
ix = np.array(np.floor(xs), dtype=int)
grid_coord = ix + unit
```

という部分について考えると、`unit.shape`は`(2**ndim, ndim)`なのですが、`grid_coord.shape`は`(npoint, 2**ndim, ndim)`とならなければなりません。
このようにブロードキャストさせるためには、`ix`と`xs`のshapeが`(npoint, 1, ndim)`となっていれば実現できます。

http://www.kamishima.net/mlmpyja/nbayes2/broadcasting.html

なお、真ん中の`1`がなければ`npoint`と`2**ndim`の要素数が一致しないというエラーが出て止まるかと思います。

### それ以降

numpyを最大限活用したおかげで、numpy配列のindexを明示的に書いているところ、あとはnumpy関数の引数のaxisを適切に変更してやれば終わりです。

* 例

```python
surround_hash = 0
for idim in range(ndim):
    surround_hash = hash[surround_hash + surround[:, :, idim]]
 #  さっきは　（中略） + surround[:, idim]]

grad_vec = sign[surround_hash % self.sign.shape[0]]
# そのまま

grad = np.sum(grad_vec * dx, axis=2)

u = fade(np.abs(dx))[:, opposite, :]
# 対面取得のための添え字入れ替えを行う次元の前に１次元増えた
weight = np.prod(u, axis=2)
# axisが1から2に変わった

val = np.sum(grad*weight, axis=1)
# axis引数が増えた
```
実際これでCPU１コアでも200x200くらいの点なら10fpsはいけます。
Pythonのくせに。


# GPU化

Chainerの付属のcupyを使います。
numpyで多くの要素を一度に処理できるように書けたら、cupyに適宜置き換えてやるだけでかなりの性能アップが見込めます。

過去にもライフゲームをGPU化したりもしました。
https://qiita.com/sage-git/items/b4c01dc1c7ebec59c0d2
この時はscipyの関数のGPU化にChainerのレイヤーを使う必要がありましたが、今回はcupyの関数のみで事足りました。

## 関数実装部分

```python
import cupy as xp
```

あとは`np.hogehoge`を`xp.hogehoge`に変えてやるだけです。

ただ、対角線上の格子点インデックスを求めるための`where`などで、cupyはnumpyほど変数の型を柔軟に処理してくれなかったので、numpyで作成したものをcupy配列に変換しました。まあ関数設定時に準備として行う処理なので問題ないかと。

## 入出力部分

データの入力や後々の処理のためにはCPUとGPUでメモリのやり取りをしなくてはなりません。今回は以下のようにしてやりました。

```
xs = xp.array(x.reshape(npoint, 1, -1), dtype=xp.float32)
```
ここで`x`はnumpy配列、つまりCPU上のデータです。これをcupy配列のコンストラクタに渡してやればGPU上に渡されている、と思います（未検証）。

一方でデータをCPUに戻すには以下のようにしました。

```python
from chainer import cuda
np.array(cuda.to_cpu(val))
```
`val`はcupy配列です。これをそのまま`numpy.array()`に渡しても怒られました。

# 異なる周波数の合成

難しいところは何もありません。
何らかのノイズ関数を`grid`とした時、

```python
class PerlinNoise(object):
    def __init__(self, amp_param=[64, 16, 4, 1], freq_param=[1, 4, 16, 64]):
        self.ap = np.array(amp_param)
        self.fp = np.array(freq_param)
    
    def __call__(self, grid, xs):
        v = 0
        for a, f in zip(self.ap, self.fp):
            v += grid(xs*f)*a
        return v
```
とすればいいです。

# 実験

以上で実装は終わりです。
その他細かい設計を詰めたり、入力データを作ったり、結果を可視化したりはありますが、ここでは省略します。
最後にまとめて実際に動かしたコードをまとめて置いていますので、ご覧ください。

とりあえず、cupyにしたことでどれだけ早くなったかを確かめたいと思います。

200x200の40kピクセルの計算を1000フレームさせた時の経過時間から、一秒間でどれくらいの処理が可能かを見積もります。
perlinノイズ以外の処理もありますので厳密ではないですが、それなりに比較できるかと思います。

* ハードウェア
  * i7-6700 CPU @ 3.40GHz
  * GTX1060 (mem 6GB)

* 設定
  * 空間2次元、時間1次元の3次元グリッド
  * 2次元の霧のようなノイズを毎時間ステップ生成
  * 0.02刻みで0～4までのグリッドを使う
  * 画像の生成まで、表示はさせない

|実装|経過時間/秒|FPS|pixel per second|
|:-----|----:|-----:|-----:|
|Numpy|55.284186|18.088|723534|
|Cupy|3.159875|316.468|12658728|

圧倒的ですね。多少CPUが並列化できてたとしても追いつけません。

なお、この測定は周波数を変えて重ねあわせとかをやっていないので、もしその処理がある場合は単純に重ねあわせた数だけ遅くなるかと思います。

理論を追うことを重視したので、たぶん最適化の余地だらけです。
ただ、それなりの性能と一般性（次元数やグリッド数などが可変）を両立できたのではと思います。

# 他の条件

https://github.com/nsf/pnoise

こちらの条件に合わせてみます。

|実装|経過時間/秒|pixel per second|
|:-----|----:|-----:|
|Numpy|7.157519|915624|
|Cupy|1.94363|3371823|

Cとかの100倍遅い・・・。
まあNumpy実装はそんなものかというところですが、Cupy実装が全然早くないのが意外です。

と思って先行ベンチマークと同様にしていた標準出力への出力をサボると、

|実装|経過時間/秒|pixel per second|
|:-----|----:|-----:|
|Numpy|3.441069|1904524|
|Cupy|0.310274|21121907|

それなりですね。
ちなみにこの数字にはimport等の準備を含んでおりません。およそ1秒程度です。

# 見た目

せっかくなので、3Dノイズで空間2D＋時間としたときの、もやもやの様子を動画にしたものを用意しました。

![perlin.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/0fd3a0cf-b73c-ce44-04bc-ac283ed7cc4a.gif)

ついでに、周波数と振幅を変えて４つほど重ねたものを作ってみました。

![perlin4.gif](https://qiita-image-store.s3.amazonaws.com/0/258019/9438caca-3087-efb6-9d60-3b173aa7a19f.gif)


# その他

$\nabla f(\vec{x})$を実装したい。
あとはFortran実装。


# ソースコード群

## CPU・GPU兼用ノイズ関数

`singlepoint`はNumpyのみです。
Numpyで実行したいときは`xp`に`np`を代入してしまえばコードをCPU・GPUで流用できます。

```python:grid.py
import numpy as np
try:
    import cupy as xp
    from chainer import cuda
except ImportError:
    xp = np
    cuda = None

def default_fade(x):
    return x*x*x*(x * (x * 6 - 15) + 10)

class Grid(object):
    def __init__(self, ndim, tick_per_dim=256, random_seed=None, fade_function = default_fade):
        np.random.seed(seed=random_seed)
        self.ndim = ndim
        self.ntick = tick_per_dim
        self.hash = xp.arange(tick_per_dim * 2, dtype=int) % tick_per_dim
        np.random.shuffle(self.hash)
        u = np.array(np.meshgrid(*([[0, 1]]*ndim))).T.reshape(-1, ndim)
        self.unit = xp.array(u)
        t = [np.where(np.all((u + ui) == 1, axis=1)) for ui in u]
        self.opposite = xp.array(t).reshape(-1)
        self.sign = xp.array(np.meshgrid(*([[-1, 0, 1]]*ndim))).T.reshape(-1, ndim)
        # self.sign = self.sign[xp.sum(xp.abs(self.sign), axis=1) == 2]
        self.sign = xp.array(self.sign[xp.count_nonzero(self.sign, axis=1) == 2])
        self.fade = fade_function

    def __call__(self, x):
        if len(x.shape) == 1:
            return self.singlepoint(x)
        return self.multipoint(x)

    def singlepoint(self, x):
        ix = np.array(np.floor(x), dtype=int)
        grid_coord = ix + self.unit
        dx = x - grid_coord
        surround = grid_coord % self.ntick
        surround_hash = 0
        for idim in range(self.ndim):
            surround_hash = self.hash[surround_hash + surround[:, idim]]
        grad_vec = self.sign[surround_hash % self.sign.shape[0]]
        grad = np.sum(grad_vec * dx, axis=1)
        u = self.fade(np.abs(dx))[self.opposite, :]
        weight = np.prod(u, axis=1)
        val = np.average(grad, weights=weight)
        return val

    def multipoint(self, x):
        npoint = x.shape[0]
        xs = xp.array(x.reshape(npoint, 1, -1), dtype=xp.float32)
        ix = xp.array(xp.floor(xs), dtype=int)
        grid_coord = ix + self.unit
        dx = xs - grid_coord
        surround = grid_coord % self.ntick
        surround_hash = 0
        for idim in range(self.ndim):
            surround_hash = self.hash[surround_hash + surround[:, :, idim]]
        grad_vec = self.sign[surround_hash % self.sign.shape[0]]
        grad = xp.sum(grad_vec * dx, axis=2)
        weight = xp.prod(xp.abs(dx)[:, self.opposite, :], axis=2)**2
        weight = weight/(xp.sum(weight, axis=1).reshape(-1, 1))
        val = xp.sum(grad*weight, axis=1)
        if cuda is not None:
            return np.array(cuda.to_cpu(val))
        return val

def main():
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import axis3d
    r = np.arange(0, 5, 0.01)
    xs_grid = np.meshgrid(r, r)
    xs = np.array(xs_grid).T.reshape(-1, 2)
    g = Grid(xs.shape[1])
    v = g(xs)

    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")
    surf = ax.plot_surface(xs_grid[0], xs_grid[1], v.reshape(xs_grid[0].shape), cmap="bwr", linewidth=0)
    fig.colorbar(surf)
    plt.show()

if __name__ == '__main__':
    main()
```

## ノイズ合成とアニメーションテスト実行

```python:perlin.py
import numpy as np

class PerlinNoise(object):
    def __init__(self, amp_param=[64, 16, 4, 1], freq_param=[1, 4, 16, 64]):
        self.ap = np.array(amp_param)
        self.fp = np.array(freq_param)

    def __call__(self, grid, xs):
        v = 0
        for a, f in zip(self.ap, self.fp):
            v += grid(xs*f)*a
        return v

def main():
    import cv2
    from grid import Grid
    dx = 0.02
    L = dx*300
    dt = np.array([0.0, 0.0, 0.1])
    r = np.arange(0, L, dx)
    xs_grid = np.meshgrid(r, r, [0.0, ])
    xs = np.array(xs_grid).T.reshape(-1, 3)

    g = Grid(xs.shape[1], xs.shape[0], random_seed=1)
    pn = PerlinNoise(freq_param=[1.], amp_param=[64.])
    # pn = PerlinNoise()

    H, W, _ = xs_grid[0].shape

    while True:
        v = pn(g, xs)
        mi = v.min()
        ma = v.max()
        img = np.array((v - mi)/(ma - mi)*255, dtype=np.uint8).reshape([H, W])
        cv2.imshow("test".format(xs[0, 2]), img)
        cv2.waitKey(1)
        xs += dt

if __name__ == '__main__':
    main()
```

## 使い方

```bash:使い方
$ ls
# grid.py perlin.py

$ python grid.py
# 地形が3Dグラフがmatplotlibで描写されます。

$ python perlin.py
# OpenCVの画像表示機能で、霧のようなノイズが動き続けます。
# 終わるにはコンソール上でCtrl-Cしてください。

```