Chapter 46

パーリンノイズ

miku
miku
2021.12.09に更新

noise() 自体の利用方法は以前の章で扱ったが、この章では noise() で利用されているアルゴリズムであるパーリンノイズの実装について学ぶ。

パーリンノイズ

noise(x, y, z);

パーリンノイズでは (x, y, z) の座標を受け取るので3次元での計算になるのだが、わかりやすいようにまずは2次元でのパーリンノイズの実装を行う。

2次元パーリンノイズ

noise(x, y);

引数として (x, y) を受け取る。これは引数名のとおり座標として扱われる。

(x, y) それぞれの値以下の最大の整数と、値より大きい最小の整数を計算する。これは floor(x), floor(x) + 1 で取得できる。この値を組み合わせて座標にすると下記の4つの点ができる。

  • (x1, y1) = (floor(x), floor(y))
  • (x2, y1) = (floor(x) + 1, floor(y))
  • (x1, y2) = (floor(x), floor(y) + 1)
  • (x2, y2) = (floor(x) + 1), floor(y) + 1)

これらの座標は整数なので、ここでは格子点と呼ぶ。

4つの格子点は座標 (x, y) を囲むような形になるだろう。2次元パーリンノイズの計算は座標 (x, y) と、それを囲っている4つの格子点で行われる。

次に、格子点ごとに固有のベクトルが付いていると考える。このベクトルは4種類あり、下記の通りとなる。

ベクトルの成分 向き
(-1, -1) 左上
(1, -1) 右上
(-1, 1) 左下
(1, 1) 右下

格子点ごとに上記のベクトルのどれかが付いている。これは格子点の座標によって決まる。

次に、各格子点から点へのベクトルを求める。

ここで各格子点と点で計算を行うのだが、計算方法は一緒なので、格子点は左上の部分にだけ注目する。今、格子点固有のベクトルと点へのベクトルの2つのベクトルがある。これらのベクトルの向きを比べて、向きが離れているほど小さい値に、逆に向きが近いほど大きい値にしたい。この計算は内積を利用する。

他の3つの格子点も同様に計算すると、4つの値ができることになる。これを v0, v1, v2, v3 とする。ノイズの戻り値としては1つの値を返したいので、たとえば4つの平均値を取るということもできるが、対象の点が左上の格子点に寄っているので、左上の v0 の影響を強くしたい。このため4つの点を2回線形補間して1つの点にする。

線形補間するには点の座標が軸ごとに 0~1 のどの位置にあるかを知らないといけないが、格子点ごとの幅と高さは 1 なので、点の小数部がそのまま割合になる。たとえば点が (12.34, 45.67) の場合は、囲まれた中では (34%, 67%) の地点にいることがわかる。

その割合を利用して2回の線形補間をすれば一つの値にまとまる。これを返り値とする。

補間関数

内積で求めた v0, v1, v2, v3 を1つの値にする際の補間は普通の線形補間を行う。ただし、点の位置を割合で求める際、x - floor(x) のように計算するが、この値を滑らかにするために下記の補間関数を適用する。

補間関数は線形補間ではなく作者が定義したsmoothstepを利用する。そのsmoothstepの多項式が 6x^5-15x^4+10x^3 になり、x^3 で括ると、x^3(6x^2-15x+10)、更にそこから 6x^2-15xx で括ると、x^3(x(6x-15)+10) となるので、あとはその式をコードに落とす。

function smootherstep(x) {
  return x * x * x * (x * (x * 6 - 15) + 10);
}

// 割合(小数部)を求めた後、更に滑らかな補間関数を適用する
const ratioX = smootherstep(x - floor(x));
const ratioY = smootherstep(x - floor(x));

内積

2つのベクトルの向きが離れているほど小さい値に、向きが近いほど大きい値にするために内積の計算を利用すると説明したが、具体的にどのような計算になるのかを確認しよう。

function dot(a, b) {
  return a.x * b.x + a.y * b.y;
}

内積は、ベクトルの成分同士をかけ合わせて、その和を返すという計算になる。

格子点に付いている固有のベクトルの成分は -11 のどちらかなので、点へのベクトルの成分の符号がそのままか逆になったものを足し合わせるという計算になるだろう。具体的には下記の通りだ。

