GPUのパストレはやっぱり速くてよい

レイトレをGPU実装するのは、再帰処理できない場合、大変。再帰的処理をコードで展開する必要があるから。

それに比べて、パストレをGPU実装するのは、楽である。確率で処理するので、上記の心配がない。ただし、確立計算に伴うノイズが絶対に生じる。デノイズ必須。

このGPUパストレの特徴

チェック床＋普通の球、ガラス球、ミラー球
デノイズなし
caustics表現なし
VSYC回避のため frameRate(960)

操作

数字キーの1~9で、積算の度合いを調整。
1ならそのまま表示
6なら32回分を積算して表示
(2^(n-1)=32)
ノイズは軽減するが、フレームレートは落ちる。

計算速度

Processing 4.4.10
MBP M1
16 Core GPU
この環境で、32回積算時に、20FPSぐらい出るので、GPU処理の恩恵を感じる。
積算無しで、デノイズするのが、一番いいと思う。
ちなみに、解像度を変更しても、ほとんどFPSには影響しません。CPUと比べて、圧倒的に有利です。

他環境

MSIのノートPC
■iGPUの場合
32回積算時に、5FPS
■GeForce 3050Ti Laptop(電力半分版)
GPU優先にしたけれど、iGPUと使用率が半々になってしまって、性能出しきれない。
iGPUに引っ張られる状態で、
32回積算時に、28FPS

デスクトップPCで後日試す。

pathtraceShader.pde

PShader pathtracing;

// 見えない積算バッファ
PGraphics accum;

// 表示専用バッファ（確定フレームだけ表示）
PGraphics displayImg;

int subSamplesPerFrame = 1; // n = 1,2,4,8...
int subSampleIndex = 0;

float animTime = 0.0;
float timeStep = 0.005;

int logicalFrameCount = 0;   // 表示したフレーム数
float logicalFPS = 0;        // 表示FPS
int lastDisplayTime = 0;     // 最後に表示した時刻（millis）

void setup() {
  //fullScreen(P2D);
  size(1280, 720, P2D);
  pixelDensity(1);
  frameRate(960);            // VSync回避してGPU100%使うため

  pathtracing = loadShader("pathtracing.txt");

  accum      = createGraphics(width, height, P2D);
  displayImg = createGraphics(width, height, P2D);

  clearAccum();
  clearDisplay();

  lastDisplayTime = millis();
}

// 積算バッファだけ消す
void clearAccum() {
  accum.beginDraw();
  accum.background(0);
  accum.endDraw();
  subSampleIndex = 0;
}

// 表示バッファだけ消す
void clearDisplay() {
  displayImg.beginDraw();
  displayImg.background(0);
  displayImg.endDraw();
}

void draw() {

  // --------- 見えない所で1サンプル積む ---------
  subSampleIndex++;

  accum.beginDraw();
  accum.shader(pathtracing);
  pathtracing.set("u_resolution", (float)width, (float)height);
  pathtracing.set("u_prevFrame", accum);
  pathtracing.set("u_sampleCount", (float)subSampleIndex);
  pathtracing.set("u_time", animTime);   // n サンプルの間は固定
  accum.rect(0, 0, width, height);
  accum.resetShader();
  accum.endDraw();

  // --------- nサンプル積み終わった瞬間だけ確定表示 ---------
  if (subSampleIndex >= subSamplesPerFrame) {

    // 確定結果を displayImg にコピー
    displayImg.beginDraw();
    displayImg.image(accum, 0, 0);
    displayImg.endDraw();

    // 表示フレームのFPSを計算
    int now = millis();
    int dt = now - lastDisplayTime;
    if (dt > 0) {
      logicalFPS = 1000.0 / dt;  // ms → fps
    }
    lastDisplayTime = now;
    logicalFrameCount++;
    
    // 次フレームへ（ライトだけここで動く）
    animTime += timeStep * subSamplesPerFrame;

    // 次のフレーム用に accum をクリア
    clearAccum();
  }

  // --------- 画面には「確定済み」だけ表示 ---------
  image(displayImg, 0, 0);

  // デバッグ表示
  fill(255);
  textAlign(LEFT, TOP);
  text(
    "n = " + subSamplesPerFrame +
    "\n subSampleIndex = " + subSampleIndex +
    //"\n animTime = " + nf(animTime, 1, 2) +
    "\n FPS = " + nf(logicalFPS, 1, 1) +
    "\n Frames = " + logicalFrameCount,
    10, 10
    );
}


