[レイマーチング] IFSとPolar Modについて雑にまとめてみる

## 概要

ずっと気になっていたIFSとPolar Mod(別の呼び方としてfoldRotateとも）を、@kaneta1992さんが詳しく解説してくれていたので色々Shadertoyでいじってみつつ、コメントを追加してまとめてみます。

**以下の記事は必読レベルのいいまとめです！**

- [魔法使いになりたい人のためのシェーダーライブコーディング入門](https://qiita.com/kaneta1992/items/21149c78159bd27e0860)



## IFS

IFSは**Iterated function system**の略で、一言で言うと「**ループによって自身の拡縮したり回転したコピーを複製するフラクタル生成システム**」ということができると思います。

[Wikipedia](https://ja.wikipedia.org/wiki/%E5%8F%8D%E5%BE%A9%E9%96%A2%E6%95%B0%E7%B3%BB)から引用すると以下のように説明されています。

>反復関数系（はんぷくかんすうけい、英: Iterated function system、IFS）はフラクタルの一種であり、一般に2次元のフラクタルの描画や計算に用いられる。IFSフラクタルは自身のいくつかのコピーの和集合から成り、各コピーは関数によって変形されている（そのため「関数系」と呼ばれる）。典型例としてはシェルピンスキーのギャスケットがある。その関数は一般に収縮写像であり、点の集合がより近くなり、形がより小さくなる。従ってIFSフラクタルは、自身の縮小コピーを（場合によっては重ね合わせて）まとめたものであり、各部を詳細に見れば、その部分もそれ自身の縮小コピーから構成されていて、これが永遠に続く。このため、フラクタルとしての自己相似性が生じる。

*3次元で言うと「[メンガーのスポンジ](https://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%B3%E3%82%AC%E3%83%BC%E3%81%AE%E3%82%B9%E3%83%9D%E3%83%B3%E3%82%B8)」が有名かな？*

### IFSはループによって作られる

IFSを適用して絵を出すのには、ループによってレイを少しずつ変化させながら絵を描いていきます。
（ちなみにレイマーチ自体もループによって作られますが、ここで言うループは意味的には**ひとつの距離関数が**ループによってフラクタル的な構造を取るという意味です）

まずはコードを見てください。
`map`関数内で行っているのがIFSの処理です。


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

float sdBox(vec3 p, vec3 b)
{
    return length(max(abs(p) - b, 0.0));
}

float map(vec3 p)
{
    vec3 z = p;
    
    float scale = 2.0;
    float sum = scale;
    float d = 1e5;
        
    for (float i = 0.; i < 4; i++)
    {
        float td = sdBox(z, vec3(0.5)) / sum;
        
        z = abs(z) - vec3(0, 2.0, 0);
        d = min(td, d);
        
        z.xy *= rot(pi * 0.25);
        z *= scale;
        sum *= scale;
    }
    
    return d;
}
```

これを実行すると以下のような絵ができます。
![IFSサンプル](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/213d763c-b03e-9e15-84fc-495f926ea6f1.png)
- [IFS with scaling ](https://www.shadertoy.com/view/wlB3Wc)

よく見ると再帰的な構造になっているのが分かるかと思います。
Cubeの上に45度回転した半分のサイズCubeの上に45度回転した半分のサイズの上に45度回転した・・・というように。
（ちなみに「**上**」が45度回転した状態になっている点に注意）

今回はループ回数が4回だったので上記のような絵になっていますが、ループ回数を増やすとどんどん小さいCubeが積み重なっていきます。



### ループごとに変化を与える


今回の例ではスケールと平行移動、そして回転すべてをループのたびに変化させています。


#### 平行移動

平行移動させている箇所は以下の部分です。

```glsl:平行移動
z = abs(z) - vec3(0, 1.0, 0);
```

`abs`のあとのベクトルを引いている箇所がそれですね。

#### 回転

次に回転です。
回転を適用しているのは以下の部分です。

```glsl:回転
z.xy *= rot(pi * 0.25);
```

今回は分かりやすいように`xy`平面に限定して回転させています。
当然この箇所を複雑に回転させることでさらに面白い形状にすることも可能です。

#### スケール

次にスケールです。

```glsl:スケール
z *= scale;
sum *= scale;
```

スケールはループのたびに`z`の値を、指定したスケール分乗算していくことで達成しています。
つまりループの回数分どんどんスケールが小さく（or 大きく）なっていく、というわけですね。
そして同時に、`sum`という変数にもループのたびに`scale`を掛けています。


#### スケールに応じて「距離」を正規化する必要がある

平行移動と回転だけなら特に問題にならないのですが、スケールを扱う際には少し考慮することがあります。
スケールを扱う場合は**レイの進む距離**に変化が生じてしまうので対策しないとなりません。

*@kaneta1992 先生のIFSの解説を見て、「なるほどなるほど。**IFS完全に理解した**」と思いながらスケールを掛けてみたらうまく動かなかったので、やっぱり実際に手を動かさないとダメですね:yum:*

ポイントはループごとに掛けられていくスケールの合計で（例では変数`sum`）、距離関数の値を正規化しているところです。

```glsl:sumで正規化
float td = sdBox(z, vec3(0.5)) / sum;
```

距離関数で得られた値を**`sum`で除算**することで正規化しています。
これは、ループのたびにレイがスケールされていくため、前回のループで計算された距離関数の距離とは大きく異なった値が返ってきてしまうためです。
なので距離関数自体を現在のスケール率に合わせて正規化する必要がある、というわけなんですね。

ちなみにこれは[@kaneta](https://twitter.com/kanetaaaaa)先生にTwitterで教えていただきました！

<blockquote class="twitter-tweet" data-lang="en"><p lang="ja" dir="ltr">ifsの最中でscaleを欠けていくと、どんどん距離場のスケールが大きくなって、最短距離からかけ離れていきます。<br>イテレーションごとに距離を正規化させてやることできれいに表示されます！</p>&mdash; かねた (@kanetaaaaa) <a href="https://twitter.com/kanetaaaaa/status/1142689342948663296?ref_src=twsrc%5Etfw">June 23, 2019</a></blockquote>




## Polar Mod

さて次はPolar Mod（またはRotateFold）について。
Polar Modとは、レイマーチングでよく利用される`mod`を利用したリピートと考え方自体は同様のものです。
ただ、リピートの仕方が`xyz`軸に対して繰り返しが表現される`mod`に対して、**極座標**的な感じで繰り返しを行う点が異なります。

ちなみに通常のrepeatに関しては以下のような計算を行うものです。

```glsl:modによるrepeat
float map(vec p)
{
    p = mod(p, 3.0) - 1.5;
}
```

これは`map`関数に渡ってきたレイの位置（`p`）に対して`mod`を適用し、繰り返しが現れるようにしています。

----------------------------------


**余談**

ちなみにShadertoyでこれを関数化して定義していたら、こうしたほうが便利だよ！　って教えてもらったのでついでに記載しておきます。

```glsl:defineによるrepeat
#define repeat(p, span) mod(p, span) - (0.5 * span)
```

`#define`にする理由は、任意の平面（例えばXZ平面とか）に限定して繰り返しさせたい場合に、関数で定義してしまうと`vec2`や`vec3`を引数に取るそれぞれの関数を定義しないとなりませんが、`#define`なら事前にプリプロセッサによって置換されるためそうした細かい配慮が必要ない、というわけです。

例えば、以下のように異なる要素数の処理についてもひとつの`#define`のみで対応できます。

```glsl:repeatの例
// ひとつの#defineだけで以下に対応できる
p.x = repeat(p.x, 3.0);

p.xz = repeat(p.xz, 3.0);

p.xyz = repeat(p.xyz, 3.0);
```

便利ですね。

----------------------------------


### Polar Mod詳細

閑話休題。

さて、ではPolar Modの詳細について。
といってもこちらは[@kaneta1992先生の記事](https://qiita.com/kaneta1992/items/21149c78159bd27e0860)を読んだほうが分かりやすいと思うので、あくまでも自分の理解のためのメモです。

```glsl:Polar-Mod
// Rotate matrix for 2 dimension
mat2 rot(float a)
{
    float c = cos(a), s =sin(a);
    return mat2(c, s, -s, c);
}

// r is splited number.
vec2 pmod(vec2 p, float r)
{
    float a = atan(p.x, p.y) + pi / r;
    // "+ pi / r" means shortcut of "+ ((pi2 / r) * 0.5)".
    // so we want to get half angles of circle splitted by r.

    float n = pi2 / r;
    a = floor(a / n) * n;
    // floor(a / n) means calculating ID.
    
    return p * rot(-a);
}

// map function. This will render the scene.
float map(vec3 p)
{
    p.xy = pmod(p.xy, 5.);
    vec3 offs = vec3(0, 2, 0);
    return length(p - offs) - 1.0;
}
```

処理自体はとても短いですね。そしてメインとなる関数は当然`pmod`関数です。


### 関数解説

それぞれの関数を見ていきましょう。

#### 回転行列を返す関数rot

最初は回転行列を返す関数`rot`です。
これ自体は取り立てて特別なところはありません。
2次元での回転を想定している、という点だけ注意が必要です。

```glsl:回転
mat2 rot(float a)
{
    float c = cos(a), s =sin(a);
    return mat2(c, s, -s, c);
}
```

普通の2次元の回転行列を生成しています。
3次元的な回転でももしかしたらいいのかもしれないですが、今回のPolar Modは「**いずれかの面**」を回転させて分割する方法です。
特に今回の例では`xy`平面に対して回転を施しています。

平面を回転することを前提にしているため、2次元での回転行列を計算しているわけなんですね。
実際に利用されている部分を見てみると以下のようになっています。

```glsl:回転の適用
// xy平面の「2次元」だけに適用している
p.xy = pmod(p.xy, 5.0);
```


#### Polar Mod本体

次に見るのは**Polar Mod**自体の関数です。
以下の部分ですね。

```glsl:PolarMod
vec2 pmod(vec2 p, float r)
{
    float a = atan(p.x, p.y) + pi / r;
    // "+ pi / r" means shortcut of "+ ((pi2 / r) * 0.5)".
    // so we want to get half angles of circle splitted by r.

    float n = pi2 / r;
    a = floor(a / n) * n;
    // floor(a / n) means calculating ID.

    return p * rot(-a);
}
```

重要なポイントを雑に書くと、

- **360度を分割数で割りイチ区画分の角度を求め、その角度の範囲内で繰り返しになるようにレイの位置を回転させる**
- **角度の計算時にイチ区画分の角度の半分だけオフセットさせる**

です。

ちなみにオフセットを加える理由ですが、@kaneta1992さんの記事から引用させてもらうと以下のように説明されています。

>オフセットを加えていない場合、空間の切れ目が画面中央に存在していますよね？
>空間の切れ目にオブジェクトが移動してしまうので、見切れた形で複製されてしまいます。


-----------------------

**余談**

少し横道にそれますが、これは通常の`mod`を利用したリピートでも同様の問題があり、補正する理由もだいたい同じなのでここで解説しておきます。

簡単な例を示すと、UVの値に対して`mod`を適用してみると以下のような繰り返しになります。

![UVのmod](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/2bef2a95-b802-4561-b1bd-63b2e9cbf59e.png)

計算式としては以下です。

```glsl:UVのmod
void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 uv = (fragCoord * 2.0 - iResolution.xy) / iResolution.xy;

    uv = mod(uv, 1.0);

    vec3 col = vec3(uv, 0);

    fragColor = vec4(col, 1.0);
}
```

図で見ると分かるように、左右対称ではなく中央からばっつりと分かれている絵になっているのが分かるかと思います。
これは`floor`の仕様として**引数の値以下の最大の整数**を返すのが理由です。

プラス方向の場合は`floor(0.5) == 0`、`floor(1.2) == 1`となりますが、マイナス方向の場合は`floor(-0.1) == -1`となります。
与えられた引数**以下**の整数なので、`-0.1`より小さい`-1`になる、というわけですね。

そしてGLSLの`mod`の[ドキュメント](https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/mod.xhtml)には以下のように書かれています。

>mod returns the value of x modulo y. This is computed as x - y * floor(x/y).

これを使って計算すると、`x`が`-0.1`のとき（つまり上の例の図で言うと、画面中央からやや左に移動したあたり）は、

```glsl:modの計算
// x = 0.1, y = 1.0の場合
float newX = -0.1 - 1.0 * floor(-0.1 / 1.0); // -0.1 * 1.0 * -1.0
// newX = 0.9
```

となり、まさに上の図のような配色なるというわけです。
画面中央すぐ左あたりを注目してもらうと（`-0.1`とか）は`0.9`となってしまっています。（赤目の色になっている＝`x`の値が`1.0`に近い）
つまり空間的には**右端**となってしまいます。そのため、レイの位置が突然右端に移動してしまうのでオフセットを与えてあげないときれいな繰り返しにならない、というわけなんですね。

そしてオフセットをする理由を図にすると以下のようなイメージ。
![modのイメージ](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/ac1c0207-39c6-46ec-cf7c-c04e51137ebe.jpeg)
上の図は、**`y`の値だけ`mod`によって繰り返している様子**を示しています。
なので矢印の高さが常に同じところを繰り返しているのが分かるかと思います。
そして大事な点は、一律`-0.5`を加えて`y`の値をオフセットさせている点です。

こうすることで、`0 ~ 1`の間になってしまっていたものを`-0.5 ~ 0.5`という範囲に変換している、というわけなんですね。
なのでオフセットを**繰り返し単位の半分**だけ与えればきれいに中心が空間の真ん中にオフセットされる、というわけです。


-----------------------

##### 回転角度の計算

閑話休題。
話を戻しましょう。

Polar Modの関数の最後から見てみると、最終的に`return p * rot(-a);`として元のレイの位置をなにがしかの角度分回転させた位置を返しています。
なのでこの`a`の算出方法を理解すれば**Polar Mod**を理解することができるはずです。

では問題の`a`の算出方法を見てみましょう。

```glsl:aの算出
float a = atan(p.x, p.y) + pi / r;
float n = pi2 / r;
a = floor(a / n) * n;
```

計算自体はそこまで複雑なことを行っていません。
まず最初に行っているのは、引数に渡ってきた`p`の位置（2次元であることに注意）を、`atan`関数を使ってそのベクトルの角度を求めます。

ちなみに`atan`関数の[ドキュメント](https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/atan.xhtml)にはこう記載されています。

>atan — return the arc-tangent of the parameters

```
genType atan(genType y,
             genType x);
```

つまりアークタンジェントの`atan`ですね。
そして返ってくる値は**ラジアン**です。
そしてその返ってきた角度に対して、さらに**180度（=`pi`）**を分割数で割った数値を足しています。
これは通常のModによるリピートと同様に、レイの位置を補正するためのオフセット処理です。


-----------------------

ちなみに通常の`mod`では**分割数の半分だけ**オフセットさせると書きましたが、今回は半分ではない？　と思った方。
上書いた最初のコードではコメントしてますが、計算自体に**半分**を示す値はありませんが、そもそも`360`度を分割数の半分で、と考えると、`pi`はすでに`360`度の半分なので、**半分にする**処理がいらないわけですね。

-----------------------

そして次の`n`の計算は空間的なインデックスの計算です。
次の行で`floor(a / n) * n;`と、`floor`関数を利用しているのがその特徴です。

`floor`関数は小数点部分を切り捨て、整数部分だけを返してくれる関数です。
これをグラフにしてみると以下のようになります。
![floorのグラフ化](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/a5c8e949-af05-53f1-2eb7-ba7835444bbe.png)
このように、**整数部分だけ**を取り出すことで、その範囲のインデックスとして利用できるわけなんですね。
（なので、作品によっては`id`っていう変数名が付けられていたりも）

そしてその値と、さらに`n`を掛けることで値を正規化しているというわけです。
あとはこの値分だけ回転させれば無事、Polar Modの完成です。

詳細については[@kaneta1992先生の記事](https://qiita.com/kaneta1992/items/21149c78159bd27e0860)を御覧ください。こちらでは図解もされているので分かりやすいと思います。

## 最後に

最後に。
IFSとPolarModを使って作ってみたのがこちら↓
![IFSとPolarModを使ったサンプル](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/5616/8d2fd967-39e4-2cfc-c498-af101e5c1bdd.gif)
- [IFS universe](https://www.shadertoy.com/view/wlSGDd)


**IFS**と**PolerMod**を使うだけでだいぶ複雑な絵になるのが分かるかと思います。
これを色々応用すれば、比較的シンプルなシーンをいっきにリッチな見た目にできるので楽しいですね。