格子点から点へのベクトル 格子点固有のベクトル 内積の計算結果
(x, y) (-1, -1) -x - y
(x, y) (1, -1) x - y
(x, y) (-1, 1) -x + y
(x, y) (1, 1) x + y

次に、格子点から点へのベクトルの成分の範囲を考える。

格子点で囲まれた範囲なので、成分が取れる範囲は -0.999...~0.999... になる。ここでは便宜上 -1~1 とする。

x, y ともに -1~1 の値になりえて、それを足し合わせるのだから、内積の結果は -2~2 になるだろう。

4つの格子点で内積の結果である -2~2 の範囲の値が計算され、それが補間されるのだが、実は補間結果は -2~2 にはならない。

左上の格子点の内積が最大値の 2 になるのは、点が対角の格子点に限りなく近づいた場合だ。この場合、右下の格子点の内積は 0 に近づくので、左下・右上がどのような内積の値であっても、補間結果は 0 に近づくことになるだろう。

補間することを加味すると、最も補間結果の値が大きくなるのは、真ん中に点を持ってきた場合だ。このとき左上の格子点から点へのベクトルの成分は (0.5, 0.5) となり、内積の最大値が 1 になる。これは他の格子点の内積の結果も同様で最大値が 1 になるので、補間結果も 1 になる。符号を考えると最小値は -1 になるので、結果、補間した値の範囲は -1~1 ということになる。

パーリンノイズの戻り値は 0~1 の範囲にしたいので、範囲に 1 を足して 0~2 に変換後、2 で割り 0~1 にする。マップ関数を利用するなら map(v, -1, 1, 0, 1) になる。

格子点固有のベクトル

これが最も厄介な仕様で、これのせいでパーリンノイズのコードがややこしくなっている節がある。本家通りの実装については後に解説するので、まずは単純な実装で固有のベクトルを求めるようにしよう。

たとえばプログラム実行中に初めて対象の格子点固有ベクトルを求める際、ランダムにベクトルを決めて、そのベクトルそのものかインデックスを記憶しておくという方法がある。

もしくはシードが指定できる疑似乱数、たとえばXorShiftなどを利用して、座標を1つの値にまとめてシードとして渡して、最初の乱数をインデックスにするというのも考えられる。

下記コード例では xy を足して下位2bitの値を見て判断する方法にしている。

// index = 格子点座標を(x + y) & 3にして下位2bitだけ取り出したもの
// x, y = 格子点から点までのベクトル
function calcDot(index, x, y) {
  if (index === 0) {
    // (-1, -1)
    return -x - y;
  } else if (index === 1) {
    // (-1, 1)
    return -x + y;
  } else if (index === 2) {
    // (1, -1)
    return x - y;
  } else {
    // (1, 1)
    return x + y;
  }
}

2次元パーリンノイズのコード例

function setup() {
  createCanvas(windowWidth, windowHeight);
  noStroke();

  const size = 10;
  for (let y = 0; y < 100; y++) {
    for (let x = 0; x < 100; x++) {
      const v = PerlinNoise.noise(x / 100, y / 100);
      fill(v * 255);
      rect(x * size, y * size, size, size);
    }
  }
}

