レイマーチングでAmbient Occlusion

## 概要

最近行われた「[GLSL Tech Night 2018](https://connpass.com/event/97565/)」で話された中で、特にレイマーチングの内容はとても面白く、レイマーチング熱が再燃してきております。
*（レイマーチングだけじゃなくて、早くレイトレとか色々なレイを飛ばせるようになりたい・・・）*

ちなみにそのときのスライド↓
[![20180926-210950.png](https://qiita-image-store.s3.amazonaws.com/0/5616/822e682c-0194-be66-a11d-cbe7479db1db.png)](https://www.slideshare.net/keimesuda/glsltech2018)
[GLSLtech2018 レイマーチングで半歩差のつく小技集](https://www.slideshare.net/keimesuda/glsltech2018)


スライドでは様々なTipsに加えて、AOを利用したGI（Global Illumination）についても考察されていました。
AO自体はレイマーチングに限った話ではありませんが、3DCGをよりよく見せるのに不可欠な要素です。

そこで今回は、レイマーチングでレンダリングした結果にAOを追加する方法を書いていこうと思います。
なお、今回実装した内容は[Distance Function](https://iquilezles.org/www/articles/distfunctions/distfunctions.htm)などでも有名なiq氏が質問サイトで回答していた内容を参考にさせていただいています。

質問サイトの回答↓
- [Raymarching Beginners' Thread :: pouët.net](http://www.pouet.net/topic.php?which=7920&page=21#c369990)

そしてそれを参考に、実際に実装して実行したレンダリング結果↓
![レイマーチングでAO](https://pbs.twimg.com/media/Dn_LO2mU8AA6QHv.jpg:large)
[Shadertoyでの実際に動くデモ](https://www.shadertoy.com/view/ltKcWc)もあります。

AOなしの方はオブジェクト同士が接している位置もパキッと別れていてCG色が強くなっているのが分かるかと思います。
一方、AOありのほうはオブジェクト同士が接している部分が暗くなっていて立体感が増していますね。

ただ、今回の実装はフェイクです。現実のものとは異なり、レイがヒットした位置の法線方向にレイを再度飛ばし直して影ができそうかを計算した結果を合成したものです。
実装にはハードコードされた定数もあり、ある程度目視での調整が入っています。


## AO実装のコンセプト

前述のように、今回実装したAOはフェイクなのでまずはその実装方針を書きます。

なお、いくつかの記事を見るとこの方法が比較的一般的なのかな、という印象でした。

**参考にした記事**

- [Raymarching Distance Fields](http://9bitscience.blogspot.com/2013/07/raymarching-distance-fields_14.html)
- [nical/GLSL-Raymarching - Github](https://github.com/nical/GLSL-Raymarching/blob/master/report/notes.mdwn)


### AO実装のフロー

まず最初に、詳細に触れる前に大まかなフローについて書いておきます。

1. 通常のレイマーチングでオブジェクトを検索
2. (1)によってオブジェクトが見つかったら、そのヒットした位置の法線方向に再びレイを飛ばす
3. (2)のレイの進行は通常のレイマーチングと異なり、「一定距離」ずつ一定回数進める（ループ処理）
4. (3)のループ処理内では、オブジェクトのヒット位置から現在のレイの位置と、その位置からのDistance Fieldの距離を比較する
5. 比較した結果と減衰率を掛けたものを現在の色として加算していく
6. 最後にそれを合成する

文章にするとなんだか微妙に分かりづらいですね；
ひとまずざっくり言ってしまうと、オブジェクト法線方向に一定距離レイを飛ばし、遮蔽対象になりそうな近くのオブジェクトがあったら影の色を加算していく、という具合です。


#### コードで理解

色々と言葉で説明するよりも先に、人によってはコードで見たほうが分かりやすいよ！　っていう人もいると思うので、今回実装したコードを先に載せておきます。

```glsl:AOの計算
// なんちゃってAOを計算する
//
// +-----------------+--------------------+
// | ro = Ray Origin | rd = Ray Direction |
// +-----------------+--------------------+
vec4 genAmbientOcclusion(vec3 ro, vec3 rd)
{
    vec4 totao = vec4(0.0);
    float sca = 1.0;

    for (int aoi = 0; aoi < 5; aoi++)
    {
        float hr = 0.01 + 0.02 * float(aoi * aoi);
        vec3 aopos = ro + rd * hr;
        float dd = distFunc(aopos);
        float ao = clamp(-(dd - hr), 0.0, 1.0);
        totao += ao * sca * vec4(1.0, 1.0, 1.0, 1.0);
        sca *= 0.75;
    }

    const float aoCoef = 0.5;
    totao.w = 1.0 - clamp(aoCoef * totao.w, 0.0, 1.0);

    return totao;
}
```

この関数を以下のように使ってAOを表現します。

```glsl:影としての計算
// ... 省略

vec4 ao = genAmbientOcclusion(pos, n);

// ... 省略（なにがしかの色計算）

col -= ao.xyz * ao.w;
```

`genAmbientOcclusion`関数の引数は、オブジェクトにヒットしたレイの位置とそのオブジェクトの法線方向です。
そしてその情報を元に計算された`ao`情報を、通常のライティングなどの「なにがしかの色計算」を行ったあとの加工として利用しています。

具体的には計算結果を減算しています。
減算している理由は、（今回の実装では）距離に応じてAOの値を加算していき、影の影響度として計算したためです。
なので、現在の色から減算してやることで影のように色が暗くなる、というわけです。

ちなみに、[参考にした記事](http://www.pouet.net/topic.php?which=7920&page=21#c369990)ではもう少し複雑な計算を行っていて、GIのように対象オブジェクトの色を反射するような計算がされています。
興味がある方はそちらの実装も見てみるといいと思います。


### オブジェクトを見つける

さてでは、実装の考え方について解説していきます。

まずは通常のレイマーチングを行ってオブジェクトを見つけるところから始めます。

レイマーチングの基本については以前自分が書いた記事「[[GLSL] レイマーチング入門 vol.1](https://qiita.com/edo_m18/items/034665d42c562da88cb6)」を参考にしてみてください。
具体的にどういう仕組で、どう実装するのかについてある程度詳細に書いています。
*（今回の記事ではレイマーチングの方法の詳細については触れません）*


まず、レイマーチングでは視点位置から画面を一枚のプレーンとして考え、そこにレイを飛ばします。
図にすると以下のような感じです。

![ドットonスクリーン.png](https://qiita-image-store.s3.amazonaws.com/0/5616/19f1fff4-af9e-8827-86e3-bc06af9a5a94.png)

こうして飛ばしたレイを、「ある距離」ずつ伸ばしていきます。
この「ある距離」とは、以下の図に示すように「最寄り」のオブジェクトと同じ距離分進めていきます。


![レイマーチイメージ.png](https://qiita-image-store.s3.amazonaws.com/0/5616/f2567f3b-4ba1-75ef-f4ae-09c060ea64a3.png)

図では緑の球体にレイがヒットしているのが分かるかと思います。
ここまでは普通のレイマーチングです。


### オブジェクト法線方向のオブジェクトを検索

レイがヒットしたオブジェクトが見つかったら、そこからAOの計算を開始します。
AOの計算ではヒットしたオブジェクトの「法線」方向にレイを再び飛ばします。
（レイマーチングの場合は法線は簡単に計算できるのでそれを利用します。詳細は[以前の記事](https://qiita.com/edo_m18/items/034665d42c562da88cb6)を参照）

このとき、通常のレイマーチングとは異なり「一定距離ずつ」レイを伸ばしていく点に注意してください。
なぜなら、新しくレイを飛ばし始める位置はすでにオブジェクトにヒットした位置です。
もし通常のレイマーチングと同様の手順でレイを伸ばしてしまうと、そのオブジェクトの表面からレイがまったく伸びていかないことになってしまうからです。


### レイの位置とDistance Fieldの距離

Distance Fieldは、一言で言えば距離関数群をまとめた関数のことです。
今回の例で言えば、床と球体、トーラスの3つのオブジェクトがシーン内にありますね。
これを実行しているのが以下の関数です。

```glsl:distance-field
float distFunc(vec3 pos)
{
    float d1 = distSphere(pos, size);
    float d2 = distPlane(pos, vec4(0.0, 1.0, 0.0, 2.0));
    float d3 = distTorus(pos, vec2(0.5, 0.2));
    return min(d3, min(d1, d2));
}
```

3つのオブジェクトの結果を合成して返しているだけですね。
これが「Distance Field」です。特に難しいところはありません。

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

英語ですが、冒頭で紹介したこちらの記事もDistance Fieldについて書かれているので参照してみてください。

- [Raymarching Distance Fields](http://9bitscience.blogspot.com/2013/07/raymarching-distance-fields_14.html)

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

続いてレイの処理です。

AO計算のためのレイの処理では、レイを「一定距離」ずつ進めていきます。
そして一定距離（ループ1回分）進めた時点でのレイの位置と、以下のふたつの距離をループごとに計算します。

1. 最初のオブジェクトとのレイのヒット位置との距離
2. AO計算用のレイの位置を用いて計算したDistance Fieldとの距離

そしてこのふたつを比較します。（つまり減算します）
仮に、AOの対象となりそうなオブジェクトが近くにない場合、この結果は常に`0`になります。
もし近くにオブジェクトがある場合はこの値はマイナスとなります。

なぜそうなるかというと、冷静になって考えると分かります。
そもそも、AO計算のためにレイを飛ばし始めた位置は、Distance Fieldによって求まったオブジェクトでした。
そこから出発して一定距離進んだ位置というのは、Distance Fieldの計算によって求まった位置から遠ざかる位置となります。

まだレイを進めていない場合のレイの位置はオブジェクト表面にありますよね。
そこでもう一度Distance Fieldの計算を行っても当然`0`になるのは自明です。

そしてそこから一定距離進んで再びDistance Fieldの計算を行えば、当然、一定距離進んだ分だけ遠ざかっているので、計算結果は進んだ距離が返ってくるわけです。

図にしてみると以下のようになります。

![AO計算イメージ.png](https://qiita-image-store.s3.amazonaws.com/0/5616/7e464bcc-e4e8-2d22-ac86-9db207dbea55.png)

ここで、`hr`はループごとに伸びていくレイの位置を表しています。
n回、ループが回った時点での位置と考えてください。

`aopos`は以下のように実装されています。

```glsl
vec3 aopos = ro + rd * hr;
```

つまり、レイの始点からレイの方向（つまりオブジェクトの法線方向）に、`hr`だけ進めた位置です。
図と一致していることを確認してください。

そう、つまり「近くにオブジェクトがない」場合、出発点から現在のレイまでの位置は、Distance Fieldの計算結果と一致しなければおかしいのです。
この事実を利用して、もし仮にその結果がマイナスになるようなことがあれば、それは「近くにオブジェクトが存在する」ということに他ならないのです。

もし近くにオブジェクトがあった場合を図解すると以下のようになります。

![AO計算イメージ2.png](https://qiita-image-store.s3.amazonaws.com/0/5616/9942a37f-5e08-1d2b-6571-72268969d483.png)

`hr`の矢印と`distFunc`の矢印を見ると、明らかに`distFunc`の矢印のほうが短いですよね。
これがマイナスになることの意味です。


AOは、近くにオブジェクトがある場合にその遮蔽効果によって影が落ちることをシミュレートするものです。
つまり、近くにオブジェクトが存在することが示唆された場合は、それを元に一定の値を加算していき、その結果を元に影を落としてやればAOらしい表現が可能になる、というわけです。

#### 減衰率を考慮

参考にさせていただいた記事の実装では、レイを進めるたびに率を`0.75`倍していきます。
これは（おそらく）減衰率的な役割を持っているものと思います。（変数名が`sca`なのでなんの略か分かりませんが・・・）

まぁ単純に考えれば、近くのオブジェクトといっても距離は様々で、AOとしての影響度は近いオブジェクトほど強く、遠くのオブジェクトほど弱くなります。
なので、遠くのものほど影響度を小さくしていくための工夫だと思います。

コードを引用すると以下になります。（ループ処理の内部）

```glsl
float hr = 0.01 + 0.02 * float(aoi * aoi); // aoiはイテレーション回数を表すint型。for (int aoi...
vec3 aopos = ro + rd * hr;
float dd = distFunc(aopos);
float ao = clamp(-(dd - hr), 0.0, 1.0);
totao += ao * sca * vec4(1.0, 1.0, 1.0, 1.0);
sca *= 0.75;
```

短いコードなので一行ずつ見ていきましょう。

```glsl
float hr = 0.01 + 0.02 * float(aoi * aoi); // aoiはイテレーション回数を表すint型。for (int aoi...
```

最初の行は、ループのイテレーション回数の2乗の距離を`0.02`倍し、`0.01`を加えたものです。
要は、イテレーション回数に応じてレイを徐々に遠くに移動させている、というわけですね。定数に関しては調整項目でしょう。

```glsl
vec3 aopos = ro + rd * hr;
```

`aopos`は、レイの始点とレイの方向へ、一定距離進めた位置の計算です。

```glsl
float dd = distFunc(aopos);
float ao = clamp(-(dd - hr), 0.0, 1.0);
```

次の2行はAOの影響度の計算部分です。
ここが肝ですね。

最初の`dd`は、前述のように、「現在のレイの位置」を元にしたDistance Fieldの計算結果です。
それを、続く行では計算済みの`hr`を減算することによってAOされるか否かの計算を行っています。
前述の通り、近くにオブジェクトがない場合は常にこの値は`0`になります。

```glsl
totao += ao * sca * vec4(1.0, 1.0, 1.0, 1.0);
sca *= 0.75;
```

最後に、計算した影響度と減衰率を掛けて影の色を計算しています。
そして次のループに備えて減衰率を0.75倍して率を減らしています。

計算の行数だけを見ればとてもシンプルに計算が行われているのが分かるかと思います。


## コード全文

最後に、Shadertoyに投稿したコードを全文載せて終わりにしたいと思います。
実際に動いているものは以下から確認できます。

- [Raymarching Ambient Occlusion - Shadertoy BETA](https://www.shadertoy.com/view/ltKcWc)


```glsl
vec3 cameraPos = vec3(0.0, -0.5, 1.5);
vec3 cameraDir = normalize(vec3(0.0, -0.3, -1.0));
vec3 lightDir = vec3(1.0, 1.0, 1.0);
float softShadow = 16.0;
float size = 0.5;

// sphereの距離関数
float distSphere(vec3 pos, float size)
{
    vec3 spPos = vec3(1.0, -1.8, -1.0);
    return length(pos - spPos) - size;
}

// Planeの距離関数
float distPlane(vec3 pos, vec4 n)
{
    return dot(pos, n.xyz) + n.w;
}

// トーラスの距離関数
float distTorus(vec3 p, vec2 t)
{
  vec2 tPos = vec2(0.0, -2.0);
  vec2 q = vec2(length(p.xy - tPos) - t.x, p.z);
  return length(q - tPos) - t.y;
}

float distFunc(vec3 pos)
{
    float d1 = distSphere(pos, size);
    float d2 = distPlane(pos, vec4(0.0, 1.0, 0.0, 2.0));
    float d3 = distTorus(pos, vec2(0.5, 0.2));
    return min(d3, min(d1, d2));
}

vec3 getNormal(vec3 pos)
{
    const float e = 0.0001;
    const vec3 dx = vec3(e, 0, 0);
    const vec3 dy = vec3(0, e, 0);
    const vec3 dz = vec3(0, 0, e);

    float d = distFunc(pos);

    return normalize(vec3(
        d - distFunc(vec3(pos - dx)),
        d - distFunc(vec3(pos - dy)),
        d - distFunc(vec3(pos - dz))
    ));
}

// 影を計算する
// オブジェクトの衝突位置からライト方向にレイを飛ばし、
// 遮蔽があったら影とする。
// また、ライトに向かうレイがなにかのオブジェクトに接近した場合は
// 最接近情報を保持し、影の影響度として利用する（ソフトシャドウ）
// +-----------------+--------------------+
// | ro = Ray Origin | rd = Ray Direction |
// +-----------------+--------------------+
float genShadow(vec3 ro, vec3 rd)
{
    // 距離関数の結果 = 距離
    float h = 0.0;

    // 現在のレイの位置
    float c = 0.001;

    // レイの最接近距離
    float r = 1.0;

    // シャドウ係数（濃さ）
    float shadowCoef = 0.5;

    // レイマーチにより影を計算する
    for (float t = 0.0; t < 50.0; t++)
    {
        h = distFunc(ro + rd * c);

        if (h < 0.001)
        {
            return shadowCoef;
        }

        // 現時点の距離関数の結果と係数を掛けたものを
        // レイの現時点での位置で割ったものを利用する
        // 計算結果のうち、もっとも小さいものを採用する
        r = min(r, h * softShadow / c);

        c += h;
    }

    return 1.0 - shadowCoef + (r * shadowCoef);
}

// なんちゃってAOを計算する
//
// +-----------------+--------------------+
// | ro = Ray Origin | rd = Ray Direction |
// +-----------------+--------------------+
vec4 genAmbientOcclusion(vec3 ro, vec3 rd)
{
    vec4 totao = vec4(0.0);
    float sca = 1.0;

    for (int aoi = 0; aoi < 5; aoi++)
    {
        float hr = 0.01 + 0.02 * float(aoi * aoi);
        vec3 aopos = ro + rd * hr;
        float dd = distFunc(aopos);
        float ao = clamp(-(dd - hr), 0.0, 1.0);
        totao += ao * sca * vec4(1.0, 1.0, 1.0, 1.0);
        sca *= 0.75;
    }

    const float aoCoef = 0.5;
    totao.w = 1.0 - clamp(aoCoef * totao.w, 0.0, 1.0);

    return totao;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord)
{
    vec2 p = (fragCoord.xy * 2.0 - iResolution.xy) / min(iResolution.x, iResolution.y);

    vec4 m = iMouse / iResolution.xxxx;
    
    vec3 col = vec3(0.0);

    vec3 cPos = cameraPos;
    vec3 cDir = normalize(vec3(m.x, -m.y, -1.0));
    vec3 cSide = normalize(cross(cDir, vec3(0, 1, 0)));
    vec3 cUp = normalize(cross(cSide, cDir));
    float targetDepth = 1.3;

    vec3 ray = normalize(cSide * p.x + cUp * p.y + cDir * targetDepth);

    float dist, depth;
    depth = 0.0;

    vec3 pos;

    const int maxsteps = 128;
    const float eps = 0.001;

    for (int i = 0; i < maxsteps; i++)
    {
        pos = cPos + ray * depth;

        dist = distFunc(pos);

        if (dist < eps)
        {
            break;
        }

        depth += dist;
    }

    float shadow = 1.0;

    if (dist < eps)
    {
        vec3 n = getNormal(pos);

        float diff = dot(n, lightDir);

        shadow = genShadow(pos + n + 0.001, lightDir);

        vec4 totao = genAmbientOcclusion(pos, n);

        float u = 1.0 - floor(mod(pos.x, 2.0));
        float v = 1.0 - floor(mod(pos.z, 2.0));
        if ((u == 1.0 && v < 1.0) || (v == 1.0 && u < 1.0))
        {
            diff *= 0.7;
        }

        col = vec3(diff) + vec3(0.1);
        col -= totao.xyz * totao.w;
    }

    fragColor = vec4(col, 1.0) * exp(-depth * 0.12);
}
```

## 参考記事

今回の実装にあたって、以下の記事も参考にさせていただきました。

- [オブジェクトにテクスチャなどを投影する - wgld.org](https://wgld.org/d/glsl/g019.html)