初めに

　この記事はBlenderのVoronoi Textureがどのような実装をされているかについて調べたことをまとめた記事です。

Voronoi Texture

Voronoi図とは？

　Voronoi図とはある空間上にいくつもの特徴点(Feature Pointまたは母点)と呼ばれる点を配置して、どの特徴点に近いかによって空間を分割したもののことを言います。一般的には最も近い特徴点までの距離(Distance)や特徴点の色(Color)、座標(Position)といった情報をテクスチャとして使用するVoronoi Textureという形でVoronoi図を見ることが多いと思います。

f:id:kinakomoti321:20220224170951p:plain — Voronoi図

　Voronoi図は上記の通り、細胞のような区分けされた不思議な模様を作り出します。

　実はこのVoronoi図の不思議な模様は自然界では度々現れ、キリンの模様やトンボの羽といった生物の模様などといったところでVoronoi図は現れます。

f:id:kinakomoti321:20220224181417j:plain — キリンの模様

f:id:kinakomoti321:20220224182154p:plain — トンボの羽の模様

　CGの分野では有機的な模様の表現によくVoronoiTextureは使われており、他にも疑似コースティクスの表現や泡の構造の生成、ステンドグラスの作成、はたまた機械的な回路の模様などといった他のテクスチャにはない独特な模様を用いることで様々な表現をすることができ、必要不可欠なプロシージャルテクスチャの１つとして使われています。

f:id:kinakomoti321:20220224183314p:plain — VoronoiTexutreでできる表現の例

f:id:kinakomoti321:20220224183306p:plain — VoronoiTexutreでできる表現の例

　BlenderでももちろんVoronoiTextureは使うことができ、シェーダーノードやジオメトリノードなどで使用することができます。このBlenderのVoronoiTextureは様々なプロパティなどを持っており、かなり多機能なものとなっています。

　BlenderはOSSなのでGitHubでソースコードが公開されており、VoronoiTextureの実装部分も見ることができます。この記事ではこのVoronoiTextureがどのように実装されているのかをソースコードを元に書いていこうと思います。

BlenderのVoronoi Texture

　ノードでは３DやF1、Eculidianといったプロパティを持っており、下にパラメーターが設定されています。これらについて各々解説します。

Demention（次元）

　ボロノイを考える空間の次元です。１から４次元まであり、それぞれ選択した空間で特徴点が配置される形になります。パラメータにおけるWやVectorというのはボロノイ図からどこの座標部分の出力を得るかという入力座標になっています。

1D ：　１次元（線）上に特徴点を配置する。この時、入力座標はWというパラメータになる
2D ：　２次元（平面）上に特徴点を配置する。入力座標はVectorで、そのうちｚパラメータは必要ないため無視される。
3D ：　３次元（立体）上に特徴点を配置する。入力座標はVector
4D ：　４次元上に特徴点を配置する。入力座標はVector(３次元ベクトル)に加えてWというパラメータで４次元座標を渡す。

f:id:kinakomoti321:20220223215353p:plain — Dimensions