const PerlinNoise = {
  noise: function (x, y) {
    // (x, y) を囲む4つの格子点
    const tl = createVector(floor(x), floor(y));
    const tr = createVector(floor(x) + 1, floor(y));
    const bl = createVector(floor(x), floor(y) + 1);
    const br = createVector(floor(x) + 1, floor(y) + 1);

    // 4つの格子点と、4つの各格子点から点までのベクトルから、4つの内積を求める
    const v0 = PerlinNoise.calcDot((tl.x + tl.y) & 3, x - tl.x, y - tl.y);
    const v1 = PerlinNoise.calcDot((tr.x + tr.y) & 3, x - tr.x, y - tr.y);
    const v2 = PerlinNoise.calcDot((bl.x + bl.y) & 3, x - bl.x, y - bl.y);
    const v3 = PerlinNoise.calcDot((br.x + br.y) & 3, x - br.x, y - br.y);

    // 割合(小数部)を求めた後、更に滑らかな補間関数を適用する
    const ratioX = PerlinNoise.smootherstep(x - floor(x));
    const ratioY = PerlinNoise.smootherstep(y - floor(y));

    // v4 = 左上と右上の補間
    // v5 = 左下と右下の補間
    // v = v4とv5の補間
    const v4 = lerp(v0, v1, ratioX);
    const v5 = lerp(v2, v3, ratioX);
    const v = lerp(v4, v5, ratioY);

     // v = -1~1なので、1を足して0~2に、2で割って0~1の範囲に変換する
    return (v + 1) / 2;
  },

  smootherstep: function (x) {
    return x * x * x * (x * (x * 6 - 15) + 10);
  },

  // index = 格子点座標を(x + y) & 3にして下位2bitだけ取り出したもの
  // x, y = 格子点から点までのベクトル
  calcDot: function (index, x, y) {
    if (index === 0) {
      return -x - y;
    } else if (index === 1) {
      return -x + y;
    } else if (index === 2) {
      return x - y;
    } else {
      return x + y;
    }
  },
};

パーリンノイズの簡易的な実装例。

パーリンノイズの動きのイメージ

コード全体を見る

const size = 5;
const ax = 20;
const ay = 20;
const w = 300;
const h = 300;
let tx, ty;
let v;
let layer;
let dx;

function setup() {
createCanvas(windowWidth, windowHeight);
stroke(240);
noFill();

v = [createVector(1, 1), createVector(0, 1), createVector(1, 0), createVector(0, -1)];

tx = 0;
ty = 0;
dx = 1;

layer = createGraphics(width, height);
}

function calcDot(index, x, y) {
if (index === 0) {
return -x - y;
} else if (index === 1) {
return -x + y;
} else if (index === 2) {
return x - y;
} else {
return x + y;
}
}

function draw() {
if (ty >= h) {
return;
}
clear();

stroke(240);

tx += size * dx;
if (tx >= w) {
tx = w - size;
dx *= -1;
ty += size;
} else if (tx < 0) {
tx = 0;
dx *= -1;
ty += size;
}

const ratioX = smootherstep(norm(tx, ax, ax + w));
const ratioY = smootherstep(norm(ty, ay, ay + h));

const v1 = v[0].dot(ratioX, ratioY);
const v2 = v[1].dot(ratioX, ratioY);
const v3 = v[2].dot(ratioX, ratioY);
const v4 = v[3].dot(ratioX, ratioY);

const ux = lerp(v1, v2, ratioX);
const bx = lerp(v3, v4, ratioX);
const value = map(lerp(ux, bx, ratioY), -1, 1, 0, 1);

layer.translate(ax, ay);
layer.noStroke();
layer.fill(value * 255);
layer.rect(tx, ty, size, size);
layer.resetMatrix();
image(layer, 0, 0);

noFill();
rect(ax, ay, w, h);

drawArrow(0, 0, v[0].x * 50, v[0].y * 50);
drawArrow(w, 0, w + v[1].x * 50, v[1].y * 50);
drawArrow(0, h, v[2].x * 50, h + v[2].y * 50);
drawArrow(w, h, w + v[3].x * 50, h + v[3].y * 50);
drawArrow(0, 0, tx, ty);
drawArrow(w, 0, tx, ty);
drawArrow(0, h, tx, ty);
drawArrow(w, h, tx, ty);
}

function smootherstep(x) {
return x * x * x * (x * (x * 6 - 15) + 10);
}

function drawArrow(x1, y1, x2, y2) {
const d = dist(x1, y1, x2, y2) - 10;
const angle = atan2(y2 - y1, x2 - x1);

stroke(240);
noFill();

translate(ax + x1, ay + y1);
rotate(angle);
line(0, 0, d, 0);
line(d, 0, d - 5, -5);
line(d, 0, d - 5, 5);
resetMatrix();
}

パーリンノイズの補間を可視化したものが上記になる。次に本家のアルゴリズム通りに実装していく。

本家パーリンノイズの特徴

noise(x, y, z);

本家パーリンノイズは引数を3つ指定する。なので今まで行ってきた計算を3次元に拡張して行わなければならない。

