LIFE LOG(MochiMochi3D)

MochiMochi3D

3D関係の記事を書きます

4K Executable Graphics (exegfx)入門

4K Executable Graphics (exegfx)とは

RevisionやSessions, TokyoDemoFestといったデモパーティではカテゴリーの一つに4K Executable Graphicsというものがあります。これは最大ファイルサイズ4KBの実行ファイルから1枚の絵を出すというもので、長いので「exegfx」と略されることが多いみたいです。

他のカテゴリーとは違いリアルタイム性は求められないので、通常GLSLCompなどで使えないぐらい重いPathtracingやSDFなどを使うことができるという特徴があります。なので、exegfxで発表される作品は部屋に飾っておきたいと思えるような絵的にすごいものが多くあります。

youtu.be

Blossom

いざexegfxを作ろうと思ってもminifyツールを使って容量減らしたりなど4Kならではの大変さもあるので、１から自分で作るのはとても大変です。

嬉しいことにそうしたごちゃごちゃを全部勝手にやってくれるexegfxに特化したフレームワーク Blossomがあります。

github.com

Blossomは先述の通りexegfxに特化したフレームワークであり、Pathtracingを念頭に置いて設計されています。そのため、1frameの処理を行うdraw, 最終的に累計したdrawから最終出力を出すpresentの２つのレンダリングモードがあり、それぞれglslで処理を書いていくというような感じの構成になっています。

また、Debug,Release,Captureの3つのビルド設定があり、Releaseビルドではややこしいminifyを自動的にやってくれるという便利な機能があります。そのため、本当にglsl書くだけに集中できるので非常に便利なフレームワークとなっています。

環境構築

先ほどのリンクからBlossomのgithubに飛んで頂いて、zipでダウンロードするなりcloneするなりでローカルに落としてください。

BlossomはVisual Studioのソリューションという形で提供されており、フォルダ内のblossom.slnがそれです。Visual Studioで開くとこんな感じの構成になっています。

draw.fragとpresent.fragが先ほど述べたdraw,presentの処理を書くところであり、基本的にこの2つのフラグメントシェーダーを触っていくだけになります。

Blossomの構成

とりあえずDebugビルドでやってみてこんな感じのデフォルトのレンダリング結果がスクリーンの左上あたりに出ているのであれば問題なくインストール出来ています。

デフォルトのビルド結果

shaderの構成

drawは後述するconfigで指定したサンプル回数だけ呼び出されるシェーダーです。draw shaderは通常のGLSLシェーダーのように解像度iResolutionと現在のフレーム数iFrameがuniform変数として渡され、最終的な出力をfragColorに入れる形で書かれます。

#version 430
layout(location = 0) out vec4 fragColor;
layout(location = 0) uniform vec4 iResolution;
layout(location = 1) uniform int iFrame;


void main()
{
//省略
    fragColor = vec4(color,1.0);
}

普通のシェーダーと違う点としてはこのdraw shaderで出力されるfragColorはaccumulateTexというテクスチャに加算されていくという点です。例えば10回サンプリングされたとすれば、10回のdraw shaderが呼び出され、その結果の総和をaccumulateTexに格納するということになります。

こうしてdraw shaderの結果を集めたをaccumulateTexはpresent shaderに送られます。present shaderはdraw shaderのサンプリングが終わった後１度だけ呼び出されるシェーダーで、accumulateTexから色を算出するポストエフェクト的な感じの立ち位置になります。デフォルトのpresent shaderはこんな感じで受け取ったaccumulateTexから平均値を出力するようなプログラムになっています。

#version 430
layout(location = 0) out vec4 fragColor;
layout(location = 0) uniform vec4 iResolution;
layout(binding = 0) uniform sampler2D accumulatorTex;

void main()
{
    // readback the buffer
    vec4 tex = texelFetch(accumulatorTex,ivec2(gl_FragCoord.xy),0);

    // divide accumulated color by the sample count
    vec3 color = tex.rgb / tex.a;

    /* perform any post-processing you like here */

    // for example, some B&W with an S-curve for harsh contrast
    //color = smoothstep(0.,1.,color.ggg);

    // present for display
    fragColor = vec4(color,1);
}

大雑把なShaderの構成

以上がBlossomにおけるシェーダーの構成です。draw shaderでは1サンプルで行う処理を書き、present shaderでは各draw shaderの色の平均値の計算+posteffectを書く、という風な感じで作品を作っていきます。パストレーシングを目的とするならばdraw shaderにパストレの1サンプルの処理を書き、present shaderでaccumulateTexから平均値を得るというような形で書くことが出来ます。

ちょっと注意としてpresent shaderは直接サンプル数を得られないので各draw shaderでalphaに+1する形でaccumulateTexのアルファにサンプル数を格納して、そのアルファからサンプル数から得ています。なのでdraw shaderにおけるfragColorのalpha値は必ず1にしておく必要があります。(特殊な事しない限り)

Build設定

BlossomではビルドにDebug, Release, Captureの3種類があります。それぞれ、Debugはとりあえずの実行、Releaseは圧縮込みの実行、Captureはレンダリング結果を画像保存する実行ファイルのビルドというようになっています。各種ビルドの解像度などのレンダリング設定は後述するconfigで設定することが出来ます。

ビルドの種類

試しにReleaseビルドをしてみましょう。ビルド項目をReleaseに変更した上でF5を押してビルドを開始してみると、なんか圧縮してそうなウィンドウが出てきて、真っ暗な画面が映ると思います。少しすれば先ほどデバッグで見た画像が全画面に映ると思います。

ビルドがReleaseになっていることを確認

圧縮してそうなウィンドウ

Releaseビルドをした後でBlossomのフォルダを見てみるとReleaseというフォルダが追加され、その中にビルドによって生成された実行ファイルが置かれています。これが本番環境で実行されるexeファイルであり、ちゃんと実行してみると先ほどの様に真っ暗な画面の後にレンダリングされた画像が現れます。

ファイルサイズを見てみると2KBと元のコードを考えると驚くほど圧縮されていることが分かります。実際のコンポではこうして生成されたexeをそのまま提出して、実際のマシンでexeを実行してその絵で判定されます。

Blossom-main/Releaseにあるexeファイル ちゃんと圧縮され2KBの容量になっている

どれぐらい圧縮されるのかはまたややこしい問題になってくるし私も良くわからないので割愛させて頂きますが、私がTDFで提出した作品だと完成ちょっと前ぐらいでようやく4Kを超えたという感じなので、でかい作品じゃない限りは結構自由にコードを書いても容量の問題はない気がします。