void keyPressed() {
  if (key >= '1' && key <= '9') {
    int k = key - '0';                  // '1'→1, '2'→2 ...
    subSamplesPerFrame = 1 << (k - 1);  // 2^(k-1) → 1,2,4,8...256
    clearAccum();
    clearDisplay();
  }
}

pathtracing.txt

#ifdef GL_ES
precision mediump float;
precision mediump int;
#endif

uniform vec2  u_resolution;
uniform float u_time;
uniform float u_sampleCount;
uniform sampler2D u_prevFrame;

const float PI = 3.14159265359;

// -------------------- 基本構造体 --------------------

struct Ray {
    vec3 ori;
    vec3 dir;
};

struct Hit {
    bool hit;
    float t;
    vec3 pos;
    vec3 normal;
    int  matId;
};

// マテリアルID
const int MAT_FLOOR  = 0;
const int MAT_DIFF   = 1;
const int MAT_GLASS  = 2;
const int MAT_MIRROR = 3;

// 球
const int NUM_SPHERES = 3;
vec3  sphereCenter[NUM_SPHERES];
float sphereRadius[NUM_SPHERES];
int   sphereMat[NUM_SPHERES];

vec3  lightPos;
vec3  lightColor = vec3(6.0);

// -------------------- 乱数 --------------------

float hash31(vec3 p) {
    p  = fract(p * 0.1031);
    p += dot(p, p.yzx + 33.33);
    return fract((p.x + p.y) * p.z);
}

float rand(inout vec3 seed) {
    float r = hash31(seed);
    seed += vec3(1.0, 1.0, 1.0);
    return r;
}

// コサイン重み付き半球サンプル
vec3 cosineSampleHemisphere(vec2 u) {
    float r = sqrt(u.x);
    float theta = 2.0 * PI * u.y;
    float x = r * cos(theta);
    float z = r * sin(theta);
    float y = sqrt(max(0.0, 1.0 - u.x));
    return vec3(x, y, z);
}

// -------------------- シーン初期化 --------------------

void initScene() {
    // 左：赤い拡散球
    sphereCenter[0] = vec3(-1.2, 0.5, 0.0);
    sphereRadius[0] = 0.5;
    sphereMat[0]    = MAT_DIFF;

    // 中央：ガラス球
    sphereCenter[1] = vec3( 0.0, 0.5, 0.0);
    sphereRadius[1] = 0.5;
    sphereMat[1]    = MAT_GLASS;

    // 右：ミラー球
    sphereCenter[2] = vec3( 1.2, 0.5, 0.0);
    sphereRadius[2] = 0.5;
    sphereMat[2]    = MAT_MIRROR;

    lightPos = vec3(2.0 * sin(u_time * 0.5), 3.0, 2.0 * cos(u_time * 0.5));
}

// -------------------- 交差計算 --------------------

bool intersectSphere(Ray ray, vec3 center, float radius, out float t) {
    vec3 oc = ray.ori - center;
    float b = dot(oc, ray.dir);
    float c = dot(oc, oc) - radius * radius;
    float h = b * b - c;
    if (h < 0.0) return false;
    h = sqrt(h);
    float t0 = -b - h;
    float t1 = -b + h;
    t = (t0 > 0.0) ? t0 : t1;
    return t > 0.0;
}

// 無限平面 y = 0
bool intersectPlane(Ray ray, float y, out float t) {
    if (abs(ray.dir.y) < 1e-4) return false;
    t = (y - ray.ori.y) / ray.dir.y;
    return t > 0.0;
}

// 背景（空）
vec3 backgroundColor(vec3 dir) {
    float t = 0.5 * (dir.y + 1.0);
    return mix(vec3(0.7, 0.8, 1.0), vec3(0.1, 0.1, 0.2), t);
}

