path tracing でレンダリング(c++) その２DOF (Depth Of Field)

先日作ったレンダラにレンズを加えてDOF（Depth Of Field）を再現してみた。
センサーとレンズの位置関係を決めるとピクセルと共役な点（集光する点）が求まる。
はじめにその共役点を求めて、絞りの各座標からそこへ向けて光線を飛ばす。

センサーの分割だけでなく絞りも分割するので、ピンホールカメラに対して計算時間がかなり増えてしまう。MacbookAirだとなかなか時間がかかる（
並列化もしてないし。）ファンも結構回る。

上の図は瞳分割4x4でに50sample/direction。下の図は瞳分割なしの50*16 sample/pixel

ところで家のデスクトップはCore i5の2.67GHｚだが、こっちのほうが遅い。
理由は多分32bitのGCCを使ったせいだろう。こんなにも速度が違うのか、初めて知った。

/*
pathtrace_basic: basic path tracing from smallpt
pathtrace_dof: introducing a camera model.
this subaperture
*/

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>
//#define drand48() (double(rand()) / RAND_MAX)

class Vec {
public:
 double x, y, z;

 // constructor
 Vec(double x_ = 0.0, double y_ = 0.0, double z_ = 0.0) {x = x_; y = y_; z = z_;} 


 //operator overloading
 Vec operator+(const Vec &b) { return Vec(x+b.x, y+b.y, z+b.z); }
 Vec operator-(const Vec &b) { return Vec(x-b.x, y-b.y, z-b.z); }
 Vec operator*(double b) { return Vec(b*x, b*y, b*z); }
 Vec operator%(const Vec &b) { return Vec(y*b.z - z*b.y, z*b.x-x*b.z, x*b.y-y*b.x);} // cross product
 Vec mult(const Vec &b) { return Vec(x*b.x, y*b.y, z*b.z); }
 Vec& norm() { return *this = *this * (1/sqrt(x*x + y*y + z*z)); }
 double dot(const Vec &b) { return x*b.x + y*b.y + z*b.z; }
};

class Ray {
public:
 Vec o, d;
 Ray(Vec o_, Vec d_) : o(o_), d(d_) {} // constructor
};

class Plane2D {
public:
 double x, y;
 int xres, yres;
 Plane2D(double x_, double y_, int xres_, int yres_): 
 x(x_), y(y_), xres(xres_), yres(yres_) {}
};

class Sensor : public Plane2D {
public:
 Vec *c; //color array
 Sensor(double x_, double y_, int xres_, int yres_): Plane2D(x_, y_, xres_, yres_) {
  c = new Vec[xres*yres];
 }
 ~Sensor(){
  delete c;
 }
};

class Lens : public Plane2D {
public:
 double f;
 Lens(double f_, double x_, double y_, int xres_, int yres_): Plane2D(x_, y_, xres_, yres_) {
  f = f_;
 } 
};

class Camera {
public:
 Sensor* snsr;
 Lens* lens;
 double a;   //distance from lens to sensor, suppose lens and snsr are parallel.
 double beta;  //magnification
 Vec pos;   //center of the lens
 Vec dir;
 Vec up;
 Vec hrz;

 //constructor
 Camera(Sensor& s, Lens& l, double a_, Vec pos_, Vec dir_, Vec up_) {
  snsr = &s; lens = &l;
  a = a_; pos = pos_; dir = dir_; up = up_;
  hrz = (up%dir).norm();
  beta = l.f / (a-l.f);
 }

 Vec conj(int xi, int yi);
 Vec lens_coord(int ui, int vi);
};

Vec Camera::conj(int xi, int yi) {
 //input coordinate on the sensor
 //returns coordinate of the conjugate position
 Vec d = hrz * (0.5 - (0.5+xi)/snsr->xres) * snsr->x + up * (0.5 - (0.5+yi)/snsr->yres) * (-snsr->y) + dir * a; //from pix to the lens' center
 return pos + d * beta;
}
Vec Camera::lens_coord(int ui, int vi) { 
 return pos + hrz*((0.5+ui)/lens->xres - 0.5)*lens->x + up*((0.5+vi)/lens->yres - 0.5) * (-lens->y);
}



