OpenGL3.3でGPUパストレーサーを実装する
2020/12/23
#c++ #glsl #opengl #raytracing
コーネルボックスをレンダリングできる比較的簡単なGPUパストレーサーを実装してみました.

この記事はレイトレ Advent Calendar 2020の23日目の記事です.

今回は最近取り組んでいたOpenGL 3.3でのフラグメントシェーダーを使ったGPUパストレーサーの実装について紹介したいと思います。本記事では純粋なパストレーシングを実装し、コーネルボックスをレンダリングできるところまでの説明を行っています。本当はBVHを組み込んでSponzaとか出すところまで紹介したかったのですが、実装が間に合わなかったのでそれについてはまた別の機会で記事を書きたいと思います。

GPUレイトレーシングと言うと最近はDirectX Ray TracingやVulkan Ray Tracingが有名ですが、ここではあえて古い(?)OpenGL 3.3での実装を行ってみました. 何故かと言うと、WebGLを使ってWebサイト上にパストレーサーのデモを配置してみたかったのと、制約の中でどこまで出来るのかに興味があったので、あえてその道を選んでみました。

今回紹介する実装はGithubで公開しているので, 興味のある方は是非ビルドして動かしてみてください。コードをこの記事の内容と合わせるにはrt-adventcalender-2020タグを参照してください。

この記事では実装上の重要なポイントだけを抜き出して説明しています。全体像はGithubのコードを見て頂くと分かりやすいかと思います。また、パストレーシングの理論や計算方法などは説明していないので、それらについては別の資料を参照して頂けると幸いです。

OpenGL 3.3でレイトレ?

当然ですがOpenGLはラスタライズにより絵を出すので、レイトレで絵を出すための自然な枠組みは存在しません. レイトレで絵を出すためにはOpenGLの頂点シェーダー→フラグメントシェーダーの流れの中にレイトレを組み込む必要があります。

そこで使うのが画面一杯に平面を表示して、フラグメントシェーダーでレイトレにより色を計算して表示するトリックです。レイマーチングを行ったことがある方にはお馴染みの方法だと思います。

つまりレイトレして3Dっぽく見えるものは実際にはただ平面に貼り付けられた画像でしかないということです。この平面を常にスクリーンを覆うように表示することであたかもレイトレで計算したシーンが目の前にあるように見せることができます。

この流れでは頂点シェーダーは頂点情報をそのままフラグメントシェーダーに渡すだけなので非常に短いコードになります。一方でフラグメントシェーダーではレイトレを実装したり表示周りを書いていったりするので、コードの大半はフラグメントシェーダーに集中します。

フラグメントシェーダーでのパストレの流れ

フラグメントシェーダーでのレイトレの流れは以下のようになります。

  1. gl_FragCoordを用いて画素に対応するレイを生成
  2. パストレーシングして色を計算(1つのサンプルに対応)
  3. サンプル蓄積用テクスチャに色を書き込む

つまりドローコールを一回呼ぶと、1つのサンプルに対応した色が計算され、それがテクスチャに加算されます。パストレーシングの結果を表示するためにはこのテクスチャの色をサンプル数で割ったものを画面に表示します。

この実装方法ではサンプル数の増加に伴って色が収束していく様子を見ることができます。

パストレーサーの実装

以降では実装上で重要と思われるポイントを抜き出して解説していきます。

レイトレで必要になるデータ型の定義

まずはレイトレで必要になるデータ型をGLSLに定義していきます。ここは人によって実装が分かれるところではあると思いますが、私は以下のようなデータ型を用意しました。

  • Ray: レイを表現する
  • IntersectInfo: レイとシーンの交差情報を表現する。交差位置, 法線, 物体IDなどを持つ。
  • Material: BRDFの種類, パラメーターを表現する。
  • Primitive: 一つの物体を表現する。物体の位置, 形状, マテリアルIDなどを持つ。

これらはGLSL中で構造体を使って表現します。

struct Ray {
    vec3 origin;    // 始点位置
    vec3 direction; // 方向ベクトル
};