// 床チェッカー
vec3 floorColor(vec3 p) {
    float scale = 2.0;
    float fx = floor(p.x * scale);
    float fz = floor(p.z * scale);
    float checker = mod(fx + fz, 2.0);
    return mix(vec3(0.9), vec3(0.1), checker);
}

// シーン全体の交差
Hit traceScene(Ray ray) {
    Hit h;
    h.hit = false;
    h.t = 1e20;
    h.matId = -1;

    float t;

    // 床
    if (intersectPlane(ray, 0.0, t)) {
        if (t < h.t) {
            h.hit = true;
            h.t = t;
            h.pos = ray.ori + ray.dir * t;
            h.normal = vec3(0.0, 1.0, 0.0);
            h.matId = MAT_FLOOR;
        }
    }

    // 球たち
    for (int i = 0; i < NUM_SPHERES; i++) {
        if (intersectSphere(ray, sphereCenter[i], sphereRadius[i], t)) {
            if (t < h.t) {
                h.hit = true;
                h.t   = t;
                h.pos = ray.ori + ray.dir * t;
                h.normal = normalize(h.pos - sphereCenter[i]);
                h.matId  = sphereMat[i];
            }
        }
    }
    return h;
}

// -------------------- シェーディング --------------------

// 拡散：床と赤球を共通処理
vec3 shadeDiffuse(Hit h, inout Ray ray, inout vec3 throughput, inout vec3 seed) {
    vec3 N = normalize(h.normal);

    vec3 baseColor;
    if (h.matId == MAT_FLOOR) {
        baseColor = floorColor(h.pos);
    } else {
        // 赤い拡散球
        baseColor = vec3(0.9, 0.2, 0.2);
    }

    // 直射光（全ての球でシャドウ）
    vec3 Ldir = normalize(lightPos - h.pos);
    float distL = length(lightPos - h.pos);

    Ray shadowRay;
    shadowRay.ori = h.pos + N * 1e-3;
    shadowRay.dir = Ldir;

    bool blocked = false;
    float tSh;

    for (int i = 0; i < NUM_SPHERES; i++) {
        if (intersectSphere(shadowRay, sphereCenter[i], sphereRadius[i], tSh) && tSh < distL) {
            blocked = true;
            break;
        }
    }

    float nl = max(dot(N, Ldir), 0.0);
    vec3 direct = vec3(0.0);
    if (!blocked && nl > 0.0) {
        float att = 1.0 / (1.0 + 0.1 * distL * distL);
        direct = baseColor * lightColor * nl * att;
    }

    vec3 ambient = baseColor * 0.05;

    // 次の間接レイ
    vec2 u = vec2(rand(seed), rand(seed));
    vec3 localDir = cosineSampleHemisphere(u);

    vec3 up = (abs(N.y) < 0.99) ? vec3(0.0, 1.0, 0.0) : vec3(1.0, 0.0, 0.0);
    vec3 T = normalize(cross(up, N));
    vec3 B = cross(N, T);

    vec3 worldDir = normalize(
        T * localDir.x +
        N * localDir.y +
        B * localDir.z
    );

    ray.ori = h.pos + N * 1e-3;
    ray.dir = worldDir;

    throughput *= baseColor;

    return direct + ambient;
}

