1. 程式人生 > >Peter Shirley-Ray Tracing The Next Week

Peter Shirley-Ray Tracing The Next Week

Peter Shirley-Ray Tracing The Next Week(2016)

原著:Peter Shirley

英文原著地址 密碼: urji

第二本書主要介紹了運動模糊,BVH(層次包圍盒),紋理貼圖,柏林噪聲,紋理對映,光照,instance,volumes,最後會渲染一張封面上的圖片。

因為機器計算能力問題,程式碼渲染的圖片解析度較小,放在The-Next-Week資料夾下,圖片使用的是原書的圖片。

目錄:

  • [x] Chapter1:Motion Blur
  • [x] Chapter2:Bounding Volume Hierarchies
  • [x] Chapter3:Solid Textures
  • [x] Chapter4:Perlin Noise
  • [x] Chapter5:Image Texture Mapping
  • [x] Chapter6:Rectangles and Lights
  • [x] Chapter7:Instances
  • [x] Chapter8:Volumes
  • [x] Chapter9:A Scene Test All New Features

Chapter1:Motion Blur

運動模糊。當你在進行ray tracing的時候,模糊反射和散焦模糊的過程中,每個畫素你需要取樣多個點,來決定最終畫素的顏色,這種效果在現實世界中是另外一種實現方法,現實世界中,相機通過控制快門的開和關,記錄下快門開閉時間內,物體運動的軌跡,通過這樣的方法實現模糊的效果。

運動模糊的意思是,現實世界中,相機快門開啟的時間間隔內,相機活著物體發生了位移,畫面最後呈現出來的畫素,是移動過程中畫素的平均值。我們可以通過隨機一條光線持續的時間,最後計算出畫素平均的顏色,這也是光線追蹤使用了很多隨機性的方法,最後的畫面接近真實世界的原因。

這種方法的基礎是當快門開機的時間段內,隨機時間點生成光線,修改之前的ray類,新增一個光線存在時間的變數。

// 增加時間資訊
class ray
{
public:
    ray(){}
    ray(const vec3& a, const vec3 & b){ A =a; B = b;}
    vec3 origin() const  { return  A;}
    vec3 direction() const { return  B;}
    vec3 point_at_parameter(float t) const { return A+t*B;}

    float time() const{ return  _time};
    vec3 A;
    vec3 B;
    // 光線的時間戳
    float _time;

};

接下來就是控制camera在時間t1和t2之間,隨機時間點生成光線,給camera類添加個float 的時間變數,記錄光線產生的時間。


class camera
{
    vec3 origin;
    vec3 u,v,w;
    vec3 horizontal;
    vec3 vertical;
    vec3 lower_left_corner;
    float len_radius;
    // 增加開始時間和結束時間
    float time0,time1;

public :
    // 建構函式增加t0,t1
    camera(vec3 lookfrom, vec3 lookat, vec3 vup, float vfov, float aspect, float aperture, float focus_dist,
    float t0,float t1)
    {
        time0 = t0;
        time1 = t1;
        len_radius = aperture/2;
        float theta = vfov*M_PI/180;
        float half_height = tan(theta/2);
        float half_width = aspect * half_height;
        origin = lookfrom;

        w = unit_vector(lookfrom - lookat);
        u = unit_vector(cross(vup, w));
        v = cross(w,u);

        lower_left_corner = origin - half_width*focus_dist*u - half_height*focus_dist*v - focus_dist*w;
        horizontal = 2*half_width*focus_dist*u;
        vertical = 2*half_height*focus_dist*v;
    }

    ray get_ray(float s,float t)
    {
        vec3 rd = len_radius * random_in_unit_disk();
        vec3 offset = u * rd.x() +v*rd.y();
        // 隨機時間戳的光線
        float time = time0 + drand48()*(time1 - time0);
        return ray(origin + offset,lower_left_corner+s*horizontal + t*vertical - origin - offset,time);
    }

};

關於如何移動物體,接下來建立一個移動球體的類,裡面儲存某個球在時間點t0從位置center0,移動到時間點t1的位置center1.這個時間間隔不用和相機快門的時間長短相同。

class moving_sphere:public  hitable
{
public:
    moving_sphere(){}
    moving_sphere(vec3 cen0,vec3 cen1,float t0,float t1,float r,material *m):
            center0(cen0),center1(cen1),time0(t0),time1(t1),radius(r),mat_ptr(m) {};

    virtual  bool hit(const ray&r,float tmin, float tmax, hit_record& rec) const;

    vec3 center(float time) const;
    vec3 center0,center1;
    float time0,time1;
    float radius;
    material *mat_ptr;
};

// 當前時間點,球心的位置
vec3 moving_sphere::center(float time) const {
    return center0 + ((time-time0)/(time1-time0))*(center1-center0);
}

重寫moving_sphere的hit函式,修改之前的center為一個時間相關的位置

bool moving_sphere::hit(const ray& r,float t_min,float t_max,hit_record & rec )const
{
    // 修改之前的center為一個時間相關的位置
    vec3 oc = r.origin() - center(r.time());
    float a = dot(r.direction(), r.direction());
    float b = dot(oc, r.direction());
    float c = dot(oc, oc) - radius*radius;
    float discriminant = b*b - a*c;
    if (discriminant > 0) {
        float temp = (-b - sqrt(discriminant))/a;
        if (temp < t_max && temp > t_min) {
            rec.t = temp;
            rec.p = r.point_at_parameter(rec.t);
            rec.normal = (rec.p - center(r.time())) / radius;
            rec.mat_ptr = mat_ptr;
            return true;
        }
        temp = (-b + sqrt(discriminant)) / a;
        if (temp < t_max && temp > t_min) {
            rec.t = temp;
            rec.p = r.point_at_parameter(rec.t);
            rec.normal = (rec.p - center(r.time())) / radius;
            rec.mat_ptr = mat_ptr;
            return true;
        }
    }
    return false;
}

最後修改上本書最後繪製的場景,小球在time=0的時候,在原來的位置,time=1的時候,移動到center+vec3(0,0.5*drand48(),0)位置,在此pre期間,光圈一直開啟。


