符号なし距離関数で作る色々なボロノイ

# ボロノイとは
プロシージャルな石畳や細胞のパターンの生成、破壊シミュレーション時のオブジェクトの分割に良く使われるノイズのことです。詳しくは[The Book of Shaders セルラーノイズ](https://thebookofshaders.com/12/?lan=jp)を参照ください

# 通常のボロノイ
## ユークリッド(p=2)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/96810d2e-208e-872b-68bc-19cce4b9df34.png)

```glsl
    d = length(r);
    d = dot2(r);
```

通常のボロノイはこんな感じですが、距離関数を変えることで色々なボロノイを作ることができます。尚、使用するのは、符号なしのディスタンスフィールドの方です。SDFは使えません。

[2次元ディスタンスフィールドの一覧はこちら](https://qiita.com/7CIT/items/fe33b9b341b9918b6c3d)

# 数学的な距離関数を使う

私たちが、普段2点間の距離を測る時に使用しているのは、三平方の定理で導かれる『ユークリッド距離』というものですが、『マンハッタン距離』、『チェビシェフ距離』、『ミンコフスキー距離』など2点間の距離の測り方にも様々な種類があります。

通常のボロノイは『ユークリッド距離』を使用していますが、この距離の測り方を『マンハッタン距離』、『チェビシェフ距離』に変えることで回路ぽいボロノイを作成できます。これはかなり使われていますね。

## マンハッタン(p=1)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/8f2a64d9-acd3-b1bf-57a8-7b4e8ac19f02.png)

```glsl
float sumabs(vec2 p){return abs(p.x)+abs(p.y);}
d = sumabs(r);
```

## チェビシェフ(p=∞)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/7c8ffc2f-8ddc-ed40-9861-e22edc541f3f.png)

```glsl
float maxabs(vec2 p){return max(abs(p.x), abs(p.y));}
d = maxabs(r);
```

## ミンコフスキー

さて、他にも数学的距離関数がないのかと言うと、ミンコフスキー距離というものがありまして、これは、『ユークリッド距離』、『マンハッタン距離』、『チェビシェフ距離』を[LPノルム](https://mathtrain.jp/lpnorm)という一つの式(一般式)で表すことができます。

以下、2次元(2D)の場合

* ミンコフスキー距離(p)

```math
\sqrt[p]{|x^p|+|y^p|}
```

* マンハッタン距離(p=1)

```math
\sqrt[1]{|x^1|+|y^1|}=|x|+|y|
```

* ユークリッド距離(p=2)

```math
\sqrt[2]{|x^2|+|y^2|}=\sqrt{x^2+y^2}
```

* チェビシェフ距離(p=∞)

```math
\sqrt[1]{|x^1|+|y^1|}=|x|+|y|=\sqrt[∞]{|x^∞|+|y^∞|}
```

各距離は、次元数pで表すことができ、p=1がマンハッタン、p=2がユークリッド、p=∞がチェビシェフ距離になります。なので、マンハッタン=L1距離、ユークリッド=L2距離と呼ばれていたりします。

|ノルム(p)|数学的距離名|
|---|---|---|
|1  |マンハッタン距離(L1距離)|
|2  |ユークリッド距離(L2距離)|
|∞  |チェビシェフ距離|
|p  |ミンコフスキー距離|

もちろん、このpには-∞や0.5なども任意に指定できるので今度は任意の次元pを指定したミンコフスキー距離を使ってみましょう。以下は面白い図形ができるp=-∞やp=0.5の場合を紹介しています。

```glsl
float lpnorm(vec2 p, float n){vec2 t=pow(abs(p),vec2(n));return pow(t.x+t.y,1./n);}
```

ところで、『チェビシェフ距離』はp=∞ですが

```glsl
d=lpnorm(r,∞);
```

とかになっていないのにお気づきでしょうか。

無限に近い数値を入れたり

```glsl
d=lpnorm(r,1e10);
```

することで近似値を得ることが出来るのですが計算負荷的にもお勧めできません。shaderでは実はもっと簡単に実現できます。max関数を利用すればいいのです。∞がmax関数なら-∞はmin関数ですね。

ディープラーニングの発展に多大な貢献をもたらした、数学的には正しくないReLU関数(=単なるmax関数ですね)のようにプログラミング言語の関数は、新しい数学の世界を開く、shaderは新しい数学的発見を見つけることができる楽しい世界です。前振り

![voronoi.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/de5dc197-82f0-c501-4bf1-5a8e10727cb7.gif)

### セル(p=3)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/d78a8c9c-bab5-132e-4762-fae588ac546a.png)

pが2以上は植物の細胞のような感じになる。pが大きくなればなるほど角ばって行きチェビシェフに近づく

### 十字(n=-∞)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/f180132a-7158-1317-eb21-87379e83a9e0.png)

```glsl
    d = minabs(r);
```

```glsl
d=lpnorm(r,-1e10);
```
でも良いですがmin使えばいいので

### アステロイド(n=0.5)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/9ff4a9f9-e926-fe6b-b304-ff3b9f1188af.png)

```glsl
    d = lpnorm(r,.5);
```

# 数学的距離関数以外

さて、本題のshaderで図形を描く時に使われる元の関数も距離関数です。
もちろん、数学的な距離関数も使われていて、以下の図形を描くのに使われています。

|ノルム(n)|図形|数学的距離名|
|---|---|---|
|-∞ |十字        |- |
|.5 |アステロイド |-  ||
|1  |正方形(45度)|マンハッタン距離(L1距離)|
|2  |円          |ユークリッド距離(L2距離)|
|∞  |正方形      |チェビシェフ距離|
|p  |-          |ミンコフスキー距離|

しかし、これだけでは図形を描くには足りないので、角丸長方形や、三角形、六角形、星型、などの距離関数が編み出されました。

当然のことながら、これら符号なし距離関数であれば、なんでも使えるので、オリジナルのボロノイを作ってみましょう。以下はその中で面白かったボロノイを紹介しています。

ただ、ここら辺の距離関数を数式で表すのはかなり困難ですね。

## 三角形
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/3ea0de52-dd46-aa90-8d5c-104626a74a48.png)

```glsl
    d = lPoly(r,3.);
```

## 星型

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/135381/0748bf5e-e6d3-cf39-dee0-7722e341cf7c.png)

```glsl
    d = lStar(r,5.);
```

# サンプルプログラム

```glsl
precision mediump float;
uniform vec2  m;       // mouse
uniform float t;       // time
uniform vec2  r;       // resolution
uniform sampler2D smp; // prev scene

const float PI = acos(-1.);
const float TAU = PI * 2.;

mat2 rot(float a){float c=cos(a),s=sin(a);return mat2(c,-s,s,c);}

// signed distance
float sd(float d,float r){return r-d;} 
float sd(float d){return 1.-d;} 

vec2 hash2( vec2 p ){
    return fract(sin(vec2(dot(p,vec2(127.1,311.7)),dot(p,vec2(269.5,183.3))))*43758.5453);
}

float lPoly(vec2 p,float n){
  float a = atan(p.x,p.y)+PI;
  float r = TAU/n;
  return cos(floor(.5+a/r)*r-a)*length(p)/cos(r*.5);
}

float lStar(vec2 p,float n){
  return min(lPoly(p,n*.5),lPoly(mod(n,2.)!=0.?vec2(-p.x,p.y):p*rot(TAU/n),n*.5));
}

float dot2(vec2 p){return dot(p,p);}
float sumabs(vec2 p){return abs(p.x)+abs(p.y);}
float maxabs(vec2 p){return max(abs(p.x), abs(p.y));}
float minabs(vec2 p){return min(abs(p.x), abs(p.y));}
float lpnorm(vec2 p, float n){vec2 t=pow(abs(p),vec2(n));return pow(t.x+t.y,1./n);}

float voronoi(vec2 x){
  vec2 n = floor(x);
  vec2 f = fract(x);
  vec2 res = vec2(8.0);
  for(int j = -1; j <= 1; j++)
  for(int i = -1; i <= 1; i++){
    vec2 g = vec2(i, j);
    vec2 r = g - f + hash2(n + g);
    float d = 0.;
    //d = length(r);
    //d = dot2(r);
    //d = sumabs(r);
    //d = sumabs(r);
    //d = maxabs(r);
    //d = minabs(r);
    //d = lpnorm(r,.5);
    //d = lpnorm(r,sin(t*.5)*4.);
    //d = lPoly(r,3.);
    d = lStar(r,5.);

    if(d < res.x){
      res.y = res.x;
      res.x = d;
    }else if(d < res.y){
      res.y = d;
    }
  }
  vec2 c = sqrt(res);
  return c.y - c.x;
}

void main(void){
    vec2 p = (gl_FragCoord.xy * 2.0 - r) / min(r.x, r.y);
    float v = voronoi(4.*p);
    float l = 1.-smoothstep(.02,.04,v);
    //l=v;
    gl_FragColor = vec4(vec3(l), 1.0);
}
```

* [サンプルプログラムを実行](http://jp.wgld.org/js4kintro/editor/#precision%20mediump%20float%3B%0Auniform%20vec2%20%20m%3B%20%20%20%20%20%20%20%2F%2F%20mouse%0Auniform%20float%20t%3B%20%20%20%20%20%20%20%2F%2F%20time%0Auniform%20vec2%20%20r%3B%20%20%20%20%20%20%20%2F%2F%20resolution%0Auniform%20sampler2D%20smp%3B%20%2F%2F%20prev%20scene%0A%0Aconst%20float%20PI%20%3D%20acos(-1.)%3B%0Aconst%20float%20TAU%20%3D%20PI%20*%202.%3B%0A%0Amat2%20rot(float%20a)%7Bfloat%20c%3Dcos(a)%2Cs%3Dsin(a)%3Breturn%20mat2(c%2C-s%2Cs%2Cc)%3B%7D%0A%0A%2F%2F%20signed%20distance%0Afloat%20sd(float%20d%2Cfloat%20r)%7Breturn%20r-d%3B%7D%20%0Afloat%20sd(float%20d)%7Breturn%201.-d%3B%7D%20%0A%0Avec2%20hash2(%20vec2%20p%20)%7B%0A%20%20%20%20return%20fract(sin(vec2(dot(p%2Cvec2(127.1%2C311.7))%2Cdot(p%2Cvec2(269.5%2C183.3))))*43758.5453)%3B%0A%7D%0A%0Afloat%20lPoly(vec2%20p%2Cfloat%20n)%7B%0A%20%20float%20a%20%3D%20atan(p.x%2Cp.y)%2BPI%3B%0A%20%20float%20r%20%3D%20TAU%2Fn%3B%0A%20%20return%20cos(floor(.5%2Ba%2Fr)*r-a)*length(p)%2Fcos(r*.5)%3B%0A%7D%0A%0Afloat%20lStar(vec2%20p%2Cfloat%20n)%7B%0A%20%20return%20min(lPoly(p%2Cn*.5)%2ClPoly(mod(n%2C2.)!%3D0.%3Fvec2(-p.x%2Cp.y)%3Ap*rot(TAU%2Fn)%2Cn*.5))%3B%0A%7D%0A%0Afloat%20dot2(vec2%20p)%7Breturn%20dot(p%2Cp)%3B%7D%0Afloat%20sumabs(vec2%20p)%7Breturn%20abs(p.x)%2Babs(p.y)%3B%7D%0Afloat%20maxabs(vec2%20p)%7Breturn%20max(abs(p.x)%2C%20abs(p.y))%3B%7D%0Afloat%20minabs(vec2%20p)%7Breturn%20min(abs(p.x)%2C%20abs(p.y))%3B%7D%0Afloat%20lpnorm(vec2%20p%2C%20float%20n)%7Bvec2%20t%3Dpow(abs(p)%2Cvec2(n))%3Breturn%20pow(t.x%2Bt.y%2C1.%2Fn)%3B%7D%0A%0Afloat%20voronoi(vec2%20x)%7B%0A%20%20vec2%20n%20%3D%20floor(x)%3B%0A%20%20vec2%20f%20%3D%20fract(x)%3B%0A%20%20vec2%20res%20%3D%20vec2(8.0)%3B%0A%20%20for(int%20j%20%3D%20-1%3B%20j%20%3C%3D%201%3B%20j%2B%2B)%0A%20%20for(int%20i%20%3D%20-1%3B%20i%20%3C%3D%201%3B%20i%2B%2B)%7B%0A%20%20%20%20vec2%20g%20%3D%20vec2(i%2C%20j)%3B%0A%20%20%20%20vec2%20r%20%3D%20g%20-%20f%20%2B%20hash2(n%20%2B%20g)%3B%0A%20%20%20%20float%20d%20%3D%200.%3B%0A%20%20%20%20%2F%2Fd%20%3D%20length(r)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20dot2(r)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20sumabs(r)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20sumabs(r)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20maxabs(r)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20minabs(r)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20lpnorm(r%2C.5)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20lpnorm(r%2Csin(t*.5)*4.)%3B%0A%20%20%20%20%2F%2Fd%20%3D%20lPoly(r%2C3.)%3B%0A%20%20%20%20d%20%3D%20lStar(r%2C5.)%3B%0A%0A%20%20%20%20if(d%20%3C%20res.x)%7B%0A%20%20%20%20%20%20res.y%20%3D%20res.x%3B%0A%20%20%20%20%20%20res.x%20%3D%20d%3B%0A%20%20%20%20%7Delse%20if(d%20%3C%20res.y)%7B%0A%20%20%20%20%20%20res.y%20%3D%20d%3B%0A%20%20%20%20%7D%0A%20%20%7D%0A%20%20vec2%20c%20%3D%20sqrt(res)%3B%0A%20%20return%20c.y%20-%20c.x%3B%0A%7D%0A%0Avoid%20main(void)%7B%0A%20%20%20%20vec2%20p%20%3D%20(gl_FragCoord.xy%20*%202.0%20-%20r)%20%2F%20min(r.x%2C%20r.y)%3B%0A%20%20%20%20float%20v%20%3D%20voronoi(4.*p)%3B%0A%20%20%20%20float%20l%20%3D%201.-smoothstep(.02%2C.04%2Cv)%3B%0A%20%20%20%20%2F%2Fl%3Dv%3B%0A%20%20%20%20gl_FragColor%20%3D%20vec4(vec3(l)%2C%201.0)%3B%0A%7D)

# 参考
* [2次元ディスタンスフィールドの一覧はこちら](https://qiita.com/7CIT/items/fe33b9b341b9918b6c3d)
* [iq先生のボロノイ解説](https://www.iquilezles.org/www/articles/voronoilines/voronoilines.htm) ※ iq先生のボロノイは完璧なボロノイを実現しますが、距離関数の変更には対応していないので注意
* [ボロノイ！シェーダーでボロノイをあそぼう！！～理論とかはいいからとにかく使い方を知りたい～](https://scrapbox.io/sayachang/%E3%83%9C%E3%83%AD%E3%83%8E%E3%82%A4%EF%BC%81%E3%82%B7%E3%82%A7%E3%83%BC%E3%83%80%E3%83%BC%E3%81%A7%E3%83%9C%E3%83%AD%E3%83%8E%E3%82%A4%E3%82%92%E3%81%82%E3%81%9D%E3%81%BC%E3%81%86%EF%BC%81%EF%BC%81%EF%BD%9E%E7%90%86%E8%AB%96%E3%81%A8%E3%81%8B%E3%81%AF%E3%81%84%E3%81%84%E3%81%8B%E3%82%89%E3%81%A8%E3%81%AB%E3%81%8B%E3%81%8F%E4%BD%BF%E3%81%84%E6%96%B9%E3%82%92%E7%9F%A5%E3%82%8A%E3%81%9F%E3%81%84%EF%BD%9E)
* [さやちゃんぐbotスクラップス Voronoi](https://scrapbox.io/sayachang/Voronoi)

# 関連
* [ボロノイズ](http://www.iquilezles.org/www/articles/voronoise/voronoise.htm) iq先生がボロノイとノイズの一般式を開発しています。