(x, y, z) には任意の実数を指定できるのだが、内部では 0~255 の範囲に丸められる。これは格子点が各軸ごとに256種類しか無いからだ。

格子点の数は 256 * 256 * 25616777216 個あり、格子点固有のベクトルは下記の16種類のどれかが付いている。しかも均等に割り振られているので、下記の固有ベクトルの1つ1つは1048576個の格子点に利用されている。

// (x, y, z)
(1,1,0), (-1,1,0), (1,-1,0), (-1,-1,0)
(1,0,1), (-1,0,1), (1,0,-1), (-1,0,-1)
(0,1,1), (0,-1,1), (0,1,-1), (0,-1,-1)

格子点固有のベクトルと、格子点から点までのベクトルの内積を取るのは同じで、点までの位置を計算するのにsmootherstepを利用して補間するのも同様である。ただし、点の周りの格子点が8個に増えているので、補間回数も3回から7回に増えている。

最後に本家の大事な特徴があり、たとえば x3.2 の場合は、x軸だけで考えると周りの格子点は 34 になるが、x255 の場合は、255256 ではなく 2550 になりループする。これは前述したとおり本家の格子点は 0~255 の256種類しかないからだ。

次に、これらの特徴を踏まえて、本家の実装方法を確認していく。

本家の格子点から固有ベクトルを取り出すまで

const permutation = [
151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225,
140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148,
247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32,
57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175,
74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122,
60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54,
65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169,
200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64,
52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212,
207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213,
119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104,
218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241,
81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, 106, 157,
184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, 205, 93,
222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180
];

本家のパーリンノイズのコードには上記の 0~255 を並び替えた256種類の値が入った配列がハードコーディングされている。これは適当にシャッフルしたものではなく均等に格子点固有のベクトルを決めるためのものだ。

前述した簡易版では、格子点固有のベクトルを求める際に、格子点固有のインデックスが必要なので、格子点の位置を足し合わせて下位2bitを取り出したものをインデックスにするという作業を行った。それに対して本家では、格子点の位置から上記の256種類の値のどれかを取り出して内積を計算するためのインデックスとして扱う。

// p = permutation.concat(permutation)
p[p[p[x] + y] + z]

格子点座標が (x, y, z) だった場合の取り出し方法が上記になる。p[] は先程の256種類の値が入った permutation[] 同士を結合したもので、512個の値に増えている。これは p[x] で参照できる値の範囲が 0~255 で更にその値に y を足しているので permutation[] のままだと配列の範囲外に参照することになるからだ。単に結合しているだけだが、うまくシャッフルされているので先程述べたとおり各固有ベクトルに均等にインデックスが割り振られるようにできている。

function grad(index, x, y, z) {
  const h = index & 15;
  const u = h < 8 ? x : y;
  let v;

  if (h < 4) {
    v = y;
  } else if (h == 12 || h == 14) {
    v = x;
  } else {
    v = z;
  }

  return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v);
}

grad() は簡易版の calcDot() と同じ役割を持つ関数である。引数の index には p[p[p[x] + y] + z] を指定する。なので index は格子点の位置を代替する値だと思っていい。残りの引数 (x, y, z) は各格子点から点までのベクトルになる。これらの引数から16種類ある固有ベクトルのうちどれか1つを返す。

ここまでで本家流の格子点から固有ベクトルを取り出すところまでを扱った。次に noise(x, y, z) を呼び出して1つの値を返すまでの全体の流れを解説する。

本家のコード全体の流れ

まず、渡した (x, y, z) から8つの格子点を求める必要がある。2Dの場合、floor(x), floor(x) + 1 のように計算したが本家では格子点がループしていることと、格子点座標は 0~255 なので負数にはならないという特徴があるので、このようには計算できず、実際には下記のように書く。

const xi = x & 255;
const yi = y & 255;
const zi = z & 255;

これは簡易版の floor(x) にあたる機能で、各軸の小さい方の格子点座標を求めるのに必要なものだ。このコードには3つの意味がある。