Future Output (特徴出力）

　最も近い特徴点の距離といった情報だけでなく他にも様々な情報を出力できるようになっています。どんな情報をボロノイ図から出力するのかというのがここで指定する部分となります。

f:id:kinakomoti321:20220223055212p:plain — Future Output

F1

　最も近い特徴点までの距離(Distance)とその特徴点の色(Color)、座標(Position)を出力とします。

f:id:kinakomoti321:20220223220946p:plain — F1

F2

　２番目に近い特徴点までの距離(Distance)とその特徴点の色(Color)、座標(Position)を出力とします。

f:id:kinakomoti321:20220223221209p:plain — F2

Smooth F1

　F1の出力をぼかしたような出力をします。このFutureではSmoothnessというぼかしの度合いを示すパラメータが追加されます。

f:id:kinakomoti321:20220223220741p:plain — Smooth F1

Distance to Edge

　ボロノイ図におけるエッジ（それぞれの特徴点の二等分線）までの距離を出力とします。

f:id:kinakomoti321:20220223221443p:plain — Distance to Edge

N-Sphere Radius

　ボロノイ図におけるボロノイセルに内接する（n次元）球の半径を出力とします。

f:id:kinakomoti321:20220223221645p:plain — N-Sphere Radius

Distance Metric（距離計量）

　今まで考えてきた一般的に言われる「距離」というのはユークリッド距離と呼ばれるものです。しかし数学では２点間の「距離」というのはユークリッド距離以外にも考えられており、別の「距離」の形式でボロノイ図を作ることも可能であるわけです。数学的な「距離」というのはDistance Metric（距離計量）と呼ばれ、F1,F2,smoothF1ではDistance Metricを指定することができます。

　それぞれ以下のように距離計量があり、それによって作られるボロノイ図は

Euclidean : ユークリッド距離
Manhattan : マンハッタン距離
Chebychev : チェビシェフ距離
Minkowski : ミンコフスキー距離　この距離計量ではexpornentというパラメータが付けられます

f:id:kinakomoti321:20220223221834p:plain — Distance Metric

f:id:kinakomoti321:20220223230521p:plain — 距離計量別のVoronoi図(F1)

Celler Noise (Worley Noise)

　Voroni図は特徴点を空間上に分布させますが、完全にランダムに設置した場合では最近傍点の探索が非常にコストがかかり、尚且つ有限の範囲でしかVoronoi図を作ることができません。

　ですが、Blenderではスケールをいくら大きくしてもVoronoi図は得られます。一般的にVoronoiTextureで使用されるVoronoi図は特徴点の配置に制限を課すことで無限に広がり、尚且つコストがかからないCeller Noise(Worley Noiseとも)と呼ばれる手法で実装されています。

　　CellerNoiseについてはThe Book of Shadersに詳しい記事があります。

thebookofshaders.com

　重要な部分だけ抜き出すと、CellerNoiseでは空間をセルという格子状に区分し、セル１つあたりに１つの特徴点をランダムに配置するという形で特徴点を空間内に分布させることでVoronoi図を生成します。実用上では整数で格子を考えます。

f:id:kinakomoti321:20220223232926p:plain — CellerNoiseの特徴点

　このようにすることでF1の探索をする際、２次元空間であれば今いる座標が含まれるセルの周囲３×３のみを探索すれば十分であり、それ以上セルに含まれる特徴点はF1の候補にならないのでどんな場所でも３×３の特徴点の探索でVoronoi図を得ることができます。（３次元とかになれば隣接するセルが3×3×3といったように増えます）

f:id:kinakomoti321:20220223233023p:plain — CellerNoiseの探索

　簡素なGLSLで書いたCellerNoiseのコードを以下に用意しました。

//ハッシュ関数
vec2 Hash2D(vec2 uv){
    vec2 st = vec2(dot(uv,vec2(134.4,314.0)),dot(uv,vec2(932.9,141.301)));
    return fract(sin(st)*39145.295039);
}

float voronoi(vec2 uv,float offset){
// uvがいるセルの番号 i_st、そのセル内での位置 f_st
    vec2 i_st = floor(uv);
    vec2 f_st = fract(uv);

// 最小距離を入れる変数
    float m_dist = 1.0;
    
//隣接する3 × 3のセルを探索
    for(int i = -1; i <= 1; i++){
        for(int j = -1; j <=1; j++){
           //現在考えるセルの番号（i_stを(0,0)として中心に考えている）
            vec2 neighbor = vec2(float(i),float(j));
          //(i_st.x + i, i_st.y + j)のセルにおける特徴点の座標を計算する
            vec2 point =  neghbor + Hash2D(i_st + neighbor);
          
         //特徴点との相対ベクトル
            vec2 diff = point - f_st;
         //特徴点までの距離
            float dist = length(diff);

         //最小距離の更新
            m_dist = min(m_dist,dist);
        }
    }
    return m_dist;
}

　処理の流れを見るとまず、uvを受け取り、その整数部分floor(uv)と小数部分fract(uv)を抜き取ります。この時、整数部分である"i_st"はuvの座標が含まれるセルの番号であり、"f_st"はセルの中での座標を示します。

// uvがいるセルの番号 i_st、そのセル内での位置 f_st
    vec2 i_st = floor(uv);
    vec2 f_st = fract(uv);

　次に最短距離を記録する変数"m_dist"を用意して、For文で探索の処理を書いていきます。まず、For文の最初にneighborというセルの番号を調べますが、あくまでi_stに対する相対的な番号なので注意してください。なので調べているセルの絶対的な番号はneighbor + i_stです。

　そして、このセル番号を用いて特徴点の座標をセル内で生成します。この座標の生成にはハッシュ関数を使用しており、セルの絶対番号neighbor + i_stを入れることによってセルの絶対番号に対応する座標を得ることができます。以上のようにしてあるセル neigbor + i_stの特徴点を作ることが出来ました。

           //現在考えるセルの番号（i_stを(0,0)として中心に考えている）
            vec2 neighbor = vec2(float(i),float(j));
          //(i_st.x + i, i_st.y + j)のセルにおける特徴点の座標を計算する
            vec2 point =  neghbor + Hash2D(i_st + neighbor);

　そして最後にやることはこの特徴点とuvの相対ベクトルを得られるので、その距離を計算し、最小距離の更新をmin関数で行っているという形です。

         //最小距離の更新
            m_dist = min(m_dist,dist);

　これがCellerNoiseの手法であり、どんなuvをもらってもちゃんとVoronoi図を返すことができ、無限に続くVoronoi Textureというのは作ることができるわけです。BlenderのVoronoiTextureの実装では特徴出力によって多少変わるものの基本的にこのCellerNoiseの処理を行っています。

f:id:kinakomoti321:20220223233133p:plain — CellerNoiseの流れ

Future Output

　Voronoi Textureのソースコードは2DのF1,3DのF1,4DのF1というように次元と特徴出力の組み合わせごとに別々に関数を作っています。ですが次元に関しては探索のfor文が増減するだけで大して変化しないので、わかりやすい2DのVoronoiに焦点を当てて調べていきます。

F1の実装

　最も近い特徴点である"F1"についての情報を出力するプロパティです。実装は以下の通りになります。

ccl_device void voronoi_f1_2d(float2 coord,
                              float exponent,
                              float randomness,
                              NodeVoronoiDistanceMetric metric,
                              float *outDistance,
                              float3 *outColor,
                              float2 *outPosition)
{
  float2 cellPosition = floor(coord);
  float2 localPosition = coord - cellPosition;

  float minDistance = 8.0f;
  float2 targetOffset = make_float2(0.0f, 0.0f);
  float2 targetPosition = make_float2(0.0f, 0.0f);
  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;
      float distanceToPoint = voronoi_distance_2d(pointPosition, localPosition, metric, exponent);
      if (distanceToPoint < minDistance) {
        targetOffset = cellOffset;
        minDistance = distanceToPoint;
        targetPosition = pointPosition;
      }
    }
  }
  *outDistance = minDistance;
  *outColor = hash_float2_to_float3(cellPosition + targetOffset);
  *outPosition = targetPosition + cellPosition;
}