struct IntersectInfo {
    float t;        // レイの始点から交差位置までの距離
    vec3 hitPos;    // 交差位置
    vec3 hitNormal; // 交差位置における法線
    vec3 dpdu;      // 交差位置における接ベクトル(tangent vector)
    vec3 dpdv;      // 交差位置における従法線(binormal vector)
    float u;        // テクスチャ座標
    float v;        // テクスチャ座標
    int primID;     // 交差した物体に対応する物体ID
};

struct Material {
    int brdf_type;  // BRDFの種類
    vec3 kd;        // Diffuse成分
    vec3 le;        // Emission成分(0以上なら光源になる)
};

struct Primitive {
    int id;               // 物体ID
    int type;             // 形状の種類
    vec3 center;          // 球の中心位置
    float radius;         // 球の半径
    vec3 leftCornerPoint; // 平面の左端の点
    vec3 up;              // 平面の上方向の長さベクトル
    vec3 right;           // 平面の右方向の長さベクトル
    int material_id;      // マテリアルID
};

Primitiveのtypeは物体の形状種別を表しています。今回は球(type=0)と平面(type=1)の2種類の形状を用意しました。center, radiusは球の形状情報, leftCornerPoint, up, rightは平面の形状情報に対応します。

この実装方法は物体の形状の種類が増えていった場合を考えると良くないです。解決策としてはShapeという物体形状を表現する構造体を別に用意し, Primitiveは対応するShapeのIDを保持するという方法が考えられます。

物体との交差計算の実装

先程定義した構造体を用いて物体との交差計算をGLSLに実装していきます。実装はCPUでやる場合とほとんど同じになります。交差計算の結果はIntersectInfoに格納し、交差したかどうかをboolで返すようにします。

bool intersectSphere(in vec3 center, in float radius, in Ray ray, out IntersectInfo info) {
    float b = dot(ray.origin - center, ray.direction);
    float len = length(ray.origin - center);
    float c = len*len - radius*radius;
    float D = b*b - c;
    if(D < 0.0) {
        return false;
    }

    float t0 = -b - sqrt(D);
    float t1 = -b + sqrt(D);
    float t = t0;
    if(t < RAY_TMIN || t > RAY_TMAX) {
        t = t1;
        if(t < RAY_TMIN || t > RAY_TMAX) { 
            return false;
        }
    }

    info.t = t;
    info.hitPos = ray.origin + t*ray.direction;

    vec3 r = info.hitPos - center;
    info.hitNormal = normalize(r);
    info.dpdu = normalize(vec3(-r.z, 0, r.x));

    float phi = atan2(r.z, r.x);
    if(phi < 0.0) {
        phi += 2.0 * PI;
    }
    float theta = acos(clamp(r.y / radius, -1.0, 1.0));
    info.dpdv = normalize(vec3(cos(phi) * r.y, -radius * sin(theta), sin(phi) * r.y));

    info.u = phi * 0.5 * PI_INV;
    info.v = theta * PI_INV;
    return true;
}

bool intersectPlane(in vec3 leftCornerPoint, in vec3 right, in vec3 up, in Ray ray, out IntersectInfo info) {
    vec3 normal = normalize(cross(right, up));
    vec3 center = leftCornerPoint + 0.5 * right + 0.5 * up;
    vec3 rightDir = normalize(right);
    float rightLength = length(right);
    vec3 upDir = normalize(up);
    float upLength = length(up);

    float t = -dot(ray.origin - center, normal) / dot(ray.direction, normal);
    if(t < RAY_TMIN || t > RAY_TMAX) {
        return false;
    }

    vec3 hitPos = ray.origin + t*ray.direction;
    float dx = dot(hitPos - leftCornerPoint, rightDir);
    float dy = dot(hitPos - leftCornerPoint, upDir);
    if(dx < 0.0 || dx > rightLength || dy < 0.0 || dy > upLength) {
        return false;
    }

    info.t = t;
    info.hitPos = hitPos;
    info.hitNormal = dot(-ray.direction, normal) > 0.0 ? normal : -normal;
    info.dpdu = rightDir;
    info.dpdv = upDir;
    info.u = dx / rightLength;
    info.v = dy / upLength;
    return true;
}

