GIパストレーサsmallptのprocessing移植

はじめに

上記サイトにsmallpt99(99行のsmall path tracing)があります。
これをProcessingに変換します。

環境

M1 Mac E2Core/P2Core
Processing 4.4.10

結果

Threads: 8 SAMPLES: 256
time(ms) = 15056
total rays = 739339043
rays / second = 4.9105944E7
5000万 ray/s
※バチクソにハマったとき

コード

オリジナル含め、壁は巨大な球の一部となっています。

300行こえています。

ray.pde

// =====================================
// smallpt風 Cornell Box Path Tracer
// - float に統一
// - Thread でマルチスレッド
// - 上下反転
// =====================================

int W = 320;
int H = 240;
int SAMPLES = 256;           // 1ピクセルのサンプル数
int THREADS;                 // CPU コア数
int MAX_DEPTH = 16;          // 再帰の最大深さ（保険）
final float S = 0.001f;      // シーン全体を 1/1000 に縮小 doubleからfloat精度にした都合で

PImage img;

// カメラ（スレッドから参照）
Ray cam;
Vec cx, cy;

// ray計算
long rayCount = 0;             // 全スレッド合計のray本数
Object rayLock = new Object(); // 合算用ロック（1スレッド1回だけ使う）


// -------------------- setup / draw --------------------
void setup() {
  size(320, 240);
  pixelDensity(1);
  img = createImage(W, H, RGB);

  THREADS = 8;//Runtime.getRuntime().availableProcessors();
  println("Threads:", THREADS, "  SAMPLES:", SAMPLES);

  cam = new Ray(
    new Vec(50f*S, 52f*S, 295.6f*S),
    new Vec(0f, -0.042612f, -1f).norm()  // 方向はそのまま
    );
  cx = new Vec(W * 0.5135f / H, 0, 0);
  cy = cx.cross(cam.d).norm().mul(0.5135f);

  println("Rendering start...");
  renderMultiThread();
  println("Rendering finished.");
}

void draw() {
  image(img, 0, 0);
}

// -------------------- ベクトル / レイ --------------------
class Vec {
  float x, y, z;
  Vec(float x_, float y_, float z_) {
    x=x_;
    y=y_;
    z=z_;
  }
  Vec() {
    this(0, 0, 0);
  }

  Vec add(Vec b) {
    return new Vec(x+b.x, y+b.y, z+b.z);
  }
  Vec sub(Vec b) {
    return new Vec(x-b.x, y-b.y, z-b.z);
  }
  Vec mul(float b) {
    return new Vec(x*b, y*b, z*b);
  }
  Vec mult(Vec b) {
    return new Vec(x*b.x, y*b.y, z*b.z);
  }

  Vec norm() {
    float n = sqrt(x*x + y*y + z*z);
    return mul(1.0f / n);
  }
  float dot(Vec b) {
    return x*b.x + y*b.y + z*b.z;
  }
  Vec cross(Vec b) {
    return new Vec(
      y*b.z - z*b.y,
      z*b.x - x*b.z,
      x*b.y - y*b.x
      );
  }
}

class Ray {
  Vec o, d;
  Ray(Vec o_, Vec d_) {
    o=o_;
    d=d_;
  }
}

enum Refl_t {
  DIFF, SPEC, REFR
}

// -------------------- 球 --------------------
class Sphere {
  float rad;
  Vec p, e, c;
  Refl_t refl;
  Sphere(float rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_) {
    rad = rad_;
    p   = p_;
    e   = e_;
    c   = c_;
    refl= refl_;
  }

  float intersect(Ray r) {
    Vec op = p.sub(r.o);
    float eps = 1e-4f;
    float b = op.dot(r.d);
    float det = b*b - op.dot(op) + rad*rad;
    if (det < 0) return 0;
    det = sqrt(det);
    float t;
    if ((t = b - det) > eps) return t;
    if ((t = b + det) > eps) return t;
    return 0;
  }
}