　多少関数名などが違うのでわかりにくいですので、いくつか関数について補足しておきます。

make_floatA : floatのA次元ベクトルを作る
hash_floatA_to_floatB : A次元ベクトルを受け取り、B次元ベクトルを返すハッシュ関数
voronoi_distance_2d : 受け取った２つのベクトルの距離を返す（詳細は距離計量の実装）

　関数さえわかれば、ほとんどCellerNoiseと同様の処理を行っていることがわかります。しかしながら、F1の出力にはF1までの距離だけではなく、追加としてF1の色や座標がありますので少々異なることを行っています。まず、最小距離を記録する"minDistance"以外にもF1の位置やセルの番号を記録できるよう

  float2 targetOffset = make_float2(0.0f, 0.0f);
  float2 targetPosition = make_float2(0.0f, 0.0f);

という変数を作っており、探索中の距離の更新ではmin関数ではなくminDistanceとの比較でif文の分岐を行い、

      if (distanceToPoint < minDistance) {
        targetOffset = cellOffset;
        minDistance = distanceToPoint;
        targetPosition = pointPosition;
      }

"targetOffset"にはセルの(相対)番号、"targetPosition"には特徴点の(相対)座標を距離と共に更新しています。これによって、最終的にはこれらにはF1のセル番号、座標が記録されており、これを元にF1の色、座標を出力しています。　

  *outDistance = minDistance;
  *outColor = hash_float2_to_float3(cellPosition + targetOffset);
  *outPosition = targetPosition + cellPosition;

カラーの出力でハッシュ関数を使用していますが、cellPosition+targetOffsetはF1のセルの絶対番号であるため、セルの番号ごとに色が出力されるようになっています。また、targetPositionはcellPositionを原点とする相対座標なのでcellPositionを足して、座標を出力しています。このようにDistance以外の出力はF1のセルや座標を元に作られています。

　また、もう一つ異なる点としてはrandomnessというパラメータが存在することです。これは特徴点の配置のランダムさを示しており、１であればハッシュによるランダムな値が特徴点に加えられますが、０になるとその値は一切使わないので特徴点は格子状に整列することになります。

 float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;

f:id:kinakomoti321:20220224184209p:plain — randomnessの変化

f:id:kinakomoti321:20220224184216p:plain — randomnessの変化

　以上がF1の実装となっています

F2の実装

　F2は２番目に近い特徴点のことを示しており、これの情報を出力するのがF2というプロパティとなります。実装は以下の通りになっています。

ccl_device void voronoi_f2_2d(float2 coord,
                              float exponent,
                              float randomness,
                              NodeVoronoiDistanceMetric metric,
                              float *outDistance,
                              float3 *outColor,
                              float2 *outPosition)
{
  float2 cellPosition = floor(coord);
  float2 localPosition = coord - cellPosition;

  float distanceF1 = 8.0f;
  float distanceF2 = 8.0f;
  float2 offsetF1 = make_float2(0.0f, 0.0f);
  float2 positionF1 = make_float2(0.0f, 0.0f);
  float2 offsetF2 = make_float2(0.0f, 0.0f);
  float2 positionF2 = make_float2(0.0f, 0.0f);
  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;
      float distanceToPoint = voronoi_distance_2d(pointPosition, localPosition, metric, exponent);
      if (distanceToPoint < distanceF1) {
        distanceF2 = distanceF1;
        distanceF1 = distanceToPoint;
        offsetF2 = offsetF1;
        offsetF1 = cellOffset;
        positionF2 = positionF1;
        positionF1 = pointPosition;
      }
      else if (distanceToPoint < distanceF2) {
        distanceF2 = distanceToPoint;
        offsetF2 = cellOffset;
        positionF2 = pointPosition;
      }
    }
  }
  *outDistance = distanceF2;
  *outColor = hash_float2_to_float3(cellPosition + offsetF2);
  *outPosition = positionF2 + cellPosition;
}

　基本的な実装はF1と同様であり、探索の際二番目に近い点を探るために少々変わっているという感じです。まず、F1用の変数とF2用の変数が用意されており、

  float distanceF1 = 8.0f;
  float distanceF2 = 8.0f;
  float2 offsetF1 = make_float2(0.0f, 0.0f);
  float2 positionF1 = make_float2(0.0f, 0.0f);
  float2 offsetF2 = make_float2(0.0f, 0.0f);
  float2 positionF2 = make_float2(0.0f, 0.0f);