シーンとの交差計算

個別物体との交差計算が実装できたので、次はシーン全体との交差計算を考えていきます。そのためにはまずシーンをGLSL中でどのように表現するかについて見ていきます。

シーンの表現

今回の構造体の定義では、シーンの表現に必要なのはPrimitiveとMaterialの配列です。これらはUniform変数として定義し、CPU側からデータを渡せるようにしました。GLSL中に配列の各要素をハードコーディングすることも考えられますが、これは拡張性がないのでやめました。

構造体の配列をUniform変数として扱う場合はUniform Buffer Object(UBO)を使うとCPU側からデータを渡すのが楽になります。具体的にはCPU側で用意した構造体の配列をそのままGPU側に渡すといった事が出来るようになります。今回はUBOを使ってPrimitive, Materialの配列を用意してみました。

以下、GLSL中でのUniform変数の定義です。

const int MAX_N_MATERIALS = 100;
layout(std140) uniform MaterialBlock {
  Material materials[MAX_N_MATERIALS];
};

const int MAX_N_PRIMITIVES = 100;
layout(std140) uniform PrimitiveBlock {
  Primitive primitives[MAX_N_PRIMITIVES];
};

CPU側からデータを送る方法については後ほど紹介します。

交差計算の実装

シーンをGLSL中で表現することができたので、シーンとの交差計算を実装していきます。今回はBVHを使わずに線形探索、つまり全ての物体との交差計算を試し、その中で最も交差位置が近いものを交差物体とするという処理を実装しました。

bool intersect_each(in Ray ray, in Primitive primitive, out IntersectInfo info) {
    switch(primitive.type) {
    // Sphere
    case 0:
        return intersectSphere(primitive.center, primitive.radius, ray, info);
    // Plane
    case 1:
        return intersectPlane(primitive.leftCornerPoint, primitive.right, primitive.up, ray, info);
    }
}

bool intersect(in Ray ray, out IntersectInfo info) {
    bool hit = false;
    info.t = RAY_TMAX;

    for(int i = 0; i < 16; ++i) {
        IntersectInfo temp;
        if(intersect_each(ray, primitives[i], temp)) {
            if(temp.t < info.t) {
                hit = true;
                temp.primID = primitives[i].id;
                info = temp;
            }
        }
    }

    return hit;
}

線形探索はintersect()が行っています。intersect_each()はPrimitiveの形状種類によって呼び出す交差計算の関数を変える補助関数です。

intersect()中ではシーンにある物体の数が16であると仮定して計算を行っています。16という数字はコーネルボックスのシーン中に存在する物体数に対応しています。拡張性を考えるとこの数字の部分はUniform変数として置き換えるのが良いでしょう。

レイの生成

次にレイを生成する部分、つまりカメラモデルの実装について見ていきます。今回はシンプルにピンホールカメラを実装してみました。

Ray rayGen(in vec2 uv, out float pdf) {
    vec3 pinholePos = camera.camPos + camera.a * camera.camForward;
    vec3 sensorPos = camera.camPos + uv.x * camera.camRight + uv.y * camera.camUp;

    Ray ray;
    ray.origin = camera.camPos;
    ray.direction = normalize(pinholePos - sensorPos);
    pdf = 1.0 / pow(dot(ray.direction, camera.camForward), 3.0);
    return ray;
}

rayGen()はスクリーンを覆う平面のテクスチャ座標を受け取り, それに対応するレイを返す関数です。pdfはレイ生成に対応する立体角確率密度です(これを考慮することでビネッティングが現れる)

camera.camPos, camera.aなどはカメラパラメーターを表現するUniform変数です。これもUBOを使ってCPU側からデータを渡せるように定義しました。

layout(std140) uniform CameraBlock {
  vec3 camPos;      // カメラの位置
  vec3 camForward;  // カメラの前方向
  vec3 camRight;    // カメラの右方向
  vec3 camUp;       // カメラの上方向
  float a;          // カメラからピンホールまでの距離
} camera;

乱数生成

パストレーシングを実装するためにはGLSL上で乱数を発生させることが必要になります。

