光线跟踪smallpt详解 (一)

smallpt详解 (一)

本文主要是为了记录我学习smallpt的过程。第一部分是关于整个程序的大致分解。第二部分主要说一下我对main函数中cx和cy的理解。第三部分主要说明光线跟踪的每一行代码所用到的数学知识。

1. 基本功能函数

1.1 随机函数

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define M_PI 3.1415

double erand48(unsigned short xsubi[3]) {
    return (double)rand() / (double)RAND_MAX;
}

erand48()函数是linux下的随机函数，可以取[0, 1]的随机浮点数，在windows下必须自己写一个类似功能的函数。

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); }   //四舍五入

这里的难点在于为什么要用pow(clamp(x), 1/ 2.2)，查了资料才知道是gamma校正方法，详情可以看
http://blog.csdn.net/candycat1992/article/details/46228771

1.2 Vec类

struct Vec
{
    double x, y, z;
    Vec(double x_ = 0, double y_ = 0, double z_ = 0) : x(x_), y(y_), z(z_) {}
    Vec operator+(const Vec& rhs) const {
        return Vec(x + rhs.x, y + rhs.y, z + rhs.z);
    }
    Vec operator-(const Vec& rhs) const {
        return Vec(x - rhs.x, y - rhs.y, z - rhs.z);
    }
    Vec operator*(double b) const {
        return Vec(x * b, y * b, z * b);
    }
    Vec mult(const Vec &rhs) const {
        return Vec(x * rhs.x, y * rhs.y, z * rhs.z);
    }
    double dot(const Vec &rhs) const { return x*rhs.x + y*rhs.y + z*rhs.z; }
    Vec cross(const Vec &rhs) {
        return Vec(y * rhs.z - z * rhs.y,
            z * rhs.x - x * rhs.z,
            x * rhs.y - y * rhs.x);
    }
    Vec& norm() { return *this = *this * (1 / sqrt(x*x + y*y + z*z)); }
}

1) 其中包括向量加减，数乘，标准化，点乘，叉乘，还有需要注意的mult函数，后来会提到。
2) 这个Vec类在之后的计算中既可以表达向量Vector，亦可以表达三维点Point，也用来保存颜色RGB的值。

1.3 Ray类

struct Ray
{
    Vec o, d;   // o is original point, d is direction of lighting
    Ray(Vec o_, Vec d_) : o(o_), d(d_) {}
};

射线方程，o是射线的起点，d代表方向

1.4 Sphere类

enum Refl_t { DIFF, SPEC, REFR};

struct Sphere
{
    double rad; // radius
    Vec position, emission, color;
    Refl_t refl;

    Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_) :
        rad(rad_), position(p_), emission(e_), color(c_), refl(refl_) {}
    double intersect(const Ray &r) const {  // return distance, 0 if nohit
        Vec op = position - r.o;
        double t, eps = 1e-4;
        double b = op.dot(r.d), det = b*b + rad*rad - op.dot(op);
        if (det < 0)    
            return 0;
        else    
            det = sqrt(det);
        return (t = b - det) > eps ? t : ((t = b + det) > eps ? t : 0);
    }
};

射线和球体/圆，求交原理图：

求交的基本原理就是将射线的参数方程代入到圆的函数中，求t的值。
1）将P(t) = O + tD 代入圆方程，会得到 t 的一元二次方程。
2）先求出Vec op，op是用球心p的坐标减去射线的起点，表示原理图中的(O - C)。
3）b = op.dot(r.d)指代 ” D * (O - C) ”
4）求det，这里要注意我们求的b和原理图中的b差了两倍，所以可以直接用

double det = b * b - op.dot(op) + rad * rad;

如果det<0说明无解，直接return 0。
否则求根号的det
5）最终的解有一个或两个，可能在 t = b - det，或者t = b + det中，选择t大于0并且两个中较小的t。

---

2. 绘制过程

2.1 放置球的位置

Sphere spheres[] = {
    Sphere(1e5, Vec(1e5 + 1, 40.8, 81.6), Vec(), Vec(.75, .25, .25), Refl_t::DIFF), //left
    Sphere(1e5, Vec(-1e5 + 99,40.8,81.6),Vec(),Vec(.25,.25,.75), Refl_t::DIFF),//Rght 
    Sphere(1e5, Vec(50,40.8, 1e5),     Vec(),Vec(.75,.75,.75), Refl_t::DIFF),//Back 
    Sphere(1e5, Vec(50,40.8,-1e5 + 170), Vec(),Vec(),           Refl_t::DIFF),//Frnt 
    Sphere(1e5, Vec(50, 1e5, 81.6),    Vec(),Vec(.75,.75,.75),Refl_t::DIFF),//Botm 
    Sphere(1e5, Vec(50,-1e5 + 81.6,81.6),Vec(),Vec(.75,.75,.75),Refl_t::DIFF),//Top 
    Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1,1,1)*.999, Refl_t::SPEC),//Mirr 
    Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(1,1,1)*.999, Refl_t::REFR),//Glas 
    Sphere(600, Vec(50,681.6 - .27,81.6),Vec(12,12,12),  Vec(), Refl_t::DIFF) //Lite 
};

1）用6个很大的球体当做平面（DIFF属性，只有漫反射），因为半径很大的话，你在近距离看起来，球面就很像一个平面。
作者这样做应该是为了避免去写平面求交，平面类等函数。
2）用1个球表示光源，就是Lite，1个Mirr球（完全反射），1个Glass球（折射和反射都有）

2.2 遍历所有的球，求交点

