シンプルなライトトレーサー（smallptベース、ピンホールカメラ、NEE）

light_tracer.cpp

#include <math.h>   // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
#include <stdio.h>  //   Modified by tyfkda, changed to light tracer.

#define USE_OPENMP
#ifdef USE_OPENMP
#include <omp.h>
#endif

struct Vec {        // Usage: time ./smallpt 5000 && xv image.ppm
  double x, y, z;                  // position, also color (r,g,b)
  explicit Vec(double xyz_=0){ x=y=z=xyz_; }
  Vec(double x_, double y_, double z_){ x=x_; y=y_; z=z_; }
  Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); }
  Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); }
  Vec operator*(double b) const { return Vec(x*b,y*b,z*b); }
  Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); }
  double squareLength() const { return x*x+y*y+z*z; }
  double length() const { return sqrt(squareLength()); }
  Vec norm() const { return *this * (1/length()); }
  double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
  Vec operator%(const Vec&b) const {return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};
struct Ray { Vec o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} };
enum Refl_t { DIFF, SPEC, REFR };  // material types, used in radiance()
struct Sphere {
  double rad;       // radius
  Vec p, e, c;      // position, emission, color
  Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive)
  Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_):
    rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
  double intersect(const Ray &r) const { // returns distance, inf if nohit
    Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
    double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
    if (det<0) return INFINITY; else det=sqrt(det);
    return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : INFINITY);
  }
};
Sphere spheres[] = {//Scene: radius, position, emission, color, material
  Sphere(1e5, Vec(-1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left
  Sphere(1e5, Vec( 1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght
  Sphere(1e5, Vec(50,40.8,-1e5),     Vec(),Vec(.75,.75,.75),DIFF),//Back
  Sphere(1e5, Vec(50,40.8, 1e5+296), Vec(),Vec(),           DIFF),//Frnt
  Sphere(1e5, Vec(50,-1e5, 81.6),    Vec(),Vec(.75,.75,.75),DIFF),//Botm
  Sphere(1e5, Vec(50, 1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top
  Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1,1,1)*.999, SPEC),//Mirr
  Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(1,1,1)*.999, REFR),//Glas
  // Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12),  Vec(), DIFF),//Lite
  Sphere(1.5, Vec(50,81.6-16.5,81.6),Vec(4,4,4)*100,  Vec(), DIFF),//Lite
};
inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }
inline bool intersect(const Ray &r, double &t, int &id, double max_dist=INFINITY){
  size_t n=sizeof(spheres)/sizeof(Sphere); double d; t=max_dist;
  for(size_t i=n;i>0;) if(d=spheres[--i].intersect(r), d<t){t=d;id=i;}
  return t<max_dist;
}
struct Random {
  unsigned short Xi[3];
  Random(unsigned long seed) { Xi[0] = seed>>32; Xi[1] = seed>>16; Xi[2] = seed; }
  inline double uniform01() { return erand48(Xi); }
  Vec uniform_unit_direction() {
    double r1=2*M_PI*uniform01(), r2=1-2*uniform01();
    return Vec(cos(r1)*sqrt(1-r2*r2), sin(r1)*sqrt(1-r2*r2), r2);
  }
  Vec cos_weighted_hemisphere_direction(const Vec& normal) {
    double r1=2*M_PI*uniform01(), r2=uniform01(), r2s=sqrt(r2);
    Vec w=normal, u=((fabs(w.x)>.1?Vec(0,1,0):Vec(1,0,0))%w).norm(), v=w%u;
    return u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2);
  }
};
inline double Gschlick(const Vec& v, const Vec& n, double k) {
  double dot = n.dot(v);
  return dot / (dot * (1 - k) + k);
};
inline double G_GGX(const Vec& v, const Vec& l, const Vec& n, double k) {
  return Gschlick(v, n, k) * Gschlick(l, n, k);
};
struct CameraScreenParam {
  const Ray* cam;
  const Vec icx;
  const Vec icy;
  int w, h;
  CameraScreenParam(const Ray* cam_, const Vec& cx, const Vec& cy, int w, int h)
    : cam(cam_), icx(cx*(1/cx.dot(cx))), icy(cy*(1/cy.dot(cy))), w(w), h(h) {}
  inline int pinhole_sensor_index(const Vec& x, const Vec& nl) const {
    Vec V = cam->o-x; double dist = V.length(); V = V * (1.0 / dist);
    if (nl.dot(V) <= 0) return -1;
    double VdotZ = -V.dot(cam->d); if (VdotZ <= 0) return -1;
    Vec sensor_pos = V * (1.0 / VdotZ);
    int ix = floor((0.5 - sensor_pos.dot(icx)) * w);
    int iy = floor((0.5 - sensor_pos.dot(icy)) * h);
    if (ix < 0 || ix >= w || iy < 0 || iy >= h) return -1;
    double t; int id; if (intersect(Ray(x, V), t, id, dist)) return -1;
    return (h-iy-1)*w+ix;
  }
};
inline Vec connect_to_camera(
    const Vec& x, const Vec& nl, const Vec& n, enum Refl_t refl, const Vec& L, const Vec& fl,
    CameraScreenParam* camera_screen, Vec* c, Random* rng
) {
  const Ray& cam = *camera_screen->cam;
  Vec V = cam.o-x; double dist = V.length(); V = V * (1.0 / dist);  // カメラ方向ベクトルと距離
  Vec nextdir; double brdf = 0; bool refracted = false;
  switch (refl) {
  case DIFF:
    brdf = M_1_PI; // Diffuse BRDF
    nextdir = rng->cos_weighted_hemisphere_direction(nl);
    break;
  case 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=-L.dot(nl), cos2t;
      if ((cos2t=1-nnt*nnt*(1-ddn*ddn))>=0){
        Vec tdir = (L*(-nnt) - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
        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 (rng->uniform01()>=Re) {
          refracted = true;
          nextdir = tdir;

          // TODO: 透過材質のBSDF
          // 無理やりBlinn-Phongを使う
          Vec L = nl * (-2 * nl.dot(nextdir)) - nextdir; // 入射方向：屈折ベクトルの法線方向以外を反転したベクトル
          Vec H = (V + L).norm(); // ハーフベクトル
          const double shininess = 100;
          double coeff = (shininess + 1) / (4 * M_PI);  // 修正Phongの正規化係数：分母が2じゃなく4なのは不明
          double HdotN = fmax(0.0, H.dot(nl));
          brdf = coeff * pow(HdotN, shininess);
          break;
        }
        // フレネル反射：下の反射に合流
      } // else 全反射：下の反射に合流
    }
    // Fallthrough
  case SPEC:
    {
      nextdir = nl*2*nl.dot(L)-L;  // 反射ベクトル

      const double roughness = 0.3, alpha = roughness * roughness, alpha2 = alpha * alpha;
      Vec H = (V + L).norm(); // ハーフベクトル
      double HdotN = fmax(0.0, H.dot(nl));
      double LdotN = fmax(0.0, L.dot(nl));
      double VdotN = fmax(0.0, V.dot(nl));
      double D = (alpha * alpha) / (M_PI * pow((HdotN * HdotN) * (alpha2 - 1) + 1, 2));
      double k = alpha / 2;
      double G_1 = G_GGX(V, L, nl, k);  // こっちだとラフネスが低い場合に暗くなりすぎる？
      double G_2 = Gschlick(V, nl, k);  // こっちだとラフネスが高い場合に光源からの輪郭部がおかしい
      double G = (G_1 - G_2) * roughness + G_2;  // ブレンドさせてみる
      double F = 1;
      brdf = (F * D * G) / (4 * LdotN * VdotN);  // GGX
    }
    break;
  }

  if (int idx = camera_screen->pinhole_sensor_index(x, nl * (refracted?-1.0:1.0)); idx >= 0) {
    double ZdotV = -cam.d.dot(V);
    double VdotN = fabs(V.dot(nl));  // 透過の場合も考慮して絶対値が必要
    double G = VdotN / (dist * dist);
    double Gsensor = pow(1.0 / ZdotV, 3);
    Vec cc = fl * (brdf * G * Gsensor);
    c[idx] = c[idx] + cc;
  }
  return nextdir;
}
void shot_light(const Sphere& light, CameraScreenParam* camera_screen, Vec* c, Random* rng) {
  Vec n = rng->uniform_unit_direction(), nl = n;
  Vec x = n * light.rad + light.p;
  double light_area = 4 * M_PI * light.rad * light.rad;
  Vec fl = light.e * (light_area * M_PI);

  enum Refl_t refl = DIFF;
  Ray r(x, n * -1);

  for (int depth=0; ; ++depth){
    Vec L = r.d * -1;  // 入射方向ベクトル
    Vec nextdir = connect_to_camera(x, nl, n, refl, L, fl, camera_screen, c, rng);

    r = Ray(x, nextdir);
    double t;                               // distance to intersection
    int id=0;                               // id of intersected object
    if (!intersect(r, t, id)) break;        // if miss, return black
    const Sphere &obj = spheres[id];        // the hit object
    x=r.o+r.d*t; n=(x-obj.p).norm(); nl=n.dot(r.d)<0?n:n*-1; Vec f=obj.c;
    double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
    if (depth>2) { if (rng->uniform01()<p) f=f*(1/p); else break; } //R.R.
    fl = fl.mult(f);
    refl = obj.refl;
  }
}
Sphere& find_light() {
  for (auto& sphere : spheres)
    if (sphere.e.x + sphere.e.y + sphere.e.z > 0)
      return sphere;
  throw "No light in the scene";
}
int main(int argc, char *argv[]){
  int w=640, h=480; long samps = argc==2 ? atol(argv[1]) : 1;
  double fov = atan(0.5135 / 2) * 2;  // deg2rad(28.8)
  Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm());
  Vec cx=Vec(w*tan(fov/2)*2/h,0,0), cy=(cx%cam.d).norm()*(tan(fov/2)*2);
  CameraScreenParam camera_screen(&cam, cx, cy, w, h);
  const Sphere& light = find_light();

#if !defined(USE_OPENMP)
  const int num_threads = 1;
  Vec *c=new Vec[w*h];
#else
  const int num_threads = 16;
	omp_set_num_threads(num_threads);
  Vec **cc=new Vec*[num_threads];
  for (int i=0; i<num_threads; ++i)
    cc[i] = new Vec[w*h];
# pragma omp parallel for schedule(dynamic, 1)
#endif
  for (long s=0; s<samps; s++){
#ifdef USE_OPENMP
    Vec *c=cc[omp_get_thread_num()];
    if (omp_get_thread_num() == 0)
#endif
    {
      long ss = (s + 1) * num_threads;
      fprintf(stderr, "\rRendering (%ld/%ld sample) %5.2f%%", ss, samps, 100.*ss/samps);
    }
    Random rng(s);
    shot_light(light, &camera_screen, c, &rng);
  }
  fprintf(stderr, "\nDone\n");

  FILE *f = fopen("image.ppm", "w");
  fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
  const double sensitivity = pow(h / (2 * tan(fov / 2)), 2);
  const double scale = sensitivity / samps;
  for (int i=0; i<w*h; i++) {
#ifdef USE_OPENMP
    Vec pixel = cc[0][i];  // スレッド別に計算した結果を合算
    for (int j=1; j < num_threads; ++j)
      pixel = pixel + cc[j][i];
#else
    const Vec& pixel = c[i];
#endif
    fprintf(f,"%d %d %d ", toInt(pixel.x*scale), toInt(pixel.y*scale), toInt(pixel.z*scale));
  }
  fclose(f);
}