次にCaptureビルドをやってみましょう。ビルド設定をCaptureに変更し、ビルドを行ってみるとDebugみたくレンダリング中の画像が見え、レンダリングが終わると勝手に終了します(デフォルトの設定だとクソデカ画面になりますがそれは解像度が高く設定されているだけなので問題ないです）。

ビルド設定

その後Blossomのフォルダを覗いてみるとCaptureというフォルダが追加され、その中に先ほどのビルドで生成されたexeファイルが置かれています。このexeを実行すれば先ほどの実行の後、レンダリング結果がpng形式で保存されます。

Captureフォルダ

Captureで生成したpng画像(解像度は落としてます)

このようにしてレンダリング画像を保存してくれる設定もあります。コンポでは結果の確認のため見本のレンダリング画像の同梱が義務付けられているので、これで生成した画像を提出フォルダに入れる感じになります。

Blossomのconfig

Blossomは画面解像度やサンプル数といったレンダリングに関わる設定はconfig.hでいじることが出来ます。config.hの中身を見てみるとこんな感じで各種設定がマクロとして定義されており、この値を変える形でレンダリングの設定を行います。

#pragma once

// like Leviathan, we have some risky tricks to shave the last few bytes
#define DESPERATE 0

// releasing at Revision? here's a handy toggle for some compo-safe config presets
#define REVISION_RULESET 0


#if _DEBUG
    #define WINDOW_FULLSCREEN 0
    #define WINDOW_AUTO_SIZE 0
    #define CANVAS_WIDTH 1280
    #define CANVAS_HEIGHT 720

    #define RENDER_MAX_TIME_MS 30000
    //#define RENDER_MIN_SAMPLES 1
    //#define RENDER_MAX_SAMPLES 256

    #define RENDER_PROGRESSIVE 1
#endif

#if RELEASE
    // party release config, based on Revision ruleset
    #if REVISION_RULESET
        #define WINDOW_FULLSCREEN 1
        #define WINDOW_AUTO_SIZE 0
        #define CANVAS_WIDTH  1920
        #define CANVAS_HEIGHT 1080

        #define RENDER_MAX_TIME_MS 30000 // Revision Rules: don't exceed 30 seconds

        #define RENDER_MAX_SAMPLES 256 // customizable, optional

        #define RENDER_PROGRESSIVE 0 // Revision Rules: no progressive rendering.
    #else
        #define WINDOW_FULLSCREEN 1
        #define WINDOW_AUTO_SIZE 1

        #define RENDER_EXACT_SAMPLES 256 // without the constraints of party rules, we can set an exact quality bar, if we want.

        #define RENDER_PROGRESSIVE 0
    #endif
#endif

#if CAPTURE
    #define WINDOW_AUTO_SIZE 0
    #define CANVAS_WIDTH 3840
    #define CANVAS_HEIGHT 2160

    #define RENDER_EXACT_SAMPLES 1024

    // which formats to save
    #define CAPTURE_SAVE_PNG 1
    #define CAPTURE_SAVE_JPG 0
    #define CAPTURE_SAVE_U8_BIN 0
    #define CAPTURE_SAVE_F32_BIN 0
#endif



// ==== housekeeping stuff, you don't need to touch this ==== //

#if CAPTURE
    #define WINDOW_FULLSCREEN 0

    #define RENDER_PROGRESSIVE 1

    #if !CAPTURE_SAVE_PNG && !CAPTURE_SAVE_JPG && !CAPTURE_SAVE_U8_BIN && !CAPTURE_SAVE_F32_BIN
        #define CAPTURE_SAVE_PNG 1
    #endif
#endif

各種定数の意味は以下の様になっています。

  • WINDOW_FULLSCREEN : レンダリング時にフルスクリーンにするか否か
  • WINDOW_AUTO_SIZE : 画面解像度を自動で設定するか
  • CANVAS_WIDTH/HEIGHT : 画面解像度の指定、WINDOW_AUTO_SIZEがfalseの時に作動
  • RENDER_MAX_TIME_MS : レンダリング時間の上限 本番環境で出来るだけサンプルを取りたいのであればこれを指定
  • RENDER_MIN_SAMPLES : 最低サンプリング数
  • RENDER_MAX_SAMPLES : 最高サンプリング数
  • RENDER_EXACT_SAMPLES : サンプリング数の指定、これが定義されていた場合RENDER_MAX_TIME/MIN_SAMPLE/MAX_SAMPLEは無視され、確実にこの値でサンプルが行われる
  • RENDER_PROGRESSIVE : レンダリング中に途中経過を表示するか、これがOFFだとサンプリングが終わるまで真っ暗で終わったら表示される形になる
  • CAPTURE_SAVE_PNG : CAPTUREにおいてPNGで保存するか否か
  • CAPTURE_SAVE_JPG : CAPTUREにおいてJPGで保存するか否か
  • CAPTURE_SAVE_U8_BIN : 8bit RGBで保存する
  • CAPTURE_SAVE_F32_BIN : float32 RGBで保存する
  • DESPERATE : 危険な方法でバイト数を減らすモード（どのような動作かはちょっとわからないです）
  • REVISION_RULESET : REVISIONのレギュレーションにするモード

例えばDebugビルドでもフルスクリーン表示にしたいというのがあれば、config.hのDebug部分に#define WINDOW_FULLSCREEN 1を追加することで実現できます。

#if _DEBUG
    #define WINDOW_FULLSCREEN 1
    #define WINDOW_AUTO_SIZE 0
    #define CANVAS_WIDTH 1280
    #define CANVAS_HEIGHT 720

    #define RENDER_MAX_TIME_MS 30000
    //#define RENDER_MIN_SAMPLES 1
    //#define RENDER_MAX_SAMPLES 256

    #define RENDER_PROGRESSIVE 1
#endif

exegfxのサンプル

UV Test

まずはUVを単に出してみるプログラムを書いてみましょう。今回はアスペクト比をyに合わせてuv.yが[-1,1]になるようなUVとして、draw shaderを以下の様に書き換えます。

これを実行してみたらいつものUVが出てきます

UV

ここでちょっとだけaccumulateの動作の確認をしてみましょう。present.fragは(コメントアウト部分除いて)以下の処理が行われています。先述の通りaccumulatorTexからRGBとサンプル数(alpha)を読み込み、平均値をcolorとして表示しています。

ここで試しにtex.aで割らずに平均値ではなくそのまま累計した値を出してみましょう。0以下じゃない限りどんどん値は増えていくので各領域でぱっきり色が分かれる形で画像として表示されるはずです。

accumulateしたRGB

こんな感じで実際accumulatorTexにはdraw.fragの結果がどんどん値が加算されていることがわかると思います。

SDFの表示

では次にSDFの球のノーマルを表示してみましょう。レイマーチングの方法とかは別段言うことはないのですが、後々それぞれのSDFにマテリアルを割り振ったりするので拡張できるようにstructを使って書いていこうと思います。まずはSDF周りの関数を以下の様に定義することとします。

//球のSDF
float sd_sphere(vec3 p){
    return length(p) - 0.5; 
}

//最も近いSDFの情報を入れる
struct SDFInfo{
    int index;
};

//シーンのSDFを構築する関数
float map(vec3 p,inout SDFInfo info){
    float d;
    d = sd_sphere(p);
    info.index = 0;
    return d;
}

//法線を求める関数
vec3 get_normal(vec3 p){
    vec2 eps = vec2(0.001,0.0);
    SDFInfo dammy;
    return normalize(vec3(
        map(p+eps.xyy,dammy)-map(p-eps.xyy,dammy),
        map(p+eps.yxy,dammy)-map(p-eps.yxy,dammy),
        map(p+eps.yyx,dammy)-map(p-eps.yyx,dammy)
    ));
}

SDFについては特に問題はないですが、シーンのSDFの情報を送れるようにSDFInfoというstructを用意し、map関数でSDFInfoをいじれるように設計しています。後の実装用にindexというプロパティを作っていますが特に今は使わないのでとりあえず0としておきます。

次はレイマーチングの処理を書いていきます。今回の演習ではレイマーチングは衝突時にtrue,非衝突時にfalseを返す「衝突判定」的な関数として定義し、衝突時にその衝突地点の情報を送れるようにSurfaceInfoというstructを用意しています。

//衝突地点の情報
struct SurfaceInfo{
    vec3 color;
    vec3 normal;
    vec3 position;
};

//レイマーチング、衝突したらtrue,しなかったらfalseを返す
#define MAX_STEP 100
bool raymarching(vec3 ro,vec3 rd,inout SurfaceInfo info){
    float dist = 0.0;
    float sum_d = 0.0;
    SDFInfo sdf_info;
    for(int i = 0; i < MAX_STEP; i++){
        dist = map(ro + rd * sum_d,sdf_info);
        if(dist < 0.001){
            info.position = ro + rd * sum_d;
            info.color = vec3(1.0); 
            info.normal = get_normal(info.position);
            return true;        
        }
        sum_d += dist;
    }

    info.color = vec3(0.0);
    info.normal = vec3(0.0);
    return false;
}

SurfaceInfoは衝突地点の情報を保持する構造体で今は法線(normal)、衝突位置(position)、色(color)を持って置くように定義しています。

以上で実装は終わったため、main関数でraymarching関数を使ってSDFの法線を出してみましょう。カメラの位置については(0,0,-3),方向はzの正方向を向くように作っています。

void main()
{
    vec2 uv = (gl_FragCoord.xy * 2.0 - iResolution.xy)/iResolution.y;
    vec3 color = vec3(0.0);

    //カメラのレイ
    vec3 cam_ori = vec3(0.0,0.0,-3.0);
    vec3 cam_dir = normalize(vec3(uv,1.0));

    SurfaceInfo info;
    if(raymarching(cam_ori,cam_dir,info)){
        //衝突時の処理    
        color = info.normal * 0.5 + 0.5;
    }
    else{
        //衝突しなかった時の処理
        color = vec3(0.0);
    }

    fragColor = vec4(color,1.0);
}

演習のレンダリング画像

SDFモデリング

これでSDFの表示ができるようになったため、簡単に部屋のようなSDFを作ってみましょう。今回は図のような感じの天井に穴の開いた部屋を作ってみましょう。まずは部屋全体のBoxを作り、そこから内部のbox,天井の穴部分のboxでくりぬく形でこの部屋を作ることが出来ます。

部屋の図

部屋のSDFの作り方

実装は先ほど言ったことをそのままやればいいので、こんな感じでmax関数でsdfをくりぬいてあげれば問題ありません。sdfの結果もある程度分けておくと見やすいので部屋のsdfはroom_dという変数でまとめて最後に球のsdfと合成する形にしています。

//球のSDF
float sd_sphere(vec3 p){
    return length(p) - 0.5; 
}

//Box SDF
float sdf_box(vec3 pos,vec3 size){
    vec3 d = abs(pos) - size;
    return min(max(d.x,max(d.y,d.z)),0.0) + length(max(d,0.0));
}

//最も近いSDFの情報を入れる
struct SDFInfo{
    int index;
};

//シーンのSDFを構築する関数
float map(vec3 p,inout SDFInfo info){
    float d;
    d = sd_sphere(p);

    float room_d;
    float outside_d = sdf_box(p,vec3(17,7,17));
    float inside_d = sdf_box(p,vec3(15,5,15));
    float hole_d = sdf_box(p - vec3(0,5,0),vec3(5,4,5));
    room_d = max(outside_d,-inside_d);
    room_d = max(room_d,-hole_d);

    info.index = 0;

    d = min(d,room_d);
    return d;
}

部屋のSDF

簡単なライティング

次は平行光源を付けて簡単なライティングをやってみましょう。main関数に色々やると見づらいのでカメラのレイを与えると色を返すrender関数を用意してそこで今後書いていきましょう。とりあえず今までのレイマーチングして色を返す処理の部分をrenderに移植しましょう。

vec3 render(vec3 ro,vec3 rd){
    vec3 color = vec3(0.0);
    SurfaceInfo info;
    if(raymarching(ro,rd,info)){
        //衝突時の処理    
        color = info.normal * 0.5 + 0.5;
    }
    else{
        //衝突しなかった時の処理
        color = vec3(0.0);
    }

    return color;
}

void main()
{
    vec2 uv = (gl_FragCoord.xy * 2.0 - iResolution.xy)/iResolution.y;
    vec3 color = vec3(0.0);

    //カメラのレイ
    vec3 cam_ori = vec3(0.0,0.0,-3.0);
    vec3 cam_dir = normalize(vec3(uv,1.0));

    color = render(cam_ori,cam_dir);

    fragColor = vec4(color,1.0);
}

ライトの方向はLIGHT_DIRと定義して、疑似環境光ありのランバート反射を書いてみましょう。今までの実装で物体の色と法線は既に取れるので次のような形で実装することが出来ます。

#define LIGHT_DIR normalize(vec3(0.5,1.0,0.0))

vec3 render(vec3 ro,vec3 rd){
    vec3 color = vec3(0.0);
    SurfaceInfo info;
    if(raymarching(ro,rd,info)){
        //衝突時の処理    
        color = info.color * (max(dot(info.normal,LIGHT_DIR),0.0) + 0.2);
    }
    else{
        //衝突しなかった時の処理
        color = vec3(0.0);
    }

    return color;
}

平行光源とLambert反射の追加

これだけだと味気ないのでシャドウを付けてみましょう。シャドウの付け方は色々ありますがここでは光源方向に物体があるか否かの0,1でハードシャドウを作ることとします。これは単純に衝突地点からライトの方向に再びレイを飛ばし、何かに当たったら影になるようにすればいいだけです。実装としてはこんな感じです。

#define LIGHT_DIR normalize(vec3(0.5,1.0,0.0))

vec3 render(vec3 ro,vec3 rd){
    vec3 color = vec3(0.0);
    SurfaceInfo info;
    if(raymarching(ro,rd,info)){
        //衝突時の処理    
        vec3 shadow_dir = LIGHT_DIR;
        vec3 shadow_ori = info.position + shadow_dir * 0.02;
        SurfaceInfo shadow_info;
        bool hit = raymarching(shadow_ori,shadow_dir,shadow_info);
        float shadow = 1.0 - float(hit);

        color = info.color * (max(dot(info.normal,LIGHT_DIR),0.0) * shadow + 0.2);
    }
    else{
        //衝突しなかった時の処理
        color = vec3(0.0);
    }

    return color;
}

シャドウ有の画像

ここでちょっと注意なのが次のレイのpositionにshadow_dirの方向にちょっとだけバイアスを与えています。positionは既にSDFの値が0.001になる場所であるため、バイアスをかけずにレイを飛ばすと最初の評価で元の場所で衝突しているものとみなされてしまいます(数値誤差による場合もあります)。これは自己衝突と呼ばれていて、バイアスを与えずに実行してみるとあらゆる場所が自己衝突して全部影扱いされるという残念なことになります。

バイアス抜きの実装　全てが自己衝突により影扱いされてしまっている

なのでこのような形で次のレイの進む方向にちょっとだけバイアスをかけてあげる必要があります。バイアスの値はレイマーチングの閾値より結構大きめに設定することをお勧めします。不十分なバイアスだと自己衝突が防ぎきれずこんな感じ輪っか上のアーティファクトが出てしまいます。

バイアスを0.001にした図 輪っか状のアーティファクトがでる

以上の実装によってこんな形で簡単なシャドウを作ることが出来ます。

マテリアリング

では次にSDFごとに色を変えられるようにマテリアルの仕組みを実装しましょう。私の実装だと各SDFにindexを設定し、衝突時最も近かったSDFのindexを基にマテリアルを設定するという形で実装しています(もっといい方法があるかもしれないです)。なのでmap関数は常に最も近かったSDFのindexを返せるように実装しておきます。

今回の演習は球のindexを0、部屋のindexを1としてmap関数で最も近いSDFのindexを返せるようにします。実装は簡単でSDFの合成の際にどっちが近いかを判定してindexを入れ替えるような形で実装できます。ということでmap関数にindexの更新をする処理を三項演算子で追加してみました。

float map(vec3 p,inout SDFInfo info){
    float d;
    d = sd_sphere(p);

    float room_d;
    float outside_d = sdf_box(p,vec3(17,7,17));
    float inside_d = sdf_box(p,vec3(15,5,15));
    float hole_d = sdf_box(p - vec3(0,5,0),vec3(5,4,5));
    room_d = max(outside_d,-inside_d);
    room_d = max(room_d,-hole_d);

    info.index = 0;
    info.index = (d < room_d) ? 1 : info.index; //Indexの判定

    d = min(d,room_d);
    return d;
}

これでSDFInfoのindexには最も近いSDFのindexが入っており、レイマーチングの衝突時にこれを基にマテリアルのカラーを変更してみましょう。

#define MAX_STEP 300
bool raymarching(vec3 ro,vec3 rd,inout SurfaceInfo info){
    float dist = 0.0;
    float sum_d = 0.0;
    SDFInfo sdf_info;
    for(int i = 0; i < MAX_STEP; i++){
        dist = map(ro + rd * sum_d,sdf_info);
        if(dist < 0.001){
            info.position = ro + rd * sum_d;
            info.normal = get_normal(info.position);
            if(sdf_info.index == 0){
                //Sphere
                info.color = vec3(0.2,0.8,0.2);
            }
            else{
                //Room
                info.color = vec3(0.8);
            }
            return true;        
        }
        sum_d += dist;
    }

    info.color = vec3(0.0);
    info.normal = vec3(0.0);
    return false;
}

これで実行してみると以下の様に球だけ緑色になり、ちゃんとマテリアルが適用されていることが分かります

マテリアル

if文で分ける方法は少々煩雑になりがちでなおかつマテリアルの追加も面倒なため、基本的にはグローバル変数としてカラーなどを配列で持って置き、indexで取得するというような実装がいいと思います(圧縮率も良くなります)。

#define NUM_MAT 2 //マテリアル数
vec3 color[NUM_MAT] = {vec3(0.2,0.8,0.2),vec3(0.8)};

//レイマーチング、衝突したらtrue,しなかったらfalseを返す
#define MAX_STEP 300
bool raymarching(vec3 ro,vec3 rd,inout SurfaceInfo info){
    float dist = 0.0;
    float sum_d = 0.0;
    SDFInfo sdf_info;
    for(int i = 0; i < MAX_STEP; i++){
        dist = map(ro + rd * sum_d,sdf_info);
        if(dist < 0.001){
            info.position = ro + rd * sum_d;
            info.normal = get_normal(info.position);
            info.color = color[sdf_info.index];
            return true;        
        }
        sum_d += dist;
    }

    info.color = vec3(0.0);
    info.normal = vec3(0.0);
    return false;
}

折角なので球を3つに増やし画像のように赤、緑、青で並べてみてください。実装例は以下にのっけて置きます。

球を三つ並べる

アンチエイリアシング

アンチエイリアシングをしてみましょう。今回は乱数にこちらのコードをお借りしてこんな感じで作ることとします。seedにちゃんと初期値を入れさえすれば後はrndを呼び出すだけで毎回異なる[0~1]の乱数を返してくれます。

www.shadertoy.com

uint seed;
uint PCGHash()
{
    seed = seed * 747796405u + 2891336453u;
    uint state = seed;
    uint word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
    return (word >> 22u) ^ word;
}

float rnd1()
{
    return PCGHash() / float(0xFFFFFFFFU);    
}

vec2 rnd2(){
    return vec2(rnd1(),rnd1());
}

これを使用してUVを1pixel内でランダムに動かせば問題なくアンチエイリアシングすることが出来ます。注意としてはフレームが異なればちゃんとシードも違うものにならなければいけないので、以下の様にiFrameを使ってシードを動かすようにしてください。

void main()
{
    seed = uint((iFrame+1) * (gl_FragCoord.x + iResolution.x * gl_FragCoord.y));

    vec2 uv = ((gl_FragCoord.xy + rnd2()) * 2.0 - iResolution.xy)/iResolution.y;
    vec3 color = vec3(0.0);

    //カメラのレイ
    vec3 cam_ori = vec3(0.0,0.0,-3.0);
    vec3 cam_dir = normalize(vec3(uv,1.0));

    color = render(cam_ori,cam_dir);

    fragColor = vec4(color,1.0);
}

アンチエイリアシング

（ちょっと光漏れがありますがこれは単にレイマーチングのイテレーションが足りないことから来ています。レイマーチングのイテレーションを増やせばいいですが、ちょっと今回の状況だとはレイマーチングと相性が悪いのでMAX_STEPを1000とかにしてもなくならないので気にしないでください）

ようやくaccumulateする意味がある実装に入ってきました。特にパストレとかではなくリアルタイム向けの手法を使い、アンチエイリアシングにaccumulateを使用する作品もあり、ここから色々作品を作ることもできます。

RTAO

次はパストレと共通の処理が多いRTAOと呼ばれる手法を実装します。RTAOはアンビエントオクルージョンの１つの手法で、衝突地点からいくつか法線方向にレイを飛ばし、指定半径内でいくつのレイが衝突するかの割合でAOを計算するという手法です。例えば図のようにレイが平面に衝突したとしたら、そこから10発のレイをランダムに衝突地点から放ちます。これらのレイは半径1の範囲で衝突するかを判定し、10発中4発衝突したらそのAOは(1.0 - 4.0 / 10.0)とするような感じで計算します。

RTAOのイメージ

RTAOの実装では衝突地点からランダムな方向をサンプリングする「方向サンプリング」の実装が必要になります。この方向サンプリングは一般的に接空間上での実装が行われます。なので今回は「接空間の基底計算」、「座標系変換」、「方向サンプリング」の3つの実装が必要になります。

接空間の基底計算の手法は色々ありますがシンプルにnormalからcrossで直交するベクトルを作る方法でやります。以下の様にnormalを与えるとtangent,binormalの基底ベクトルを算出してくれるtangentSpaceBasisという関数を作ります（どうゆうものなのかは割愛させて頂きます)。