GLSLで乱数生成というと以下のようなコードが有名です。

float rand(vec2 co){
    return fract(sin(dot(co.xy ,vec2(12.9898,78.233))) * 43758.5453);
}

これ以外にもvec3を受け取って乱数を生成するコードなどもあります。これらの関数はGLSLで手っ取り早く乱数を扱いたい時には便利なのですが、パストレーシングでは高次元の乱数が必要になるので扱いにくいです。

今回の実装ではこれらの便利な関数を使わずに、高速でそれなりに質の良い乱数を生成してくれるXORShift32をGLSL上で実装して乱数を生成することにしました。

乱数生成器の状態をいかにして保持するか?

XORShift32で乱数を生成するためには乱数生成器の状態を保持しておくことが必要になります。これは今回のように1つのドローコールが1つのサンプルに対応するという実装では悩ましい問題です。単純にGLSL上でグローバル変数として乱数生成器の状態を定義してしまうと、ドローコールを呼ぶ度に以前の乱数生成器の状態が破棄されてしまいます。

ここでは解決策として乱数生成器の状態を保持するuintテクスチャを用意し、各ドローコールでそのテクスチャから状態を取得し、計算が終わったら逆にテクスチャに状態を書き込むことで状態を保持することにしました。

GLSL中でuintテクスチャは以下のようなUniform変数で定義します。

uniform usampler2D stateTexture

CPU側では以下のようにしてuintテクスチャを用意しておきます。CPU側でuint32の乱数を生成し、その値でテクスチャを初期化しておきます。

    // setup RNG state texture
    glGenTextures(1, &stateTexture);
    glBindTexture(GL_TEXTURE_2D, stateTexture);
    std::vector<uint32_t> seed(width * height);
    std::random_device rnd_dev;
    std::mt19937 mt(rnd_dev());
    std::uniform_int_distribution<uint32_t> dist(
        1, std::numeric_limits<uint32_t>::max());
    for (unsigned int i = 0; i < seed.size(); ++i) {
      seed[i] = dist(mt);
    }
    glTexImage2D(GL_TEXTURE_2D, 0, GL_R32UI, width, height, 0, GL_RED_INTEGER,
                 GL_UNSIGNED_INT, seed.data());
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glBindTexture(GL_TEXTURE_2D, 0)

乱数の生成

乱数生成器の状態を保持するstateTextureを用いてXORShift32による乱数生成を実装していくと以下のようになります。

struct XORShift32_state {
    uint a;
};

// 乱数生成器の状態
XORShift32_state RNG_STATE;

// XORShift32の実装
uint xorshift32(inout XORShift32_state state) {
    uint x = state.a;
    x ^= x << 13u;
    x ^= x >> 17u;
    x ^= x << 5u;
    state.a = x;
    return x;
}

// [0, 1]の乱数を返す
float random() {
    return float(xorshift32(RNG_STATE)) * 2.3283064e-10;
}

// シード値をセットする
void setSeed(in vec2 uv) {
    RNG_STATE.a = texture(stateTexture, uv).x;
}

パストレの実装で使うのはrandom()setSeed()のみです。setSeed()はパストレの計算を行う前に呼ぶ必要があります(これについては後でコード例で見ていきます)

BRDFサンプリングの実装

乱数が生成できるようになったので、次にパストレで必要なBRDFサンプリング(物体に交差した後の次のレイの方向をBRDFに比例するような重点的サンプリングで決めること)を実装していきます。

方向サンプリングは交差位置のローカル座標系(x方向: 接ベクトル, y方向: 法線, z方向: 従法線)で行うことにします。こうすることで物体のワールド座標系での位置や回転を考えなくて良くなるので、実装がシンプルになります。

// フレネル反射率
float fresnel(in vec3 v, in float n1, in float n2) {
    float F0 = pow((n1 - n2) / (n1 + n2), 2.0);
    return F0 + (1.0 - F0) * pow(1.0 - abs(v.y), 5.0);
}