hitable *random_scene() {
    int n = 500;
    hitable **list = new hitable *[n + 1];
    list[0] = new sphere(vec3(0, -700, 0), 700, new lambertian(vec3(0.5, 0.5, 0.5)));
    int i = 1;
    for (int a = -11; a < 11; a++) {
        for (int b = -11; b < 11; b++) {
            float choose_mat = drand48();
            vec3 center(a + 0.9 * drand48(), 0.2, b + 0.9 * drand48());
            if ((center - vec3(4, 0.2, 0)).length() > 0.9) {
                if (choose_mat < 0.8) {  // diffuse
                    // 運動模糊的小球
                    list[i++] = new moving_sphere(center, center + vec3(0, 0.5 * drand48(), 0), 0.0, 1.0, 0.2,
                                                  new lambertian(vec3(drand48() * drand48(), drand48() * drand48(),
                                                                      drand48() * drand48())));
                } else if (choose_mat < 0.95) { // metal
                    list[i++] = new sphere(center, 0.2,
                                           new metal(vec3(0.5 * (1 + drand48()), 0.5 * (1 + drand48()),
                                                          0.5 * (1 + drand48())), 0.5 * drand48()));
                } else {  // glass
                    list[i++] = new sphere(center, 0.2, new dielectric(1.5));
                }
            }
        }
    }

    list[i++] = new sphere(vec3(0, 1, 0), 1.0, new dielectric(2.5));
    list[i++] = new sphere(vec3(-4, 1, 0), 1.0, new lambertian(vec3(0.4, 0.2, 0.1)));
    list[i++] = new sphere(vec3(4, 1, 0), 1.0, new metal(vec3(1, 1, 1), 0.0));

    return new hitable_list(list, i);
}

camera類的get_ray函式返回了一條隨機時間t在t0和t1之間時間點的光線,這個時間t被用在moving_sphere中,決定了center球心的位置。在迴圈取樣ns的位置,不停的get_ray,不停的和隨機時間t位置的球求交,這樣就形成了動態模糊的效果。

最後渲染出來達到的效果如下:

Chapter2:Bounding Volume Hierarchies

層次包圍盒
第二章,是比較重要的一部分,層次包圍盒的出現,可以使我們的程式碼“跑的更快“,主要是通過重構hitable類,新增rectangles和boxes。

之前寫的ray tracing的複雜度是線性的,有多少調光線多少個物體,複雜度是線性相關。我們可能同時發出來幾百萬的光線,但其實這個過程我們可以通過二分查詢的思想來進行。這個過程分為2個關鍵的部分

  • 1)劃分空間
  • 2)劃分物體物件

關鍵的思想是使用bounding volume(包圍盒),包圍盒就是一個普通的立方體,這個立方體將物體完全包裹著。舉個簡單的例子,現在有10個物體,你用一個bounding sphere將他們包住,如果光線沒有射到這個包圍球,那麼肯定沒有射到這10個物體,如果光線射到了包圍球,再進行後面的判斷,虛擬碼如下:

if(ray hit bounding object)
    return whether ray hit bounded objects  // 是否擊中包圍內部的物體
else
    return false    

還有個關鍵的點是,如何劃分物體形成子集。實際上我們不是直接劃分螢幕活著volume的,每個物體都有一個bounding volume,而且bounding volume可以重疊。建立一個bounding volume的層級關係。舉個例子,我們將物體的總集分為紅藍2個子集,分別用一個bounding volume包圍起來,就有了下面的這張圖:

紅色和藍色都在紫色的包圍盒內,他們發生了重疊,就有了右邊的樹型結構,紅藍分別是紫的左右孩子,虛擬碼如下:

if(hit purple) // 紫色
    hit0 = hits blue enclosed objects
    hit1 = hits red enclosed objects
    if(hit0 or hit1)
        return true and info of closer hit         //返回hit的資訊
else
    return false

為了更好的效能,一個好的bounding volume結構是很有必要的,須要方便劃分,有要儘可能少的計算量,axis-aligned bounding boxes(AAABB)包圍盒就是一種很好的結構,我們只需要知道是否hit到了物體,不需要知道hit到的點,和法線。

很多人用一種叫“slab”的方法,這是一種基於n個緯度的AABB,就是從n個軸上取n個區間表示。3<x<5 , x in (3,5)這樣表示更加簡潔。
2D的時候,x,y2個區間可以現成一個矩形。

如果判斷一條光線是否射中一個區間,須要先判斷光線是否擊中分界線。在2d平面內,邊界是2條線,而這條ray有2個引數t0和t1,就可以在平面內確定一條光線;如果是在3d空間,邊界是2個面,假設為x=x0和x=x1(這是x方向上的2個面),對於時刻t,有個關於p(t)的函式

p(t)= A + tB

這個公式是適用與xyz三個座標系,比如

x(t) = Ax + t*Bx

當t0時刻,射線擊中平面的位置 x=x0 ,即

x0 = Ax +t0*Bx

求出來

t0 = (x0 - Ax) / Bx

同理