// -------------------- シーン --------------------
Sphere[] spheres = {
  new Sphere(1e5f*S, new Vec((1e5f+1f)*S, 40.8f*S, 81.6f*S), new Vec(), new Vec(.75f, .25f, .25f), Refl_t.DIFF),
  new Sphere(1e5f*S, new Vec((-1e5f+99f)*S, 40.8f*S, 81.6f*S), new Vec(), new Vec(.25f, .25f, .75f), Refl_t.DIFF),
  new Sphere(1e5f*S, new Vec(50f*S, 40.8f*S, 1e5f*S), new Vec(), new Vec(.75f, .75f, .75f), Refl_t.DIFF),
  new Sphere(1e5f*S, new Vec(50f*S, 40.8f*S, (-1e5f+170f)*S), new Vec(), new Vec(), Refl_t.DIFF),
  new Sphere(1e5f*S, new Vec(50f*S, 1e5f*S, 81.6f*S), new Vec(), new Vec(.75f, .75f, .75f), Refl_t.DIFF),
  new Sphere(1e5f*S, new Vec(50f*S, (-1e5f+81.6f)*S, 81.6f*S), new Vec(), new Vec(.75f, .75f, .75f), Refl_t.DIFF),
  new Sphere(16.5f*S, new Vec(27f*S, 16.5f*S, 47f*S), new Vec(), new Vec(1, 1, 1).mul(.999f), Refl_t.SPEC),
  new Sphere(16.5f*S, new Vec(73f*S, 16.5f*S, 78f*S), new Vec(), new Vec(1, 1, 1).mul(.999f), Refl_t.REFR),
  new Sphere(600f*S, new Vec(50f*S, (681.6f-.27f)*S, 81.6f*S), new Vec(12, 12, 12), new Vec(), Refl_t.DIFF)
};

// -------------------- 共通ユーティリティ --------------------
float clamp(float x) {
  return x<0?0:x>1?1:x;
}
int   toInt(float x) {
  return int(pow(clamp(x), 1.0f/2.2f)*255.0f+0.5f);
}

// 最近ヒット
int intersectScene(Ray r, float[] tOut) {
  float inf = 1e20f;
  float t = inf;
  int id = -1;
  for (int i=0; i<spheres.length; i++) {
    float d = spheres[i].intersect(r);
    if (d!=0 && d<t) {
      t = d;
      id= i;
    }
  }
  tOut[0] = t;
  return id;
}

// -------------------- マルチスレッドワーカー --------------------
class RenderWorker extends Thread {
  int id;
  int seed;  // スレッド専用 RNG 状態
  long localRayCount = 0;   // このスレッドが使ったrayの本数

  RenderWorker(int id_) {
    id = id_;
    seed = 1234567 + id * 12345;   // 適当な初期値
  }

  // 線形合同法 RNG（0〜1）
  float rnd() {
    seed = seed * 1664525 + 1013904223;
    int bits = (seed >>> 9) & 0x007FFFFF;    // 上位 23bit
    return bits / (float)(1<<23);
  }