// コサインに比例した半球内での方向サンプリング
vec3 sampleCosineHemisphere(in float u, in float v, out float pdf) {
    float theta = 0.5 * acos(clamp(1.0 - 2.0 * u, -1.0, 1.0));
    float phi = 2.0 * PI * v;
    float y = cos(theta);
    pdf = y * PI_INV;
    return vec3(cos(phi) * sin(theta), y, sin(phi) * sin(theta));
}

// BRDFサンプリング
vec3 sampleBRDF(in vec3 wo, out vec3 wi, in int brdf_type, in vec3 kd, out float pdf) {
    switch(brdf_type) {
    // lambert
    case 0:
        wi = sampleCosineHemisphere(random(), random(), pdf);
        return kd * PI_INV;
        break;

    // mirror
    case 1:
        pdf = 1.0;
        wi = reflect(-wo, vec3(0, 1, 0));
        return kd / abs(wi.y);
        break;

    // glass
    case 2:
        pdf = 1.0;

        // set appropriate normal and ior
        vec3 n = vec3(0, 1, 0);
        float ior1 = 1.0;
        float ior2 = 1.5;
        if(wo.y < 0.0) {
            n = vec3(0, -1, 0);
            ior1 = 1.5;
            ior2 = 1.0;
        }
        float eta = ior1 / ior2;

        // fresnel
        float fr = fresnel(wo, ior1, ior2);

        // reflection
        if(random() < fr) {
            wi = reflect(-wo, n);
        }
        // refract
        else {
            wi = refract(-wo, n, eta);
            // total reflection
            if(wi == vec3(0)) {
                wi = reflect(-wo, n);
            }
        }

        return kd / abs(wi.y);
        break;
    }
}

sampleBRDFがBRDFサンプリングを行う関数です。wiにサンプリングした方向, pdfに立体角確率密度が格納され、RGBとしてBRDFの値を返すようにしています。

sampleBRDFの中ではbrdf_typeによってサンプリングの処理を切り替えます。brdf_type = 0はLambertian BRDFに対応します。この中ではsampleHemisphere()を使って重点的サンプリングによる方向生成を行っています。sampleHemisphere()の引数にrandom()を使うことでランダムな方向生成が行われます。

brdf_type=1は完全鏡面, brdf_type=2はガラスに対応します。

パストレーシングの実装

以上を持ってパストレーシングを実装する準備が整ったので、その実装について見ていきます。

// 純粋なパストレーシングで受け取ったレイの方向から来る放射輝度を計算する
vec3 computeRadiance(in Ray ray_in) {
    Ray ray = ray_in;

    const float russian_roulette_prob = 0.99; // ロシアンルーレットの確率
    vec3 color = vec3(0); // 放射輝度
    vec3 throughput = vec3(1); // スループット
    for(int i = 0; i < MAX_DEPTH; ++i) {
        // ロシアンルーレット
        if(random() >= russian_roulette_prob) {
            break;
        }
        throughput /= russian_roulette_prob;

        IntersectInfo info;
        // レイがシーン中の物体と交差した場合
        if(intersect(ray, info)) {
            Primitive hitPrimitive = primitives[info.primID];
            Material hitMaterial = materials[hitPrimitive.material_id];
            vec3 wo = -ray.direction;

            // レイの方向をローカル座標系に変換
            vec3 wo_local = worldToLocal(wo, info.dpdu, info.hitNormal, info.dpdv);

            // 光源に当たった場合
            if(any(greaterThan(hitMaterial.le, vec3(0)))) {
                color += throughput * hitMaterial.le;
                break;
            }

            // BRDFサンプリングにより次のレイの方向を生成
            float pdf;
            vec3 wi_local;
            vec3 brdf = sampleBRDF(wo_local, wi_local, hitMaterial.brdf_type, hitMaterial.kd, pdf);
            // prevent NaN
            if(pdf == 0.0) {
                break;
            }

            // サンプリングした方向をワールド座標系に変換
            vec3 wi = localToWorld(wi_local, info.dpdu, info.hitNormal, info.dpdv);

            // スループットを更新
            float cos_term = abs(wi_local.y);
            throughput *= brdf * cos_term / pdf;

            // 次のレイをセット
            ray = Ray(info.hitPos, wi);
        }
        // レイが空に飛んでいった場合
        else {
            color += throughput * vec3(0);
            break;
        }
    }

    return color;
}