t1 = (x1 - Ax) / Bx  (當x = x1時

剛才聚的例子是1緯空間的,xy分開計算,但是要知道2緯空間,是否擊中物體,就要計算,x空間和y空間擊中的物體是否發生重疊,就想下圖中藍色和綠色表示x,y空間擊中物體的2個平面,4個平面重疊的部分。

虛擬碼如下:

    compute(tx0,tx1);
    compute(ty0,ty1);
    return overlap?((tx0,tx1),(ty0,ty1))

三維的時候就再加上z空間的判斷。

注意事項:

  • 對於求解出來的tx0和tx1,構成的區間可能是(7,3)這樣的形式,那麼就須要對tx0和tx1做下翻轉,轉成(3,7)

    tx0 = min((x0 - Ax)/Bx,(x1 - Ax)/Bx);
    tx1 = max((x0 - Ax)/Bx,(x1 - Ax)/Bx);

  • 如果除數是0,既Bx=0,或者分子是0,既(x0-Ax)=0或(x1-Ax)=0那麼求解出來的答案,求出來可能無意義,分子是0,表示只有一個解,等於光線是擦邊,不好界定是射中了還是沒有射中。

對於bvh在判斷重疊的方法,在一維平面內原理就是比較2個區間,看2個區間是否重疊,比如區間(d,D)和區間(e,E),計算出來的重疊區間為(f,F),如果滿足

// 計算是否重疊
bool overlap(d,D,e,E,f,F)
{
    f = max(d,e);
    F = min(D,e);
    return f<F;
}

我自己總結了下就是 左大右小(左區間區max,右區間取min,比較2個值,如果左<右,為真)發生重疊。

// aabb包圍盒
class aabb
{
public:
    aabb(){}
    aabb(const vec3 a,const vec3 &b)
    {
        _min = a;_max = b;
    }

    vec3 min()const{ return  _min};
    vec3 max()const{ return  _max};

    bool hit(const ray& r,float tmin,float tmax)const
    {
        for(int a =0;a<3;a++)
        {
            float invD = 1.0f/r.direction()[a];
            float t0 = (min()[a] - r.direction()[a]) * invD;
            float t1 = (max()[a] - r.direction()[a]) * invD;

            if(invD<0.0f)
                std::swap(t0,t1);
            tmin = t0>tmin?t0:tmin;
            tmax = t1<tmax?t1:tmax;
            if(tmax <= tmin)
                return false;
        }
        return true;
    }
    vec3 _min;
    vec3 _max;
};

對於hitable的類,須要加一個bounding_box的虛擬函式,方便派生類實現

class hitable
{
public:
    virtual bool hit(const ray& r,float t_min,float t_max,hit_record & rec)const =0;
    virtual bool bounding_box(float t0,float t1,aabb & box)const =0;
};

之前寫的球類實現bounding_box的函式,球的boundingbox很簡單,就是球心加半徑。

bool sphere::bounding_box(float t0, float t1, aabb &box) const {
    box = aabb(center - vec3(radius, radius, radius), center + vec3(radius, radius, radius));
    return true;
}

對於動態運動的球,對t0時刻的box和t1時刻的box,取一個更大的boundingbox

aabb moving_sphere::surrounding_box(aabb &box0, aabb &box1) const {
    vec3 small(fmin(box0.min().x(), box1.min().x()),
               fmin(box0.min().y(), box1.min().y()),
               fmin(box0.min().z(), box1.min().z()));
    vec3 big(fmax(box0.max().x(), box1.max().x()),
             fmax(box0.max().y(), box1.max().y()),
             fmax(box0.max().z(), box1.max().z()));
    return aabb(small,big);
}

hitable類也要加點東西,因為BVH涉及左右子樹,所以以連結串列的形式,新增左右子樹。

class bvh_node:public hitable
{
public:
    bvh_node(){}
    bvh_node(hitable **l,int n,float time0,float time1);
    virtual bool hit(const ray&r,float tmin,float tmax,hit_record &rec)const;
    virtual bool bounding_box(float t0,float t1,aabb &box) const;

    hitable *left;
    hitable *right;
    aabb box;
};

對於左右子樹進行遞迴操作,直到射到葉子節點,擊中重疊的部分,擊中的資料用引用rec傳出去。


bool bvh_node::hit(const ray &r, float tmin, float tmax, hit_record &rec) const {
    if(box.hit(r,tmin,tmax))
    {
        hit_record left_rec,right_rec;
        bool hit_left = left->hit(r,tmin,tmax,left_rec);
        bool hit_right = right->hit(r,tmin,tmax,right_rec);
        if(hit_left && hit_right)           // 擊中重疊部分
        {
            if(left_rec.t<right_rec.t)
                rec = left_rec;             // 擊中左子樹
            else
                rec = right_rec;            // 擊中右子樹
            return true;
        } else if(hit_left)
        {
            rec = left_rec;
            return  true;
        } else if(hit_right)
        {
            rec = right_rec;
            return true;
        } else
            return false;
    } else
        return false;                       // 未擊中任何物體
}

這種bvh的結構,bvh_node節點記錄了擊中子類的record資訊,而且是一種二分的結構。如果bounding box劃分的合理是很高效的,最完美的是滿二叉樹的情況。
boundingbox的hitrecord的資料,使用qsort函式重寫compare函式來進行排序。

int box_x_compare(const void *a,const void *b)
{
    aabb box_left,box_right;
    hitable *ah = *(hitable**)a;
    hitable *bh = *(hitable**)b;
    if(!ah->bounding_box(0,0,box_left) || !bh->bounding_box(0,0,box_right))
        std::cerr <<"No bounding box in bvh_node constructor\n";
    if(box_left.min().x() - box_right.min().x()<0.0)
        return  -1;
    else
        return 1;
}

Chapter3:Solid Textures

固體紋理。紋理在圖形學中表示為一個面片上有關聯的不同的顏色,這裡須要完善一個texture的類,表現物體表面的紋理。


class texture {
public:
    virtual vec3 value(float u, float v, const vec3 &p) const = 0;
};

class constant_texture : public texture {
public:
    constant_texture() {}

    constant_texture(vec3 c) : color(c) {}

    virtual vec3 value(float u, float v, const vec3 &p) const {
        return color;
    }

    vec3 color;
};

這樣就可以使用texture類,通過uv取樣紋理顏色的方法,替換之前寫的vec3 的color。比如把之前的lambertain材質重寫,使用新的texture來表現顏色。


class lambertian : public material {
public:
    lambertian(texture *a) : albedo(a) {}
    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const  {
        vec3 target = rec.p + rec.normal + random_in_unit_sphere();
        scattered = ray(rec.p, target-rec.p);
        attenuation = albedo->value(0,0,rec.p);
        return true;
    }

    texture *albedo;    
};

新建一個lambertain材質,同時新建一個checker_texture(棋盤紋理),繼承普通的紋理,不過包含2個指標分別指向棋盤的間隔顏色。


// 棋盤紋理
class checker_texture:public texture
{
public:
    checker_texture(){}
    checker_texture(texture *t0,texture *t1):even(t0),odd(t1){}
    virtual vec3 value (float u,float v, const vec3 &p)const {
        float sines = sin(10*p.x())*sin(10*p.y())*sin(10*p.z());
        if(sines<0)
            return odd->value(u,v,p);
        else
            return  even->value(u,v,p);
    }
    // 棋盤紋理的間隔顏色
    texture *odd;
    texture *even;
};

更新main函式中的vec3的color,使用新的texture紋理,注意lambertain材質的建構函式改成 texture的指標了,之前是一個v3的物件。

      // 棋盤紋理
    texture *checker = new checker_texture(new constant_texture(vec3(0.2, 0.3, 0.1)),
                                           new constant_texture(vec3(0.9, 0.9, 0.9)));
    list[0] = new sphere(vec3(0, -700, 0), 700, new lambertian(checker));

最後得到的圖案就是把底部大球的紋理,改成了棋盤紋理,效果如下:

Chapter4:Perlin Noise

Perlin噪聲 ( Perlin noise )指由Ken Perlin發明的自然噪聲生成演算法 。

柏林噪聲有2個關鍵的部分,第一是輸入相同的3D點,總能返回相同的隨機值,第二是使用一些hack的方法,達到快速近似的效果。

noise函式,通過傳入一個三維空間的點,返回一個float型別的噪聲值。


#include "vec3.h"

class perlin {
public:
    float noise(const vec3 &p) const {
        float u = p.x() - floor(p.x());
        float v = p.y() - floor((p.y()));
        float z = p.z() - floor(p.z());
        int i = int(4 * p.x()) & 255;
        int j = int(4 * p.y()) & 255;
        int k = int(4 * p.z()) & 255;
        return ranfloat[perm_x[i] ^ perm_y[j] ^ perm_z[k]];
    }

    static float *ranfloat;
    static int *perm_x;
    static int *perm_y;
    static int *perm_z;
};

static float *perlin_generate() {
    float *p = new float[256];
    for (int i = 0; i < 256; i++) {
        p[i] = drand48();
    }
    return p;
}

// 改變序列函式
void permute(int *p, int n) {
    for (int i = n - 1; i > 0; i--) {
        int target = int(drand48() * (i + 1));
        int tmp = p[i];
        p[i] = p[target];
        p[target] = tmp;
    }
}

static int *perlin_generate_perm() {
    int *p = new int[256];
    for (int i = 0; i < 256; i++) {
        p[i] = i;
    }
    permute(p, 256);
    return p;
}


float *perlin::ranfloat = perlin_generate();
int *perlin::perm_x = perlin_generate_perm();
int *perlin::perm_y = perlin_generate_perm();
int *perlin::perm_z = perlin_generate_perm();

在texture標頭檔案中,新增生成噪聲紋理的程式碼,通過在0-1之間選取float,建立噪聲紋理

// 噪聲紋理
class noise_texture:public texture{
public:
    noise_texture(){}
    noise_texture(float sc):scale(sc){}
    virtual vec3 value(float u,float v,const vec3& p)const
    {
        return vec3(1,1,1)*0.5*(1+sin(scale*p.x())+ 5*noise.noise(p));
    }
    perlin noise;
    float scale;
};

在lambertian的球上應用噪聲紋理

hitable *two_perlin_spheres()
{
    texture *pertext = new noise_texture();
    hitable **list = new hitable*[2];
    list[0] = new sphere(vec3(0,-1000,0),1000,new lambertian(pertext));
    list[1] = new sphere(vec3(0,2,0),2,new lambertian(pertext));
    return new hitable_list(list,2);
}

得到的效果如下

再使紋理變得平滑一些,使用線性插值的方法:

inline float trilinear_interp(float cp[2][2][2], float u, float v, float w) {
    float accum = 0;
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; ++j) {
            for (int k = 0; k < 2; ++k) {
                accum += (i * u + (1 - i) * (1 - u)) * (j * v + (1 - j) * (1 - v)) * (k * w + (1 - k) * (1 - w)) *
                         cp[i][j][k];
            }
        }
    }
    return accum;
}