探索の際、F1が発見された場合と、F2が発見された場合で更新を行う形になっています。F1が発見された場合にはF1だった特徴点はF2になるため、F1の情報をF2に移し、F1を更新します。一方、F1ではなかったもののF2が発見された場合にはF2のみの更新を行います。

      if (distanceToPoint < distanceF1) {
        distanceF2 = distanceF1;
        distanceF1 = distanceToPoint;
        offsetF2 = offsetF1;
        offsetF1 = cellOffset;
        positionF2 = positionF1;
        positionF1 = pointPosition;
      }
      else if (distanceToPoint < distanceF2) {
        distanceF2 = distanceToPoint;
        offsetF2 = cellOffset;
        positionF2 = pointPosition;
      }

　以上の実装によりF2は探索することができ、出力としてはF2の情報を用いてF1と同様にカラーや座標を出力する形となっています。

Smooth F1の実装

　Smooth F1とはF1の結果をSmoothnessというパラメータで滑らかな感じにぼかすものです。この実装は実はドキュメントができるぐらいには背景の技術がかなり複雑ですが、実装そのものは単純なためこの記事では簡単に流れだけ説明しようと思います。詳しい原理の説明は別の記事にて書く予定です。

　実装は以下のようになっています。

ccl_device void voronoi_smooth_f1_2d(float2 coord,
                                     float smoothness,
                                     float exponent,
                                     float randomness,
                                     NodeVoronoiDistanceMetric metric,
                                     float *outDistance,
                                     float3 *outColor,
                                     float2 *outPosition)
{
  float2 cellPosition = floor(coord);
  float2 localPosition = coord - cellPosition;

  float smoothDistance = 8.0f;
  float3 smoothColor = make_float3(0.0f, 0.0f, 0.0f);
  float2 smoothPosition = make_float2(0.0f, 0.0f);
  for (int j = -2; j <= 2; j++) {
    for (int i = -2; i <= 2; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;
      float distanceToPoint = voronoi_distance_2d(pointPosition, localPosition, metric, exponent);
      float h = smoothstep(
          0.0f, 1.0f, 0.5f + 0.5f * (smoothDistance - distanceToPoint) / smoothness);
      float correctionFactor = smoothness * h * (1.0f - h);
      smoothDistance = mix(smoothDistance, distanceToPoint, h) - correctionFactor;
      correctionFactor /= 1.0f + 3.0f * smoothness;
      float3 cellColor = hash_float2_to_float3(cellPosition + cellOffset);
      smoothColor = mix(smoothColor, cellColor, h) - correctionFactor;
      smoothPosition = mix(smoothPosition, pointPosition, h) - correctionFactor;
    }
  }
  *outDistance = smoothDistance;
  *outColor = smoothColor;
  *outPosition = cellPosition + smoothPosition;
}

　F1の実装と比べるとかなり異なる部分が多く、順々に説明していきます。まず、最初に最終的な出力用の変数として

  float smoothDistance = 8.0f;
  float3 smoothColor = make_float3(0.0f, 0.0f, 0.0f);
  float2 smoothPosition = make_float2(0.0f, 0.0f);

というのを用意しています。SmoothF1ではぼかしのような絵を出すため、F1だけではなく他の特徴点の距離、色、座標などを用いて補間を行います。従って、カラーなどの出力を最後ではなく探索中において逐次計算していくこととなり、そのための変数を用意しているわけです。