void tangentSpaceBasis(vec3 normal,inout vec3 tangent,inout vec3 binormal){
    vec3 d = vec3(0,1,0);
    if(abs(normal.y) > 0.99) d = vec3(0,0,1);
    tangent = normalize(cross(normal,d));
    binormal = normalize(cross(tangent,normal));
}

次は座標系の変換を行う関数を作ります。world座標から与えられた基底の座標系の座標へと変換するworldToLocalとその逆を行うlocalToWorldは以下の様に書かれます。これは単純に線形代数であるような座標変換の式をそのままやっただけのものです。

vec3 worldToLocal(vec3 tangent,vec3 normal,vec3 binormal,vec3 world){
    return vec3(dot(world,tangent),dot(world,normal),dot(world,binormal));
}

vec3 localToWorld(vec3 tangent,vec3 normal,vec3 binormal,vec3 local){
    return tangent * local.x + binormal * local.z + normal * local.y;
}

こうして接空間の変換を作ることが出来ました。では、次は方向サンプリングですがあらゆる方向に一様にサンプリングすることが出来る半球一様サンプリングを用います。半球一様サンプリングでは乱数2つを受け取り、方向の偏角\thetaと方位角\phiを算出し、そこから方向を出すというようなサンプリングを行います。詳しい内容についてはここでは省略しますが以下のように実装することが出来ます。

