3点形式でパストレーサーを動かしてみる

struct Vertex {
  Vec p, n, e, c;  // 交点、法線、放射輝度、アルベド
  Refl_t refl;     // 反射タイプ (DIFF, SPEC, REFR)
  double pdf;      // 反射レイの確率密度
  Vertex(Vec p_, Vec n_, Vec e_, Vec c_, Refl_t refl_, double pdf_=1)
    : p(p_), n(n_), e(e_), c(c_), refl(refl_), pdf(pdf_) {}
};
void traceScene(Ray r, unsigned short *Xi, vector<Vertex>* path, int maxDepth=100){
  for (int depth = 0; depth < maxDepth; ) {
    double t; int id=0;
    if (!intersect(r, t, id)) break;
    const Sphere &obj = spheres[id];        // the hit object
    Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1;
    if (obj.e.x + obj.e.y + obj.e.z > 0) {  // 光源に到達
      path->emplace_back(x, nl, obj.e, obj.c, obj.refl);
      return;
    }
    double maxrefl = fmax(obj.c.x, fmax(obj.c.y, obj.c.z)), pdf=1;
    if (++depth>5||!maxrefl) { if (erand48(Xi)<maxrefl) pdf=maxrefl; else break; } //R.R
    if (obj.refl == DIFF) {                     // Ideal DIFFUSE reflection
      Vec d = cosine_weighted_random_hemisphere(nl, Xi);
      r = Ray(x, d);
      pdf *= nl.dot(d) / M_PI;  // 反射方向の確率密度を反映
    } else {
      Vec wi = reflect(r.d, n);
      if (obj.refl == REFR) {
        bool into = n.dot(nl)>0;            // Ray from outside going in?
        double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
        if ((cos2t=1-nnt*nnt*(1-ddn*ddn))>=0) { // Total internal reflection
          Vec tdir = refract(r.d, nl, nnt);
          double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
          double Re=R0+(1-R0)*c*c*c*c*c;
          if (erand48(Xi)>=Re) wi = tdir;
        }
      }
      r = Ray(x, wi);
    }
    path->emplace_back(x, nl, obj.e, obj.c, obj.refl, pdf);
  }
}

Vec radiance(const Ray &r, unsigned short *Xi){
  vector<Vertex> path;
  traceScene(r, Xi, &path);
  if (path.empty())
    return Vec(0);

  constexpr int LIGHT_ID = numSpheres - 1;
  const Sphere &light = spheres[LIGHT_ID];
  // 光源の表面をランダムにサンプリング
  Vec y_n = random_unit_direction(Xi);
  Vec y = light.p + y_n * light.rad;
  double light_area = 4 * M_PI * light.rad * light.rad;
  double pA_light = 1 / light_area; // 光源上の一様点サンプリング面積PDF

  Vec color(0);
  Vec throughput(1); // パスの累積スループット比 W_i = f(x_0...x_i) / p(x_0...x_i)

  for (size_t i = 0; ; ++i) {
    const Vertex& vi = path[i];

    // --- 1. 光源に到達した場合 ---
    if (vi.e.x + vi.e.y + vi.e.z > 0) {
      // ここで処理されるのは視点から直接またはスペキュラ反射で光源に辿り着く場合のみ（MISなし）。
      color = color + throughput.mult(vi.e);
      break; // 光源に到達したため以降の延伸・接続処理を終了
    }

    // --- 2. 光源上の点 y をサンプリングして接続する直接光計算 (Next Event Estimation) ---
    // (拡散反射面 DIFF の場合)
    if (vi.refl == DIFF) {
      Vec wi = y - vi.p;
      double dist = wi.length();
      wi = wi * (1 / dist);

      double cos_out = vi.n.dot(wi);
      double cos_in = -y_n.dot(wi);
      if (cos_out > 0 && cos_in > 0 &&
          is_light_visible(Ray(vi.p, wi), dist)) {
        // 接続エッジ (vi.p, y) の幾何因子 G(vi.p, y)
        double G = (cos_out * cos_in) / (dist * dist);
        double bsdf = 1 / M_PI;

        // 経路積分形式における接続パス (x_0...x_m, y) の寄与評価:
        // f(x)/p(x) = throughput * bsdf * G * L_e(y) / pA_light
        Vec L = vi.c.mult(light.e) * (bsdf * G / pA_light);
        // この接続点 y を BSDF サンプリングで選んだ場合の面積 PDF
        double pA_bsdf1 = vi.pdf * (cos_in / (dist * dist));
        double weight = mis_weight(pA_light, pA_bsdf1);
        color = color + throughput.mult(L) * weight;
      }
    }

    if (i + 1 >= path.size())  // 次の頂点がない：終了
      break;

    // --- 3. 視点サブパスを次の頂点へ延伸するための累積スループットの更新 ---
    // (経路積分形式 Path Integral Formulation による明示的な幾何因子 G と 面積 PDF pA の計算)
    if (vi.refl == DIFF) {
      const Vertex& vj = path[i + 1];
      Vec wi = vj.p - vi.p;
      double dist = wi.length();
      wi = wi * (1 / dist);

      double cos_out = vi.n.dot(wi);
      double cos_in = -vj.n.dot(wi);

      // エッジ (vi, vj) の幾何因子 G(vi, vj)
      double G = (cos_out * cos_in) / (dist * dist);
      // 頂点 v における BSDF f_r = 1 / pi
      double bsdf = 1 / M_PI;

      // 方向 PDF p_w(wi) = (cos_out / pi) * pdf
      // 角度積分 d_w から面積測度 dA への変換による面積 PDF pA(vj) の計算:
      double pA = vi.pdf * (cos_in / (dist * dist));

      // 経路積分形式におけるスループットの更新比: (bsdf * G) / pA
      throughput = throughput.mult(vi.c * (bsdf * G / pA));

      // 次の頂点が光源にたどり着いた場合：MIS
      if (vj.e.x + vj.e.y + vj.e.z > 0) {
        Vec wi = vj.p - vi.p;
        double dist = wi.length();
        wi = wi * (1 / dist);

        // double cos_out = vi.n.dot(wi);
        double cos_in = -vj.n.dot(wi);
        // BSDF サンプリングにおける面積 PDF: pA_bsdf(v) = (cos_out / pi * vi.pdf) * (cos_in / dist^2)
        double pA_bsdf2 = vi.pdf * (cos_in / (dist * dist));
        double weight = mis_weight(pA_bsdf2, pA_light);
        color = color + throughput.mult(vj.e) * weight;
        break; // 光源に到達したため以降の延伸・接続処理を終了
      }
    } else {
      // Specular / Refractive (Delta BSDF)
      // デルタ関数の評価により (bsdf * G) / pA = vi.c / vi.pdf
      throughput = throughput.mult(vi.c * (1 / vi.pdf));
    }
  }

  return color;
}