效果如下:

為了達到更好的平滑效果,使用hermite cubic方法去做平滑。

    float noise(const vec3 &p) const {
        float u = p.x() - floor(p.x());
        float v = p.y() - floor((p.y()));
        float w = p.z() - floor(p.z());
        // hermite cubic 方法平滑
        u = u*u*(3-2*u);
        v = v*v*(3-2*v);
        w = w*w*(3-2*w);
        int i = floor(p.x());
        int j = floor(p.y());
        int k = floor(p.z());

效果如下:

同時縮放輸入的點p來使噪聲變化的更快。

class noise_texture:public texture{
public:
    noise_texture(){}
    noise_texture(float sc):scale(sc){}
    virtual vec3 value(float u,float v,const vec3& p)const
    {
        return vec3(1,1,1)*noise.noise(scale * p);
    }
    perlin noise;
    float scale;
};

得到的效果如下:

現在仍然能看到網格,因為這種模式下,max和min總是收到具體的xyz值影響,然後Ken就用了一個trick,使用隨機的vectors替代原來的floats,通過點乘的方法改變格子上的max和min值。

// perlin 插值
inline float perlin_interp(vec3 c[2][2][2],float u,float v, float w)
{
    float uu =u*u*(3-2*u);
    float vv = v*v*(3-2*v);
    float ww = w*w*(3-2*w);
    float accum = 0;
    for (int i = 0; i < 2; ++i) {
        for (int j = 0; j < 2; ++j) {
            for (int k = 0; k < 2; ++k) {
                vec3 weight_v(u-i,v-j,w-k);
                accum += (i*uu + (1-i)*(1-uu))*
                        (j*vv +(1-j)*(1-vv))*
                        (k*ww +(1-k)*(1-ww))*dot(c[i][j][k],weight_v);
            }

        }
    }

    return accum;
}

再新增turb擾動的噪聲,達到更加自然的效果

    // 噪聲擾動
    float turb(const vec3& p, int depth=7) const {
        float accum = 0;
        vec3 temp_p = p;
        float weight = 1.0;
        for (int i = 0; i < depth; i++) {
            accum += weight*noise(temp_p);
            weight *= 0.5;
            temp_p *= 2;
        }
        return fabs(accum);
    }

並在texture的noise紋理中應用

// 噪聲紋理
class noise_texture:public texture{
public:
    noise_texture(){}
    noise_texture(float sc):scale(sc){}
    virtual vec3 value(float u,float v,const vec3& p)const
    {
        // 加縮放和擾動後
        return vec3(1,1,1)*0.5*(1 + sin(scale*p.x() + 5*noise.turb(scale*p))) ;
    }
    perlin noise;
    float scale;
};

最終達到的效果如下

補充下Perlin Noise的擴充套件閱讀,Building Up Perlin Noise

Chapter5:Image Texture Mapping

紋理對映,通過讀取一張圖片,使用uv對映的方法,直接將一張圖片的紋理繪製在物體表面。

直接的方法是縮放uv,uv是[0,1]之間的float。而畫素肯定大於這個區間,所以需要進行縮放,用(i,j)表示當前畫素,nx和ny表示紋理的解析度,所以對於任意畫素(i,j)位置,對應的uv座標就是

u = i / (nx - 1)
v = j / (ny - 1)

這種是對於平面座標的uv對映,如果是一個球體的話,使用極座標可以更方便的表示對映關係

u = phi / (2*Pi)
v = theta / Pi

通過hitpoint的xyz,可以計算出theta 和phi,對於單位球體,他們之間的關係如下

x = cos(phi)cos(theta)
y = sin(phi)cos(theta)
z = sin(theta)

然後math.h中提供了atan2()方法,可以計算反三角函式

phi = atan2(y,x)

atan2返回的值是在(-Pi,Pi)之間的

theta = asin(z)

theta值在(-Pi/2,Pi/2)之間。

最終就在hit 檔案中寫了一個獲取球體uv的函式

void get_sphere_uv(const vec3& p, float& u, float& v) {
    float phi = atan2(p.z(), p.x());
    float theta = asin(p.y());
    u = 1-(phi + M_PI) / (2*M_PI);
    v = (theta + M_PI/2) / M_PI;
}

以及使用stb_image從圖片讀取rgb的標頭檔案


class image_texture : public texture {
public:
    image_texture() {}
    image_texture(unsigned char *pixels, int A, int B) : data(pixels), nx(A), ny(B) {}
    virtual vec3 value(float u, float v, const vec3& p) const;
    unsigned char *data;
    int nx, ny;
};

vec3 image_texture::value(float u, float v, const vec3& p) const {
    int i = (1- u)*nx;
    int j = (1-v)*ny-0.001;
    if (i < 0) i = 0;
    if (j < 0) j = 0;
    if (i > nx-1) i = nx-1;
    if (j > ny-1) j = ny-1;
    float r = int(data[3*i + 3*nx*j]  ) / 255.0;
    float g = int(data[3*i + 3*nx*j+1]) / 255.0;
    float b = int(data[3*i + 3*nx*j+2]) / 255.0;
    return vec3(r, g, b);
}

須要注意main函式中

// 須要先宣告巨集,不然stb_image 會報錯找不到圖片格式
#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"

Chapter6:Rectangles and Lights

矩形和光照。如何做一個自發光的材質,首先須要在hit_record裡面加一個 emitted的方法。比如說背景如果是純黑的話,就相當於光線來了的時候,他不反射任何光線。

// 自發光材質
class diffuse_light:public material
{
public:
    diffuse_light(texture *a):emit(a){}
    virtual bool scatter(const ray& r_in,const hit_record &rec,vec3 & attenuation,ray& scattered)const {
        return false;
    }
    virtual vec3 emitted(float u,float v,const vec3 &p)const {
        return emit->value(u,v,p);
    }
    texture *emit;
};

材質類須要加個emitted的虛擬函式,預設return的是黑色。方便子類重寫

class material  {
public:
    // 散射虛擬函式
    // 引數:r_in 入射的光線, rec hit的記錄, attenuation v3的衰減,scattered 散射後的光線
    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const = 0;
    // 非自發光材質,預設返回黑色
    virtual vec3 emitted(float u,float v,const vec3 &p)const {
        return vec3(0,0,0);
};

接下來寫一個rect的類,用來表示空間中的矩形。

以xy平面為例,在z=k的情況下,用2條直線,滿足x=x0,x=x1,y=y0,y=y1可以得到一個區域。

當判斷ray是否擊中這個rectangle的時候,ray的表示式為:

p(t) = a + t*b

在xy平面上,等價於:

z(t) = az + t*bz

解這個關於t的方程,當z=k的時候

t = (k - az) / bz

在知道t之後,

x = ax + t * bx
y = ay + t * by

如果滿足 x在區間[x0,x1],y在[y0,y1]上的話,ray就擊中了這個rect。

程式碼如下:

// xy平面的矩形
class xy_rect: public hitable  {
public:
    xy_rect() {}
    xy_rect(float _x0, float _x1, float _y0, float _y1, float _k, material *mat) : x0(_x0), x1(_x1), y0(_y0), y1(_y1), k(_k), mp(mat) {};
    virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
    virtual bool bounding_box(float t0, float t1, aabb& box) const {
        box =  aabb(vec3(x0,y0, k-0.0001), vec3(x1, y1, k+0.0001));
        return true; }
    material  *mp;
    float x0, x1, y0, y1, k;
};

具體實現如下:

// 是否擊中,形參傳了hit_record的引用。
bool xy_rect::hit(const ray& r, float t0, float t1, hit_record& rec) const {
    float t = (k-r.origin().z()) / r.direction().z();
    if (t < t0 || t > t1)
        return false;
    float x = r.origin().x() + t*r.direction().x();
    float y = r.origin().y() + t*r.direction().y();
    if (x < x0 || x > x1 || y < y0 || y > y1)
        return false;
    rec.u = (x-x0)/(x1-x0);
    rec.v = (y-y0)/(y1-y0);
    rec.t = t;
    rec.mat_ptr = mp;
    rec.p = r.point_at_parameter(t);
    rec.normal = vec3(0, 0, 1);
    return true;
}

在場景中放個rect做為光源

// 帶rect和光源的場景
hitable *simple_light()
{
    texture *pertext = new noise_texture(4);
    texture *checker = new checker_texture(new constant_texture(vec3(0.2, 0.3, 0.1)),
                                           new constant_texture(vec3(0.9, 0.9, 0.9)));
    hitable **list = new hitable*[4];
    list[0] = new sphere(vec3(0,-700,0),700,new lambertian(checker));
    list[1] = new sphere(vec3(0,2,0),2,new lambertian(pertext));
    list[2] = new sphere(vec3(0,7,0),2,new diffuse_light(new constant_texture(vec3(4,4,4))));
    list[3] = new xy_rect(3,5,1,3,-2,new diffuse_light(new constant_texture(vec3(4,4,4))));
    return new hitable_list(list,4);
}

得到如下的圖片

接下來補全yz平面和xz平面的程式碼

class xz_rect: public hitable  {
public:
    xz_rect() {}
    xz_rect(float _x0, float _x1, float _z0, float _z1, float _k, material *mat) : x0(_x0), x1(_x1), z0(_z0), z1(_z1), k(_k), mp(mat) {};
    virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
    virtual bool bounding_box(float t0, float t1, aabb& box) const {
        box =  aabb(vec3(x0,k-0.0001,z0), vec3(x1, k+0.0001, z1));
        return true; }
    material  *mp;
    float x0, x1, z0, z1, k;
};

class yz_rect: public hitable  {
public:
    yz_rect() {}
    yz_rect(float _y0, float _y1, float _z0, float _z1, float _k, material *mat) : y0(_y0), y1(_y1), z0(_z0), z1(_z1), k(_k), mp(mat) {};
    virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
    virtual bool bounding_box(float t0, float t1, aabb& box) const {
        box =  aabb(vec3(k-0.0001, y0, z0), vec3(k+0.0001, y1, z1));
        return true; }
    material  *mp;
    float y0, y1, z0, z1, k;
};

hit方法:


bool xz_rect::hit(const ray& r, float t0, float t1, hit_record& rec) const {
    float t = (k-r.origin().y()) / r.direction().y();
    if (t < t0 || t > t1)
        return false;
    float x = r.origin().x() + t*r.direction().x();
    float z = r.origin().z() + t*r.direction().z();
    if (x < x0 || x > x1 || z < z0 || z > z1)
        return false;
    rec.u = (x-x0)/(x1-x0);
    rec.v = (z-z0)/(z1-z0);
    rec.t = t;
    rec.mat_ptr = mp;
    rec.p = r.point_at_parameter(t);
    rec.normal = vec3(0, 1, 0);
    return true;
}

bool yz_rect::hit(const ray& r, float t0, float t1, hit_record& rec) const {
    float t = (k-r.origin().x()) / r.direction().x();
    if (t < t0 || t > t1)
        return false;
    float y = r.origin().y() + t*r.direction().y();
    float z = r.origin().z() + t*r.direction().z();
    if (y < y0 || y > y1 || z < z0 || z > z1)
        return false;
    rec.u = (y-y0)/(y1-y0);
    rec.v = (z-z0)/(z1-z0);
    rec.t = t;
    rec.mat_ptr = mp;
    rec.p = r.point_at_parameter(t);
    rec.normal = vec3(1, 0, 0);
    return true;
}

再在場景中放5個牆,一個燈,做個經典的cornell box。

// cornell_box經典場景
hitable *cornell_box() {
    hitable **list = new hitable*[8];
    int i = 0;
    material *red = new lambertian( new constant_texture(vec3(0.65, 0.05, 0.05)) );
    material *white = new lambertian( new constant_texture(vec3(0.73, 0.73, 0.73)) );
    material *green = new lambertian( new constant_texture(vec3(0.12, 0.45, 0.15)) );
    material *light = new diffuse_light( new constant_texture(vec3(15, 15, 15)) );
    list[i++] = new flip_normals(new yz_rect(0, 555, 0, 555, 555, green));
    list[i++] = new yz_rect(0, 555, 0, 555, 0, red);
    list[i++] = new xz_rect(213, 343, 227, 332, 554, light);
    list[i++] = new flip_normals(new xz_rect(0, 555, 0, 555, 555, white));
    list[i++] = new xz_rect(0, 555, 0, 555, 0, white);
    list[i++] = new flip_normals(new xy_rect(0, 555, 0, 555, 555, white));
    return new hitable_list(list,i);
}

camera的引數做一些調整

    vec3 lookfrom(278,278,-800);
    vec3 lookat(278, 278, 0);
    float dist_to_focus = 10.0;
    float aperture = 0.1;
    float vfov = 40.0;
    camera cam(lookfrom, lookat, vec3(0, 1, 0), vfov, float(nx) / float(ny), aperture, dist_to_focus, 0.0, 1.0);

會發現渲染出來的有幾個面是黑色的,是因為法向量的問題。

須要翻轉法向量

// 翻轉法向量
class flip_normals : public hitable {
public:
    flip_normals(hitable *p) : ptr(p) {}
    virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
        if (ptr->hit(r, t_min, t_max, rec)) {
            rec.normal = -rec.normal;
            return true;
        }
        else
            return false;
    }
    virtual bool bounding_box(float t0, float t1, aabb& box) const {
        return ptr->bounding_box(t0, t1, box);
    }
    hitable *ptr;
};

最後渲染出來的圖片長這樣:

注意

開始的時候渲染出來一片黑,排查了很久,是color裡面的tmin設定的問題,原來設定是0,原始碼中是0.001.

當tmin設0的時候會導致,遍歷hitlist時候,ray的t求解出來是0,hit的時候全走了else,導致遞迴到50層的時候,最後return的是0,* attenuation結果還是0。距離越遠,散射用到random_in_unit_sphere生成的ray誤差越大,就像上面的圖一樣。所以cornel 距離5,600的時候,場景中的lambert就全黑了。

    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const  {
        vec3 target = rec.p + rec.normal + random_in_unit_sphere();
        scattered = ray(rec.p, target-rec.p, r_in.time());
        attenuation = albedo->value(rec.u, rec.v, rec.p);
        return true;
    }

Chapter7:Instances

上一章渲染的結果是cornel box,但其實還不是完整的,完整版的在空間中還會有2個有輕微偏移的立方體。所以首先寫一個box的類,用之前的rect來實現一個立方體,box類繼承hitable,實現hit和bounding_box的虛方法。

class box: public hitable  {
    public:
        box() {}
        box(const vec3& p0, const vec3& p1, material *ptr);
        virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
        virtual bool bounding_box(float t0, float t1, aabb& box) const {
               box =  aabb(pmin, pmax);
               return true; }
        vec3 pmin, pmax;
        hitable *list_ptr;
};

box::box(const vec3& p0, const vec3& p1, material *ptr) {
    pmin = p0;
    pmax = p1;
    hitable **list = new hitable*[6];
    list[0] = new xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p1.z(), ptr);
    list[1] = new flip_normals(new xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p0.z(), ptr));
    list[2] = new xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p1.y(), ptr);
    list[3] = new flip_normals(new xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p0.y(), ptr));
    list[4] = new yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p1.x(), ptr);
    list[5] = new flip_normals(new yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p0.x(), ptr));
    list_ptr = new hitable_list(list,6);
}