vec3 hemisphereSampling(vec2 uv){
    float theta = acos(uv.x);
    float phi = 2.0 * PI * uv.y;
    return vec3(sin(theta) * cos(phi),cos(theta),sin(theta) * sin(phi));
}

ここで使用してるPIは単純に円周率です。マクロとして数値を定義してもいいですが、短縮の関係でacosの数値で定義した方がいいらしいです。

const float PI = acos(-1);

また、RTAOでは衝突したとしても閾値の値以下の距離でないと衝突したとみなしません。なのでレイマーチング部分で衝突地点までの距離を返すように機能を追加しましょう。SurfaceInfoに衝突地点までの距離を入れるray_distという変数を用意し、レイマーチング中のsum_dを返すように実装します。衝突しなかった場合は一応未定義動作を防ぐため、適当な大きな値を入れておきます。

bool raymarching(vec3 ro,vec3 rd,inout SurfaceInfo info){
    float dist = 0.0;
    float sum_d = 0.0;
    SDFInfo sdf_info;
    for(int i = 0; i < MAX_STEP; i++){
        dist = map(ro + rd * sum_d,sdf_info);
        if(dist < 0.001){
            info.ray_dist = sum_d;
            info.position = ro + rd * sum_d;
            info.normal = get_normal(info.position);
            info.color = color[sdf_info.index];
            return true;        
        }
        sum_d += dist;
    }
    info.ray_dist = 10000.0;
    info.color = vec3(0.0);
    info.normal = vec3(0.0);
    return false;
}

以上でRTAOの準備が整いました。アルゴリズムは衝突→接空間で半球一様サンプリング→ワールド空間に方向を戻す→衝突点から一定範囲内で衝突したらAOに+1→繰り返しというような流れで、最終的にサンプリング数で割ることでAOの値を得ます。閾値を1.0としてこのアルゴリズムを実装してみるとこのような感じになります。

#define RTAO_NUM 16
vec3 render(vec3 ro,vec3 rd){
    vec3 color = vec3(0.0);
    SurfaceInfo info;
    if(raymarching(ro,rd,info)){
        //衝突時の処理    
        vec3 shadow_dir = LIGHT_DIR;
        vec3 shadow_ori = info.position + shadow_dir * 0.02;
        SurfaceInfo shadow_info;
        bool hit = raymarching(shadow_ori,shadow_dir,shadow_info);
        float shadow = 1.0 - float(hit);

        vec3 tangent,binormal;
        tangentSpaceBasis(info.normal,tangent,binormal);

        float RTAO = 0.0;
        for(int i = 0; i < RTAO_NUM; i++){
            vec3 dir = hemisphereSampling(rnd2());
            dir = localToWorld(tangent,info.normal,binormal,dir);
            vec3 ori = info.position + dir * 0.02;
            SurfaceInfo rtao_info;
            bool hit = raymarching(ori,dir,rtao_info);

            if(rtao_info.ray_dist < 1.0){
                RTAO += float(hit);
            }
        }
        RTAO = (1.0 - RTAO / RTAO_NUM);

        color = info.color * (max(dot(info.normal,LIGHT_DIR),0.0) * shadow + 0.2) * RTAO;
    }
    else{
        //衝突しなかった時の処理
        color = vec3(0.0);
    }

    return color;
}

