簡単なRayMarching（レイマーチング）を実装してみた。

「ShaderToyでカメラのテストをやってみた。」に続き、半分くらい本格的なRayMarching（レイマーチング）を使用した描画を実装してみました。面白かった！というか、数学の重要性が切実に感じられた気がします。まだ複雑なモデルは実装してみませんでしたが、時間があればどんどんと複雑なものまで実装してレンダリングをやってみたいですね。

分析

レイマーチングってなに？とする方は「[GLSL] レイマーチング入門 vol.1」をご覧ください。

追加された共通関数

/// @brief Rotate position (pos) with radian on counter-clockwise.
vec3 RotateY(vec3 pos, float radian)
{
    return vec3(
        cos(radian) * pos.x + sin(radian) * pos.z, 
        pos.y,
        -sin(radian) * pos.x + cos(radian) * pos.z);
}

RotateY関数は原点を軸として、Y軸・反時計回りでposを回転します。一つ軸だけで回転しますので簡単に三角関数を使って実装しました。
実はこの関数は「ShaderToyで線が引きたくて実装・メモを書いてみた。」のRotate関数とかぶっちゃいます。

構造体ら

/// @brief Directional light structure.
struct DLight { vec3 mDir; vec3 mCol; };

/// @brief Plane structure.
struct DPlane { vec3 mCenterPos; vec3 mNormal; vec3 mCol; };

/// @brief Sphere structure.
struct DSphere { vec3 mPos; vec3 mCol; float mRadius; };

レイマーチングを使用して任意モデルをレンダリングするにはSDF（Signed Distance Function）を使って、値によってレンダリングフェイズに転換したりして描画します。ですので、SDFに必要となるポイント、モデルごとの要素を宣言しました。
DLightはDirectional Light（直射日光）を構造化したものです。直射日光は環境にすべて平行になりますので、mDirだけで大丈夫です。
DPlaneは $n_x(x - a_x) + n_y(y - a_y) + n_z(z - a_z) = 0$ の数式モデルに従うため、中心となる点・法線を宣言しました。DPlaneの長所は法線をを算出するのが定数時間に終わる点です。
DSphereは中心となる点・直径だけを宣言します。これで光線が衝突したところに適した法線を求められます。

mainImage

基本的なUVの座標・カメラの位置およびrayDirの算出方法は同じです。
スクリーンのUVの座標はスクリーン領域にて中心が $[0, 0$ ]となります。スクリーンの右上の端が必ず $[0.5, 0.5$ ]であるわけではありません。

DPlane far;
far.mCenterPos = camera + 10.0 * forward;
far.mNormal = -forward;

f:id:neuliliilli:20190525225144p:plain — http://help.agi.com/AGIComponentsJava/html/GraphicsCameraViewFrustum.htm

普通のリアルタイムレンダリングのカメラが持っているfarを実装します。near部分は仮想のスクリーンが代わりに行ってくれます。放った光線がスクリーン領域のfar面を超えると、その部分はモデルを描画せずにそのまま初期値を描画します。
スクリーン領域（スペース）でのfar面となりますので、normalは-forwardになります。10.0が境界値となります。

// Setting up light
DLight dirLight;
dirLight.mDir = normalize(vec3(0, 1, 1));
dirLight.mCol = vec3(1, 1, 1);

// Setting up background color.
vec3 result = vec3(0);
result = vec3(0., 56. / 255., 91. / 255.);

直射日光を設定します。日光の色は白（当然ですけど）、そして方向は右上で、zにちょっと出してる方向になります。
resultはスクリーンのテクセル（Texel）に最終に描画する予定の色を指します。初期値として青色を付けます。

// Setting up structure (model)
DPlane plane;
plane.mCenterPos = vec3(0, -1, 0);
plane.mNormal = vec3(0, 1, 0);
plane.mCol = vec3(1, .5, 0);

// Setting up sphere1, sphere2...

DSphere sphere3;
sphere3.mPos = RotateY(vec3(0, 0, 1), 1.17);
sphere3.mCol = vec3(1, 0, 0);
sphere3.mRadius = 0.125;

描画するモデルの値を設定します。Planeは下の床部分になります。Sphereは３つ設定して、丸型で配置するようにRotateYを使用しました。

Shooting Ray・検出・描画の部分

// Raycasting & marching & rendering.
for (int i = 0; i < 64 && sdfPlane(camera, far) > 0.5f; ++i)
{
    float shortestDist = sdfPlane(camera, far);
    shortestDist = min(shortestDist, sdfSphere(camera, sphere));
    if (shortestDist < 0.01)
    {
        vec3 normal = GetNormalSphere(camera, rayDir, sphere);
        vec3 qLambert = clamp(dot(dirLight.mDir, normal) * sphere.mCol * dirLight.mCol, 0., 1.) * 0.75 + 0.25;
        result = qLambert;
        break;
    }
    
    // ...
    
    shortestDist = min(shortestDist, sdfPlane(camera, plane));
    if (shortestDist < 0.01)
    {
        vec3 normal = GetNormalPlane(camera, rayDir, plane);
        vec3 qLambert = clamp(dot(dirLight.mDir, normal) * plane.mCol * dirLight.mCol, 0., 1.) * 0.75 + 0.25f;
        result = qLambert;
        break;
    }
    
    camera += shortestDist * rayDir;
}