コードが何をやっているかはコメントを見てもらえればと思います。CPUで実装する場合でもそうですが、純粋なパストレーシングの実装は周辺パーツをちゃんと用意してあげればこのように50行程度のコードで書くことができます。

コード中に出てくるlocalToWorld, worldToLocalは方向ベクトルをワールド座標系とローカル座標系の間で相互変換する関数ですが、実装は以下のようになっています。

vec3 worldToLocal(in vec3 v, in vec3 lx, in vec3 ly, in vec3 lz) {
    return vec3(dot(v, lx), dot(v, ly), dot(v, lz));
}
vec3 localToWorld(in vec3 v, in vec3 lx, in vec3 ly, in vec3 lz) {
    return vec3(dot(v, vec3(lx.x, ly.x, lz.x)), dot(v, vec3(lx.y, ly.y, lz.y)), dot(v, vec3(lx.z, ly.z, lz.z)));
}

サンプルの蓄積と画面表示

computeRadiance()は1つのサンプルに対応する放射輝度を計算するだけです。パストレーシングした結果を表示するためにはcomputeRadiance()で返ってくる放射輝度を蓄積し、それをサンプル数で割ったものを画面に表示する必要があります。

サンプルの蓄積

サンプルの蓄積はfloatのRGBテクスチャを用意して、そこに加算していくことで実現します。

uniform sampler2D accumTexture;

CPU側では以下のようにして蓄積用テクスチャを用意します。

    // setup accumulate texture
    glGenTextures(1, &accumTexture);
    glBindTexture(GL_TEXTURE_2D, accumTexture);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB32F, width, height, 0, GL_RGB,
                 GL_FLOAT, 0);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glBindTexture(GL_TEXTURE_2D, 0);

これを用いて、いよいよフラグメントシェーダーのmain()中の処理を書いていきます。

layout (location = 0) out vec3 color;
layout (location = 1) out uint state;

void main() {
    // 乱数生成器の状態をセット
    setSeed(texCoord);

    // 最初のレイを生成
    vec2 uv = (2.0*(gl_FragCoord.xy + vec2(random(), random())) - resolution) * resolutionYInv;
    uv.y = -uv.y;
    float pdf;
    Ray ray = rayGen(uv, pdf);
    float cos_term = dot(camera.camForward, ray.direction);

    // 放射輝度をサンプル蓄積用テクスチャに加算
    vec3 radiance = computeRadiance(ray) / pdf;
    color = texture(accumTexture, texCoord).xyz + radiance * cos_term;

    // 乱数生成器の状態をテクスチャに書き込み
    state = RNG_STATE.a;
}

処理の流れはコメントを見てもらえれば分かると思います。main()中ではサンプル蓄積用テクスチャ(colorに対応)と、乱数生成器の状態を保持するテクスチャ(stateに対応)への2つの書き込みを行っていることに注意してください。つまり Multi Render Target(MRT) を利用しています。

CPU側ではこれを行うためには Frame Buffer Object(FBO) のセットアップが必要です。以下のようにすれば良いです。

    // setup accumulate FBO
    glGenFramebuffers(1, &accumFBO);
    glBindFramebuffer(GL_FRAMEBUFFER, accumFBO);
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D,
                           accumTexture, 0);
    glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT1, GL_TEXTURE_2D,
                           stateTexture, 0);
    GLuint attachments[2] = {GL_COLOR_ATTACHMENT0, GL_COLOR_ATTACHMENT1};
    glDrawBuffers(2, attachments);
    glBindFramebuffer(GL_FRAMEBUFFER, 0)

画面表示

accumTextureに蓄積されている放射輝度をサンプル数で割ったものを画面に表示します。上で実装した放射輝度計算&サンプル蓄積用のフラグメントシェーダーとは別に、画面表示用のフラグメントシェーダーを用意し、そこに処理を書いていきます。

#version 330 core

uniform float samplesInv;
uniform sampler2D accumTexture;

in vec2 texCoord;
out vec4 fragColor;