bool box::hit(const ray& r, float t0, float t1, hit_record& rec) const {
    return list_ptr->hit(r, t0, t1, rec);
}

新建2個box的物件

    list[i++] = new box(vec3(130,0,65),vec3(295,165,230),white);
    list[i++] = new box(vec3(265,0,295),vec3(430,330,460),white);

渲染出來的影象如下

當前從側面看的話,這2個box的關係如下:

但其實目前的做法和真正cornelbox中是不一樣的,我們是在空間中擺放了2個不同位置的box,但其實第二個box是可以通過transform屬性,通過第一個box來表示出來的,一般相同形狀的模型都可以用instance的方法模擬出來。在hitable.h中實現,translate,繼承hitable,同樣實現hit和boundingbox的虛擬函式,這2個虛擬函式都用到了translate這個類中的一個成員變數vec3 的offset表示偏移量。


// 用於instance的移動
class translate : public hitable {
public:
    translate(hitable *p, const vec3& displacement) : ptr(p), offset(displacement) {}
    virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const;
    virtual bool bounding_box(float t0, float t1, aabb& box) const;
    hitable *ptr;
    vec3 offset;    // vec3的偏移
};

bool translate::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
    ray moved_r(r.origin() - offset, r.direction(), r.time());
    if (ptr->hit(moved_r, t_min, t_max, rec)) {
        rec.p += offset;
        return true;
    }
    else
        return false;
}

