MineCraftの地形。矩形だけども美しい(^-^)
その３ パーリンノイズとフラクタル
～ フラクタブルな地形を作る基盤 ～
① ノイズのお話
オープンワールドを採用しているゲームでは壮大な地形が用意されています。山あり川あり海ありと、地球上のどこかに本当に存在しているかのようなその仮想世界は本当に美しく、ゲームを忘れただひたすらに駆け回りたくなってしまいます。そういう地形の中にはもちろん開発者が丹念に作り上げているものもありますが、例えばマインクラフトのワールドなどは開発者が最初から地形を用意しているわけではなく、完全に自動生成されています。マイクラをやった事のある人なら、その地形がとてもプログラムが作り上げ得たとは思えない程自然っぽくてパターンぽい所が見受けられず変化に富んでいる事を知っているはずです：
MineCraftの地形。矩形だけども美しい(^-^)
上のスクリーンショットのような地形は、基本的には特定の座標にある地面を盛り上げたり凹ませて出来る（マイクラは洞窟や穴開きなどがあるので実際はもっと複雑な調整をしていると思います）訳ですが、単純にランダムに凸凹させても決して上のような絵にはなりません。「単純にランダムに」というのは特定の座標の地形の高さを乱数で決めてしまうという事です。試しに1次元なグラフでそれを可視化してみましょう：
これは横軸が座標、縦軸がその座標での地形の高さを表しています。確かに凸凹してはいるのですが、お世辞にも普段良く目にする自然の地形のようには見えませんよね。自然の地形は完全ランダムでは無いという事です。
自然を観察すると、まず大きな山や谷があります。その山肌をズームするとつるんとしているわけではなく、より小さな山と谷が折り重なって構成されています。その小さな山も似たように凸凹しています。このように大局的に見ても微視的に見てもその形状に似たような傾向がある図形を「フラクタル」と言います。そう、自然には至る所にフラクタルが溢れているんです。
大きな山と谷の中にさらに細かい山と谷…、こういう波のような形状を表現する一つの方法に「三角関数による波の合成」があります。様々な周波数の三角関数を重ね合わせると複雑な波を作り出す事が出来ます。下のグラフは大きな周波数と振幅（波の高さ）を持ったcos波と、それより周波数も振幅も小さいcos波をいくつか重ね合わせてみたものです：
最初の完全ランダムな凸凹よりは山と谷がしっかり識別できます。大きな谷の中にも小さな山と谷があって、フラクタルな図形になっているように見えますよね。ただ…これが自然かと言うと何か違う気がします。違和感の理由はあまりにも機械的な「繰り返し感」のためです。実際の自然はフラクタルではあるのですが、ここまで無機質な周期性がある訳でもありません。ありとあらゆる周波数の三角関数を重ねれば周期性を感じにくくはできますが、マイクラのような地形をそれで作るには何千何万という数のパラメータを調整する必要があり、多分実用的ではありません。
適度にランダムで自然な山や谷があり、且つ繰り返しを感じさせない地形を作るにはどうしたら良いのでしょうか？この問題に光を灯した一つのアイデアがパーリンノイズです。
② 2つのパーリンノイズ
パーリンノイズは1985年にKen PerlinがSIGGRAPHで発表した適度なランダム性がありつつもフラクタルな性質を持つノイズ生成法です（論文はこちら(PDF)）。このパーリンノイズの応用性の高さからPerlinは1987年にアカデミー科学技術賞を受賞しています。さらにPerlinは自身のパーリンノイズの欠点を改善した「改良パーリンノイズ」を2002年に発表しています（論文はこちら(PDF)）。ここでは以後「パーリンノイズ」と言うと1985年に発表された最初の物、「改良パーリンノイズ」は2002年の改良版を指すものとします。
③ パーリンノイズの作り方
実はパーリンノイズも先の三角関数のように周波数を持っています。三角関数と違うのは波の高さに繰り返し感がありません。実際のパーリンノイズ（1次のパーリンノイズ）をグラフ化するとこんな感じになっています：
このパーリンノイズの周波数は1です。グラフとX軸との交点を見ると1,2,3...と整数部分で必ず交差しているのが分かると思います。上下の波の高さに何かパターン的なものは感じられませんよね。この周期性がありつつも繰り返し感が無い性質が自然物を表現するのに非常に適しているんです。ではパーリンノイズの作り方を見ていきましょう。
まずX軸方向を整数間隔に区切ります。各区画を「格子」と呼ぶ事にします。上のグラフがそうであるように、その格子の境目ではノイズの大きさ（上のグラフの高さ）は必ず0になるというルールを敷きます。次に「ウェーブレット関数（Wavelet）」という特別な関数を格子の境目上に並べます。ウェーブレット関数の典型を下に示します：
ウェーブレット関数と言うのは、ある有限な区間にだけ作用する波の事です。パーリンノイズとは、格子上に存在するウェーブレットが合成されたノイズなんです。パーリンノイズで使用するウェーブレット関数は以下の条件を満たす必要があります：
上のグラフはこの4つの条件を満たしているのが分かると思います。ウェーブレット関数の形を決めるのは原点での傾きaだけです。このウェーブレット関数の原点を格子の境界上に来るようにし、各境界で適当な傾きaを決めると、境界別にウェーブレット関数が並ぶことになります：
青いのはx=0をその中心に、オレンジのはx=1をその中心にしたウェーブレット関数です。各中心点での傾きが違うので波の高さが異なっているのがわかりますよね。このように境界ごとにウェーブレット関数を並べたら、次に任意のxでのノイズ値を計算します。やり方は簡単で、青い方の値W0とオレンジの値W1とを0～1の間で線形補間してあげるだけです。式に書くと、
となります。t=0の時はW0(0)の値が、t=1の時にはW1(1)の値が採用され、0<t<1の間は線形に値がブレンドされます。実際この式で上の0～1の間を補間するとこんなグラフが出来あがります：
緑色のが補間したグラフであり、これが求めたいパーリンノイズその物です(^-^)。パーリンノイズを作るには「各境界に設けた適当な傾きを持つウェーブレット関数を使って矩形内の値を補間する」というだけなんです。本質はとても簡単ですね。
ではそのウェーブレット関数をどうやって作るのか？ここがポイントとなります。
④ ウェーブレット関数の作り方
ウェーブレット関数は先に挙げた条件さえ満たせば実は何でも良いので、ここではPerlinが用いたウェーブレット関数の作り方を紹介します。
まず用意するのは原点(x=0)で1、x=1及びx=-1で0になる次のような3次関数C(t)です：
絶対値がポイントで、これによって原点を中心に鏡写しのようなグラフになります。このグラフの主たる目的は、まず左右対称である事、そして「x=±1の所で緩やかにゼロになる」事にあります。これ、ウェブレット関数の条件の一つですよね。でもこのグラフだけだと原点で0を通るとか、原点での傾きがaであるなどその他の条件は満たしていません。そこで次に傾きがaで原点を通る直線L(t)をC(t)に掛け算します。すると…
黒いグラフが傾きaの直線です。青い3次曲線C(t)にこの直線の値を掛け算すると紫色のグラフで示したウェブレット関数W(t)が得られます。直線は原点を通りそこで値が反転するので、掛け算すると上手い事原点で反転が起きてウェブレット関数の残りの条件を全部作ってくれるわけです。素晴らしい(^-^)
3次関数C(t)はtに対して常に同じ値を返すのに対し、直線L(t)は傾きaによって値が変わります。傾きが大きい程掛けられる値も大きくなるので、傾きaが大きい程ウェブレット関数の波は高くなります：
左がa=0.7、右が2.0です。また、波の山谷は傾きの符号で入れ替わります：
格子の境界での直線の傾きaの大きさと符号により色々なウェブレット関数が作れます。これはつまり、「格子の境界毎に適当な傾きaを乱数で振れば格子上に様々なウェブレット関数が乗る」ということになります。各境界での様々なウェーブレット関数が矩形内で補間される事で、様々な凸凹がありかつ特定の周波数も維持しているノイズ値となります。これがパーリンノイズです！
⑤ 2次のパーリンノイズ
ここまでは矩形が1次元であるパーリンノイズ（1次のパーリンノイズ）で作り方の基本を説明してきました。同じような考え方をXYの2次元の世界に拡張する事が出来ます。パーリンノイズは多次元なほど面白いんです(^-^)
1次元で矩形は数直線上の線分でしたが、2次元だと正方形になります。そして矩形の境目（整数座標）は「角」という事になります。1次元の時と同じで矩形の角にはウェーブレット関数を乗せますが、次元数が増えているので乗せるのも2次元のウェーブレット関数です。その性質は1次元とほぼ変わりませんが、次元が増えたのでより一般的な表現になります：
1次の時と少し違う部分がありますので説明します。ウェーブレッド関数の原点での値がゼロであるというのは同じですが、そこでは傾きではなくて勾配(ax, ay)がある、としています。勾配（gradient）というのは傾きを一般化したようなものです。簡単な説明をしますね。中学校で習う原点を通る直線の式には傾きの項があります：
上式のaが傾きで、単純に特定のx座標を与えると対応する値が算出できました。XYZの3Dの世界になると、傾き具合というのは直線では無くて平面で表現しなければなりません。原点を通る平面の式はこんな感じになります：
axとayは平面のどこかに立った時のX軸方向の傾きとY軸方向の傾き（勾配）になります。先程の直線と同じように、座標(x,y)での値は上式にその座標を放り込めば計算できます。この式、見方を変えると座標(x, y)と勾配(ax, ay)の「内積」になっている事がポイントです。一般にもっと次元の高い世界での平面の式も各成分とその傾きを掛けて全部足す式になります。
次に「周囲の境界矩形の所でなだらかにゼロになる」という条件は次の図のような状態を言っています：
上の図は中心点に2次元ウェーブレット関数の原点を置いています。赤い所は値がプラスの所、青い所はマイナスの所で、黒くなる程値がゼロに近づいていきます（灰色の線は等高線）。この図の領域の端の方はどこもなだらかにゼロになっているのがわかると思います。これは1次のパーリンノイズのウェーブレット関数がx=-1と1の所でなだらかにゼロになっていたのと一緒です。2次のパーリンノイズではこんな感じのウェーブレット関数を使用します。「何か難しくなった～」と思うかもしれませんが、大したことはありません。1次のウェーブレット関数の考え方が大分使えます。
今矩形のある角に注目します。例えば上の図の中心点だと思って下さい。そこからの相対座標を(u, v)とします。uとvの範囲をそれぞれ-1～1の間とすると、上図の範囲を(u, v)で指し示す事ができます。ここでX軸Y軸それぞれの3次関数Cx(t)及びCy(t)を定義し、u成分及びv成分に該当する値Cx(u)とCy(v)を求めます：
求めたCx(u)とCy(v)を掛け算すると、中心点で値が1で周囲の矩形に行くほど値がゼロに向かってなだらかに減少する山になります：
同心円ではなくてちょっと矩形っぽい形になります
式で書くと、
このように3次関数同士を掛け算した形になります。
次に中心点を通って適当な勾配a=(ax, ay)を持つ平面を作り、その平面の値を上の山C(x, y)に掛け算します。これは原点で値をゼロにしたり山と谷を作って積分を0にしたり（＝原点を対称に勾配に沿った山谷をつくるため）という残りの条件を満たすようにするためです。1次のウェーブレット関数を作る時も同じような事をしましたよね。勾配がaである平面は、
で表せて、式の形としては勾配と座標との内積になっているのでした。傾きaxとayをランダムで与えてL(u, v)を定め上の山C(x, y)に掛け算すると2次のウェーブレット関数が出来上がります：
さて、矩形の4角に適当な勾配を設定したウェーブレット関数を作ったら、後は座標(u,
v)での値を線形補間で算出する事になります。補間の方法ですが、まず左下角(W00)と右下角(10)の対、および左上角(W01)と右上角(W11)の対でそれぞれウェーブレット関数の値を線形補間します。具体的には、
こんな感じで線形補間し座標u上の上下端にある2つの値を算出します：
最後にY軸方向も線形補間すると目出度く座標(u, v)での2次のパーリンノイズの値が算出されます。
実際に上のアルゴリズムで出力したパーリンノイズがこちらです：
白い所が値が大きく、黒くなる程値が小さい事を表しています。矩形単位の周波数を保った細かい凹凸があり、繰り返しを感じさせない良い感じなノイズになってます(^-^)
この多次元へ拡張した方法は、もっと次元が高くなっても基本一緒です。矩形からの相対座標(u0, u1, u2,...)から3次関数C(t)を使ってなだらかな山をまずは作り、それに原点で適当な勾配（a=(a0, a1, a2,...)）を持った平面（超平面）を掛け算するとN次元のウェーブレット関数が出来上がります。矩形内の相対座標(u0, u1, u2,...)に対応した値の補間の方法もまずはu0軸で補間、次にu1軸で補間、というのを繰り返していくだけです。
⑥ 改良パーリンノイズ
ここまでは1985年にPerlinが発表した元祖パーリンノイズについてのお話でした。2002年にPerlinは自らのパーリンノイズの欠点を補う「改良パーリンノイズ」を公表しています。
元祖パーリンノイズからの改良点は2つあります。一つはウェーブレット関数を作る時に使ってきた3次関数を次の5次関数に変更しています：
この関数と3次関数とを比較したグラフがこちら：
オレンジがこれまでの3次関数、青が改良した5次関数です。5次関数の方が傾きがゼロになる所がよりなだらかになっています。見た目にはちょっとしか違わないように見えますが、5次関数の方がt=0の所での「連続性」が高くなっています。上の3次関数は2階微分するとt<0の範囲では-6-12t、t>=0の範囲では-6-12tとなってt=0の所で不連続となります。これは曲線の傾きの変化量がt=0の所で不連続である事を表します。接線と密接にかかわっているのは法線です。つまり上の3次関数はt=0の所で法線の変化量が不連続になっているんです。法線が不連続に変化するという事はライティングした時に不連続な所が出るという事で、見た目のクオリティに諸に影響してしまいます。
一方5次関数は2階微分で-120t^3+180t^2+60tとなり、t=0で値がゼロになります。よって法線の変化の仕方もスムーズになり、t=0でのライティングのアーティファクトを防ぐことができます。
もう一つの改良はパフォーマンスの向上（アルゴリズムの最適化）のための妥協です。パーリンノイズの勘所であり唯一のパラメータは「勾配」です。この勾配に適当な値を振る事で繰り返し感の無いノイズが出来るのでした。勾配は以下の平面を定義するのに使われます：
この右辺は内積の形をしていて勾配と座標の「掛け算」になっています。一つのウェーブレット関数には一つの平面が必要で、2次のパーリンノイズの場合は四つ角があるので4回この内積計算を行う必要があります。3次の場合は立方体の角なので8回。この掛け算のコストを減らすために、Perlinは大胆にも勾配を「矩形の中心から辺へ伸びるベクトル」のみに減らす事を提唱しました（正確には3次のパーリンノイズに対してです）。
2次のパーリンノイズの場合、矩形は正方形です。その中心から辺へ向かうベクトルというのは次の4パターンしかありません：
この単位長な4つの勾配しか使わないとすると、上の平面の式からは±xと±yのいずれかしか出てこない事になります。よって、ランダムで上の4つの勾配を選んだら、内積を計算するのではなくて、その勾配に該当する±xか±yのどれかを返すだけで良くなります。ちなみに、Perlinの原著では3次のパーリンノイズで説明されていて、その場合の勾配ベクトルは以下の図のようになります：
ただですねぇ…実際2Dバージョンで上の4つの勾配を用いてパーリンノイズを描画してみると、正直使い物になりません。理由は方向が単純過ぎるためです。実際に描画した結果をご覧ください：
左が通常のパーリンノイズ、真ん中が勾配が4つだけの場合の改良パーリンノイズ、そして右は8方向に増やしてみたバージョンです。通常のパーリンノイズは左のように一つの角に対して任意の勾配方向を設定します。勾配ベクトルの長さを正規化すればそれは円状になるわけです。これにより波の方向が様々に重なり良い感じのノイズが作られます。一方改良パーリンノイズの妥協を2Dの場合に適用した真ん中の場合だとあまりに波のバリエーションが単調過ぎて残念なノイズになっています。そこで、右のように勾配を8方向に増やすと多少改善が見られました。でもやっぱり右のでもまだ人工的な連なりを感じる箇所があります。
実験結果から鑑みると、この勾配を妥協する改良パーリンノイズは2Dの場合だとあまり良くないのかもしれません。んー (-_-;
⑦ パーリンノイズでフラクタル
パーリンノイズは大小の周波数のノイズを重ね合わせてフラクタル化した時に真の真価が発揮されます。パーリンノイズの周波数は矩形のサイズで調整ができます。もちろん様々な周波数のを重ね合わせても良いのですが、典型的には一番大きな矩形の周波数を1とした時に、その半分の1/2、さらに半分の1/4、1/8…とサイズを半々にして周波数を倍に上げていきます。
また個々の周波数で生成するパーリンノイズにはスケールを掛けてから重ね合わせるのが一般的です。例えば自然にある大きな山の高さに対してその表面のもっと小さい山の高さはうんと小さいはずです。山の形が相似形なら、フラクタル図形の縮小スケールがそのまま山の高さのスケールになるはずです。つまり周波数をω(i)=2^iとした時（iは整数）、重ね合わせるノイズのスケールをその逆数（1/ω）にすると相似形なフラクタル図形になります。実際に単純な2次のパーリンノイズをフラクタル化してみたのがこちら：
とっても良い感じ(^-^)
雲のようでもあり、山岳地帯を上から眺めたような感じでもあり、岩肌のようでもあり、何となく自然に見受けられるような見た目ですよね。ちなみに、Photoshopを触ったことがある方ならば「あれ？これどっかで見たことある」と思うかもしれません。そうPhotoshopのフィルタの定番である「雲模様1」とそっくりです。ただPhotoshopの方は幅高をさ256単位にすると上下左右でループする模様を作ってくれますが、これはそういうループテクスチャにはなっていません。でもパーリンノイズから作られるのをループテクスチャ化するのは割と簡単で、左右及び上下端の矩形の勾配を一緒にしてあげるだけで出来てしまいます。
⑧ 矩形にランダムだけど固定的な勾配を与えるには
パーリンノイズの要は「勾配をランダムに振る」です。これは単純に考えれば、勾配の値を決める時に疑似乱数から値をゲットするだけのような気がしてしまいますが、良く考えるとそれはあまり上手くいかない事がわかります。疑似乱数は最初にランダムシード（乱数の種）を与え、後は乱数を引き続けます。もしあるパーリンノイズのデータ群を一度だけ取得するような作りであれば何とかなりそうですが、多くの場合パーリンノイズは「座標(x, y)でのノイズ値を取得する」という使い方をします。座標(x, y)が所属する矩形の角ではいつでも同じ勾配設定になっていないといけません。つまり「座標(x, y)に対応した固定的な勾配」が得られなければならない訳で、その時の乱数を何らかの方法で覚えておくか再現できないといけないんです。
2Dパーリンノイズの場合を例に考えてみます。座標(x, y)が所属する矩形の境界（角）は4つあります。単純の為座標(x, y)の整数部分を矩形の左下の境界座標としましょう：
座標(x, y)が与えられたらその境界座標がわかります。各境界座標にはランダムな勾配を振る必要がありますが、ノイズ計算時には同じ境界座標はいつも同じ勾配でないといけません。これを実現するには「整数座標(i, j)に対して勾配(ax, ay)を返すハッシュ」が必要になるという事です。
一番単純なのはすべての境界座標に対して最初にランダムな勾配を与えてしまい、その値をキャッシュする方法です。しかしこれは境界座標分のバッファを要するのでメモリを多量に食います。また新しい境界座標に対応できませんし、別のパーリンノイズを生成する度にバッファが必要になります。限定的な範囲では可能とは思いますがちょっとコストが高そうです。
この問題の解決方法としてPerlinが提案しているのは「ハッシュテーブル」です。最初に0～255までの数値が一つずつ格納された整数配列を用意し、それをシャッフルします（256個という整数配列のサイズはPerlin曰く「十分な大きさ」だそうです）。そして境界座標から次のように勾配値を算出します：
/** 勾配成分を取得
* i, j: 境界整数座標
* elem: 勾配の軸番号（x軸->0、y軸->1など）
*/
float getGrad( int i, int j, int elem ) {
// hashは0～255がシャッフルされた配列
int hashX = hash[ i % 256 ];
int hashXY = hash[ ( hashX + j ) % 256 ];
int hashGrad = hash[ ( hashXY + elem ) % 256 ];
return ( hashGrad / 255.0f ) * 2.0f - 1.0f; // -1～1の範囲を返す
}
まず整数成分iを0～255に丸めてそのハッシュ値を取得します。ここで1回目のシャッフルをしているという事です。得られたシャッフル値に成分jを足し、0～255に丸めて2回目のシャッフルを行います。最後に軸の番号を使って3回目のシャッフルをしてそれを-１～1の範囲に変換した値を返します。簡易的な方法ではありますが割と上手くいきます。ただ、早い段階で周期性が出ますので注意が必要です。例えばgetGrad(0, 0, 0)とgetGrad(256, 0, 0)は最初のhashXが同じ値になるため1回目のシャッフル効果が無くなります。つまり256セル単位でX軸方向に全く同じ勾配の並び方が出現する事になります。
多分実用の範囲内で困る事はあまりないとは思いますが「アルゴリズムとしてそれはどうかなぁ」とちょっと悩みまして、Webをあれこれ調べてみた所、ハッシュテーブルに取って代われそうな手法を見つけました。それが「Xorshift疑似乱数」です。Xorshift疑似乱数は乱数生成方法の一つでGeorge Marsagliaが2003年に発表しました。次のようにXOR演算とビットシフトだけで乱数を発生させます：
/** 引数の整数をランダムなハッシュ値に変換
* xorshift 32bit版
*/
unsigned xorshift32( unsigned v ) {
v = v ^ ( v << 13 );
v = v ^ ( v >> 17 );
v = v ^ ( v << 15 );
return v;
}
引数の任意の整数値に対して32bitの適当な疑似乱数を返します。加減乗除が無く3行で終わるため高速です。「これ偏らないの？」と思うかもしれませんが、意外や意外結構良い乱数列になるそうです。Google ChromeがJavaScriptエンジンで使われている乱数をこれに採用し直したという実績もある方法だそうです。
実際にこの乱数を使って勾配を振る方法の一つが以下の実装です。疑似乱数が返す値は32bitの整数ですが、適当な範囲に丸めて使った方が良さそうです：
/** 引数の整数座標に対応するランダムな傾き値を取得
* -1.0f～1.0fの範囲
*/
float getGrad( int x, int y, int elem, int seed = 0 ) {
unsigned r = xorshift32( x + seed );
r = xorshift32( r + y );
r = xorshift32( r + elem ) % 10000;
return (float)r / 5000 - 1.0f;
}
Xorshift乱数を使う利点はハッシュバッファを用意しなくて良い事と、動的に計算するのでseedを変えるだけで全く別のパーリンノイズがすぐに得られるという点にあります（ハッシュテーブルを用いる場合は都度ハッシュテーブルを作り直す必要がある）。
パーリンノイズとフラクタルは大理石や雲などを表現する手続きテクスチャ（プロシージャルテクスチャ）や、地形や波などの手続きモデル（プロシージャルモデル）の基盤を支えています。ここはアルゴリズムを紹介する編なので、適用例は別の編に譲るとします。パーリンノイズ、自由に出力して適材適所で使うと表現の幅が広がりそうです。
