ランダム地形生成 Part2～フラクタルブラウン運動

#パーリンノイズだけで出来ることは限られています
前回話したパーリンノイズ一見自然的ですが、でも実際使って見たら効果がそんなに理想的ではありません。
実際の効果は話より画像を見た方が分かりやすいと思いますので
以下のノイズとそのノイズを使った地形をご覧ください。

単純なパーリンノイズだと見た目はやはり。。。とてもノイズ的ですね。
もちろん、波やテクスチャなどとして使うのなら、これても結構ですが
しかし、リアル的な地形を目標としたら、これだけでは満足できません。
#パーリンノイズとリアル地形の差
では、一体パーリンノイズとリアル地形との差はどうこでしょう。

リアルの地形とパーリンノイズの効果を比べたら、その差がすぐ分かると思います。
リアルの方が丘や山や坂などありまして、それぞれの間、高さの落差が結構あると見えます。
もう一方のパーリンノイズ地形ではそうではありません。
高さの差がほぼないと言うか、みんな一律丘にしか見えません。

これ以外、リアル地形は全体から見て高さの変動がありだけではなく（ピークと山腹との高さの差か）
局部から見ても高さの変動があります（同じ山腹だとしても、高さがちょっとだけ差があります）。
つまり、似ているパターン（高さの波動）が尺度と関わらずに存在していることです。

しかし、パーリンノイズで生成した地形は全体から見て高さの波動がありますが
局部から見ると、相当スムーズです。
残念ながら、地形生成においてスムーズが必ずいいこととは限らません。
#フラクタルとフラクタルブラウン運動
上記の特徴（尺度と関わらずに似ているパターンが存在すること）はフラクタルと言われていて、自然界によくある現象です。

そして、それを利用し、真似する方法はフラクタルブラウン運動です。
もちろん、コンピューティングリソースが有限である限り、これはあくまでもシミュレーションに過ぎません。
フラクタルブラウン運動の原理はいくつ頻度（ frequency）と振幅（ amplitude）が違う八度
（octave、ここではノイズのことを指していますので、以下はノイズと言います）を重ねてフラクタルと似た結果を作り出すことです。
この過程で使われているノイズの数が多いほど、結果ももっとリアルに見えますが、計算するにかかる時間も長くなります。
#パーリンノイズ+フラクタルブラウン運動
では、パーリンノイズにおいた頻度と振幅とはどういうことでしょう。
フラクタルブラウン運動にある頻度とは一定範囲内にノイズの波動回数です。
ざっくり言うと、ハイトマップの白と黒の変化密度。
パーリンノイズの場合、それを影響するのはコントロールポイントの数です。
つまり、頻度の調整はコントロールポイントの数の調整です。

振幅とはノイズの変化の激しさです。
同じくハイトマップを例として説明すると
白い部分がどれほど真っ白に近いか、黒い部分がどれほど真っ黒に近いかとのことです。
こちらの調整は結構シンプルで、ノイズをサンプルした結果をスケールすれば大丈夫です。

重ねられるノイズ毎の頻度と振幅は勝手に設定しても構いませんが、効果が比較的にリアル的なやり方は頻度を二倍つづ増やして
振幅を底が0～1の数字となっている指数関数を基に変化させることです。
指数関数の底となる数字が一致性（persistent）と呼ばれています。
この数字が低くほど、最終結果もスムーズになります。

数式で表示すると、第i番のノイズの頻度と振幅は以下となります。

何故頻度を増やす同時に振幅を減らすかと言うと、上の樹の画像を見れば分かると思います。
ブランチが分枝するほど、密度もどんどん高くなります。
代わりに、ブランチのサイズがどんどん小さくなっていきます。
上記の数式はまさにこの現象をシミュレーションすることです。

#実装例
実装例は以下となります。（綺麗に見えるように、前回のコードをちょっとクラス化しました）

NoiseGenerator.cs

public class ComplexNoise2D
{
    private IEnumerable<Noise2D> m_Noises = null;
    private double m_TotalAmplitude = 0.0;

    internal ComplexNoise2D(IEnumerable<Noise2D> noises)
    {
        m_Noises = noises;
            
        foreach (var n in noises)
        {
            m_TotalAmplitude += n.Amplitude;
        }
    }

    public double EvalAt(double x, double y)
    {
        double result = 0.0;
        foreach (var n in m_Noises)
        {
            result += n.EvalAt(x, y);
        }

	    //　未処理の生成結果が-m_TotalAmplitude～m_TotalAmplitudeになっているため
	    //　ここでその結果を正規化しておきます
        return (result + m_TotalAmplitude) / (m_TotalAmplitude * 2);  　
    }
}