1つ目は格子点座標が 255 より大きい値にならないように 255 と論理積を取っていること。2つ目は小数点以下を切り捨てること。これは論理積を取ると自動で切り捨てれられる。最後は負の値ならマイナスを外すことだ。2進数の一番左のビットが 1 なら負数として扱うが、255 と論理積を取ることで負数ではなくなる。

計算した (xi, yi, zi) それぞれの値に 1 を足すか足さないかで、8種類の格子点を求めることができるが、ループ対応ができていない。なので、もし +1 をすることで格子点座標が 256 になったら 0 に戻す必要がある。

function inc(num) {
  return (num + 1) % PerlinNoise.permutation.length;
}
(xi, yi, zi)
(inc(xi), yi, zi)
(xi, inc(yi), zi)
(xi, yi, inc(zi))
(inc(xi), inc(yi), zi)
(xi, inc(yi), inc(zi))
(inc(xi), yi, inc(zi))
(inc(xi), inc(yi), inc(zi))

その処理を関数 inc() としてまとめると、8つの座標は上記のとおりに計算できる。

const p000 = p[p[p[xi] + yi] + zi];
const p100 = p[p[p[inc(xi)] + yi] + zi];
const p010 = p[p[p[xi] + inc(yi)] + zi];
const p001 = p[p[p[xi] + yi] + inc(zi)];
const p110 = p[p[p[inc(xi)] + inc(yi)] + zi];
const p011 = p[p[p[xi] + inc(yi)] + inc(zi)];
const p101 = p[p[p[inc(xi)] + yi] + inc(zi)];
const p111 = p[p[p[inc(xi)] + inc(yi)] + inc(zi)];

格子点から格子点固有のインデックスを求める。

const xf = x - floor(x);
const yf = y - floor(y);
const zf = z - floor(z);

実数から整数部を引くことで小数部だけを取り出す。

const u = PerlinNoise.smootherstep(xf);
const v = PerlinNoise.smootherstep(yf);
const w = PerlinNoise.smootherstep(zf);

0以上1未満の値に対して補間をかける。

inc: function (num) {
  num++;
  if (PerlinNoise.n > 0) num %= PerlinNoise.permutation.length;

  return num;
},
  
const aaa = p[p[p[xi] + yi] + zi];
const aba = p[p[p[xi] + inc(yi)] + zi];
const aab = p[p[p[xi] + yi] + inc(zi)];
const abb = p[p[p[xi] + inc(yi)] + inc(zi)];
const baa = p[p[p[inc(xi)] + yi] + zi];
const bba = p[p[p[inc(xi)] + inc(yi)] + zi];
const bab = p[p[p[inc(xi)] + yi] + inc(zi)];
const bbb = p[p[p[inc(xi)] + inc(yi)] + inc(zi)];

(x, y, z) に対する固有のベクトル番号を得る。このように取得すると、順列テーブルと兼ね合いで均等に勾配ベクトルが決まる。

3次元の場合、(x, y, z) を囲む格子点は8点あり、(x, y, z) 各値以下の最小の整数は (xi, yi, zi) である。あとはそれぞれに 1 を足した (xi + 1, yi + 1, zi + 1) が必要だが、順列テーブルの最大値は 255 なので、たとえば xi255 の場合は xi + 10 にしてループにする。

const xf = x - floor(x);
const yf = y - floor(y);
const zf = z - floor(z);

const u = PerlinNoise.smootherstep(xf);
const v = PerlinNoise.smootherstep(yf);
const w = PerlinNoise.smootherstep(zf);

(xf, yf, zf)(x, y, z) の小数部を取り出したものだ。これらの値は別に必要になるので変数に入れている。更にこの値を滑らかにするために smootherstep() を適用したものを (u, v, w) とする。

const g000 = grad(p000, xf, yf, zf);
const g100 = grad(p100, xf - 1, yf, zf);
const g010 = grad(p010, xf, yf - 1, zf);
const g001 = grad(p001, xf, yf, zf - 1);
const g110 = grad(p110, xf - 1, yf - 1, zf);
const g011 = grad(p011, xf, yf - 1, zf - 1);
const g101 = grad(p101, xf - 1, yf, zf - 1);
const g111 = grad(p111, xf - 1, yf - 1, zf - 1);

