この記事はTUT Advent Calendar 2022 123日目の記事です。

はじめに

初めまして、sh00t(.er)と申します。 12月となり、今年もデスクトップPCが暖房器具として実用的になる時期になりました。 (前年も同じこと言ってる)

PCの排気で暖まるにはGPUに負荷をかけるのが一番ですが、レイマーチングというGPUに負荷をかけるのにうってつけの方法があるので紹介します。

レイマーチング

レイマーチングとは

コンピュータで画を出すには普通「ラスタライズ」と呼ばれる手法で行っています。これは3D空間上に頂点と辺で構成されたポリゴンを配置し、その後カメラなどの情報から座標変換されたポリゴンを画面に表示するために塗る必要があるピクセルを決定する処理などのいくつかの工程を経て画面が描写されます。ラスタライズ法ではポリゴン数が多いほど綺麗に3Dモデルが表示されます。

映画などの映像製作ではより綺麗な映像を作るために「レイトレーシング」と呼ばれる方法が用いられることがあります。現実では光源からの光が様々な経路を辿り最終的に目に入ることで物体を見ています。レイトレーシングでは現実の逆で、カメラ(目にあたる)から光線(レイ)を飛ばして追跡(トレース)して、空間にある物体やライトの情報からカメラに入る光をシミュレートする方法です。レイは全てのピクセルの数飛ばすので、1920*1080の解像度なら2,073,600回行われます。

「レイマーチング」はレイトレーシングの種類で、カメラから飛ばしたレイを進めていく手法です。その中でも空間内の物体への最短距離の分だけ進めていく方法を「スフィアトレーシング」といいます。今回はスフィアトレーシングを用いて描写してみます。

色々書きましたがなんかいい感じの絵を作れる手法程度の理解でいいと思います。

レイマーチング入門

フラグメントシェーダにレイマーチングのコードを記述していきます。フラグメントシェーダは簡単に言うと画面上の各画素に対して色を決定する処理です。 ShaderToyというフラグメントシェーダを書いて実行できるサイトがあるので利用します。右上のNewから新規作成。最初は下画像と同じ画面になっていると思います。最終的な色はfragColorというvec4型の変数に(R,G,B,A)の順番で入れる必要があります。 800*450の各画素についてこのコードが実行されています。これを下のコードに変更してみましょう。

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    // Normalized pixel coordinates (from 0 to 1)
    vec2 uv = fragCoord/iResolution.xy;

    vec3 col = vec3(uv, 0.0);

    // Output to screen
    fragColor = vec4(col,1.0);
}

