2013年7月31日水曜日

path tracing でレンダリング(c++)


4D Light Fieldのデータを取得するためにレンダリングソフトを作った。
参考にしたのはhttp://www.kevinbeason.com/smallpt/ 。（ただかなり見づらい）
あとRussian Rouletteについてはこの本を参考にした。

単に4D Light Fieldを得るだけならBlenderとかで視点を変えた画像を取ればいいけど、 Heterodyned Light Fieldを取得したい場合にも対応するにはシミュレーションソフトを自作しておいたほうがよいと思った。

画像はレンダリング結果。512x512で2000 samples/pixel。


 少し上記参考ソースとの相違点の説明。
  •  カメラは画角と場所を変えられるように新しくClassを作成した。 
  • ピクセルの同じ場所から光線を発射（２ｘ２分割しなかった） 
  • 光線の反射はすべてロシアンルーレットで決定 
  • Cornell Boxらしきものも自分で定義しなおした。 
  • Seedの要らないdrand48()を使用。（これはWindowsでは動作しなかった。）
一応自分で整理しなおしたので理解した。
何かあったら連絡ください。

次回はLytroのようなRefocusingをやってみたい。



#include 
#include 
#include 

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 Camera {
public:
 double w, h; // length
 int wres, hres; //resolution. pixel pithc = w/wres
 double a; //distance from pinhole to sensor
 Vec pos; //position @pinhole
 Vec dir; //direction
 Camera(double w_, double h_, int wres_, int hres, Vec pos_, Vec dir_, double a_);
 void move(Vec pos_);
};

Camera::Camera(double w_, double h_, int wres_, int hres_, Vec pos_, Vec dir_, double a_) {
 w = w_; h = h_; wres = wres_; hres = hres_;
 pos = pos_; dir = dir_;
 a = a_;
}

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(-27, -50+16.5, -4), Vec(),     Vec(1,1,1)*.999,  SPEC),//Mirr 
  Sphere(16.5,   Vec( 22, -50+16.5, 15), Vec(),     Vec(1,1,1)*.999,  REFR),//Glas 
  Sphere(300,    Vec(0.0, 349., -0.0),   Vec(3,3,3), 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-5 && 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[]){
 int samps = argc==2 ? atoi(argv[1]): 1; // # samples
   int wres = 512, hres = wres;
   Vec *c = new Vec[wres*hres];
   Camera cam(15.1, 15.1, wres, hres, Vec(3, 2.0, 155.0), Vec(0.0, 0.0, -1.0).norm(), 15);
  
   Vec cx = Vec(-cam.w, 0.0, 0.0), cy = (cx%cam.dir).norm()*cam.h;
   Vec r;
   int i=0;
//#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP
   for (int y=0; y < cam.hres; y++){                       // Loop over image rows
      fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps,100.*y/(cam.hres-1));
     for (int x=0; x < cam.wres; x++) {   // Loop cols
        r = Vec();
         for (int s=0; s < samps; s++){
             Vec d = cx*(0.5 - (double)x/cam.wres) + cy*(0.5 - (double)y/cam.hres) + cam.dir*cam.a;
             r = r + radiance(Ray(cam.pos, d.norm()), 0)*(1./samps);
          }
          i = (cam.hres-1-y)*cam.wres + x;
           c[i] = Vec(clamp(r.x),clamp(r.y),clamp(r.z));
        }
   }
   FILE *f = fopen("image.ppm", "w");         // Write image to PPM file.
   fprintf(f, "P3\n%d %d\n%d\n", cam.wres, cam.hres, 255);
   for (int i=0; i < cam.wres*cam.hres; i++) {
         fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
    }
}
ラベル: , , 時刻:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)