public class Noise2D
{
    private Vector<double>[,] m_Gradients;
    private int m_X;
    private int m_Y;
    private double m_Amplitude = 1.0;
    public double Amplitude
    {
        get
        {
            return m_Amplitude;
        }
    }

    internal Noise2D(Vector<double>[,] normalised_gradients, int x, int y, double amplitude)
    {
        m_Gradients = normalised_gradients;
        m_X = x;
        m_Y = y;
        m_Amplitude = amplitude;
    }

    public double EvalAt(double x, double y)
    {
        x = x * m_X;
        y = y * m_Y;

        int i = (int)x;
        int j = (int)y;

        var g00 = m_Gradients[i, j];
        var g10 = m_Gradients[i + 1, j];
        var g01 = m_Gradients[i, j + 1];
        var g11 = m_Gradients[i + 1, j + 1];

        double u = x - i;
        double v = y - j;

        double n00 = g00.DotProduct(Vector<double>.Build.DenseOfArray(new double[] { u, v }));
        double n10 = g10.DotProduct(Vector<double>.Build.DenseOfArray(new double[] { u - 1, v }));
        double n01 = g01.DotProduct(Vector<double>.Build.DenseOfArray(new double[] { u, v - 1 }));
        double n11 = g11.DotProduct(Vector<double>.Build.DenseOfArray(new double[] { u - 1, v - 1 }));

        double nx0 = Interpolate(n00, n10, u, F);
        double nx1 = Interpolate(n01, n11, u, F);

        double nxy = Math.Round(Interpolate(nx0, nx1, v, F), 5, MidpointRounding.AwayFromZero);

        return nxy * m_Amplitude;
   }

    private double Interpolate(double from, double to, double percent, Func<double, double> f)
    {
        double f_out = f(percent);

        return from * (1 - f_out) + to * f_out;
    }

    private double F(double t)
    {
        return (double)(6.0 * Math.Pow(t, 5.0) - 15.0 * Math.Pow(t, 4.0) + 10 * Math.Pow(t, 3.0));
    }
}

public class NoiseGenerator
{
    private IList<Vector<double>> m_Gradients = null;

    public Noise2D Generate2D(int x, int y, double amplitude = 1.0)
    {
	    // 一番右の列と一番下の行にあるコントロールポイントと重ねるポジションのノイズ値を
	    // 計算するために、xとy両方もコントロールポイントをもう一つ入れておきます
	    x = x + 1;
	    y = y + 1;

        Vector<double>[,] gradients = new Vector<double>[x, y];
        for (var i = 0; i < x; ++i)
            for (var j = 0; j < y; ++j)
            {
                gradients[i, j] = CalculateRandomGradientFor(i, j);
            }

        Noise2D noise = new Noise2D(gradients, x, y, amplitude);

        return noise;
    }

    public ComplexNoise2D GenerateComplex2D(int iterations, double persistence, int base_x, int base_y)
    {
        if (iterations < 1)
            throw new Exception("You need to at least iterate once!");

        Noise2D[] noises = new Noise2D[iterations];

        for (int i = 0; i < iterations; ++i)
        {
            int frequency = (int)(Math.Pow(2.0, i) + 0.5);  // 四捨五入をする
            double amplitude = (double)Math.Pow(persistence, i);

            int x = frequency * base_x;
            int y = frequency * base_y;

            noises[i] = Generate2D(x, y, amplitude);
        }

        return new ComplexNoise2D(noises);
    }

    private Vector<double> CalculateRandomGradientFor(int x, int y)
    {
	    // ランダムの勾配を生成する            
    }
}

Program.cs

class Program
{
    static void Main(string[] args)
    {
        int iterations = 6;
        double persistence = 0.5;

        int base_ctrl_points_x = 2;
        int base_ctrl_points_y = 2;

        int height_map_x = 512;
        int height_map_y = 512;

        var generator = new NoiseGenerator();
        var noise = generator.GenerateComplex2D(iterations, persistence, base_ctrl_points_x, base_ctrl_points_y);

        double[] height_map = new double[height_map_x * height_map_y];

        for (int j = 0; j < height_map_y; ++j)
        {
            double normalised_pos_y = (double)j / (height_map_y - 1);

            for (int i = 0; i < height_map_x; ++i)
            {
                double normalised_pos_x = (double)i / (height_map_x - 1);
                int height_map_idx = j * output_y + i;

                height_map[height_map_idx] = noise.EvalAt(normalised_pos_x, normalised_pos_y);
            }
        }

        // ハイトマップを使う・出力...
    }
}

以下はフラクタルブラウン運動を加えたパーリンノイズで生成したハイトマップと地形です。