コードの通りfragColorには(col, 0.0)という変数が、colには(uv, 0.0)が格納されています。 'uv'変数はvec2型で'uv.x'、uv.yとアクセスできるので、 fragColor‘には(uv.x, uv.y, 0.0, 1.0)が格納されていることになります。 [f:id:sh00t:20221213001314p:plain] 実行すると左下:黒、右下:赤、左上:緑、右上:黄色の画面が出力されます。 mainImageの最初で処理中の画素を表すfragCoordを解像度iResolution.xyで割った値を変数uvに入れています。左下が原点で0なので、uvには現在処理中の座標が左下が0、右と上に行くほど1に近づく値が入っています。fragColorは(R, G, B, A)なので、uv.xが1に近いほど赤く、uv.y`が1に近いほど緑となります。

uvの正規化

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

    vec3 col = vec3(uv, 0.0);

    fragColor = vec4(col,1.0);
}

変わっているのはuvに代入しているところだけですが、これを実行すると中央から左下:黒、右下:赤、左上:緑、右上:黄色の画像が出力されます。このコードではuvの値を中心(0,0)にして正規化して下の画像のように値が入っています。これからの処理で扱いやすいので正規化しています。

次からいよいよレイマーチングに入ります。

レイマーチングで球の描写

// 球の距離関数
float sdSphere(vec3 p, float r)
{
    return length(p) - r;
}

// シーンの距離関数
float map(vec3 p )
{
    float d = sdSphere(p, 1.0);
    return d;
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (2.0 * fragCoord - iResolution.xy) / min(iResolution.x, iResolution.y);
    
    // レイを飛ばす位置のz座標
    float screenZ = -10.0;
    // カメラの位置
    vec3 cameraOrigin = vec3(0.0, 0.0, -20.0);
    // 飛ばすレイの方向
    vec3 rayDirection = normalize(vec3(uv, screenZ) - cameraOrigin);
    
    // 最短距離
    float dist = 0.0;
    float t = 0.0;
    // マーチングループ
    for(int i = 0; i < 32; i++)
    {
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        dist = map(rayPosition);
        
        t += dist;
        if(dist < 0.0001)
        {
            break;
        }
        
    }
    
    vec3 col = vec3(0.0);
    if(dist < 0.0001)
    {
        col += 1.0;
    }
    
    fragColor = vec4(col,1.0);
}

画面の中心に白い円が表示されますが、処理上は球に当たって色が塗られています。まず、カメラの位置とレイを飛ばす方向を決めなければなりません。カメラの位置はcameraOrigin、例を飛ばす方向はrayDirectionに格納されています。 screenZはuvによる仮想的なスクリーンの位置のz座標です。スクリーンの位置(uv.x, uv.y, screenZ)からカメラの座標を引いて正規化することで、レイを飛ばす方向ベクトルを求めています。今回はカメラの座標を(0, 0, -20)としているので、カメラからz座標で10離れた場所に仮想的なスクリーンがあり、現在処理しているシェーダの座標に対応してuvの値が変化してレイの方向も変化します。下画像はxz平面のカメラとスクリーンと物体の位置関係です。 for文の中ではレイを進める処理をしています。 vec3 rayPosition = cameraOrigin + rayDirection * t;の部分ですが、これはレイの先端の座標を求めています。 dist = map(rayPosition);ではレイの先端から最も近い物体の面までの距離を求めています。座標pから半径rの球の表面までの距離はlength(p) - rで求まります。やっていることは $\sqrt{x^{2}+y^{2}+z^{2}}-r=0$ と同じです。最短距離を求めたらレイの方向にかけるtを求めた最短距離分進めます。もし最短距離distが0.0001より小さかったらレイと物体が衝突しているとしてループを抜けます。ループを抜けた後の処理ですが、もしレイと物体が衝突していたら現在処理中のピクセルの色を白にします。

次は立体的に表示するために簡単にライティングしてみます。

ライティング

#define EPSILON 0.0001

// 球の距離関数
float sdSphere(vec3 p, float r)
{
    return length(p) - r;
}

// シーンの距離関数
float map(vec3 p )
{
    float d = sdSphere(p, 1.0);
    return d;
}

// 法線を求める関数
vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        map(vec3(p.x + EPSILON, p.y, p.z)) - map(vec3(p.x - EPSILON, p.y, p.z)),
        map(vec3(p.x, p.y + EPSILON, p.z)) - map(vec3(p.x, p.y - EPSILON, p.z)),
        map(vec3(p.x, p.y, p.z  + EPSILON)) - map(vec3(p.x, p.y, p.z - EPSILON))
    ));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (2.0 * fragCoord - iResolution.xy) / min(iResolution.x, iResolution.y);
    
    // レイを飛ばす位置のz座標
    float screenZ = -10.0;
    // カメラの位置
    vec3 cameraOrigin = vec3(0.0, 0.0, -20.0);
    // 飛ばすレイの方向
    vec3 rayDirection = normalize(vec3(uv, screenZ) - cameraOrigin);
    
    // 最短距離
    float dist = 0.0;
    float t = 0.0;
    // マーチングループ
    for(int i = 0; i < 32; i++)
    {
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        dist = map(rayPosition);
        
        t += dist;
        if(dist < 0.0001)
        {
            break;
        }
        
    }
    
    vec3 col = vec3(0.0);
    if(dist < 0.0001)
    {
        // 物体表面のレイ
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        // 物体表面の法線を求める
        vec3 normal = estimateNormal(rayPosition);
        // Lambert拡散反射
        vec3 lightDir = normalize(vec3(1.0, 1.0, -1.0));
        float l = dot(normal, lightDir);
        col += clamp(l, 0.0, 1.0);
        col += 0.1;
    }
    
    fragColor = vec4(col,1.0);
}

前回と比べてestimateNormal関数の追加とif(dist < 0.0001){...}にライティング処理が追加されています。 estimateNormal関数ですが、x,y,z座標を微小量だけずらして物体との最短距離を求めています。呼び出す際は物体表面のレイの座標を渡すので偏微分を行っていることになります。自分もあまり詳しくないのですが距離関数は陰関数なため偏微分をすることで法線が求まります。なぜ陰関数の偏微分で法線が求められるのかは下のサイトなどで解説されています。 manabitimes.jp 法線を求めることでライティングの計算ができるようになります。今回はLambert拡散反射のライティングをします。Lambert拡散反射は光が強く当たっている表面ほど明るく光が当たっていない表面ほど暗くなる光の計算方法です。具体的には光源と法線のなす角が小さいほど光が強く当たっているとみなします。コード上ではライトのある方向lightDirを決め、法線との内積を取っています。内積は同じ向き程1に近づき、反対向きで-1となるので内積で明るさが計算できます。最後に環境光として全体に光が当たっているとするために0.1を足しています。環境光は光源からの光が様々な物体と反射して間接光として届く光で、まじめに計算するのは難しいので一律に明るさを足しています。ライティングにより立体的に球を表示することができました！

Boxの表示

今までは球だけでしたが、次はBoxの表示をしてみます。

#define EPSILON 0.001

// 球の距離関数
float sdSphere(vec3 p, float r)
{
    return length(p) - r;
}

// 箱の距離関数
float sdBox( vec3 p, vec3 b)
{
    vec3 q = abs(p) - b;
    return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}

// シーンの距離関数
float map(vec3 p )
{
    float d = sdSphere(p + vec3(3.0, -1.0, -8.0), 2.0);
    // min関数で複数の物体の最短距離を求める
    d = min(d, sdBox(p + vec3(0.0, 1., -3.0), vec3(0.5, 0.5, 6.0)));
    d = min(d, sdBox(p + vec3(-3.0, 0.0, -2.0), vec3(0.3, 2.0, 2.0)));
    return d;
}

// 法線を求める関数
vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        map(vec3(p.x + EPSILON, p.y, p.z)) - map(vec3(p.x - EPSILON, p.y, p.z)),
        map(vec3(p.x, p.y + EPSILON, p.z)) - map(vec3(p.x, p.y - EPSILON, p.z)),
        map(vec3(p.x, p.y, p.z  + EPSILON)) - map(vec3(p.x, p.y, p.z - EPSILON))
    ));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (2.0 * fragCoord - iResolution.xy) / min(iResolution.x, iResolution.y);
    
    // レイを飛ばす位置のz座標
    float screenZ = -10.0;
    // カメラの位置
    vec3 cameraOrigin = vec3(0.0, 0.0, -20.0);
    // 飛ばすレイの方向
    vec3 rayDirection = normalize(vec3(uv, screenZ) - cameraOrigin);
    
    // 最短距離
    float dist = 0.0;
    float t = 0.0;
    // マーチングループ
    for(int i = 0; i < 400; i++)
    {
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        dist = map(rayPosition);
        
        t += dist;
        if(dist < 0.0001)
        {
            break;
        }
        
    }
    
    vec3 col = vec3(0.0);
    if(dist < 0.0001)
    {
        // 物体表面のレイ
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        // 物体表面の法線を求める
        vec3 normal = estimateNormal(rayPosition);
        vec3 lightDir = normalize(vec3(2., 2.0, -4.));
        float l = dot(normal, lightDir);
        col += clamp(l, 0.0, 1.0);
        col += 0.1;
    }
    
    fragColor = vec4(col,1.0);
}

前のコードからの変更点は、sdBox関数の追加とmap関数の変更です。実行すると、左上に球、下と右に箱が表示されます。箱の距離関数は

float sdBox( vec3 p, vec3 b)
{
    vec3 q = abs(p) - b;
    return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}

となっています。この距離関数はiq氏のサイトInigo Quilez :: computer graphics, mathematics, shaders, fractals, demoscene and moreからそのまま拝借しました。作者の方の解説もあるので載せておきます。

youtu.be

XY平面のみで考えると、幅2、高さ1のBoxは $\sqrt{\max\left(\operatorname{abs}\left(x\right)-2,0\ \right)^{2}+\max\left(\operatorname{abs}\left(y\right)-1,\ 0\right)}+\min\left(\max\left(\operatorname{abs}\left(x\right)-2,\ \operatorname{abs}\left(y\right)-1\right),\ 0\right)=0$ と表せます。 map関数では各物体への最短距離をmin関数で比較して最短距離を採用しています。距離関数にレイの先端の座標pを渡す時に足しているのは表示する座標をずらすためです。

float map(vec3 p )
{
    float d = sdSphere(p + vec3(3.0, -1.0, -8.0), 2.0);
    // min関数で複数の物体の最短距離を求める
    d = min(d, sdBox(p + vec3(0.0, 1., -3.0), vec3(0.5, 0.5, 6.0)));
    d = min(d, sdBox(p + vec3(-3.0, 0.0, -2.0), vec3(0.3, 2.0, 2.0)));
    return d;
}

次は形状を組み合わせてTUTを表示します。

形状の合成

#define EPSILON 0.001

// 箱の距離関数
float sdBox( vec3 p, vec3 b)
{
    vec3 q = abs(p) - b;
    return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}

float sdT(vec3 p)
{
    // Boxの合成
    float d1 = sdBox(p + vec3(0.0, -4.0, 0.0), vec3(5.0, 1.0, 1.0));
    float d2 = sdBox(p, vec3(1.0, 5.0, 1.0));
    return min(d1, d2);
}

float sdLink( vec3 p, float le, float r1, float r2 )
{
  vec3 q = vec3( p.x, max(abs(p.y)-le,0.0), p.z );
  return length(vec2(length(q.xy)-r1,q.z)) - r2;
}

float sdU(vec3 p)
{
    float d1 = sdLink(p + vec3(0.0, -2.5, 0.0), 4.0, 3.2, 1.0);
    // 上部分を切り取る
    float d2 = sdBox(p + vec3(0.0, -9.0, 0.0), vec3(5.0, 4.0, 1.0));
    return max(d1, -d2);
}

// シーンの距離関数
float map(vec3 p )
{
    float d1 = sdT(p + vec3(10.0, 0.0, 0.0));
    float d2 = sdU(p + vec3(0.0, 0.0, 0.0));
    float d3 = sdT(p + vec3(-10.0, 0.0, 0.0));
    return min(d1, min(d2, d3));
}

// 法線を求める関数
vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        map(vec3(p.x + EPSILON, p.y, p.z)) - map(vec3(p.x - EPSILON, p.y, p.z)),
        map(vec3(p.x, p.y + EPSILON, p.z)) - map(vec3(p.x, p.y - EPSILON, p.z)),
        map(vec3(p.x, p.y, p.z  + EPSILON)) - map(vec3(p.x, p.y, p.z - EPSILON))
    ));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (2.0 * fragCoord - iResolution.xy) / min(iResolution.x, iResolution.y);
    
    // レイを飛ばす位置のz座標
    float screenZ = -18.0;
    // カメラの位置
    vec3 cameraOrigin = vec3(0.0, 0.0, -20.0);
    // 飛ばすレイの方向
    vec3 rayDirection = normalize(vec3(uv, screenZ) - cameraOrigin);
    
    // 最短距離
    float dist = 0.0;
    float t = 0.0;
    // マーチングループ
    for(int i = 0; i < 400; i++)
    {
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        dist = map(rayPosition);
        
        t += dist;
        if(dist < 0.0001)
        {
            break;
        }
        
    }
    
    vec3 col = vec3(0.0);
    if(dist < 0.0001)
    {
        // 物体表面のレイ
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        // 物体表面の法線を求める
        vec3 normal = estimateNormal(rayPosition);
        vec3 lightDir = normalize(vec3(2., 2.0, -4.));
        float l = dot(normal, lightDir);
        col += clamp(l, 0.0, 1.0);
        col += 0.1;
    }
    
    fragColor = vec4(col,1.0);
}

実行するとTUTが表示されます。変更点は距離関数の部分のみです。

float sdT(vec3 p)
{
    // Boxの合成
    float d1 = sdBox(p + vec3(0.0, -4.0, 0.0), vec3(5.0, 1.0, 1.0));
    float d2 = sdBox(p, vec3(1.0, 5.0, 1.0));
    return min(d1, d2);
}

sdTで(5.0,1.0,1.0)の横に長いBoxと(1.0,5.0,1.0)の縦に長いBoxをmin関数で合成しています。

float sdLink( vec3 p, float le, float r1, float r2 )
{
  vec3 q = vec3( p.x, max(abs(p.y)-le,0.0), p.z );
  return length(vec2(length(q.xy)-r1,q.z)) - r2;
}

float sdU(vec3 p)
{
    float d1 = sdLink(p + vec3(0.0, -2.5, 0.0), 4.0, 3.2, 1.0);
    // 上部分を切り取る
    float d2 = sdBox(p + vec3(0.0, -9.0, 0.0), vec3(5.0, 4.0, 1.0));
    return max(d1, -d2);
}

sdLinkはドーナツを引き伸ばした形状の立体の距離関数です。 sdUでは、sdLinkの最短距離とsdBoxの最短距離をマイナスの値の大きい方を返しています。最短距離d1のある形状から最短距離d2の形状を引いた形状を作りたい場合、max(d1,-d2)とすることでd2の形状が切り取られた形になります。これによりsdLinkの形状の上部分がsdBoxによって切り取られた形状になります。 map関数ではレイのx座標をずらして左から各TUT形状の距離関数を求めて、一番小さい距離関数を返しています。

modを使って同じ形状の複製をしてみます。

物体の複製

#define EPSILON 0.001

// 箱の距離関数
float sdBox( vec3 p, vec3 b)
{
    vec3 q = abs(p) - b;
    return length(max(q,0.0)) + min(max(q.x,max(q.y,q.z)),0.0);
}

float sdT(vec3 p)
{
    // Boxの合成
    float d1 = sdBox(p + vec3(0.0, -4.0, 0.0), vec3(5.0, 1.0, 1.0));
    float d2 = sdBox(p, vec3(1.0, 5.0, 1.0));
    return min(d1, d2);
}

float sdLink( vec3 p, float le, float r1, float r2 )
{
  vec3 q = vec3( p.x, max(abs(p.y)-le,0.0), p.z );
  return length(vec2(length(q.xy)-r1,q.z)) - r2;
}

float sdU(vec3 p)
{
    float d1 = sdLink(p + vec3(0.0, -2.5, 0.0), 4.0, 3.2, 1.0);
    // 上部分を切り取る
    float d2 = sdBox(p + vec3(0.0, -9.0, 0.0), vec3(9.0, 4.0, 1.1));
    return max(d1, -d2);
}

// シーンの距離関数
float map(vec3 p )
{
    float c = 70.0;
    // 繰り返す
    vec3 q = mod(p, c) - 0.5 * c;
    float d1 = sdT(q + vec3(20.0, 0.0, 0.0));
    float d2 = sdU(q + vec3(-20.0, 0.0, 0.0));
    return min(d1, d2);
}

// 法線を求める関数
vec3 estimateNormal(vec3 p) {
    return normalize(vec3(
        map(vec3(p.x + EPSILON, p.y, p.z)) - map(vec3(p.x - EPSILON, p.y, p.z)),
        map(vec3(p.x, p.y + EPSILON, p.z)) - map(vec3(p.x, p.y - EPSILON, p.z)),
        map(vec3(p.x, p.y, p.z  + EPSILON)) - map(vec3(p.x, p.y, p.z - EPSILON))
    ));
}

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 uv = (2.0 * fragCoord - iResolution.xy) / min(iResolution.x, iResolution.y);
    
    // レイを飛ばす位置のz座標
    float screenZ = 0.0;
    // カメラの位置
    vec3 cameraOrigin = vec3(sin(iTime * 0.2) * 10.0, cos(iTime * 0.2) * 10.0, -20.0);
    // 飛ばすレイの方向
    vec3 rayDirection = normalize(vec3(uv, screenZ) - cameraOrigin);
    
    // 最短距離
    float dist = 0.0;
    float t = 0.0;
    // マーチングループ
    for(int i = 0; i < 400; i++)
    {
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        dist = map(rayPosition);
        
        t += dist;
        if(dist < 0.0001)
        {
            break;
        }
        
    }
    
    vec3 col = vec3(0.0);
    if(dist < 0.0001)
    {
        // 物体表面のレイ
        vec3 rayPosition = cameraOrigin + rayDirection * t;
        // 物体表面の法線を求める
        vec3 normal = estimateNormal(rayPosition);
        vec3 lightDir = normalize(vec3(2., 2.0, -4.));
        float l = dot(normal, lightDir);
        col += clamp(l, 0.0, 1.0);
        col += 0.1;
    }
    
    fragColor = vec4(col,1.0);
}

実行するとTとUが置くまで並べられた映像になります。ポリゴンでTUの形状を作成して描写する方法でこの数を描写することは負荷が厳しいと思うので、簡単に複製できるのはレイマーチングの利点かもしれません。オブジェクトの複製はmap関数内のvec3 q = mod(p, c) - 0.5 * c;で行っています。modを利用することでqに入る値が繰り返される仕組みです。cは繰り返しの周期で、0.5-cは繰り返しの中心を半周期分ずらしています。ちなみに、q.xy = mod(p.xy, c) - 0.5 * c;のようにすると、xy平面の繰り返しとなります。

float c = 20.0;
// xy平面での複製
vec3 q = p;
q.z += -2000.0;
q.xy = mod(p.xy, c) - 0.5 * c;

TUT(^^)

おわりに

レイマーチングは短いコードで良い感じの絵を出力できるので皆さんもやってみましょう！

最後に自分の作った作品をいくつか載せます。画像の10MB制限に収めるためにかなり画質を荒くしていますが、実際はもっと綺麗な映像です。いずれも先程までのコードのノウハウ+αで作ることができます。特にmodによる複製は便利です。

レイマーチングではないですがシェーダ―芸では下のようなものも作成できます。

明日のアドカレ記事はどどどど氏のの「自宅サーバは良いぞ」です。

参考リンク

Inigo Quilez :: computer graphics, mathematics, shaders, fractals, demoscene and more レイマーチングやシェーダに関する幅広いノウハウが載っているサイトです。自分が使う基本的な距離関数はほとんどここから拝借しています。

[GLSL] レイマーチング入門 vol.1 - Qiita

ゼロから始めるレイマーチング - Qiita

sh00tの日記

レイマーチングことはじめ