先日作ったレンダラにレンズを加えて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));
}
}
.jpeg)
.jpeg)