　次に探索部分ですが、今までとは違い5×5の更に広い範囲で探索を行います。これはsmoothnessが高くなると他の特徴点の影響が強くなることから１つ跨いだセルからの影響も考慮しなくてはならないことになるためです。実際にこの探索範囲を3×3にするとその分の影響がすっぱり消されてしまうため途切れる部分が出てきてしまいす。実際に3×3で済ませた画像を見るとわかりづらいですがうっすら切れ目のようなものが見えており、5×5でないとスムーズな画像が出ないことがわかります。

  for (int j = -2; j <= 2; j++) {
    for (int i = -2; i <= 2; i++) {

f:id:kinakomoti321:20220224185535p:plain — 3×3の場合と5×5の場合

f:id:kinakomoti321:20220224185539p:plain — 3×3の場合と5×5の場合

　そして肝心のDistanceなどの計算ですが、Smoothnessを用いてまずhというパラメータを生成しています。

      float h = smoothstep(
          0.0f, 1.0f, 0.5f + 0.5f * (smoothDistance - distanceToPoint) / smoothness);

このhは以後の補間に関するパラメータとなります。smoothstepを使用している理由は単にSmoothnessを動かしたときhの動きも滑らかになるようにしているだけです。本質的には

     float h = clamp(0.5f + 0.5f * (smoothDistance - distanceToPoint) / smoothness,0.0f,1.0f);

とみても大丈夫なはずです。

　hの計算後、さらに補間用のファクターとして以下のようなものを計算します。

      float correctionFactor = smoothness * h * (1.0f - h);

これを行った後に最終的にsmoothDistanceには

      smoothDistance = mix(smoothDistance, distanceToPoint, h) - correctionFactor;

という形の計算が行われます。この式は実はmin関数とほとんど同じ働きをしており、min関数を滑らかにしたような出力を行います。そのため、ここで言うことはCellerNoiseで言う

         //最小距離の更新
            m_dist = min(m_dist,dist);

に相当する処理になります。ただし、ここで滑らかなmin関数を使用することでいわゆるぼかしのような作用が生まれているわけです。

　カラー、座標についても同様なことを行いますが、補間のファクター部分に対して

      correctionFactor /= 1.0f + 3.0f * smoothness;

f:id:kinakomoti321:20220224184637p:plain — ファクターの処理を行った場合と行わなかった場合

f:id:kinakomoti321:20220224184645p:plain — ファクターの処理を行った場合と行わなかった場合

という計算を行った上で処理を行っています。おそらくこれは距離と同様の補間をしてしまうとぼかしが強すぎてしまうために、このような調整がなされているのだと考えられます。補間そのものは距離と同様に行います。

      smoothColor = mix(smoothColor, cellColor, h) - correctionFactor;
      smoothPosition = mix(smoothPosition, pointPosition, h) - correctionFactor;

以上のような形でSmoothF1の実装がなされています。

Distance to Edgeの実装

　Voronoi図において、それぞれ領域が区分されていますがこの際の境界線となる部分をエッジと呼びます。Distance to Edgeはそのエッジまでの最小距離について出力するものです。

　エッジはF1が切り替わる境界線であり、ある特徴点２つから同じ距離にあるものです。つまり、エッジは特徴点の二等分線によって構築されており、実際に出るエッジとはこの二等分線らの中から最も近い部分だけ取り出したものとなります。

f:id:kinakomoti321:20220224020834p:plain — エッジと特徴点の関係

　ある特徴点 $a,b$ が作る二等分線から座標 $x$ までの距離 $distance$ は幾何的な関係性から、以下のような関係式にあります。なお $x$ は $a$ 側にあり、 $\theta$ は $b-a$ と $x - (a+b) / 2$ が成す角として定義しています。

$\displaystyle{distance = |x - \frac{a+b}{2}| \cos{\theta} = (x - \frac{a+b}{2}) \cdot (\frac{b - a}{|b - a|}) }$

f:id:kinakomoti321:20220223234720p:plain — 二等分線までの距離の計算

　これによって $a,b$ が作るエッジまでの距離を計算できるため、 $a$ とその周辺の特徴点が作るエッジまでのそれぞれの距離を計算することでエッジまでの最短距離というのが求めることができるというわけです。

　ただし、これはあくまで距離計量がユークリッドの場合の話をしています。他の距離計量においてはそもそも二等分線が直線でなかったりとセルのエッジがかなり複雑となってしまうため、ユークリッド以外の距離計量でのDistance to Edgeは実装されていません。なので、Voronoi TextureノードではDistance to Edgeを選択すると距離計量の欄がなく、ユークリッドのみに固定されることとなります。

　では実装の方を見ていきます。実際の実装は以下のようになっています。

ccl_device void voronoi_distance_to_edge_2d(float2 coord, float randomness, float *outDistance)
{
  float2 cellPosition = floor(coord);
  float2 localPosition = coord - cellPosition;

  float2 vectorToClosest = make_float2(0.0f, 0.0f);
  float minDistance = 8.0f;
  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 vectorToPoint = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness -
                             localPosition;
      float distanceToPoint = dot(vectorToPoint, vectorToPoint);
      if (distanceToPoint < minDistance) {
        minDistance = distanceToPoint;
        vectorToClosest = vectorToPoint;
      }
    }
  }

  minDistance = 8.0f;
  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 vectorToPoint = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness -
                             localPosition;
      float2 perpendicularToEdge = vectorToPoint - vectorToClosest;
      if (dot(perpendicularToEdge, perpendicularToEdge) > 0.0001f) {
        float distanceToEdge = dot((vectorToClosest + vectorToPoint) / 2.0f,
                                   normalize(perpendicularToEdge));
        minDistance = min(minDistance, distanceToEdge);
      }
    }
  }
  *outDistance = minDistance;
}

　この実装を見ると２回ほどfor文があり、探索が２度行われています。それぞれどうゆうものかを順にみていきます。

　まず、最初のfor文ですがこれは単にF1の探索です。この際、距離の計算はdotで行っていますがこれは単にF1までの距離そのものはあとで使わないので重い処理であるsqrtを避けただけと考えられます。また、ここでは特徴点の座標はlocalPositionからの相対座標であることに注意してください。

  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 vectorToPoint = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness -
                             localPosition;
      float distanceToPoint = dot(vectorToPoint, vectorToPoint);
      if (distanceToPoint < minDistance) {
        minDistance = distanceToPoint;
        vectorToClosest = vectorToPoint;
      }
    }
  }