if、break分岐が入ります。率直にGPUで動くコードでそういう分岐のコードを入れるのはGPUの構造上良くないです。ですけど仕方ありませんからifを入れて読みやすくします。
forを使用して、64回ループ・光線がfar面に到達するまでshortestDistを $t$ にし、rayDir方向へ進行させます。shortestDistの初期値は現在の光線の位置からfarまでの最短距離です。
このコードでは、shortestDistが一定数値以下の場合にデフォルト値で光線を進ませるようにしませんでした。これがないと、光線が進まなくなり、farに到達せずに済むことがあります。これは後ほどコードを追加します。

shortestDist = min(shortestDist, sdfSphere(camera, sphere));
if (shortestDist < 0.01)
{
    vec3 normal = GetNormalSphere(camera, rayDir, sphere);
    vec3 qLambert = clamp(dot(dirLight.mDir, normal) * sphere.mCol * dirLight.mCol, 0., 1.) * 0.75 + 0.25;
    result = qLambert;
    break;
}

一つのモデルごとに上のようなコードを通します。SDF数値の比較・条件が当たるときにレンダリング（イテレーションの終了）をします。
shortestDistはモデルを通す度に最小の値を持たなければなりません。この値を使用して、出来るだけ効率的に光線を進ませます。sdfSphere、sdfPlane関数は別途で説明します。
shortestDistが一定数値より下だったら、光線が今通しているモデルのあるフラグメントに当たってたと仮定します。そして当たった部分の法線を算出します。これを使用して拡散モデルの色を求め、記録してからイテレーションを終えます。
- 拡散モデルは3/4Quarter Lambertを使用しました。（任意のカスタムモデルです。普通はHalf Lambertを使いますが）

SDF関数・法線算出関数

/// @brief Get signed-distance funciton weight value of plane.
float sdfPlane(vec3 pos, DPlane plane)
{
    // Calculus Early Transcendentals 7th p823.
    float denom = length(plane.mNormal);
    float numem = abs(dot(plane.mNormal, pos - plane.mCenterPos));
    return numem / denom;
}

/// @brief Get signed-distance function weight value of sphere.
/// 0 is isosurface, and negative value is in sphere.
float sdfSphere(vec3 pos, DSphere sphere)
{
    return distance(pos, sphere.mPos) - sphere.mRadius;
}

sdfPlaneはDPlane（面）のSDF数値を、sdfSphereはDSphere（球体）のSDF数値を返します。sdfSphereの部分はLINKを参照しました。sdfPlaneは面とポイントの距離を求める数式を実装しました。
sdfSphereは負の数値を持てます。だとしても $\le 0$ になればレンダリング範囲に入るのでただ0にすることとは同じ効果を持ちます。
このsdf関数で進行する最短値を求めます。

/// @brief Get normal vector of plane.
vec3 GetNormalPlane(vec3 pos, vec3 dir, DPlane plane) { return plane.mNormal; }

/// @brief Get proper normal vector from intersection of ray with sphere,
/// when ray is intersected with sphere.
vec3 GetNormalSphere(vec3 pos, vec3 dir, DSphere sphere)
{
    float longHypot  = distance(pos, sphere.mPos);
    float shortDist2 = pow(DistLineToPoint(pos, dir, sphere.mPos), 2.0);
    
    float t = sqrt(pow(longHypot, 2.0) - shortDist2);
    float minusT = sqrt(pow(sphere.mRadius, 2.0) - shortDist2);
    t -= minusT;
    
    return normalize((pos + t * dir) - sphere.mPos);
}

DPlaneの法線は、DPlane自体に法線の情報があるからこれをリターンすることで済みます。

f:id:neuliliilli:20190526003042p:plain

DSphereの法線を求むことはちょっとトリッキーです。上の図のように、tを最終値を求めて、posをdir方向へ移動させた新しい位置を、球体の中心点からの方向を計算して返します。
実はこの方法ではなくて、もっと汎用性があるノーマル取り方があります。（下）

法線の汎用的算出方法

iquilezles.org

vec3 getNormal(vec3 pos, float size)
{
    float ep = 0.0001;
    return normalize(vec3(
            dist_func(pos, size) - dist_func(vec3(pos.x - ep, pos.y, pos.z), size),
            dist_func(pos, size) - dist_func(vec3(pos.x, pos.y - ep, pos.z), size),
            dist_func(pos, size) - dist_func(vec3(pos.x, pos.y, pos.z - ep), size)
        ));
}