enum Refl_type { DIFF, SPEC, REFR};

class Sphere {
public:
 double rad;
 Vec pos, emi, color;
 Refl_type refl;

 //constructor
 Sphere(double rad_, Vec pos_, Vec emi_, Vec color_, Refl_type refl_):
  rad(rad_), pos(pos_), emi(emi_), color(color_), refl(refl_) {}

 double intersect(const Ray &r) {
  //returns distance, 0 if missed
  Vec op = pos - r.o;
  double t, eps=1e-4;
  double b = op.dot(r.d);
  double det = b*b - op.dot(op) + rad*rad;
  if (det < 0) {
   return 0;
  } else { 
   det = sqrt(det);
  }

  return (t = b-det) > eps ? t : ((t = b + det) > eps ? t : 0);
 }
};

Sphere spheres[] = {//Scene: 
 //     radius, position,               emission,    color,     material 
   Sphere(1e5,    Vec(-1e5-50, 0.0, 0.0), Vec(),     Vec(.75,.25,.25), DIFF),//Left 
  Sphere(1e5,    Vec( 1e5+50, 0.0, 0.0), Vec(),     Vec(.25,.25,.75), DIFF),//Rght 
  Sphere(1e5,    Vec(0.0, 0.0, -1e5-50), Vec(),     Vec(.75,.75,.75), DIFF),//Back 
  //Sphere(1e5,    Vec(0.0, 0.0,  1e5+1000),Vec(),     Vec(),            DIFF),//Frnt 
  Sphere(1e5,    Vec(0.0, -1e5-50, 0.0), Vec(),     Vec(.75,.75,.75), DIFF),//Botm 
  Sphere(1e5,    Vec(0.0,  1e5+50, 0.0), Vec(),     Vec(.75,.75,.0), DIFF),//Top 
  Sphere(16.5,   Vec(-20, -50+16.5, -50), Vec(),     Vec(.99,0.0,.0),  DIFF),//
  Sphere(16.5,   Vec(  0, -50+16.5, 0), Vec(),     Vec(1,1,1)*.999,  REFR),//Glas 
 Sphere(16.5,   Vec( 20, -50+16.5, 50), Vec(),     Vec(.99, .99, .0),  DIFF),//
  Sphere(300,    Vec(0.0, 349., 50.0),   Vec(6,6,6), 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) {
 //loop over the spheres,
 //and store the closest distance and the id.
 int n = sizeof(spheres) / sizeof(Sphere);
 double d, inf = t = 1e20;

 for (int i = 0; i < n; i++) {
  if ((d=spheres[i].intersect(r))>1e-4 && d < t) {
   t = d;
   id = i;
  }
 }
 return t<inf;
}

Vec radiance(const Ray &r_, int depth_) {
 double t; //distance to intersection
 int id = 0;
 Ray r = r_;
 int depth = depth_;

 Vec color(0,0,0); //accumulated color
 Vec reflect(1,1,1); //accumulated reflectance

 while(1) {
  if (!intersect(r, t, id)) //if missed
   return color;

  const Sphere &obj = spheres[id];  // hit object
  Vec x = r.o + r.d * t;     // intersection point
  Vec n = (x-obj.pos).norm();   //surface normal
  Vec nl = n.dot(r.d) < 0 ? n : n*(-1); //dot becomes negative if ray is goin into the sphere.
  Vec oc = obj.color;     //object color (reflectance).

  color = color + reflect.mult(obj.emi);  //accumulate color if any emission
  //printf("%f, %f, %f", color.x, color.y, color.z);

  double p = oc.x > oc.y && oc.x > oc.z ? oc.x : oc.y > oc.z ? oc.y : oc.z; // max object reflectance
  
  //russian roulet
  if (drand48() < p) {
   oc = oc * (1/p); //normalize by maximum reflectnce
  } else {
   return color;
  }

  reflect = reflect.mult(oc); //accumulate reflectance

  if (obj.refl == DIFF) { //lambert diffusion
   double r1 = 2 * M_PI * drand48();
   double r2 = drand48(), r2s = sqrt(r2);

   // new axis @ a intersection point
   Vec w = nl;
   Vec u = ((fabs(w.x) > 0.1 ? Vec(0,1,0): Vec(1,0,0))%w).norm();
   Vec v = w % u;

   //random sample in a semisphere
   Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();

   r = Ray(x, d);
   continue;


  } else if (obj.refl == SPEC) { //specular reflection (like mirror)
   r = Ray(x, r.d - n*2*n.dot(r.d));
   continue;
  }

  //dielectric Refraction
  Ray reflRay(x, r.d - n*2*n.dot(r.d)); //reflection ray
  bool into = n.dot(nl) > 0;     // true if ray is goin into the object
  double ni = 1.0, nr = 1.5;
  double nir = into ? ni/nr : nr/ni;
  double ddn = r.d.dot(nl);
  double cos2t = 1 - nir*nir*(1-ddn*ddn);

  if (cos2t < 0) { //total reflection
   r = reflRay;
   continue;
  }
  Vec refr_dir = (r.d * nir - n*((into?1:-1) * (ddn*nir + sqrt(cos2t)))).norm(); //refracted ray
  double R0 = (nr-ni)*(nr-ni)/((nr+ni)*(nr+ni));  //normal reflectance
  double co = 1 - (into ? -ddn : refr_dir.dot(n)); 

  double Re = R0 + (1 - R0) * co * co * co * co * co; //sheered reflectance
  double Tr = 1 - Re;

  if (drand48() < Re) { 
   r = reflRay;  //reflected
  } else {
   r = Ray(x, refr_dir); //refracted
  }
  continue;
 }
}

int main(int argc, char *argv[]){
 srand(time(NULL));
 clock_t sclock, eclock;

 int samps = argc==2 ? atoi(argv[1]): 1; // # samples
   int xres = 256, yres = xres;
   int ures = 1, vres = ures;
   Sensor snsr(100., 100., xres, yres);
   double fl = 100;
   Lens lens(fl, fl/1.4, fl/1.4, ures, vres);
   Camera cam(snsr, lens, 2*fl, Vec(0.0, -0.1, 50+fl*2), Vec(0., 0., -1.), Vec(0.,1.,0.));
   int uv = cam.lens->xres * cam.lens->yres;
   Vec r;

//#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP
   sclock = clock(); 
   for (int y=0; y < cam.snsr->yres; y++){                       // Loop over image rows
      fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps,100.*y/(cam.snsr->yres-1));
     for (int x=0; x < cam.snsr->xres; x++) {   // Loop over img cols
      Vec pconj = cam.conj(x,y);
      r = Vec(); //radiance

      for (int v=0; v<cam.lens->yres; v++) {
       for (int u=0; u<cam.lens->xres; u++) {
        Vec p = cam.lens_coord(u,v);
        Vec d = (pconj - p).norm();//lens to the conj point
        
        for(int s=0; s<samps; s++) {
         r = r + radiance(Ray(p,d), 0) * (1./(samps*uv));
        }
       }
      }
      int i = (cam.snsr->xres-1-y)*cam.snsr->xres + x;
           cam.snsr->c[i] = Vec(clamp(r.x),clamp(r.y),clamp(r.z));
        }
   }
   eclock = clock();
   printf("time: %.2f1[s]\n",(double)(eclock-sclock)/CLOCKS_PER_SEC);

   FILE *f = fopen("image.ppm", "w");         // Write image to PPM file.
   fprintf(f, "P3\n%d %d\n%d\n", cam.snsr->xres, cam.snsr->yres, 255);
   for (int i=0; i< cam.snsr->xres * cam.snsr->yres; i++) {
     fprintf(f,"%d %d %d ", toInt(cam.snsr->c[i].x), toInt(cam.snsr->c[i].y), toInt(cam.snsr->c[i].z));
    }
}