// ガラス：確率的な反射 / 屈折
vec3 shadeGlass(Hit h, inout Ray ray, inout vec3 throughput, inout vec3 seed) {
    vec3 N = normalize(h.normal);
    float ior = 1.5;

    bool entering = dot(ray.dir, N) < 0.0;
    vec3 n = entering ? N : -N;
    float eta = entering ? (1.0 / ior) : (ior / 1.0);

    float cosI = clamp(-dot(ray.dir, n), 0.0, 1.0);
    float R0 = pow((1.0 - ior) / (1.0 + ior), 2.0);
    float F  = R0 + (1.0 - R0) * pow(1.0 - cosI, 5.0);

    float pReflect = clamp(F, 0.05, 0.95);
    float r = rand(seed);

    vec3 glassTint = vec3(0.98, 0.99, 1.0);

    if (r < pReflect) {
        // 反射
        vec3 reflDir = reflect(ray.dir, n);
        ray.ori = h.pos + n * 1e-3;
        ray.dir = normalize(reflDir);
    } else {
        // 屈折
        vec3 refrDir = refract(ray.dir, n, eta);
        if (length(refrDir) == 0.0) {
            // 全反射
            vec3 reflDir = reflect(ray.dir, n);
            ray.ori = h.pos + n * 1e-3;
            ray.dir = normalize(reflDir);
        } else {
            ray.ori = h.pos - n * 1e-3;
            ray.dir = normalize(refrDir);
        }
        throughput *= glassTint;
    }
    return vec3(0.0);
}

// ミラー：完全鏡面（レイの向きを変えるだけ）
vec3 shadeMirror(Hit h, inout Ray ray, inout vec3 throughput, inout vec3 seed) {
    vec3 N = normalize(h.normal);
    vec3 reflDir = reflect(ray.dir, N);
    ray.ori = h.pos + N * 1e-3;
    ray.dir = normalize(reflDir);
    // ミラー自体は発光しない
    return vec3(0.0);
}

// -------------------- main --------------------

void main() {
    initScene();

    vec2 fragCoord = gl_FragCoord.xy;
    vec2 uv = (fragCoord / u_resolution) * 2.0 - 1.0;
    uv.x *= u_resolution.x / u_resolution.y;

// ===== カメラ：y軸周りを回転しながら中央を見る =====
vec3 camTarget = vec3(0.0, 0.5, 0.0);  // 真ん中のガラス球あたりを見る

float camRadius = 3.0;   // 中心からの距離
float camHeight = 1.5;   // カメラの高さ
float camSpeed  = 0.2;   // 回転速度

// u_time は Processing 側の animTime から来ている
vec3 camPos = vec3(
    camRadius * sin(u_time * camSpeed),
    camHeight,
    camRadius * cos(u_time * camSpeed)
);

vec3 camDir   = normalize(camTarget - camPos);
vec3 up       = vec3(0.0, 1.0, 0.0);
vec3 camRight = normalize(cross(camDir, up));
vec3 camUp    = cross(camRight, camDir);

float fov   = radians(45.0);
float focal = 1.0 / tan(fov * 0.5);

vec3 dir = normalize(
    camRight * uv.x +
    camUp    * uv.y +
    camDir   * focal
);


    Ray ray;
    ray.ori = camPos;
    ray.dir = dir;

    vec3 seed = vec3(fragCoord, u_sampleCount);

    const int MAX_BOUNCES = 6;

    vec3 L = vec3(0.0);
    vec3 throughput = vec3(1.0);

    for (int bounce = 0; bounce < MAX_BOUNCES; bounce++) {
        Hit h = traceScene(ray);

        if (!h.hit) {
            L += throughput * backgroundColor(ray.dir);
            break;
        }

        if (h.matId == MAT_FLOOR || h.matId == MAT_DIFF) {
            vec3 direct = shadeDiffuse(h, ray, throughput, seed);
            L += throughput * direct;
        } else if (h.matId == MAT_GLASS) {
            vec3 direct = shadeGlass(h, ray, throughput, seed);
            L += throughput * direct;
        } else if (h.matId == MAT_MIRROR) {
            vec3 direct = shadeMirror(h, ray, throughput, seed);
            L += throughput * direct;
        }

        // ロシアンルーレット
        if (bounce > 2) {
            float p = max(throughput.r, max(throughput.g, throughput.b));
            p = clamp(p, 0.05, 1.0);
            if (rand(seed) > p) break;
            throughput /= p;
        }
    }

    // 累積
    vec2 texCoord = fragCoord / u_resolution;
    vec3 prev = texture2D(u_prevFrame, texCoord).rgb;
    float n = max(u_sampleCount, 1.0);
    vec3 accumColor = (prev * (n - 1.0) + L) / n;

    gl_FragColor = vec4(accumColor, 1.0);
}