bool translate::bounding_box(float t0, float t1, aabb& box) const {
    if (ptr->bounding_box(t0, t1, box)) {
        box = aabb(box.min() + offset, box.max()+offset);
        return true;
    }
    else
        return false;
}

這樣就可以實現平移操作了,對於旋轉,座標點在三維空間中繞z軸旋轉的示意圖如下,z座標保持不變,x,y座標旋轉,旋轉角度為theta。

繞z軸旋轉時,xy的座標變化如下

x' = cos(theta) * x - sin(theta) * y
y' = sin(theta) * x + cos(theta) * y

同理繞y軸旋轉時,xz如下

x' = cos(theta) * x + sin(theta) * z
z' = -sin(theta) * x + cos(theta) * z

繞x軸旋轉,yz如下

y' = cos(theta) * y - sin(theta) * z
z' = sin(theta) * y + cos(theta) * z

相比平移,旋轉還需要考慮的一個問題就是,當面發生轉動的時候,面的法向量也是會發生轉動的,面法向量發生轉動,散射的出射光線也會發生改變。實現rotate_y繼承hitable,也是實現hit和bounding_box的虛擬函式,加入2個新的成員變數sin_theta和cos_theta用於角度計算。

class rotate_y : public hitable {
public:
    rotate_y(hitable *p, float angle);
    virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const;
    virtual bool bounding_box(float t0, float t1, aabb& box) const {
        box = bbox; return hasbox;}
    hitable *ptr;
    float sin_theta;
    float cos_theta;
    bool hasbox;
    aabb bbox;
};