　VectorToClosestにはF1の座標が入っており、これは次の探索に使用されます。次の探索は最も近いエッジの計算となります。minDistanceをまたリセットし、再び同じ探索を行います。今回の探索で行うことはそれぞれの特徴点とF1が作るエッジの最短距離を求めることです。

  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 vectorToPoint = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness -
                             localPosition;
      float2 perpendicularToEdge = vectorToPoint - vectorToClosest;
      if (dot(perpendicularToEdge, perpendicularToEdge) > 0.0001f) {
        float distanceToEdge = dot((vectorToClosest + vectorToPoint) / 2.0f,
                                   normalize(perpendicularToEdge));
        minDistance = min(minDistance, distanceToEdge);
      }
    }
  }

　まず、先ほどと同様に特徴点の座標を生成します。

      float2 cellOffset = make_float2(i, j);
      float2 vectorToPoint = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness -
                             localPosition;

　この次にやることはさっき求めたF1と今の特徴点が作るエッジ（二等分線）からの距離を計算し、最短距離を更新することです。実装ではif文で分岐が行われていますが、これはF1とF1自身の計算を無視するためです。if文の中はF1以外の特徴点の場合に処理されるもので、やってることはさっきの式をそのまま行っており、最短距離の更新を行っています。最終的にこれを出力することでエッジまでの最短距離を求めることができるわけです。

      if (dot(perpendicularToEdge, perpendicularToEdge) > 0.0001f) {
        float distanceToEdge = dot((vectorToClosest + vectorToPoint) / 2.0f,
                                   normalize(perpendicularToEdge));
        minDistance = min(minDistance, distanceToEdge);
      }

　以上がDistance to Edgeの処理の流れとなります。　

N-Sphere Radiusの実装

　N-SphereRadiusはボロノイセルに対して内接する特徴点を中心とするn次元球（二次元であれば円）の半径を出力するものです。その半径はF1からエッジの最短距離に相当します。実際の実装では、F1とF1に最も近い特徴点を探索し、その２つの点の距離の半分を求めることでこれを求めます。なので実装そのものはかなり単純です。

ccl_device void voronoi_n_sphere_radius_2d(float2 coord, float randomness, float *outRadius)
{
  float2 cellPosition = floor(coord);
  float2 localPosition = coord - cellPosition;

  float2 closestPoint = make_float2(0.0f, 0.0f);
  float2 closestPointOffset = make_float2(0.0f, 0.0f);
  float minDistance = 8.0f;
  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;
      float distanceToPoint = distance(pointPosition, localPosition);
      if (distanceToPoint < minDistance) {
        minDistance = distanceToPoint;
        closestPoint = pointPosition;
        closestPointOffset = cellOffset;
      }
    }
  }

  minDistance = 8.0f;
  float2 closestPointToClosestPoint = make_float2(0.0f, 0.0f);
  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      if (i == 0 && j == 0) {
        continue;
      }
      float2 cellOffset = make_float2(i, j) + closestPointOffset;
      float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;
      float distanceToPoint = distance(closestPoint, pointPosition);
      if (distanceToPoint < minDistance) {
        minDistance = distanceToPoint;
        closestPointToClosestPoint = pointPosition;
      }
    }
  }
  *outRadius = distance(closestPointToClosestPoint, closestPoint) / 2.0f;
}

　こちらも探索が二度行われているため、それぞれ見ていきます。まず、最初の探索は例のごとくF1を探索しており、F1の情報を得ておきます。

  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      float2 cellOffset = make_float2(i, j);
      float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;
      float distanceToPoint = distance(pointPosition, localPosition);
      if (distanceToPoint < minDistance) {
        minDistance = distanceToPoint;
        closestPoint = pointPosition;
        closestPointOffset = cellOffset;
      }
    }
  }

　その次の探索ではF1に最も近い特徴点を調べるために”F1のセルを中心に”３×３の探索を行います。どのようにしてやっているかというと単にcellOffsetの値にF1のセルのオフセットを加えることで実装しています。また、F1のセルは無視するため i = 0,j = 0の際に処理を行わないようになっています。

  for (int j = -1; j <= 1; j++) {
    for (int i = -1; i <= 1; i++) {
      if (i == 0 && j == 0) {
        continue;
      }
      float2 cellOffset = make_float2(i, j) + closestPointOffset;
      float2 pointPosition = cellOffset +
                             hash_float2_to_float2(cellPosition + cellOffset) * randomness;
      float distanceToPoint = distance(closestPoint, pointPosition);
      if (distanceToPoint < minDistance) {
        minDistance = distanceToPoint;
        closestPointToClosestPoint = pointPosition;
      }
    }
  }

　この探索によりF1とF1に最も近い特徴点の２つを見つけられたので、その２点間の距離を求めて半分の値を出力することで内接円の半径を求めることができます。

  *outRadius = distance(closestPointToClosestPoint, closestPoint) / 2.0f;

　以上がN-sphere Radiusの実装となります。

Distance Metricの実装

　上記の通り、エッジは特徴点同士の「垂直二等分線」で構成されていることはどんな距離計量でも共通することです。ユークリッド以外の距離計量で定義される「垂直二等分線」はユークリッドのように直線ではなく、折曲がった直線、はたまた曲線であったりします。そのため、距離計量によってVoronoi図が変化するというわけです。

