光源・視点双方からの経路を接続するパストレーサー（経路積分）

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, double initialPdfOmega, int maxDepth){
  double pdfOmega = initialPdfOmega;
  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;
    double pdfArea = isinf(pdfOmega) ? 1.0 : pdfOmega * -nl.dot(r.d) / (t * t);
    path->emplace_back(x, nl, obj.e, obj.c, obj.refl, pdfArea);
    if (obj.e.x + obj.e.y + obj.e.z > 0) {  // 光源に到達した場合は終了
      return;
    }
    ++depth;  // ロシアンルーレットなし
    if (obj.refl == DIFF) {                     // Ideal DIFFUSE reflection
      Vec d = cosine_weighted_random_hemisphere(nl, Xi);
      r = Ray(x, d);
      pdfOmega = 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);
      pdfOmega = INFINITY;  // デルタ関数の反射面を示す
    }
  }
}

void generateEyeSubpath(const Ray& r, unsigned short *Xi, vector<Vertex>* path, int n){
  traceScene(r, Xi, path, INFINITY, n);
}
void generateLightSubpath(unsigned short *Xi, vector<Vertex>* path, int n){
  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 d = cosine_weighted_random_hemisphere(y_n, Xi);
  double pdfOmega = y_n.dot(d) / M_PI;

  // 光源の始点を追加
  path->emplace_back(y, y_n, light.e, Vec(1), light.refl, pA_light);  // light.c の代わりにVec(1)

  if (n > 1)
    traceScene(Ray(y, d), Xi, path, pdfOmega, n - 1);
}

Vec updateThroughput(const Vec& throughput, const vector<Vertex>& path, int i, bool is_light) {
  if (i + 1 >= path.size())  // 次の頂点がない：終了
    return throughput;

  const Vertex& vi = path[i];
  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 = (is_light && i == 0) ? 1 : 1 / M_PI;  // 光源部分経路の[0]は光源上の頂点からの放射で反射ではないため、特別扱いする必要がある

    // 経路積分形式におけるスループットの更新比: bsdf * G
    return throughput.mult(vi.c * (bsdf * G));
  } else {
    // Specular / Refractive (Delta BSDF)
    return throughput.mult(vi.c);
  }
}

bool connectVertices(const Vertex& ys, int s, const Vertex& zt, Vec* out) {
  Vec wi = ys.p - zt.p;
  double dist = wi.length();
  wi = wi * (1 / dist);

  double cos_out = zt.n.dot(wi);
  double cos_in = -ys.n.dot(wi);
  if (cos_out <= 0 || cos_in <= 0 ||
      !is_visible(Ray(zt.p, wi), dist))
    return false;

  double G = (cos_out * cos_in) / (dist * dist);

  Vec y_brdf = (s == 0 ? Vec(1) : ys.c * (1 / M_PI));
  Vec z_brdf = zt.c * (1 / M_PI);
  *out = y_brdf.mult(z_brdf) * G;
  return true;
}

Vec radiance(const Ray &r, unsigned short *Xi){
  // 部分経路同士の接続だけでは L.*(DD).*E 以外が扱えないので、一定確率でナイーブなパストレを行う
  constexpr double pnaive=0.005;
  if (erand48(Xi) < pnaive) {  // ナイーブなパストレ：直接光を描画するため
    vector<Vertex> path;
    generateEyeSubpath(r, Xi, &path, 100);
    const Vertex& zt = path.back();
    Vec color = zt.e * (1 / pnaive);
    if (color.x + color.y + color.z > 0) {  // 光源にたどり着いた場合にのみ寄与が得られる
      for (size_t t = 0; t < path.size(); ++t) {
        if (t + 1 < path.size() && path[t].refl == DIFF && path[t + 1].refl == DIFF)
          return Vec(0);  // DIFFが2回続いているため部分経路接続で生成可能、なのでこちらでは除外
        color = updateThroughput(color, path, t, false) * (1 / path[t].pdf);
      }
    }
    return color;
  }

  // 乱数により経路長を選択
  constexpr double p = 0.4;
  size_t k = 2 + (size_t)fmin(floor(log(erand48(Xi)) / log(1 - p)), 100);  // 合計頂点数：幾何分布、短い経路長ほど発生しやすくする
  size_t ne = irand(1, k, Xi);  // 視点部分経路頂点数：1 <= ne <= k - 1
  size_t nl = k - ne;           // 光源部分経路頂点数：1 <= nl <= k - 1

  vector<Vertex> eyeSubpath;
  generateEyeSubpath(r, Xi, &eyeSubpath, ne);
  vector<Vertex> lightSubpath;
  generateLightSubpath(Xi, &lightSubpath, nl);
  if (nl != lightSubpath.size() ||
      ne != eyeSubpath.size())  // 望みの長さの経路が得られなかった
    return Vec(0);

  // 部分経路の終端頂点同士をつなげる(拡散反射面 DIFF 同士の場合のみ)
  const Vertex& ys = lightSubpath.back();  // 光源側終端頂点
  const Vertex& zt = eyeSubpath.back();    // 視点側終端頂点
  Vec color(0);
  if (Vec connect; ys.refl == DIFF && zt.refl == DIFF &&
      connectVertices(ys, nl - 1, zt, &connect)) {
    int valid = 1;
    // 光源部分経路のスループットと確率密度を計算
    const Vec& Le = lightSubpath[0].e;
    Vec lightThroughput(Le);
    double pdfLight = 1.0;
    for (size_t s = 0; s < nl; ++s) {
      lightThroughput = updateThroughput(lightThroughput, lightSubpath, s, true);
      pdfLight *= lightSubpath[s].pdf;
      if (s + 1 < nl && lightSubpath[s].refl == DIFF && lightSubpath[s + 1].refl == DIFF)
        ++valid;
    }

    // 視点部分経路のスループットと確率密度を計算
    Vec eyeThroughput(1);
    double pdfEye = 1.0;
    for (size_t t = 0; t < ne; ++t) {
      eyeThroughput = updateThroughput(eyeThroughput, eyeSubpath, t, false);
      pdfEye *= eyeSubpath[t].pdf;
      if (t + 1 < ne && eyeSubpath[t].refl == DIFF && eyeSubpath[t + 1].refl == DIFF)
        ++valid;
    }

    double pdfPathLen = pow(1 - p, k - 2) * p * (1 - pnaive);
    double validRate = (double)valid / (k - 1);
    color = connect.mult(lightThroughput).mult(eyeThroughput) * (1 / (pdfPathLen * pdfLight * pdfEye * validRate));
  }
  return color;
}