void main() {
  vec3 color = texture(accumTexture, texCoord).xyz * samplesInv;
  fragColor = vec4(pow(color, vec3(0.4545)), 1.0);
}

samplesInvはサンプル数の逆数を格納するUniform変数です。accumTextureの値にsamplesInvをかけ合わせることでパストレーシングで計算した放射輝度の値が得られます。それをfragColorに出力することで画面表示を行っています。この際にpow(color, vec3(0.4545))と計算している部分はガンマ補正です。

CPU側での処理

ここまでの説明でGPU側でのパストレーシングの実装と画面表示の方法を見てきました。ここからは途中で飛ばしたUBOを用いてCPU側からデータを送る方法や、四角形の描画方法などについて見ていきます。

UBOを用いてUniform変数をセットアップする

まずはCPU側でデータを用意します。

GLSLでの構造体と全く同じような構造体をCPU側でも用意します。

struct alignas(16) Primitive {
  int id;                                 // 4
  int type;                               // 8
  alignas(16) glm::vec3 center;           // 24
  float radius;                           // 28
  alignas(16) glm::vec3 leftCornerPoint;  // 44
  alignas(16) glm::vec3 up;               // 60
  alignas(16) glm::vec3 right;            // 76
  int material_id;                        // 80
};

struct alignas(16) Material {
  int brdf_type;             // 4
  alignas(16) glm::vec3 kd;  // 20
  alignas(16) glm::vec3 le;  // 36
};

struct alignas(16) CameraBlock {
  alignas(16) glm::vec3 camPos;
  alignas(16) glm::vec3 camForward;
  alignas(16) glm::vec3 camRight;
  alignas(16) glm::vec3 camUp;
  float a;

  CameraBlock(const glm::vec3& camPos, const glm::vec3& camForward,
              const glm::vec3& camRight, const glm::vec3& camUp)
      : camPos(camPos),
        camForward(camForward),
        camRight(camRight),
        camUp(camUp) {}
};

ここでalignasを用いて構造体や構造体中の変数を16バイト境界にアラインメントしていることに注意してください。これはGLSLの方でUniform変数を定義する際にlayout(std140)を指定したため必要になっています。例えばvec3std140では16バイト境界にアラインメントされることになっています。詳細については公式ドキュメントを参照してください。

以上で定義した構造体を用いて、CPU側でシーンデータを用意します。

std::vector<Primitive> primitives; // これにPrimitiveを追加していく
std::vector<Material> materials;   // これにMaterialを追加していく
CameraBlock params;                // これにカメラパラメーターを突っ込む

次にUBOを作成します。

    // setup UBO
    glGenBuffers(1, &cameraUBO);
    glBindBuffer(GL_UNIFORM_BUFFER, cameraUBO);
    glBufferData(GL_UNIFORM_BUFFER, sizeof(CameraBlock), &params,
                 GL_DYNAMIC_DRAW);
    glBindBuffer(GL_UNIFORM_BUFFER, 0);

    glGenBuffers(1, &materialUBO);
    glBindBuffer(GL_UNIFORM_BUFFER, materialUBO);
    glBufferData(GL_UNIFORM_BUFFER, sizeof(Material) * 100,
                 materials.data(), GL_STATIC_DRAW);
    glBindBuffer(GL_UNIFORM_BUFFER, 0);

    glGenBuffers(1, &primitiveUBO);
    glBindBuffer(GL_UNIFORM_BUFFER, primitiveUBO);
    glBufferData(GL_UNIFORM_BUFFER, sizeof(Primitive) * 100,
                 primitives.data(), GL_STATIC_DRAW);
    glBindBuffer(GL_UNIFORM_BUFFER, 0);

UBOの作成はVBOやEBOなどの他のバッファの作成とほとんど同じ流れです。UBOを作成した後はglBindBufferBaseを用いて作成したUBOをGLSL側のバインディングポイントに対応させます。

    glBindBufferBase(GL_UNIFORM_BUFFER, 1, cameraUBO);
    glBindBufferBase(GL_UNIFORM_BUFFER, 2, materialUBO);
    glBindBufferBase(GL_UNIFORM_BUFFER, 3, primitiveUBO);