RTAO

こうしてRTAOを作ることが出来ました。レイマーチングでは他のAOの手法はいくつもありますが、ここではパストレでも共通の処理があるという理由のためRTAOを実装しました。これさえ実装できればもはやあまり追加することなくパストレーシングをすることが出来るようになります。

パストレーシング

では折角なので次はパストレーシングのコードを書いてみましょう。ここではあまり詳しくパストレーシングの理論については深く話さないで、簡単にアルゴリズムだけ説明して実装をお見せするというような感じでやっていきます。意外と実装はそこまで複雑ではないので安心してください。

そもそもパストレーシングとは何ぞやという人もいるかもしれません。パストレーシングとは現実的にあり得る光の経路(パス)を探索することで現実と同じように光の計算を行う手法です。PBRに基づいて反射などを計算すればパストレーシングはとてもリアルな画像を出してくれます。

パストレーシングのアルゴリズムは比較的簡単です。流れとしてはこのような感じです。

  1. 衝突判定
  2. if 光源に衝突or何にも衝突しない 2.1光源の寄与を加算して終了
  3. 衝突位置で次に進むレイの方向をサンプリング
  4. 反射率の計算
  5. 次のレイをセット
  6. 繰り返し

パストレのイメージ

実装を見ながら説明した方が分かりやすいと思いますので、実装を見ていきましょう。

vec3 cosineSampling(vec2 uv,inout float pdf){
    float theta = acos(1.0 - 2.0f * uv.x) * 0.5;
    float phi = 2.0 * PI * uv.y;
    pdf = cos(theta) / PI;
    return vec3(sin(theta) * cos(phi),cos(theta),sin(theta) * sin(phi));
}

vec3 IBL(vec3 dir){
    return vec3(1.0);
}

#define MAX_DEPTH 10
vec3 render(vec3 ro,vec3 rd){
    vec3 LTE = vec3(0.0); //最終結果を入れる場所
    vec3 throughput = vec3(1.0); //反射率を入れる場所
    
    vec3 ray_ori = ro;
    vec3 ray_dir = rd;

    for(int i = 0; i < MAX_DEPTH; i++){
        //衝突判定
        SurfaceInfo info;   
        if(!raymarching(ray_ori,ray_dir,info)){
            //衝突しなかった場合   
            LTE += throughput * IBL(ray_dir);   
            break;
        }

        //衝突した場合
        vec3 normal = info.normal;
        vec3 tangent,binormal;
        tangentSpaceBasis(normal,tangent,binormal);
        
        vec3 local_wo = worldToLocal(tangent,normal,binormal,-ray_dir);

        //方向サンプリング
        float pdf;
        vec3 local_wi = cosineSampling(rnd2(),pdf);

        vec3 wi = localToWorld(tangent,normal,binormal,local_wi);

        //BSDFの計算
        vec3 bsdf = info.color / PI; //Lambert
        float cosine = dot(wi,normal);

        //throughputの更新
        throughput *= bsdf * cosine / pdf;

        //レイの更新
        ray_dir = wi;
        ray_ori = info.position + ray_dir * 0.01;
    }

    return LTE;
}

まずは変数の部分から、LTEは最終的な結果を入れる場所でthrouputは反射率を入れる変数です。各反射において反射率を計算して、throughput にかけていくことで最終的な反射率を求めることができます。次に飛ばすレイはray_ori,ray_dirで定義しており、最初はカメラからのレイを代入しています

 vec3 LTE = vec3(0.0); //最終結果を入れる場所
    vec3 throughput = vec3(1.0); //反射率を入れる場所
    
    vec3 ray_ori = ro; 
    vec3 ray_dir = rd;

では、for文の中身を見ていきましょう。最初は衝突判定から始まります。

     SurfaceInfo info;   
        if(!raymarching(ray_ori,ray_dir,info)){
            //衝突しなかった場合   
            LTE += throughput * IBL(ray_dir);   
            break;
        }

ここで衝突しなかった場合空にレイは飛んでいくものとして考えます。その場合というのは通常のレンダリング同じようにIBLから光が来ているということになりますので、IBLから光の値を取り今までの経路の反射率が入っているthrouputをかけてあげた値がそのパスの明るさとなります。ここではそれをLTEに入れているというわけですね。これ以上反射はすることはないのでこの処理が終わり次第breakして終了とします。

今回IBLは適当に全方面から真っ白なvec3(1.0,1.0,1.0)の光が来ているものとして定義しています。

vec3 IBL(vec3 dir){
    return vec3(1.0);
}

次は物体に衝突した場合となりますが、まずは後々の計算のため接空間の基底を生成してrayが入ってきた方向ray_dirを接空間座標に変換します(local_wo)。ちなみに一般的にシェーディングでは衝突地点を原点として方向を考えるので、ここではray_dirを反転した-ray_dirを接空間に変換していることに注意してください。

        //衝突した場合
        vec3 normal = info.normal;
        vec3 tangent,binormal;
        tangentSpaceBasis(normal,tangent,binormal);
        
        vec3 local_wo = worldToLocal(tangent,normal,binormal,-ray_dir);

ここから先は理論的にちょっと難しい話になってきます。PBRの分野では光の反射率というのはBSDFと呼ばれる関数によって表されており、光が入射してきた方向ベクトル \omega_iと出ていく方向\omega_oとすると一般的に f(\omega_i,\omega_o)という関数として表現することが出来ます

例えば、私たちが良く使っていたLambertはAlbedo(カラー) \rhoに対して以下のようなBSDFとして定義されます。

\displaystyle{f_{Lambert}(\omega_i,\omega_o) = \frac{\rho}{\pi}}

ちょっとややこしい話ですが、ある方向から来た\omega_iから来た光が\omega_o方向に飛んで行った時に光の強さ(輝度)の減衰率という意味での「反射率」はBSDFに \cos \theta_iをかけた値として扱われます。理論的な話は長くなるのでそうゆうものだと考えてください。リアルタイムで使われるLambertに\cosが存在するのはこの話から来ています。

\displaystyle{反射率 = f(\omega_i,\omega_o) \cos{\theta_i}}

パストレーシングではランダムに次に進む方向をサンプリングするわけですが、取る方向がなんでもいいというわけではありません。BSDFはものによっては0になるような方向を持つものがあったりします。そんな方向を取ってしまっては反射率は既に０であり、例え光源に当たっても全く光が存在しないものになりますのでそのパスの寄与は0という値になります。

それを防ぐため、一般的には方向サンプリングはBSDFが高い値を取る部分を中心にサンプリングする「重点的サンプリング」という手法を使います。これはそれぞれBSDFに対応したサンプリング方法があり、Lambertでは以下のようなコサイン重点的サンプリングというものを使用します。実装ではこのように実装されており、接空間上で次に飛ばすレイの方向wiを求めています。

vec3 cosineSampling(vec2 uv,inout float pdf){
    float theta = acos(1.0 - 2.0f * uv.x) * 0.5;
    float phi = 2.0 * PI * uv.y;
    pdf = cos(theta) / PI;
    return vec3(sin(theta) * cos(phi),cos(theta),sin(theta) * sin(phi));
}
        //方向サンプリング
        float pdf;
        vec3 local_wi = cosineSampling(rnd2(),pdf);

        vec3 wi = localToWorld(tangent,normal,binormal,local_wi);

方向サンプリングの関数を見てみるとPDFという新たな数値が返されていることがわかると思います。これは方向サンプリングをした時、その方向がサンプリングされる確率のようなものを示す確率密度関数(Probability Density Function)と呼ばれるものです。

方向サンプリングは元々サンプリングされる方向がこのPDFの分布に従うように設計して作られるものであり、PDFの値が1に近い方向ほど多く取られるようになっています(なのでPDFはBSDFの式を基に作られています)。サンプリング方法それぞれで固有のPDFがあり、コサイン重点サンプリングのPDFはこのようになっています。

\displaystyle{PDF = \frac{\cos{\theta_i}}{\pi}}

なぜこんなものが必要かというとパストレーシングの理論的な背景に「モンテカルロ法」というものがあるためです。この方法は関数の評価値をその値が取られるPDFで割ってあげることでうまいこと期待値が積分値に一致するようにしており、パストレも同様にPDFで割るという作業が必要になります。

