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

概要

こちらを読んで理解した範囲でまとめ、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}$のベクトル

\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$として１つの格子内に注目し、式を簡単にします。

勾配

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

\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}$で作られる直方体の体積となります。

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

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$

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

まとめ

これらをまとめると、

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. ノイズ関数の合成

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

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

たとえば、

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しています。

import numpy as np

グリッドの準備

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

格子点

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

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

unit = np.array(np.meshgrid(*([[0, 1]]*ndim))).T.reshape(-1, ndim)

meshgridの中身はmeshgrid([0, 1], [0, 1], ・・・ )と展開されます。

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

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ベクトルは足すと全成分が１になります。これを探します。

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

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

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

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

(array([5]),)

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

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

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つ合わせてシャッフルします。

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

勾配ベクトル

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

sign = np.array(np.meshgrid(*([[-1, 0, 1]]*ndim))).T.reshape(-1, ndim)

ndim = 3だと

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のベクトルのみとしています。

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

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以降だと

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$

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}$

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}$

dx = x - grid_coord

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

勾配とその内積

• ハッシュ値

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

surround = grid_coord % ntick

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

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}$

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です。

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関数を適用するにあたり、体積の計算をするために絶対値をとっておきます。

u = fade(np.abs(dx))

これから、重みとして使うのは反対側の格子点での直方体ですので、添え字も並べ替えておきます。

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

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

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

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

まとめ

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

いよいよ結論です。

val = np.sum(grad*weight)
val = (val + 1.0)/2.0

スカラー値のvalが得られました。これが$f(\vec{x})$の値になります。
最後の１足して２で割る行はあってもなくてもいいです。

お疲れ様でした。

複数点を一度に

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

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

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

ポイント

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

入力ベクトル

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

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)となっていれば実現できます。

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

それ以降

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

• 例

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の関数のみで事足りました。

関数実装部分

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に戻すには以下のようにしました。

from chainer import cuda
np.array(cuda.to_cpu(val))

valはcupy配列です。これをそのままnumpy.array()に渡しても怒られました。

異なる周波数の合成

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

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が並列化できてたとしても追いつけません。

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

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

他の条件

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

実装	経過時間/秒	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＋時間としたときの、もやもやの様子を動画にしたものを用意しました。

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

その他

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

ソースコード群

CPU・GPU兼用ノイズ関数

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

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()

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

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()

使い方

使い方

$ ls
# grid.py perlin.py

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

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