  Vec radiance(Ray r, int depth) {
    localRayCount++;
    if (depth > MAX_DEPTH) return new Vec();

    float[] tArr = new float[1];
    int idHit = intersectScene(r, tArr);
    if (idHit == -1) return new Vec();  // miss

    Sphere obj = spheres[idHit];
    Vec x  = r.o.add(r.d.mul(tArr[0]));
    Vec n  = x.sub(obj.p).norm();
    Vec nl = (n.dot(r.d)<0) ? n : n.mul(-1);
    Vec f  = obj.c;

    float p = max(f.x, max(f.y, f.z));

    depth++;
    if (depth > 5) {
      if (rnd() < p) {
        f = f.mul(1.0f/p);
      } else {
        return obj.e;
      }
    }

    if (obj.refl == Refl_t.DIFF) {
      float r1  = TWO_PI * rnd();
      float r2  = rnd();
      float r2s = sqrt(r2);

      Vec w = nl;
      Vec u = ((abs(w.x) > 0.1f ? new Vec(0, 1, 0) : new Vec(1, 0, 0)).cross(w)).norm();
      Vec v = w.cross(u);

      Vec d = u.mul(cos(r1)*r2s)
        .add(v.mul(sin(r1)*r2s))
        .add(w.mul(sqrt(1.0f-r2)))
        .norm();

      return obj.e.add(f.mult(radiance(new Ray(x, d), depth)));
    }

    if (obj.refl == Refl_t.SPEC) {
      Vec reflDir = r.d.sub(n.mul(2.0f * n.dot(r.d)));
      return obj.e.add(f.mult(radiance(new Ray(x, reflDir), depth)));
    }

    // 屈折
    Ray reflRay = new Ray(x, r.d.sub(n.mul(2.0f * n.dot(r.d))));
    boolean into = n.dot(nl) > 0.0f;

    float nc = 1.0f;
    float nt = 1.5f;
    float nnt = into ? nc/nt : nt/nc;
    float ddn = r.d.dot(nl);
    float cos2t = 1.0f - nnt*nnt*(1.0f - ddn*ddn);

    if (cos2t < 0) {
      return obj.e.add(f.mult(radiance(reflRay, depth)));
    }

    Vec tdir = r.d.mul(nnt)
      .sub(n.mul((into?1.0f:-1.0f)*(ddn*nnt+sqrt(cos2t))))
      .norm();

    float a  = nt-nc;
    float b  = nt+nc;
    float R0 = (a*a)/(b*b);
    float c  = 1.0f - (into ? -ddn : tdir.dot(n));
    float Re = R0 + (1.0f - R0)*pow(c, 5.0f);
    float Tr = 1.0f - Re;
    float P  = 0.25f + 0.5f*Re;
    float RP = Re / P;
    float TP = Tr / (1.0f - P);

    if (depth > 2) {
      if (rnd() < P) {
        return obj.e.add(f.mult(radiance(reflRay, depth).mul(RP)));
      } else {
        return obj.e.add(f.mult(radiance(new Ray(x, tdir), depth).mul(TP)));
      }
    }

    return obj.e.add(
      f.mult(
      radiance(reflRay, depth).mul(Re)
      .add(radiance(new Ray(x, tdir), depth).mul(Tr))
      )
      );
  }

  public void run() {
    for (int y = id; y < H; y += THREADS) {
      int row = (H-1-y) * W;   // 上下反転
      for (int x = 0; x < W; x++) {
        Vec col = new Vec();

        for (int sy=0; sy<2; sy++) {
          for (int sx=0; sx<2; sx++) {
            Vec sub = new Vec();
            for (int s=0; s<SAMPLES; s++) {
              float r1 = 2.0f * rnd();
              float dx = (r1 < 1.0f) ? sqrt(r1)-1.0f : 1.0f - sqrt(2.0f-r1);
              float r2 = 2.0f * rnd();
              float dy = (r2 < 1.0f) ? sqrt(r2)-1.0f : 1.0f - sqrt(2.0f-r2);

              float u = (((sx+0.5f+dx)/2.0f + x) / (float)W - 0.5f);
              float v = (((sy+0.5f+dy)/2.0f + y) / (float)H - 0.5f);

              Vec d = cx.mul(u)
                .add(cy.mul(v))
                .add(cam.d)
                .norm();

              Vec sample = radiance(new Ray(cam.o.add(d.mul(140.0f*S)), d), 0);
              sub = sub.add(sample.mul(1.0f/SAMPLES));
            }
            col = col.add(new Vec(
              clamp(sub.x),
              clamp(sub.y),
              clamp(sub.z)
              ).mul(0.25f));
          }
        }

        int idx = row + x;
        img.pixels[idx] = color(
          toInt(col.x),
          toInt(col.y),
          toInt(col.z)
          );
      }
    }
    // このスレッド分のrayを合計に足す
    synchronized (rayLock) {
      rayCount += localRayCount;
    }
  }
}

// -------------------- レンダリング起動 --------------------
void renderMultiThread() {
  img.loadPixels();

  RenderWorker[] ws = new RenderWorker[THREADS];
  for (int i=0; i<THREADS; i++) {
    ws[i] = new RenderWorker(i);
    ws[i].start();
  }
  for (int i=0; i<THREADS; i++) {
    try {
      ws[i].join();
    }
    catch(Exception e) {
    }
  }

  img.updatePixels();
  
  // 結果表示（起動からの経過時間でOKという指定）
  long dt = millis();
  println("time(ms)      = " + dt);
  println("total rays    = " + rayCount);
  println("rays / second = " + (rayCount * 1000.0 / dt));
}