實現hit和rotate_y的建構函式

rotate_y::rotate_y(hitable *p, float angle) : ptr(p) {
    float radians = (M_PI / 180.) * angle;
    sin_theta = sin(radians);
    cos_theta = cos(radians);
    hasbox = ptr->bounding_box(0, 1, bbox);
    vec3 min(FLT_MAX, FLT_MAX, FLT_MAX);
    vec3 max(-FLT_MAX, -FLT_MAX, -FLT_MAX);
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            for (int k = 0; k < 2; k++) {
                float x = i*bbox.max().x() + (1-i)*bbox.min().x();
                float y = j*bbox.max().y() + (1-j)*bbox.min().y();
                float z = k*bbox.max().z() + (1-k)*bbox.min().z();
                float newx = cos_theta*x + sin_theta*z;
                float newz = -sin_theta*x + cos_theta*z;
                vec3 tester(newx, y, newz);
                // 旋轉之後重新計算bounding box
                for ( int c = 0; c < 3; c++ )
                {
                    if ( tester[c] > max[c] )
                        max[c] = tester[c];
                    if ( tester[c] < min[c] )
                        min[c] = tester[c];
                }
            }
        }
    }
    bbox = aabb(min, max);
}

bool rotate_y::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
    vec3 origin = r.origin();
    vec3 direction = r.direction();
    origin[0] = cos_theta*r.origin()[0] - sin_theta*r.origin()[2];
    origin[2] =  sin_theta*r.origin()[0] + cos_theta*r.origin()[2];
    direction[0] = cos_theta*r.direction()[0] - sin_theta*r.direction()[2];
    direction[2] = sin_theta*r.direction()[0] + cos_theta*r.direction()[2];
    ray rotated_r(origin, direction, r.time());
    if (ptr->hit(rotated_r, t_min, t_max, rec)) {
        vec3 p = rec.p;
        vec3 normal = rec.normal;
        // normal 也做相應的旋轉,因為是繞y軸,所以改p[0]和p[2]
        p[0] = cos_theta*rec.p[0] + sin_theta*rec.p[2];
        p[2] = -sin_theta*rec.p[0] + cos_theta*rec.p[2];
        normal[0] = cos_theta*rec.normal[0] + sin_theta*rec.normal[2];
        normal[2] = -sin_theta*rec.normal[0] + cos_theta*rec.normal[2];
        rec.p = p;
        rec.normal = normal;
        return true;
    }
    else
        return false;
}