PDFの話はちゃんとやろうとなるととても大変なのでここではあまり話しませんが、実装上ではサンプリングしたPDFで割ってあげた値を反射率と見なして計算すれば問題ありません。なのでこんな感じに反射率を計算すれば問題ありません。こうして求められた反射率をthrouputにかけてあげているわけですね。

\displaystyle{反射率 = \frac{f(\omega_i,\omega_o) \cos{\theta_i}}{pdf}} 
        //BSDFの計算
        vec3 bsdf = info.color / PI; //Lambert
        float cosine = dot(wi,normal);

        //throughputの更新
        throughput *= bsdf * cosine / pdf;

そして最後に次に進むレイを設定して繰り返し何回も反射していくということになります

        //レイの更新
        ray_dir = wi;
        ray_ori = info.position + ray_dir * 0.01;

以上でパストレのアルゴリズムは終わりで、これで実行してみると次のようにとてもリアルな画像が得られます。

パストレーシング

未だとIBLのみの光源を考えていましたが、折角なのでSDFも光源にできるように書いていきましょう。実装は非常に簡単でマテリアルにライトの明るさとしてEmissionというパラメーターを追加しましょう。

//衝突地点の情報
struct SurfaceInfo{
    float ray_dist;
    vec3 color;
    vec3 emission; //ライトの明るさ
    vec3 normal;
    vec3 position;
};

#define NUM_MAT 4
vec3 color[NUM_MAT] = {vec3(0.2,0.8,0.2),vec3(0.8),vec3(0.8,0.2,0.2),vec3(0.2,0.2,0.8)};
vec3 emission[NUM_MAT] = {vec3(1.0),vec3(0.0),vec3(0.0),vec3(0.0)}; //ライトの明るさ

//レイマーチング、衝突したらtrue,しなかったらfalseを返す
#define MAX_STEP 300
bool raymarching(vec3 ro,vec3 rd,inout SurfaceInfo info){
    float dist = 0.0;
    float sum_d = 0.0;
    SDFInfo sdf_info;
    for(int i = 0; i < MAX_STEP; i++){
        dist = map(ro + rd * sum_d,sdf_info);
        if(dist < 0.001){
            info.ray_dist = sum_d;
            info.position = ro + rd * sum_d;
            info.normal = get_normal(info.position);
            info.color = color[sdf_info.index];
            info.emission = emission[sdf_info.index]; //ライトの明るさを返す
            return true;        
        }
        sum_d += dist;
    }
    info.ray_dist = 10000.0;
    info.color = vec3(0.0);
    info.normal = vec3(0.0);
    return false;
}

こうしてマテリアルにEmissionというパラメーターを追加できました。パストレーシングのコードに戻り、Emissionという値が0以上であればその物体は光源として扱うようにしましょう。光源だった場合の処理はIBLと同様でthrouputと光源のemissionの値をかけたものをLTEに加算するだけです。

#define MAX_DEPTH 10
vec3 render(vec3 ro,vec3 rd){
    vec3 LTE = vec3(0.0); //最終結果を入れる場所
    vec3 throughput = vec3(1.0); //反射率を入れる場所
    
    vec3 ray_ori = ro;
    vec3 ray_dir = rd;

    for(int i = 0; i < MAX_DEPTH; i++){
        //衝突判定
        SurfaceInfo info;   
        if(!raymarching(ray_ori,ray_dir,info)){
            //衝突しなかった場合   
            LTE += throughput * IBL(ray_dir);   
            break;
        }

        if(length(info.emission) > 0.0){
                        //光源に衝突した場合
            LTE += throughput * info.emission;
            break;
        }

        //衝突した場合
        vec3 normal = info.normal;
        vec3 tangent,binormal;
        tangentSpaceBasis(normal,tangent,binormal);
        
        vec3 local_wo = worldToLocal(tangent,normal,binormal,-ray_dir);

        //方向サンプリング
        float pdf;
        vec3 local_wi = cosineSampling(rnd2(),pdf);

        vec3 wi = localToWorld(tangent,normal,binormal,local_wi);

        //BSDFの計算
        vec3 bsdf = info.color / PI; //Lambert
        float cosine = dot(wi,normal);

        //throughputの更新
        throughput *= bsdf * cosine / pdf;

        //レイの更新
        ray_dir = wi;
        ray_ori = info.position + ray_dir * 0.01;
    }

    return LTE;
}

SDFを光源にしてみる

ちょっと分かりずらかったのでIBLを消して、球の明るさをもう少し上げてみましょう

vec3 IBL(vec3 dir){
    return vec3(0.0);
}
vec3 emission[NUM_MAT] = {vec3(10.0),vec3(0.0),vec3(0.0),vec3(0.0)}; //ライトの明るさ

このような感じでパストレでは面光源も簡単に作ることが出来ます。

以上で簡単なパストレの実装が終わりました。試しにこのコードでRealeaseビルドしてみて、どれぐらいの容量になっているか見てみると、ようやく2Kを超えた程度でまだまだ余裕がありますね。

exeの容量

こんな感じでexegf作ることが出来ました。minimalなパストレのコードですがSDFをちょっと変えたりするだけでも結構いい感じの作品を作ることが出来るので是非挑戦してみてください。好きなSDFを表示するだけでも楽しいですよ

作品例

中級者向け

ここからはちょっと発展的な手法について簡単に書いておきます。

ロシアンルーレット

throughputが小さくなってしまったパスというのは例え次に光源に当たっても大した寄与を持たないことが多いです。そうしたパスは正直計算する意義があまりないので何とか棄却したいところですが、愚直にthrouputが小さくなったらもう計算しないみたいな手法ではそれはそれで問題がありそうです。

そうしたパスの棄却をうまいこと確率的に選択することでアンバイアスに軽量化することが出来る手法としてロシアンルーレットというものがあります。理論的に説明するとちょっとややこしい話がありますが、アルゴリズムとしてはすごい単純で各depthにおいて何らかの確率 pを設定し、乱数がそれ以下だったら継続、それ以上なら棄却するというような確率的に棄却の選択を行うだけです。継続だった場合は補正項としてその確率でthrouputを割ってあげる必要があります。

        float rossian_p = clamp(max(max(throughput.x,throughput.y),throughput.z),0.0,1.0);
        if(rossian_p < rnd1()){
            break;
        }
        throughput /= rossian_p;

このロシアンルーレットの良いこととして、確率 pはどんな時でも任意のものを使ってよいというところです。throughputが小さい物はある程度棄却したいと考えているため、throuputのRGBの内最大値を選択確率 pにすることでthrouputが全体的に小さい時にかなりの確率で棄却することが可能となります。そしてこの棄却をしたとしても平均的にはちゃんとロシアンルーレットなしの場合の結果と一致することが知られており、絵的にも問題がないという保証がされています。

vec3 render(vec3 ro,vec3 rd){
    vec3 LTE = vec3(0.0); //最終結果を入れる場所
    vec3 throughput = vec3(1.0); //反射率を入れる場所
    
    vec3 ray_ori = ro;
    vec3 ray_dir = rd;

    for(int i = 0; i < MAX_DEPTH; i++){
    
        float rossian_p = clamp(max(max(throughput.x,throughput.y),throughput.z),0.0,1.0);
        if(rossian_p < rnd1()){
            break;
        }
        throughput /= rossian_p;

        //衝突判定
        SurfaceInfo info;   
        if(!raymarching(ray_ori,ray_dir,info)){
            //衝突しなかった場合   
            LTE += throughput * IBL(ray_dir);   
            break;
        }
        
        if(length(info.emission) > 0.0){
                        //光源に衝突した場合
            LTE += throughput * info.emission;
            break;
        }

        //衝突した場合
        vec3 normal = info.normal;
        vec3 tangent,binormal;
        tangentSpaceBasis(normal,tangent,binormal);
        
        vec3 local_wo = worldToLocal(tangent,normal,binormal,-ray_dir);

        //方向サンプリング
        float pdf;
        vec3 local_wi = cosineSampling(rnd2(),pdf);

        vec3 wi = localToWorld(tangent,normal,binormal,local_wi);

        //BSDFの計算
        vec3 bsdf = info.color / PI; //Lambert
        float cosine = dot(wi,normal);

        //throughputの更新
        throughput *= bsdf * cosine / pdf;

        //レイの更新
        ray_dir = wi;
        ray_ori = info.position + ray_dir * 0.01;
    }

    return LTE;
}