さらに各シェーダーのUniform変数を先程指定したバインディングポイントに対応させます。

// programはシェーダープログラム
// block_nameはUniform変数名が入る
// indexはバインディングポイント
const GLuint index = glGetUniformBlockIndex(program, block_name.c_str());
glUniformBlockBinding(program, index, binding_number);

このようにUBOを用いると構造体や構造体の配列などのデータをglBufferDataを使って簡単にGPU側に送ることができます。UBOを使わない場合はprimitives[0].type = 0などのように構造体の各要素、配列の各要素を一つ一つ指定していかなければならないのでとても不便です。

また、複数のシェーダーで同じUniform変数を扱う場合には、一度UBOを作ってしまえばバインディングポイントを共有することで簡単に同じUniform変数を扱うことができます。

描画

ここではパストレーシングを行うシェーダーを呼び出して画面表示を行う流れについて見ていきます。

頂点情報の用意

四角形の頂点情報を用意します。 頂点情報として頂点座標とテクスチャ座標を用意しておきます。VBOとEBO作ってデータ送ってそれをVAOにリンクして・・・という流れです。

    // setup VAO
    glGenVertexArrays(1, &VAO);
    glBindVertexArray(VAO);

    // setup VBO;
    // positions and texture coords
    GLfloat vertices[] = {-1.0f, -1.0f, 0.0f, 0.0f, 0.0f, 1.0f, -1.0f,
                          0.0f,  1.0f,  0.0f, 1.0f, 1.0f, 0.0f, 1.0f,
                          1.0f,  -1.0f, 1.0f, 0.0f, 0.0f, 1.0f};
    glGenBuffers(1, &VBO);
    glBindBuffer(GL_ARRAY_BUFFER, VBO);
    glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW);

    // setup EBO;
    GLuint indices[] = {0, 1, 2, 2, 3, 0};
    glGenBuffers(1, &EBO);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, EBO);
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(indices), indices,
                 GL_STATIC_DRAW);

    // position attribute
    glEnableVertexAttribArray(0);
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 5 * sizeof(GLfloat),
                          (GLvoid*)0);
    glEnableVertexAttribArray(1);
    glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 5 * sizeof(GLfloat),
                          (GLvoid*)(3 * sizeof(float)));

    // unbind VAO, VBO, EBO
    glBindVertexArray(0);
    glBindBuffer(GL_ARRAY_BUFFER, 0);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0)

描画

長くなりましたが、最後にCPU側での描画周りです。四角形の頂点情報を持つ Vertex Array Object(VAO) をバインドし、まずは放射輝度の計算&サンプル蓄積用のシェーダーによる描画、次にその結果を画面表示するシェーダーを呼び出します。

// 放射輝度の計算とサンプル蓄積
// Frame Buffer Objectをバインドしてからシェーダーを呼び出す
glBindFramebuffer(GL_FRAMEBUFFER, accumFBO);
glUseProgram(pt_shader);

glBindVertexArray(VAO);
glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, 0);
glBindVertexArray(0);

glUseProgram(0);
glBindFramebuffer(GL_FRAMEBUFFER, 0);

// サンプル数の更新
samples++;

// 画面表示
// サンプル数の逆数の更新
glUseProgram(output_shader);
const GLuint location = glGetUniformLocation(output_shader, "samplesInv");
glUniform1f(location, 1.0f / samples);

glBindVertexArray(VAO);
glDrawElements(GL_TRIANGLES, 6, GL_UNSIGNED_INT, 0);
glBindVertexArray(0);

glUseProgram(0);

結果

動かすと以下のようになります。GUIもつけてインタラクティブな操作ができるようにしました。

動作している様子は動画でyoutubeにアップロードしてみました。

https://youtu.be/-dmQk2q3FTo

CPUで計算させる場合と比べて圧倒的に計算が早くて感動しました。じゅわっと収束する様子が見ていて気持ち良いです。

最後に

ここまでお読みくださりありがとうございました。何か気になる点や、ここはこうした方が良いよ!という点があったら是非コメント欄やDMなどで教えて頂けると嬉しいです。

参考資料