inline bool intersect(const Ray &r, double &t, int &id) {
    double n = sizeof(spheres) / sizeof(Sphere), d, inf = t = 1e20;
    for (int i = int(n); i--;) {
        if ((d = spheres[i].intersect(r)) && d < t) 
        {
            t = d;
            id = i;
        }
    }
    return t < inf;
}

此光线射出去，在所有的球体中求交点。
求出距离camera最近的交点，这就是待会要绘制在屏幕上的主要的点。

2.3 main函数

1）camera的位置是在（50， 52， 295.6）， 往z轴的负方向看。

int main(int argc, char *argv[])
{
    int w = 1024, h = 768, samps = argc == 2 ? atoi(argv[1]) / 4 : 10; // # samples 
    Ray cam(Vec(50, 52, 295.6), Vec(0, -0.042612, -1).norm()); // cam pos, dir 
    Vec cx = Vec(w*.5135 / h), cy = (cx.cross(cam.d)).norm()*.5135, r, *c = new Vec[w*h];

2）遍历每个像素点，用随机采样的方式求得要射出的光线的方向d。

for (int y = 0; y<h; y++) {                       // Loop over image rows
        fprintf(stderr, "\rRendering (%d spp) %5.2f%%", samps * 4, 100.*y / (h - 1));
        for (unsigned short x = 0, Xi[3] = { 0,0,y*y*y }; x<w; x++)   // Loop cols 
            for (int sy = 0, i = (h - y - 1)*w + x; sy<2; sy++)     // 2x2 subpixel rows 
                for (int sx = 0; sx<2; sx++, r = Vec()) {        // 2x2 subpixel cols 
                    for (int s = 0; s<samps; s++) {
                        double r1 = 2 * erand48(Xi), dx = r1<1 ? sqrt(r1) - 1 : 1 - sqrt(2 - r1);
                        double r2 = 2 * erand48(Xi), dy = r2<1 ? sqrt(r2) - 1 : 1 - sqrt(2 - r2);
                        Vec d = cx*(((sx + .5 + dx) / 2 + x) / w - .5) +
                            cy*(((sy + .5 + dy) / 2 + y) / h - .5) + cam.d;
                        r = r + radiance(Ray(cam.o + d * 140, d.norm()), 0, Xi)*(1. / samps);
                    } // Camera rays are pushed ^^^^^ forward to start in interior 
                    c[i] = c[i] + Vec(clamp(r.x), clamp(r.y), clamp(r.z))*.25;
                }
    }
    FILE *f = fopen("image.ppm", "w");         // Write image to PPM file. 
    fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
    for (int i = 0; i<w*h; i++)
        fprintf(f, "%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));

2.4光线跟踪递归函数

radiance函数在第二部分介绍。
以下是大致的分块。

1)判断是否相交，求交点，求表面法向

Vec radiance(const Ray &r, int depth, unsigned short *Xi) {
    double t;                               // distance to intersection 
    int id = 0;                               // id of intersected object 
    if (!intersect(r, t, id)) 
        return Vec(); // if miss, return black 
    const Sphere &obj = spheres[id];        // the hit object 
    Vec x = r.o + r.d*t, n = (x - obj.position).norm(); // calculate vector n,球面法向量
    Vec nl = n.dot(r.d) < 0 ? n : n*-1, f = obj.color;
    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||!p) 
        if (erand48(Xi)<p) 
            f = f*(1 / p); 
        else  
            return obj.emission; 

2）漫反射（DIFF）
如果材质是漫反射，那么就随机生成一个方向进行漫反射。

    if (obj.refl == DIFF) {                  // Ideal DIFFUSE reflection 
        double r1 = 2 * M_PI*erand48(Xi), r2 = erand48(Xi), r2s = sqrt(r2);
        Vec w = nl, u = ((fabs(w.x)>.1 ? Vec(0, 1) : Vec(1)).cross(w)).norm(), v = w.cross(u);  //w，v，u为正交基
        Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1 - r2)).norm();
        return obj.emission + f.mult(radiance(Ray(x, d), depth, Xi));
    }

3）镜面反射（材质为SPEC）
计算镜面反射的方向，然后继续递归

else if (obj.refl == SPEC)            // Ideal SPECULAR reflection 
        return obj.emission + f.mult(radiance(Ray(x, r.d - n * 2 * n.dot(r.d)), depth, Xi));

4）反射和折射（材质为REFR）
玻璃材质，有一部分光进行反射，有一部分光进行折射。
这里用到了轮盘赌方法。

Ray reflRay(x, r.d - n * 2 * n.dot(r.d));     // Ideal dielectric REFRACTION 由平行四边形的方法求得反射光的direction
    bool into = n.dot(nl)>0;                // Ray from outside going in? 
    double nc = 1, nt = 1.5, nnt = into ? nc / nt : nt / nc, ddn = r.d.dot(nl), cos2t;
    if ((cos2t = 1 - nnt*nnt*(1 - ddn*ddn))<0)    // Total internal reflection 
        return obj.emission + f.mult(radiance(reflRay, depth, Xi));
    Vec tdir = (r.d*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, Tr = 1 - Re, P = .25 + .5*Re, RP = Re / P, TP = Tr / (1 - P);
    return obj.emission + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette 
        radiance(reflRay, depth, Xi)*RP : radiance(Ray(x, tdir), depth, Xi)*TP) :
        radiance(reflRay, depth, Xi)*Re + radiance(Ray(x, tdir), depth, Xi)*Tr);

---

参考链接：
[1]https://drive.google.com/file/d/0B8g97JkuSSBwUENiWTJXeGtTOHFmSm51UC01YWtCZw/view

[2]http://www.kevinbeason.com/smallpt/#moreinfo