這時候就可以把main中的box引用上translate和roate了,記得先旋轉再平移

        list[i++] = new translate(new rotate_y(new box(vec3(0, 0, 0), vec3(165, 165, 165), white), -18), vec3(130,0,65));
        list[i++] = new translate(new rotate_y(new box(vec3(0, 0, 0), vec3(165, 330, 165), white),  15), vec3(265,0,295));

Chapter8:Volumes

第八章是比較激動人心的一張,Volumes,常見的體渲染包括 煙、霧等。volumes的另一個特性是,可以在內部發生散射,就像光線會在稠密的霧中發生散射一樣。體渲染通常的做法是,在體的內部,放很多隨機性的面,來實現散射的效果。比如一束煙可以表示為,在這束煙的內部任意point位置,都可以存在一個面,面的集合實現了煙的物理效果(感覺翻譯的不好,大概就是這個意思0.0)。

對於一個常量密度的volume,一條ray通過其中的時候,在volume中傳播的時候也會發生散射,光線在volume中能傳播多遠,也是由volume的密度決定的,密度越高,傳播的效率越低,光線傳播的距離也越短。

當光線穿過volume的時候,volume中的任意位置都可以發生散射,

probability = C * dL

這裡的C是一個volume的係數,表示這個volume視覺上的濃密程度,dL是任意可以發生散射的一小段距離。對於一個Constant volume,我們只需要調節 密度 C和Boundary包圍盒。這裡再寫一個constant_medium繼承hitable。
在材質裡面新增各向異性的材質

// 各向異性材質
class isotropic : public material {
public:
    isotropic(texture *a) : albedo(a) {}
    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const  {
        scattered = ray(rec.p, random_in_unit_sphere());
        attenuation = albedo->value(rec.u, rec.v, rec.p);
        return true;
    }
    texture *albedo;
};
// 體,恆量介質
class constant_medium : public hitable  {
public:
    constant_medium(hitable *b, float d, texture *a) : boundary(b), density(d) { phase_function = new isotropic(a); }
    virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const;
    virtual bool bounding_box(float t0, float t1, aabb& box) const {
        return boundary->bounding_box(t0, t1, box); }
    hitable *boundary;
    float density;
    // 材質為各項異性的材質
    material *phase_function;
};

bool constant_medium::hit(const ray &r, float t_min, float t_max, hit_record &rec) const {
    hit_record rec1, rec2;
    if (boundary->hit(r, -FLT_MAX, FLT_MAX, rec1)) {
        if (boundary->hit(r, rec1.t+0.0001, FLT_MAX, rec2)) {
            if (rec1.t < t_min)
                rec1.t = t_min;
            if (rec2.t > t_max)
                rec2.t = t_max;
            if (rec1.t >= rec2.t)
                return false;
            if (rec1.t < 0)
                rec1.t = 0;
            float distance_inside_boundary = (rec2.t - rec1.t)*r.direction().length();
            float hit_distance = -(1/density)*log(drand48());
            if (hit_distance < distance_inside_boundary) {
                rec.t = rec1.t + hit_distance / r.direction().length();
                rec.p = r.point_at_parameter(rec.t);
                rec.normal = vec3(1,0,0);  // arbitrary
                rec.mat_ptr = phase_function;
                return true;
            }
        }
    }
    return false;
}

使用main中的cornell_smoke 渲染出來的場景如下:

Chapter9:A Scene Testing All New Features

最後一張是渲染運用第二本書上的知識點,渲染出一張和封面一樣的圖片。

我渲染的場景和原書略有不同,引數如下:

解析度1000x1000 sample 100

hitable *final() {
    int nb = 10;
    hitable **list = new hitable*[3000];
    material *white = new lambertian( new constant_texture(vec3(0.73, 0.73, 0.73)) );
    material *ground = new lambertian( new constant_texture(vec3(0.48, 0.83, 0.53)) );
    int b = 0;
    int l = 0;
    for (int i = 0; i < nb; i++) {
        for (int j = 0; j < nb; j++) {
            float w = 100;
            float x0 = i*w;
            float z0 = j*w;
            float y0 = 0;
            float x1 = x0 + w;
            float y1 = 100*(drand48()+0.01);
            float z1 = z0 + w;
            cout << "("<<x0<<","<<y0<<","<<z0<<") ("<<x1<<","<<y1<<","<<z1<<")"<<endl;
            list[l++] = new box(vec3(x0, y0, z0), vec3(x1, y1, z1), ground);
        }
    }
    material *light = new diffuse_light( new constant_texture(vec3(7, 7, 7)) );
    list[l++] = new xz_rect(123, 423, 147, 412, 554, light);
    vec3 center(400, 400, 200);
    list[l++] = new moving_sphere(center, center+vec3(30, 0, 0), 0, 1, 50,
                                  new lambertian(new constant_texture(vec3(0.7, 0.3, 0.1))));
    list[l++] = new sphere(vec3(260, 150, 45), 50, new dielectric(1.5));
    list[l++] = new sphere(vec3(0, 150, 145), 50, new metal(vec3(0.8, 0.8, 0.9), 10.0));
    hitable *boundary = new sphere(vec3(360, 150, 145), 70, new dielectric(1.5));
    list[l++] = boundary;
    list[l++] = new constant_medium(boundary, 0.2, new constant_texture(vec3(0.2, 0.4, 0.9)));
    boundary = new sphere(vec3(0, 0, 0), 5000, new dielectric(1.5));
    list[l++] = new constant_medium(boundary, 0.0001, new constant_texture(vec3(1.0, 1.0, 1.0)));
    texture *pertext = new noise_texture(0.1);
    list[l++] =  new sphere(vec3(220,280, 300), 80, new lambertian( pertext ));
    int ns = 1000;
    for (int j = 0; j < ns; j++) {
        list[l++] = new sphere(vec3(165*drand48()-100, 165*drand48()+270, 165*drand48()+395),10 , white);
    }
    cout<< "len(l) = " << l << endl;
    return new hitable_list(list,l);
}