　こうした距離計量は実装上どのように行われているか見ていきます。まず、各VoronoiTextureの処理では二点間の距離を求める関数としてdistanceToPointという距離の関数が別途に用意されており、２つの座標と共にmetricやexponentという変数を入れています。

 　float distanceToPoint = voronoi_distance_2d(pointPosition, localPosition, metric, exponent);

　この距離関数は以下のように定義されています。この実装を見るとmetricという引数に距離計量が指定されており、それに応じてif文で距離の計算を分けているという感じになっています。

ccl_device float voronoi_distance_2d(float2 a,
                                     float2 b,
                                     NodeVoronoiDistanceMetric metric,
                                     float exponent)
{
  if (metric == NODE_VORONOI_EUCLIDEAN) {
    return distance(a, b);
  }
  else if (metric == NODE_VORONOI_MANHATTAN) {
    return fabsf(a.x - b.x) + fabsf(a.y - b.y);
  }
  else if (metric == NODE_VORONOI_CHEBYCHEV) {
    return max(fabsf(a.x - b.x), fabsf(a.y - b.y));
  }
  else if (metric == NODE_VORONOI_MINKOWSKI) {
    return powf(powf(fabsf(a.x - b.x), exponent) + powf(fabsf(a.y - b.y), exponent),
                1.0f / exponent);
  }
  else {
    return 0.0f;
  }
}

　また、各次元それぞれで1d,2d,3d,4dでの距離関数が作られています。2~4dに関してはベクトルの次元が増えるだけなので特筆すべき点はありませんが、1dに関しては距離計量をどう取っても変わらないため、絶対値を取るだけで済ましています。

ccl_device float voronoi_distance_1d(float a,
                                     float b,
                                     NodeVoronoiDistanceMetric metric,
                                     float exponent)
{
  return fabsf(b - a);
}

ユークリッド距離

　 $n$ 次元空間における $\vec{a} = (a_1,a_2, ... , a_n)$ , $\vec{b} = (b_1,b_2,... , b_n)$ の２点間のユークリッド距離はよく知られているように以下のように定義されており、ユークリッド距離の等高線グラフを出すと円形の模様が出ます。

$\displaystyle{Euclidian = \sqrt{\sum_{n}^{i} |a_i - b_i|^2} }$

f:id:kinakomoti321:20220224020353p:plain — ユークリッド距離

　実装においては単にdistanceという関数で済まされていますが、中身は別の場所で

float distance(float a, float b)
{
  return abs(a - b);
}

float distance(vector2 a, vector2 b)
{
  return length(a - b);
}

float distance(vector4 a, vector4 b)
{
  return length(a - b);
}

というようにちゃんとユークリッド距離を返すようになっています。

マンハッタン距離

　マンハッタン距離は直角にしか動けない場合における最短距離を示します。よく言われるのが囲碁の目のようなマス目を上下左右だけで移動した際にかかる手数というイメージです。定義としては以下のように各成分の絶対値の総和を取るだけであり、その等高線グラフはユークリッド距離とは異なりひし形のような形状を示します。

$\displaystyle{Manhattan = \sum_{n}^{i} |a_i - b_i| }$

f:id:kinakomoti321:20220224014734p:plain — マンハッタン距離

　実装の方は以下のようになります。fabsfは単なるabs関数であるため、上の式をそのまま行っているだけなのがわかります。

    return fabsf(a.x - b.x) + fabsf(a.y - b.y);

チェビシェフ距離

　チェビシェフ距離はマンハッタンと似ていますが、少し違う点として横軸と縦軸の距離の内、長い方の軸の距離を距離として扱います。定義としては以下のようにMax関数を用いてそれぞれの成分の絶対値の内、最も大きい値を求めます。等高線のグラフにすると四角形のような模様が表れます。

$\displaystyle{Chebychev = \mathrm{max}(|a_1 - b_1|,|a_2 - b_2|,...,|a_n - b_n|)}$

f:id:kinakomoti321:20220224015627p:plain — チェビシェフ距離

　実装の方ではこの定義をそのまま実装している形となっています。

    return max(fabsf(a.x - b.x), fabsf(a.y - b.y));

ミンコフスキー距離

　ミンコフスキー距離は距離計量を一般化したものです。 $n$ 次元空間において、ミンコフスキー距離は以下のように定義されています。

$\displaystyle{Minkowski = (\sum_{n}^{i} |a_i - b_i|^p)^{1/p}}$

f:id:kinakomoti321:20220224005028p:plain — ミンコフスキー距離の等高線

　このようにミンコフスキー距離では次数 $p$ というパラメータを持ち、これによって距離計量の取り方が変わってきます。先述の通り、これは距離計量を一般化したものであり、上の等高線のように $p$ の値によって他の距離計量に一致することとなります。p = 1の際はマンハッタン距離、p = 2の時はユークリッド距離、さらに大きくしていくとチェビシェフ距離に近づいていきます。

　実装ではexponentという変数はミンコフスキー距離の次数 $p$ であり、以上の定義をそのまま計算している形となります。

    powf(powf(fabsf(a.x - b.x), exponent) + powf(fabsf(a.y - b.y), exponent),
                1.0f / exponent);