コード量もかなり少なく済む割には（シーンによりますが）結構な高速化を見込める手法であり、通常のオフラインレンダリングでも使われている手法です。exegfxでもバイト数に見合う効果があるのではないかなと思います。

GGXの実装

なかなかLambertだけだとあまりマテリアルに味気なさを感じると思います。実際Specularのマテリアルを入れると劇的に絵が変わるので是非ともDiffuseだけではなくSpecularもマテリアルに組み込みたいところです。BSDF周りの部分を関数にしてまとめた方が後々拡張するとき便利なのでこんな感じでBSDFという関数でまとめちゃいます。

vec3 BSDF(vec3 wo,inout vec3 wi,SurfaceInfo info){
    float pdf;
    //Lambert
    wi = cosineSampling(rnd2(),pdf);
    vec3 bsdf = info.color / PI; 

    //Cosine Term
    float cosine = wi.y;

    return bsdf * cosine / pdf;
}
     vec3 local_wo = worldToLocal(tangent,normal,binormal,-ray_dir);
        vec3 local_wi;

        //throughputの更新
        throughput *= BSDF(local_wo,local_wi,info);

        vec3 wi = localToWorld(tangent,normal,binormal,local_wi);

SpecularのBSDFはデモシーンだとBlim-Phongがありますが、少々パラメーターが扱いづらいことや(許容するなら別に問題ではないですけど)現実に即していないということもありますし、一般的にCGで用いられるGGXを使用した方が個人的にはいい気がします。パストレのフレームワークだとそこまでコード量が多くなくて済みますし、もし複雑なマテリアルとかをちゃんとPBR的にやりたいのであればGGXが基本なのでこちらを使った方がいい気がします。GGXに関する話は私の記事を参照ください。

kinakomoti321.hatenablog.com

GGXは接空間表現であれば三角関数を使う必要はなく、こんな感じで各D,Gを書くことが出来ます(粗さは\alpha、異方性は考えないものとする)。この式は[Heitz 2018]Sampling the GGX Distribution of Visible Normalsに載っているものです。

float GGX_Lambda(vec3 v,float alpha) {
    float delta = 1.0f + (alpha * alpha * v.x * v.x + alpha * alpha * v.z * v.z) / (v.y * v.y);
    return (-1.0 + sqrt(delta)) / 2.0f;
}

float GGX_D(vec3 wm,float alpha) {
    float term1 = wm.x * wm.x / (alpha * alpha) + wm.z * wm.z / (alpha * alpha) + wm.y * wm.y;
    float term2 = PI * alpha * alpha * term1 * term1;
    return 1.0f / term2;
}

float GGX_G1(vec3 w,float alpha) {
    return 1.0f / (1.0f + GGX_Lambda(w,alpha));
}

float GGX_G2_HeightCorrelated(vec3 wi, vec3 wo,float alpha) {
    return 1.0f / (1.0f + GGX_Lambda(wi,alpha) + GGX_Lambda(wo,alpha));
}

GGXの重点的サンプリングはちょっと特殊で出射方向ベクトル \omega_iを直接サンプリングするのではなく、ハーフベクトル\omega_mをサンプリングします。ハーフベクトルのサンプリング方法は色々ありますが、一番シンプルな手法としては[Walter et al. 2007]Microfacet Models for Refraction through Rough Surfacesで発表されたDに対するサンプリングがあります。ここではWalter Samplingと呼称することとします。

Walter Sampling

vec3 ggx_halfsampling(vec2 uv,float alpha){
    float theta = atan(alpha * sqrt(uv.x) / sqrt(max(1.0 - uv.x,0.0)));
    float phi = 2.0 * PI * uv.y;
    return vec3(sin(theta) * cos(phi),cos(theta),sin(theta) * sin(phi));
}

こうして得られたハーフベクトルに対して、ハーフベクトルを法線とする反射ベクトルを\omega_iとしてサンプリングします。

 //GGX
    vec3 wm = ggx_halfsampling(rnd2(),alpha);
    wi = reflect(-wo,wm);

注意点として、wmは半球中のどの方向も取り得るのでその反射ベクトルである\omega_iは地面にめり込む方向\omega_i.y &lt; 0.0に行くことがあります。その場合は反射できなかったものとしてもう計算を打ち切るのが良いとされています。私の実装ではBSDFを0にしてロシアンルーレットで弾いてもらうようにしています。

 if(wi.y < 0.0){
        return vec3(0.0);
    }

PDFに関しては以下の様に求めることが出来ます。

PDF

 //Walter Sampling
    pdf = D * wm.y / (4.0 * dot(wm,wo));

また、このサンプリングより良いとされる可視法線分布(VNDF)サンプリングというものがあります。詳しい理論としては置いておいて、今のところ一番コードがシンプルになる手法でやるとこんな感じの処理です。([Dupuy, Benyoub 2023] Sampling Visible GGX Normals with Spherical Caps)

vec3 sampleVisibleNormal(vec2 uv, vec3 wo,float alpha) {
    vec3 strech_wo = normalize(vec3(wo.x * alpha, wo.y, wo.z * alpha));
    float phi = 2.0f * PI * uv.x;
    float z = fma((1.0f - uv.y), (1.0f + strech_wo.y), -strech_wo.y);
    float sinTheta = sqrt(clamp(1.0f - z * z, 0.0f, 1.0f));
    float x = sinTheta * cos(phi);
    float y = sinTheta * sin(phi);
    vec3 c = vec3(x, z, y);
    vec3 h = c + strech_wo;

    vec3 wm = normalize(vec3(h.x * alpha, h.y, h.z * alpha));
    return wm;
}

VNDFサンプリングにおけるPDFはこのような形で書くことが出来ます。

\displaystyle{PDF = \frac{D(\omega_m) |\omega_i \cdot \omega_m| G_1(\omega_i,\omega_m) }{4|\omega_i \cdot n| |\omega_m \cdot \omega_o|} } 
 //Visible Normal Sampling
    pdf = 0.25f * GGX_D(wm,alpha) * dot(wo, wm) / (dot(wm, wo) * abs(wo.y)* (1.0f + GGX_Lambda(wo,alpha)));

全体的な実装はこんな感じになっています。

vec3 BSDF(vec3 wo,inout vec3 wi,SurfaceInfo info){
    float pdf;
    //Lambert
    //wi = cosineSampling(rnd2(),pdf);
    //vec3 bsdf = info.color / PI; 

    float alpha = clamp(info.roughness * info.roughness,0.001,1.0);

    //GGX
    //WalterSampling
    //vec3 wm = ggx_halfsampling(rnd2(),alpha);

    //VisibleNoraml Sampling
    vec3 wm = sampleVisibleNormal(rnd2(),wo,alpha);
    wi = reflect(-wo,wm);

    if(wi.y < 0.0){
        return vec3(0.0);
    }

    float D = GGX_D(wm,alpha);
    float G = GGX_G1(wo,alpha) * GGX_G1(wi,alpha);
    vec3 F = shlickFresnel(info.color,dot(wm,wo));
    
    vec3 bsdf = D * G * F / (4.0 * wo.y * wi.y);
    
    //Walter Sampling
    //pdf = D * wm.y / (4.0 * dot(wm,wo));

    //Visible Normal Sampling
    pdf = 0.25f * GGX_D(wm,alpha) * dot(wo, wm) / (dot(wm, wo) * abs(wo.y)* (1.0f + GGX_Lambda(wo,alpha)));

    //Cosine Term
    float cosine = wi.y;

    return bsdf * cosine / pdf;
}

GGXの実装

Specular付のマテリアル実装

Diffuse、Specularのマテリアルはこうして実装できました。次はDiffuseとSpecularを組み合わせた現実的なBSDFを作ってみましょう。

一般的に現実の物体に当たった光はそのまま反射する光と内部に入って行く光に分裂します。この際、反射する割合を表す数値がフレネルFと呼ばれています。こうしてフレネルの割合で反射した光はいわゆるSpecular,内部に入った光はDiffuseとしてレンダリングの分野で扱われています。すなわち、DiffuseとSpecularを合わせたBSDFというのは以下の様に書くことが出来ます。