grad() に、格子点の位置を代替するインデックス (pXXX) と、位置を割合にして滑らかにした (xf, yf, zf) を渡す。xf - 1 のように 1 を引いている箇所があるのは、点から格子点までのベクトルを表すからだ。これで grad() は内積の結果を返す。

const x0 = lerp(g000, g100, u);
const x1 = lerp(g010, g110, u);
const x2 = lerp(g001, g101, u);
const x3 = lerp(g011, g111, u);
const y0 = lerp(x0, x1, v);
const y1 = lerp(x2, x3, v);
const z0 = lerp(y0, y1, w);

return (z0 + 1) / 2;

8つの格子点ごとに内積の結果が得られたので、点の位置まで1つの値になるように線形補間を7回行う。線形補間するには割合を表す t が必要だが、これは先程計算した (u, v, w) を利用する。

本家パーリンノイズの実装例

function setup() {
  createCanvas(windowWidth, windowHeight);
  stroke(240);
  noFill();

  PerlinNoise.init();

  let prev = PerlinNoise.noise(0, 0, 0, 4, 0) * height;
  for (let i = 1; i < width; i++) {
    const v = PerlinNoise.noise(i / 100, 0, 0, 4, 0.5) * height;
    line(i - 1, prev, i, v);

    prev = v;
  }
}

const PerlinNoise = {
  permutation: [
    151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
    129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180,
  ],
  p: [],

  init: function () {
    PerlinNoise.p = PerlinNoise.permutation.concat(PerlinNoise.permutation);
  },

  smootherstep: function (x) {
    return x * x * x * (x * (x * 6 - 15) + 10);
  },

  inc: function (num) {
    num++;
    if (PerlinNoise.n > 0) num %= PerlinNoise.permutation.length;

    return num;
  },

  grad: function (code, x, y, z) {
    const h = code & 15;
    const u = h < 8 ? x : y;
    let v;

    if (h < 4) {
      v = y;
    } else if (h == 12 || h == 14) {
      v = x;
    } else {
      v = z;
    }

    return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v);
  },

  noise: function (x, y = 0, z = 0) {
    const xi = x & 255;
    const yi = y & 255;
    const zi = z & 255;

    const xf = x - floor(x);
    const yf = y - floor(y);
    const zf = z - floor(z);

    const u = PerlinNoise.smootherstep(xf);
    const v = PerlinNoise.smootherstep(yf);
    const w = PerlinNoise.smootherstep(zf);

    const p = PerlinNoise.p;
    const inc = PerlinNoise.inc;
    const grad = PerlinNoise.grad;
    const slerp = PerlinNoise.slerp;

    const p000 = p[p[p[xi] + yi] + zi];
    const p100 = p[p[p[inc(xi)] + yi] + zi];
    const p010 = p[p[p[xi] + inc(yi)] + zi];
    const p001 = p[p[p[xi] + yi] + inc(zi)];
    const p110 = p[p[p[inc(xi)] + inc(yi)] + zi];
    const p011 = p[p[p[xi] + inc(yi)] + inc(zi)];
    const p101 = p[p[p[inc(xi)] + yi] + inc(zi)];
    const p111 = p[p[p[inc(xi)] + inc(yi)] + inc(zi)];

    const g000 = grad(p000, xf, yf, zf);
    const g100 = grad(p100, xf - 1, yf, zf);
    const g010 = grad(p010, xf, yf - 1, zf);
    const g001 = grad(p001, xf, yf, zf - 1);
    const g110 = grad(p110, xf - 1, yf - 1, zf);
    const g011 = grad(p011, xf, yf - 1, zf - 1);
    const g101 = grad(p101, xf - 1, yf, zf - 1);
    const g111 = grad(p111, xf - 1, yf - 1, zf - 1);

    const x0 = lerp(g000, g100, u);
    const x1 = lerp(g010, g110, u);
    const x2 = lerp(g001, g101, u);
    const x3 = lerp(g011, g111, u);
    const y0 = lerp(x0, x1, v);
    const y1 = lerp(x2, x3, v);
    const z0 = lerp(y0, y1, w);

    return (z0 + 1) / 2;
  },
};