次元について

　基本的に次元が変わったところで処理はそこまで変わることはありません。一例として３次元のF1の関数についてみてみますと

ccl_device void voronoi_f1_3d(float3 coord,
                              float exponent,
                              float randomness,
                              NodeVoronoiDistanceMetric metric,
                              float *outDistance,
                              float3 *outColor,
                              float3 *outPosition)
{
  float3 cellPosition = floor(coord);
  float3 localPosition = coord - cellPosition;

  float minDistance = 8.0f;
  float3 targetOffset = make_float3(0.0f, 0.0f, 0.0f);
  float3 targetPosition = make_float3(0.0f, 0.0f, 0.0f);
  for (int k = -1; k <= 1; k++) {
    for (int j = -1; j <= 1; j++) {
      for (int i = -1; i <= 1; i++) {
        float3 cellOffset = make_float3(i, j, k);
        float3 pointPosition = cellOffset +
                               hash_float3_to_float3(cellPosition + cellOffset) * randomness;
        float distanceToPoint = voronoi_distance_3d(
            pointPosition, localPosition, metric, exponent);
        if (distanceToPoint < minDistance) {
          targetOffset = cellOffset;
          minDistance = distanceToPoint;
          targetPosition = pointPosition;
        }
      }
    }
  }
  *outDistance = minDistance;
  *outColor = hash_float3_to_float3(cellPosition + targetOffset);
  *outPosition = targetPosition + cellPosition;
}

２次元のF1との大きな差異は探索のfor文の数ぐらいとなっています。

  for (int k = -1; k <= 1; k++) {
    for (int j = -1; j <= 1; j++) {
      for (int i = -1; i <= 1; i++) {

基本的にn次元のボロノイはこのように探索をn重for文で行うだけなので、３次元、４次元にしたければfor文を追加して、それに伴い座標とかのベクトルの次元も異なるのでそこさえ合わせれば次元を変更することができます。

おわりに

　VoronoiTextureの中身ってどうなってるのかなと疑問に思い調べたので、勉強したことをまとめるためにこの記事を書きました。実装を見てみると意外とわかりやすく、知らなかった手法とかを幅広く知ることができてとてもいい勉強になりました。Blenderは広く使われているだけあって結構幅広い技術が使われていて、尚且つ簡潔なコードになっているのでBlenderの実装を元にノイズとかを勉強していくのもありなんじゃないかと感じました。

　勉強がてら今回の記事で書いたVoronoiTextureをGLSLでひとまとめに書いたものをGitHubに上げましたので良ければご覧ください。

github.com

その他Voronoi（おまけ）

　Blenderには実装されてないその他Voronoi図の表現について書いていきます

Crackle

　BlenderのVoronoiTextureは2.8から2.9への移行で結構変わっていたようで、昔はF1とかの他に「Crackle」というFutureOutputがあったらしいです。これは簡単なDistanceToEdgeであり、ユークリッド距離意外にも適用することができました。なのでManhattan距離とかでのエッジを作りたい場合によく使われていたそうです。　

f:id:kinakomoti321:20220321041140p:plain — Crackle

　これの実装はかなり簡単でF2のDistanceをF1のDistanceで引き、係数をかけるだけという処理で作ることができます。

float voronoi_Crackle_2d(vec2 uv,int MetricMode,float exponent,float randomness,float scale){
    float F1;
    float F2;
    vec3 col;
    vec2 pos;
    voronoi_f1_2d(uv,exponent,randomness,MetricMode,F1,col,pos);
    voronoi_f2_2d(uv,exponent,randomness,MetricMode,F2,col,pos);

    return clamp((F2 - F1)*scale,0.0,1.0);
}

　ちなみにBlenderのシェーダーノードで再現するとこんな感じになります。

f:id:kinakomoti321:20220321041159p:plain — Crackleをシェーダーノードで再現

　Crackleはユークリッド距離以外にも使えるDistanceToEdgeみたいな感じですが、別段エッジまでの距離をちゃんと出すものでもなんでもなく、それっぽい感じの出力をするものというだけです。所々、色が変化なかったりしきれいなエッジを作ることができなかったりします。おそらくですが、こうした理由から2.9の移行の際削除され、より正確なエッジを出せるDistanceToEdgeを実装し置き換えたのだと思われます。（割とほかの距離計量でのエッジも欲しい時は結構あるのでなくさなくてもよかったのでは感がある）

Weightあり Voronoi

　WeightありVoronoiは各セル $c$ に対して固有のウェイト $w_c$ を考えてあげて、距離の計算の際にこのウェイトを引くことで作られるVoronoi図のことです。これを生成するには単に距離の計算時にウェイトを引いてあげるだけでできるので実装上ではこんな感じに距離を計算する部分を書き換えてあげればできます。

float distanceToPoint = voronoi_distance_2d(pointPosition,localPosition,MetricMode,expornent) - Hash_2D_to_1D(cellPosition + cellOffset) * weight;

f:id:kinakomoti321:20220321043101p:plain — WeightありVoronoi

　正直使い道はあんまり思いつかないです。

BlenderのVoronoiの実装を見る

初めに