\displaystyle{ f(\omega_i,\omega_o) = (1.0 - F) * f_{diffuse}(\omega_i,\omega_o) + F * f_{specular}(\omega_i,\omega_o) }

このようなシェーディングモデルはPhongの反射モデルと呼ばれています。DisneyBRDFを参考に1なら金属,0なら非金属になるmetallicというパラメーターを追加して以下のようなモデルを考えることにします。

\displaystyle{ f(\omega_i,\omega_o) = (1.0 - F) * f_{diffuse}(\omega_i,\omega_o) * (1.0 - metallic) + F * f_{specular}(\omega_i,\omega_o) }

光の反射において金属と非金属は大きく異なっており、金属ではDiffuse光は一切存在しないということになっています。また、F0の値も異なっており、非金属ではおおかた0.04程度ですが、金属であれば0.8といった高い値を自由に取ることが出来ます。なのでDisneyBRDFだとこのような形でF0を定義しています。

F0 = mix(vec3(0.04),color,metallic)

BSDFはこのように考えることが出来ましたが、問題は重点的サンプリングです。LambertとGGXそれぞれの重点的サンプリングは既に分かっていますが、それを組み合わせたBSDFに対してはどのような重点的サンプリングを考えるべきでしょうか？例えばLambertのサンプリング方法を採用したとすると、うまいことDiffuseの部分はサンプリングできますが一方でSpecularは全く効率的にサンプリングすることはできません。だからと言ってGGXを取ったとしても、今度は逆にDiffuseがうまく出なくなってしまいます。つまるところどちらのサンプリング方法も重要であるというわけです。

このような複数のサンプリング方法がある時、それらをうまいこと組み合わせる手法としてMultiple Importance Sampling(MIS)というものがあり、その中でOne SampleモデルというのがBSDFではよく用いられています。詳しい説明についてはyumchawiz氏の記事に詳しく書かれています。

blog.teastat.uk

詳細は割愛するとして、やり方だけ説明すると2つのサンプリングにそれぞれ確率 c_d,c_sを考え、それに従って乱数でどちらのサンプリングを使うかをランダムに選択します。この時、確率については自由に決めて問題がありません。今回は Metallicが0の時はお互い半々になるようにして、Metallicが1の時はGGXのみを取るようにウェイトを考えて確率を決めます。

 float dif_weight = (1.0 - info.metallic);
    float spec_weight = 1.0;
    float sum_weight = dif_weight + spec_weight;

    float cd = dif_weight / sum_weight;
    float cs = spec_weight / sum_weight;

選択したサンプリング手法に従って方向サンプリングをして出射方向ベクトル\omega_iを得ますが、この時に選択されなかったサンプリングでもし\omega_iを得た時のPDFも同時に計算する必要があります。そのため、\omega_iからPDFを逆算できる関数を用意する必要があります。

float pdfLambert(vec3 wi){
    return wi.y / PI;
}

float pdfGGX(vec3 wo,vec3 wm,float alpha){
    return 0.25f * GGX_D(wm,alpha) * dot(wo, wm) / (dot(wm, wo) * abs(wo.y)* (1.0f + GGX_Lambda(wo,alpha)));
}

これを用いてサンプリングの選択部分の処理はこのような形で書くことが出来ます。

 vec3 wm;
    float pdf_diffuse;
    float pdf_specular;
    if(rnd1() < cd){
                //Diffuse選択
        wi = cosineSampling(rnd2(),pdf_diffuse);    
        wm = normalize(wi + wo);    
        pdf_specular = pdfGGX(wo,wm,alpha);
    }
    else{
                //Specular選択
        wm = sampleVisibleNormal(rnd2(),wo,alpha);
        wi = reflect(-wo,wm);
        pdf_specular = pdfGGX(wo,wm,alpha);
        pdf_diffuse = pdfLambert(wi);
    }

何でわざわざ選択されなかった方のPDFを計算したかというとOne SampleモデルではPDFが2つのPDFの合成として定義されているためです。このような合成を考えることでうまいこと確率的サンプリングをアンバイアス、かつ効率的にできるようにしています（バランスヒューリスティックというウェイトの付け方)。

\displaystyle{PDF = c_d * pdf_{diffuse} + c_s * pdf_{specular} }

こうしてサンプリングをすることが出来ました。その後は単純にBSDFの評価をすることで複雑なBSDFをパストレで扱うことが出来るようになります。

vec3 BSDF(vec3 wo,inout vec3 wi,SurfaceInfo info){

    float alpha = clamp(info.roughness * info.roughness,0.001,1.0);
    vec3 F0 = mix(vec3(0.04),info.color,info.metallic); 

    float dif_weight = (1.0 - info.metallic);
    float spec_weight = 1.0;
    float sum_weight = dif_weight + spec_weight;

    float cd = dif_weight / sum_weight;
    float cs = spec_weight / sum_weight;

    vec3 wm;
    float pdf_diffuse;
    float pdf_specular;
    if(rnd1() < cd){
        wi = cosineSampling(rnd2(),pdf_diffuse);    
        wm = normalize(wi + wo);    
        pdf_specular = pdfGGX(wo,wm,alpha);
    }
    else{
        wm = sampleVisibleNormal(rnd2(),wo,alpha);
        wi = reflect(-wo,wm);
        pdf_specular = pdfGGX(wo,wm,alpha);
        pdf_diffuse = pdfLambert(wi);
    }
    
    if(wi.y < 0.0){
        return vec3(0.0);
    }
    
    float pdf = cd * pdf_diffuse + cs * pdf_specular;
    
    vec3 F = shlickFresnel(F0,dot(wi,wm));
    float D = GGX_D(wm,alpha);
    float G = GGX_G2_HeightCorrelated(wi,wo,alpha); 

    vec3 lambert = info.color / PI;
    
    vec3 bsdf = lambert * (vec3(1.0 ) - F) * (1.0 - info.metallic) + F * D * G / (4.0 * wo.y * wi.y);
    
    float cosine = wi.y;

    return bsdf * cosine / pdf;
}

ここまでやればつるつるした陶器などのかなり様々な質感を作ることが出来るようになります。ここまでのコードでも2.47KBになるのでまだまだ色々実装することが出来ます。

実装例

おわり

記事を読んでいただきありがとうございました。この記事で紹介している手法はあくまでパストレを愚直に書いたものなので必ずしもexegfxに適した書き方&手法ではありません。正直絵が良くてコードが少なければ正義なのでいい手法があったら教えてください。

私がexegfxを知ったのはSESSIONが初めで、Combined Graphics Compoにそうゆう部門があるんだ～という感じであまり詳しく調べてなかったのですが、0b5vr氏の20230419_fillercubeを見て「え！これが4Kで生成できるの！」とかなり衝撃を受けました。そしてその作品がパストレでやっていると聞き、普段やってることが使えて面白そうなカテゴリーだと感じいつかやってみたいな～と思っていました。

www.shadertoy.com

そんな中TDF16msが開催ということになったのでよしじゃあ今回はexegfxを作ってみようと思い立ち、exegfxの作品を提出させて頂きました。割と記事にあることをそのままやってSDFをこねこねして作りました。

TDF16msで提出した作品 「Bon appetite」

exegfxをやってみて思ったのですが、他のカテゴリーに比べるとかなりやりやすいな～と思いました（多分私がレイトレやってるのもあると思いますが）。Demoとかのカテゴリーではグラフィックス以外にも音楽とかアニメーションとかの要素が求められますが、exegfxは純粋なグラフィックス勝負という感じなため一つのことに集中して制作できるというのはとてもいいなと感じました。パストレさえできれば割とプリミティブなシーンでもかなりいい感じの絵になりますし、SDFをぐにゃぐにゃするだけでも中々見ることが出来ない面白い画像を作ることが出来ます。環境設定も楽だし

なので、個人的には初心者にとって敷居が低いカテゴリーなんじゃないかなと思っています。もちろんパストレが必須というわけではなく良くやるリアルタイムな手法で絵を作って、DOFやアンチエイリアシングを綺麗にするだけでも綺麗な作品を作ることが出来ます。普段ならやらない重いシェーディングとかSDFを試すこともできますし、結構楽しい分野だと思います。

Revisionでは結構exegfxの作品数は多いですが、日本だとexegfxをやっている方が少ないので是非この記事を機会に始める人が出てほしいな～と思っています。