パストレーサーで直接光計算（重点的サンプリング、合成pdf）

Vec radiance(const Ray &r, int depth, unsigned short *Xi,int E=1){
  ...
  // 拡散反射面の箇所
    // Loop over any lights
    Vec e;
    for (int i=0; i<numSpheres; i++){
      const Sphere &s = spheres[i];
      if (s.e.x<=0 && s.e.y<=0 && s.e.z<=0) continue; // skip non-lights

      Vec sw=s.p-x, su=((fabs(sw.x)>.1?Vec(0,1):Vec(1))%sw).norm(), sv=sw%su;
      double cos_a_max = sqrt(1-s.rad*s.rad/(x-s.p).dot(x-s.p));
      double eps1 = erand48(Xi), eps2 = erand48(Xi);
      double cos_a = 1-eps1+eps1*cos_a_max;
      double sin_a = sqrt(1-cos_a*cos_a);
      double phi = 2*M_PI*eps2;
      Vec l = su*cos(phi)*sin_a + sv*sin(phi)*sin_a + sw*cos_a;
      l.norm();
      if (intersect(Ray(x,l), t, id) && id==i){  // shadow ray
        double omega = 2*M_PI*(1-cos_a_max);
        e = e + f.mult(s.e*l.dot(nl)*omega)*M_1_PI;  // 1/pi for brdf
      }
    }

    return obj.e*E+e+f.mult(radiance(Ray(x,d),depth,Xi,0));

struct Rectangle {
  Vec p, a1, a2, n; double area;
  Vec e, c;         // emission, color
  Rectangle(const Vec& p_, const Vec& a1_, const Vec& a2_, const Vec& e_, const Vec& c_)
    : p(p_), a1(a1_), a2(a2_), e(e_), c(c_) {
    n = a1 % a2; area = n.length(); n = n * (1.0 / area);
  }
  double intersect(const Ray &ray) const { // returns distance, 0 if nohit
    const double eps=1e-4;
    Vec r = ray.o - p, alpha = ray.d % a2, beta = r % a1;
    double inv = 1.0 / alpha.dot(a1), t = beta.dot(a2) * inv, u, v;  // tはnanになり得る: 0.0 * inf = nan
    if (t >= eps) {
      if ((u = alpha.dot(r) * inv) >= 0 && u <= 1 &&
          (v = beta.dot(ray.d) * inv) >= 0 && v <= 1) {
        return t; }
    }
    return INFINITY;
  }
  Vec random_sample(unsigned short *Xi) const {
    // TODO: 正方形じゃない場合の一様サンプル方法
    double r1 = erand48(Xi), r2 = erand48(Xi);
    return p + a1 * r1 + a2 * r2;
  }
};
const Rectangle light(Vec(50-16,81.5,81.6-16), Vec(32,0,0), Vec(0,0,32), Vec(12,12,12), Vec());

Vec radiance(const Ray &r, int depth, unsigned short *Xi){
  double t;                               // distance to intersection
  Vec x, n, f, e; Refl_t refl;
  if (int id; intersect(r, t, id)) {
    const Sphere &obj = spheres[id];        // the hit object
    x=r.o+r.d*t; n=(x-obj.p).norm(); f=obj.c; e=obj.e; refl=obj.refl; }
  if (double t2 = light.intersect(r); t2 < t) {  // 光源の交差判定を追加
    t=t2; x=r.o+r.d*t; n=light.n; f=light.c; e=light.e; refl=DIFF; }
  if (!isfinite(t)) return Vec(); // if miss, return black
  Vec nl=n; if (n.dot(r.d)>=0) {nl=n*-1; e=Vec(0);}
  double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
  if (++depth>5) { if (erand48(Xi)<p) f=f*(1/p); else return e; } //R.R.
  if (refl == DIFF){                  // Ideal DIFFUSE reflection
    const double direct_prob = 0.5;  // 直接光のパスを選ぶ確率
    Vec nextdir;
    // サンプルする方向の選択
    if (erand48(Xi) < direct_prob) {  // 光源上の点を一様サンプリングし、その方向を候補とする
      Vec lp = light.random_sample(Xi); nextdir = (lp - x).norm();
    } else {  // 反射特性に従い、半球状をコサインに比例してサンプリング
      double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
      Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1,0):Vec(1,0,0))%w).norm(), v=w%u;
      nextdir = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
    }
    Ray nextray(x, nextdir);
    // 選んだ方向が「選ばれる確率密度」を逆算
    double pdf_light = 0;  // 光源サンプルのpdf
    double pdf_brdf = 0;   // 反射面のpdf
    if (double t3 = light.intersect(nextray); t3 < INFINITY) {
      double distance_squared = t3 * t3;
      double cosine = fabs(nextdir.dot(light.n));
      pdf_light = distance_squared / (cosine * light.area);
    }
    pdf_brdf = fmax(nextdir.dot(nl), 0.0) * M_1_PI;
    double pdf_value = direct_prob * pdf_light + (1 - direct_prob) * pdf_brdf;  // 合成pdf
    if (pdf_value <= 0.0)  return e;
    double scattering_pdf = fmax(nextdir.dot(nl), 0.0) * M_1_PI;
    return e + f.mult(radiance(nextray, depth, Xi)) * (scattering_pdf / pdf_value);
  // 以降、鏡面反射と透過材質の処理（元ソースに合流）
  // obj.~は適宜置き換えること
  } else ...

#include <time.h>

// Rectangle, lightなどは上のソースから

double random(double min, double max) {
  return (rand() * (max - min) / (RAND_MAX + 1.0)) + min;
}
Vec random_in_unit_sphere() {
  for (;;) {
    Vec point(random(-1.0, +1.0), random(-1.0, +1.0), random(-1.0, +1.0));
    if (point.squareLength() < 1.0)
      return point;
  }
}
Vec random_unit_vector() {
  return random_in_unit_sphere().norm();
}
int main(){
  srand(time(NULL));
  rand();
  Vec x(random(1.0, 99.0), random(0.0, 80), random(1.0, 99.0));
  printf("x(%.2f, %.2f, %.2f), light.area=%g\n", x.x, x.y, x.z, light.area);

  const long N = 100000000;
  double sum = 0;
  int n = 1;
  Vec nl(0, 1, 0);
  for (long i = 0; i < N; ++i) {
    Vec nextdir = random_unit_vector();
    if (nl.dot(nextdir) <= 0)
      nextdir = nextdir * -1;
    double pdf_value = 1.0 / (2 * M_PI);

    Ray nextray(x, nextdir);
    if (double t3 = light.intersect(nextray); t3 < INFINITY) {
      double distance_squared = t3 * t3;
      double cosine = fabs(nextdir.dot(light.n));
      double pdf_light = distance_squared / (cosine * light.area);
      sum += pdf_light / pdf_value;
    }

    if ((i + 1) == n) {
      printf("%10ld: result: %g\n", i+1, sum / n);
      n *= 10;
    }
  }
}