前言

  • 本系列总结自Ray Tracing: The Next Week,前篇为rayTracingInOneWeekend – 爱莉希雅 – 博客园 (cnblogs.com)。需要一点点的图形学基础,不懂的同学可以移至GAMES101-现代计算机图形学入门-闫令琪_哔哩哔哩_bilibili
  • 这个版本的光线追踪非常适合初学者,难度既不大但也感受到光追的魅力

总结

动态模糊

  • 动态模糊:在生活中,摄像机快开打开的时间间隔中,摄像机和物体都在移动,拍出来的图像肯定是糊的。大多数特效都可以通过对每个像素进行多重采样实现,动态模糊也属于这个范畴。

  • 方法:我们可以随机的方法在不同时间发射多条射线来模拟快门的打开。也就是说,在快门打开时,随时间变化随机生成光线,并同时射出与模型相交。一般而言,我们让摄像机和物体同时运动,并让每条射线都拥有自己的时间点

    只要物体在那个时间处于它自己正确的位置,我们便能得出光线在那个时间的精确平均值

  • 实现。让漫反射材质的球动起来:

    1. 让每条射线拥有自己的时间点

      class ray{    public:    ray() {}    ray(const vec3& origin, const vec3& direction, double time = 0.0)        : orig(origin), dir(direction), tm( time )        {}    vec3 get_origin() const { return orig; }    vec3 get_direction() const { return dir; }    vec3 at(double t) const    {        return orig + dir * t;    }    private:    vec3 orig;    vec3 dir;    double tm;};
    2. 让camera在[time1, time2]中随机生成光线

      class camera{    public:        camera(vec3 lookfrom, vec3 lookat, vec3 up_vector, double boc, double aspect, double aperture, double focus_dist, double t0 = 0.0, double t1 = 0.0)    {        origin = lookfrom;        lens_radius = aperture / 2;        time0 = t0;        time1 = t1;        auto theta = degrees_to_radians(boc);        auto half_height = tan(theta / 2);        auto half_width = aspect * half_height;        w = unit_vector(lookfrom - lookat);        u = unit_vector(cross(up_vector, w));        v = cross(w, u);        lower_left_corner = origin - u * half_width * focus_dist - v * half_width * focus_dist - w * focus_dist;        horizontal = u * 2 * half_width * focus_dist;        vertical = v * 2 * half_height * focus_dist;    }    ray get_ray(double s, double t)    {        vec3 rd = random_in_unit_disk() * lens_radius;        vec3 offset = u * rd.x() + v * rd.y();        return ray(origin + offset, lower_left_corner + horizontal * s + vertical * t - origin - offset, random_double(time0, time1));    }    private:    vec3 lower_left_corner;    vec3 horizontal;    vec3 vertical;    vec3 origin;    vec3 u, v, w;    double lens_radius;    double time0, time1;};
    3. 让sphere class的球心在[time0, time1]内从center0运动至center1。超出这个时间范围,球心依然在动,也就是说,线性插值time既可小于0.0,亦可大于1.0。这两个时间变量和快门的开关时刻无需一一对应

      class moving_sphere : public hittable{    public:    moving_sphere() = default;    moving_sphere(vec3 cen0, vec3 cen1, double t0, double t1, double r, shared_ptr m_ptr)        : center0(cen0), center1(cen1), time0(t0), time1(t1), radius(r), mat_ptr(m_ptr)        {}    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const;    vec3 center(double time) const;    private:    vec3 center0, center1;    double time0, time1;    double radius;    shared_ptrmat_ptr;};inline vec3 moving_sphere::center(double time) const{    //线性插值    return center0 + ( (center1 - center0) * (time - time0) / (time1 - time0) );}//与sphere的变化只在centerinline bool moving_sphere::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const{    vec3 oc = r.get_origin() - center(r.get_time());    auto a = r.get_direction().length_squared();    auto half_b = dot(r.get_direction(), oc);    auto c = oc.length_squared() - radius * radius;    auto delta = half_b * half_b - a * c;    if (delta > 0)    {        auto root = sqrt(delta);        auto temp = (-half_b - root) / a;        if (temp > t_min && temp  t_min && temp < t_max)        {            rec.t = temp;            rec.p = r.at(rec.t);            vec3 outward_normal = (rec.p - center(r.get_time())) / radius;            rec.set_face_normal(r, outward_normal);            rec.mat_ptr = mat_ptr;            return true;        }    }    return false;}
    4. 确保moving_sphere的材质在运算散射时,散射光线与入射光线存在的时间点相同

      //漫反射材质class lambertian : public material{    public:    lambertian( const vec3& a) : albedo(a) {}    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const    {        vec3 scatter_direction = rec.normal + random_unit_vector();        scattered = ray(rec.p, scatter_direction, r_in.get_time());        attenuation = albedo;        return true;    }    private:    vec3 albedo;};
    5. 快门在time0打开,在time1关闭,每个球心从位置L运动到L + (0,r/2, 0),r为[0,1]的随机数。跳动肯定时变换y轴

      hittable_list random_scene(){    hittable_list world;    world.add(make_shared(vec3(0, -1000, 0), 1000, make_shared(vec3(0.5, 0.5, 0.5))));    int i = 1;    for (int a = -10; a < 10; ++a)    {        for (int b = -10; b  0.9)            {                if (choose_mat < 0.8)                {                    //diffuse                    auto albedo = vec3::random() * vec3::random();                    world.add(make_shared(center, center + vec3(0, random_double(0, 0.5), 0), 0.0, 1.0, 0.2, make_shared(albedo)));                }                else if (choose_mat < 0.95)                {                    //metal                    auto albedo = vec3::random(0.5, 1);                    auto fuzz = random_double(0, 0.5);                    world.add(make_shared(center, 0.2, make_shared(albedo, fuzz)));                }                else                {                    //glass                    world.add(make_shared(center, 0.2, make_shared(1.5)));                }            }        }    }    world.add(make_shared(vec3(0, 1, 0), 1.0, make_shared(1.5)));    world.add(        make_shared(vec3(-4, 1, 0), 1.0, make_shared(vec3(0.4, 0.2, 0.1))));    world.add(        make_shared(vec3(4, 1, 0), 1.0, make_shared(vec3(0.7, 0.6, 0.5), 0.0)));    return world;}void output_image(){    //...    auto world = random_scene();    vec3 lookfrom(13, 2, 3);    vec3 lookat(0, 0, 0);    vec3 up_vector(0, 1, 0);    auto dist_to_focus = 10.0;    auto aperture = 0.0;    const auto aspect_ratio = double(image_width) / image_height;    //...}
  • 输出:

  • 总结:用随机的方法在不同时间发射多条射线来模拟快门的打开。让摄像机和物体同时运动,并让每条射线都拥有自己的时间点

层次包围盒

  • 在这之前,光线的求交运算会耗费大量时间,且计算所消耗的时间与场景中的物体数量线性相关。使用遍历反复查找同一个模型会有许多多余的计算,因此本节将使用一种树型结构来进行加速,即层次包围盒(Bounding Volume Hierarchy)

  • 层次包围盒又属于空间数据结构,空间数据结构有以下优点

    1. 易于将二维和三维的概念扩展到高维
    2. 可用于加速查询。如裁减算法、相交测试、光线追踪、碰撞检测等。访问速度很快。时间复杂度可以从 O(n)提升到 O(log n)
  • 层次包围盒的思想:用体积略大而几何特征简单的包围盒来近似地描述复杂的几何对象,从而只需对包围盒重叠的对象进行进一步的相交测试。此外,通过构造树状层次结构,可以越来越逼近对象的几何模型,直到几乎完全获得对象的几何特征。不再以空间作为划分依据,而是从对象的角度考虑

    层次包围盒的代码大致如下:

    if (ray hits bounding object)    return whether ray hits bounded objectselse    return false
  • 层次包围盒的缺点:对于复杂的模型,效果并不好

  • 举个例子,以下场景以层次树结构进行组织,包含一个根节点(root)、一些内部节点(internal nodes),以及一些叶子节点 (leaves)。顶部的节点是根,其无父节点。叶子节点(leaf node)包含需渲染的实际几何体,且其没有子节点

  • AABB包围盒求交

    • 运用层次包围盒的前提是 如何判断光线是否与物体相交,这里就用到了AABB包围盒:由三对平面的交集构成,且任意一对平面都与x轴,y轴、z轴垂直

      如图:

    • 思想: 光线若不会与一个物体相交,那就没有必要去遍历该物体的所有三角形面。因此利用一个包围盒包住该物体,在与该物体计算求交之前先判断光线是否与包围盒相交,若包围盒与光线没有交点,那么显然不会与物体有交点

    • 举个例子:我们以2维的AABB包围盒为例

      如下图,左边的图计算x轴上 最晚进来的交点t\(min\) 和 最先出去的t\(max\);中间的图计算y轴上 最先进来的交点t\(min\)和最先出去的t\(max\);右边的图则综合前两张图的计算结果,计算出最终的结果。最终结果是如何计算的呢?主要是以下几点:

      1. 只有当光线进入了所有平面,才是真正地进入了包围盒。\(\large t_{enter} = max(tmin)\)
      2. 只有当光线离开了任一平面,就算是真正地离开了包围盒。\(\large t_{exit} = min(t_{max})\)
      3. 当且仅当,\(\large t_{enter} 0\)成立时,才有交点。
      4. 若t小于0,说明盒子在光线背后;若\(\large t_{exit} >= 0, t_{enter} < 0\),说明光线起点在盒子内部
    • 为什么使用AABB?因为求t时,只需使用某一轴的信息,不必使用整个坐标,与点乘相比更容易

    • 伪代码实现:

      compute(tx0, tx1)compute(ty0, ty1)compute(tz0, tz1)return overlap?( (tx0, tx1), (ty0, ty1), (tz0, tz1))
  • 实现:

    1. 若光线向x轴负方向射出并与包围盒相交,很显然此时\(t_{enter} > t_{exit}\),因此我们需要交换一下

      而对于Bx = 0来说,此时答案是无意义的,为NaN,因此对于任何NaN结果来说,我们都返回false;对于(x0 – Ax)和(x1 – Ax)两分子为0来说,此时只有一个解,相当于是光线擦边,认定光线射中或没射中都可以,我们认定它没有射中,即return false

      class aabb{    private:    vec3 m_min;    vec3 m_max;    public:    aabb() = default;    aabb(const vec3& a, const vec3& b) { m_min = a; m_max = b; }    vec3 get_min() const { return m_min; }    vec3 get_max() const { return m_max; }    bool hit(const ray& r, double tmin, double tmax) const    {        for (int i = 0; i < 3; ++i)        {            //若光线向x轴负方向射出            //库函数考虑了异常、NaN,速度要慢些            //以下考虑了分母为0的情况            auto t0 = ffmin((m_min[i] - r.get_origin()[i]) / r.get_direction()[i],                            (m_max[i] - r.get_origin()[i]) / r.get_direction()[i]);            auto t1 = ffmax((m_min[i] - r.get_origin()[i]) / r.get_direction()[i],                            (m_max[i] - r.get_origin()[i]) / r.get_direction()[i]);            //只有当光线进入了所有平面,才是真正地进入了包围盒            tmin = ffmax(t0, tmin);            //只有当光线离开了任一平面,就算是真正地离开了包围盒            tmax = ffmin(t1, tmax);            if (tmax <= tmin)            {                return false;            }        }        return true;    }};
    2. 每个图元构建aabb的方式。在hittable class中加入bounding_box()来计算包裹hittable class的包围盒,再构建一个层次包围盒的树形结构。在层次树中,所有图元都处于叶子节点,又因为物体会动,所以还需接受time1和time2

      class hittable{    public:    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const = 0;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const = 0;};//计算两个包围盒的包围盒aabb surrounding_box(aabb box0, aabb box1){    vec3 small(ffmin(box0.get_min().x(), box1.get_min().x()),               ffmin(box0.get_min().y(), box1.get_min().y()),               ffmin(box0.get_min().z(), box1.get_min().z()));    vec3 big(ffmax(box0.get_max().x(), box1.get_max().x()),             ffmax(box0.get_max().y(), box1.get_max().y()),             ffmax(box0.get_max().z(), box1.get_max().z()));    return aabb(small, big);}//sphere构建包围盒很简单bool sphere::bounding_box(double t0, double t1, aabb& output_box) const{    output_box = aabb(center - vec3(radius, radius, radius), center + vec3(radius, radius, radius));    return true;}//moving_sphere构建包围盒,先求球体在t0时的包围盒,再求t1时的包围盒,再计算这两个包围盒的包围盒bool moving_sphere::bounding_box(double t0, double t1, aabb& output_box) const{    aabb box0(center(t0) - vec3(radius, radius, radius), center(t0) + vec3(radius, radius, radius));    aabb box1(center(t1) - vec3(radius, radius, radius), center(t1) + vec3(radius, radius, radius));    output_box = surrounding_box(box0, box1);    return true;}//hittable_list构建包围盒,在运行时计算(构造函数中进行运算亦可),因为包围盒的计算一般只有在BVH构造时才调用inline bool hittable_list::bounding_box(double t0, double t1, aabb& output_box) const{    if (objects.empty() == true)    {        return false;    }    aabb temp_box;    bool first_box = true;    for (const auto& object : objects)    {        //图元无法构建包围盒        if (object->bounding_box(t0, t1, temp_box) == false)        {            return false;        }        //是否是第一个包围盒,不是则进行构造两个包围盒的包围盒        output_box = first_box ? temp_box : surrounding_box(output_box, temp_box);        first_box = false;    }    return true;}
    3. 构建层次包围盒

      1. 构建层次包围盒的分割原则。每次分割时沿轴把物体列表分成两半:

        1. 随机选取任一一轴进行分割
        2. 使用sort()对图元进行排序,也就是对图元任一轴的坐标大小进行排序
        3. 根据排序对半分,每个子树分一半的图元
      2. BVH也应为hittable一员。两种设计方案,我选择前一种:

        1. 为树和树的节点设计两个不同的class

        2. 一个class加上指针

          //constantAndTool.hinline int random_int(int min, int max){    return static_cast(random_double(min, max + 1));}//hittable.h//设计bvh树节点class bvh_node : public hittable{    private:    shared_ptrleft;    shared_ptrright;    aabb box;    public:    bvh_node() = default;    bvh_node(hittable_list& list, double time0, double time1)        : bvh_node(list.get_objects(), 0, list.get_objects().size(), time0, time1)        {}    bvh_node(vector<shared_ptr>& objects, size_t start, size_t end, double time0, double time1);    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const;};//若光线打中包围盒,则继续判断的子节点inline bool bvh_node::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const{    //终止条件    if (box.hit(r, t_min, t_max) == false)    {        return false;    }    //左右中    bool hit_left = left->hit(r, t_min, t_max, rec);    bool hit_right = right->hit(r, t_min, hit_left == true ? rec.t : t_max, rec);    return hit_left || hit_right;}//std::sort的compare函数//先判断是哪个轴,再为compare函数赋值inline bool box_compare(const shared_ptr a, const shared_ptrb, int axis){    aabb box_a;    aabb box_b;    //两图元中任一图元无法构建包围盒    if (a->bounding_box(0, 0, box_a) == false || b->bounding_box(0, 0, box_b) == false)    {        cerr << "No bounding box in bvh_node constructor." << endl;    }    return box_a.get_min().e[axis] < box_b.get_min().e[axis];}inline bool box_x_compare(const shared_ptra, const shared_ptr b){    return box_compare(a, b, 0);}inline bool box_y_compare(const shared_ptra, const shared_ptr b){    return box_compare(a, b, 1);}inline bool box_z_compare(const shared_ptra, const shared_ptr b){    return box_compare(a, b, 2);}//构建bvh//当数组传入的元素个数只剩下两个时,在左右子树各放一个;当只有一个元素时,为了省去检查空指针这步,在左右子树都放一个inline bvh_node::bvh_node(const vector<shared_ptr>& objects, size_t start, size_t end, double time0, double time1){    auto temp_objects = objects;    int axis = random_int(0, 2);    auto comparator = (axis == 0) ? box_x_compare : (axis == 1) ? box_y_compare : box_z_compare;    //当前还需分割多少图元    size_t object_span = end - start;    if (object_span == 1)    {        left = right = temp_objects[start];    }    else if (object_span == 2)    {        if (comparator(temp_objects[start], temp_objects[start + 1]) == true)        {            left = temp_objects[start];            right = temp_objects[start + 1];        }        else        {            left = temp_objects[start + 1];            right = temp_objects[start];        }    }    else    {        std::sort(temp_objects.begin() + start, temp_objects.begin() + end, comparator);        auto mid = start + object_span / 2;        left = make_shared(temp_objects, start, mid, time0, time1);        right = make_shared(temp_objects, mid, end, time0, time1);    }    aabb box_left, box_right;    //检验图元是否有aabb    if (left->bounding_box(time0, time1, box_left) == false || right->bounding_box(time0, time1, box_right) == false)    {        cerr << "No bounding box in bvh_node constructor." << endl;    }    box = surrounding_box(box_left, box_right);}
    4. 在random_scene()中进行运用

      hittable_list random_scene(){    //...    return static_cast(make_shared(world, 0.0, 1.0));}
  • 输出:结果图和上一小节一样但是速度更快

  • 总结:

纹理

  • 纹理贴图:通过将投影方程(projector function)运用于空间中的点 ,从而得到一组称为参数空间值(parameterspacevalues)的关于纹理的数值

    纹理贴图常常意味着一个将颜色赋予图元表面的过程,这个过程可以是纹理生成代码,亦可以是一张图片,亦可以是前两者的结合

  • 棋盘格上的球体

    1. texture class。两种设计思路:

      1. 把静态rgb颜色和贴图写成两class,以此来区分
      2. 把静态rgb颜色和贴图写成同一类别class,以此将任何颜色弄成纹理
      class texture{    public:    virtual vec3 value(double u, double v, const vec3& p) const = 0;};class constant_texture : public texture{    public:    constant_texture() = default;    constant_texture( vec3 c ) : color(c) {}    virtual vec3 value(double u, double v, const vec3& p) const override    {        return color;    }    private:    vec3 color;};
    2. 更新hit_record,让其存储相交点的uv信息

      struct hit_record{    //...    double t;    double u;    double v;    bool front_face;//法线朝向//...};
    3. 更新漫反射材质,对其进行纹理贴图

      将vec3改为纹理指针

      class lambertian : public material{    public:    lambertian( shared_ptr a ) : albedo(a) {}    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const    {        vec3 scatter_direction = rec.normal + random_unit_vector();        scattered = ray(rec.p, scatter_direction, r_in.get_time());        attenuation = albedo->value( rec.u, rec.v, rec.p);        return true;    }    private:    shared_ptr albedo;};
    4. 漫反射材质呈现棋盘格纹理

      运用sine和cosine函数周期性的变化,在三个维度都乘上周期函数

      class checker_texture : public texture{    public:    checker_texture() = default;    checker_texture( shared_ptr t0, shared_ptr t1 ) : even(t0), odd(t1){}    //use the sine and cosine functions to change periodically to create a checkerboard texture    virtual vec3 value(double u, double v, const vec3& p) const override    {        auto sine = sin(10 * p.x()) * sin(10 * p.y()) * sin(10 * p.z());        if (sine value(u, v, p);        }        else        {            return even->value(u, v, p);        }    }    private:    //The odd and even pointer points to a static texture    shared_ptr odd;    shared_ptr even;};
    5. 更新random_scene()

      hittable_list random_scene(){    //...    auto checkerboard = make_shared(make_shared(vec3(0.2, 0.3, 0.1)),make_shared(vec3(0.9, 0.9, 0.9)));    world.add(make_shared(vec3(0, -1000, 0), 1000, make_shared(checkerboard)));    //...    if (choose_mat < 0.8)    {        //diffuse        auto albedo = vec3::random() * vec3::random();        world.add(make_shared(center, center + vec3(0, random_double(0, 0.5), 0), 0.0, 1.0, 0.2, make_shared(make_shared(albedo))));    }    //...    world.add(make_shared(vec3(-4, 1, 0), 1.0, make_shared(make_shared(vec3(0.4, 0.2, 0.1)))));}
  • 输出:

  • 棋盘格球体

    1. 更换场景

      hittable_list two_spheres(){    hittable_list objects;    auto checkerboard = make_shared        (        make_shared(vec3(0.2, 0.3, 0.1)),        make_shared(vec3(0.9, 0.9, 0.9))    );    objects.add(make_shared(vec3(0, -10, 0), 10, make_shared(checkerboard)));    objects.add(make_shared(vec3(0, 10, 0), 10, make_shared(checkerboard)));    return static_cast(make_shared(objects, 0.0, 1.0));}void output_image(){    //...    auto world = two_spheres();    vec3 lookfrom(13, 2, 3);    vec3 lookat(0, 0, 0);    vec3 up_vector(0, 1, 0);    auto dist_to_focus = 10.0;    auto aperture = 0.0;    const auto aspect_ratio = double(image_width) / image_height;    camera cam(lookfrom, lookat, up_vector, 20, aspect_ratio, aperture, dist_to_focus, 0.0, 1.0);}
  • 输出:

  • 总结:让材质接受纹理(vec3),根据记录的纹理参数(hit_record)进行散射运算

柏林噪声概念

  • 什么是柏林噪声?

    大多数程序贴图纹理(计算机算法生成的)都是基于噪声函数的,而柏林噪声就属于噪声函数,这些程序贴图纹理旨在创建用于纹理映射的自然元素(大自然中的一些纹理)的真实表面或三维物体创建的纹理图像

  • 为什么使用柏林噪声?

    • 普通噪声毫无规律可言,而想要模拟一些逼真的随机效果需要在随机中又有规律

      以下是普通噪声:

    • 柏林噪声可以生成伪随机数,但其并不是完全随机,对其加以一些物体的规律的限制,可以用于模拟大理石的纹理、云、火焰等(这些自然现象都是随机的,但随机中又有规律)

      噪声函数的目的:在三维空间中,提供一种高效地实现、可重复、伪随机的信号,此信号大部分能量集中在一个空间频率附近,而视觉上各向同性的(统计上不随旋转变化)

      以下是柏林噪声:创建某种类似于完全的随机信号,通过低通滤波后,滤除所有的空间高频率而变得模糊

  • 柏林噪声的特点

    1. 所有视角细节都是相同的大小
    2. 可复现性。对于每个输入位置,提供一个可重复的伪随机值
    3. 实现简单快捷
    4. 有一个已知的范围(通常是[-1,1])
    5. 空间频率带宽有限
    6. 对于重复表现的不明显
    7. 空间频率在平移下是不变的
  • 如何实现柏林噪声?

    三个步骤:

    1. 初始化

      在3维空间的每个整数(x,y,z)位置产生一个可重复的伪随机值,这运用哈希函数实现,且哈希函数基于一个包含随机数在[0,255]的整数排列表

      排列表(Permutation Table):乱序存放一系列 索引

      梯度表(Gradient Table):存放一系列随机梯度值

    2. 建立采样空间,对于噪声图上的像素,找到它在单位矩形/单位立方体上对应的一点,并求出它的参考点的坐标

      1. 对于一维柏林噪声,采样空间为一个一维的坐标轴,轴上整数坐标位置均有一个点,参考点为两侧最近的整数点
      2. 对于二维柏林噪声,采样空间为一个二维坐标系,坐标系中横纵坐标为整数的地方均有一点,参考点为组成包围该点的单位正方体的四个点
      3. 三维同理,参考点为组成包围该点的单位立方体的八个点
    3. 对于不同类型的噪声,采样点在不同空间中, 根据伪随机整数的梯度来计算 采样空间中小数位置到参考点的距离向量和梯度向量的点积,最后进行插值

      一维噪声需一次插值,二维需三次插值,三维需七次插值,呈指数增长,时间复杂度是\(O(2^N)\)

  • 一维柏林噪声

    1. 初始化

      1. 决定一个点的梯度,我们需要结合哈希函数,以点坐标作为参数,以所得值作为索引,在排列表中取得对应值。也就是说,将点坐标转换为排列表的索引,最后取得索引对应的梯度表的值

        //以下是二维的哈希函数int perm = {...};int grad = {...};void hash( int[] gradient, int x, int y){    //找到对应值在排列表中的索引    int permIdx[] = new int[2];    permIdx[0] = FindIndex(x);    permIdx[1] = FindIndex(y);        //通过排列表的索引值查找梯度表的索引    int gradIdx[] = new int[2];    gradIdx[0] = perm[permIdx[0]];    gradIdx[0] = perm[permIdx[1]];        //通过梯度表的索引找到对应的梯度值    gradient[0] = grad[gradIdx[0]];    gradient[1] = grad[gradIdx[1]];}
      2. 在一维中,各个整点的梯度可以用斜率来表现,因此我们可以不使用梯度表,直接将排列表对应索引取得的值作为梯度

      3. 在经典噪声算法中,使用的排列表是由数字0到255散列组成的

        int permutation[] = { 151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43,172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181, 199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, 205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180 };
    2. 采样、插值、计算梯度向量和其偏移量之间的点积

      采样:一维柏林噪声函数的参数 输入点 是坐标轴上的一点,其参考点就是它左右两侧最近的整数点。因此,我们需要找到这两个点,并求出这两个点对应的梯度及相对输入点的方向向量

      插值:对梯度值进行插值,公式为\(\large x = x1 + t(x2 – x1)\),不如线性插值效果不太好,产生的噪声曲线很尖锐。我们使用\(\large t = 3t^2 – 2t^3\)来替换线性插值公式里的\(t\)

      float perlin( float x ){    //假设x1为x左边的索引值且x - x1的方向向量为正,x2为右边的索引值    int x1 = floor(x);    int x2 = x1 + 1;    //方向向量    float vec1 = x - x1;    float vec2 = x - x2;        //求梯度值    float grad1 = perm[x1 % 255] * 2.0 - 255.0;    float grad2 = perm[x2 % 255] * 2.0 - 255.0;        //插值    float t = 3 * pow( vec1, 2 ) - 2 * pow(vec1, 3);    float product1 = grad1 * vec1;    float product2 = grad2 * vec2;        return product1 + t * (product2 - product1);}
  • 三维柏林噪声

    1. 将空间(排列表)中的每一点(i,j,k)(坐标为整数)赋值为0,并赋予可以通过散列得到的伪随机斜率
    2. 定义坐标(x,y,z)为一个整数加上小数部分,\((x,y,z) = (i + u, j + v, k + w)\),找到围绕这个点的单位立方体的8个参考点,\((i,j,k),(i+1,j,k),·····(i+1, j+1, k+1)\)
    3. 采样、七次线性插值

    如下图所示:

  • 柏林噪声的缺点:

    1. 插值函数每个维度都是\(\large t = 3t^2 – 2t^3\),其二次导数含有非0值。如此在使用噪声函数的导数会导致失真的视觉效果
    2. 从(i,j,k)得到斜率,需要在一组256个伪随机斜率矢量中查找,这组矢量散布在3D立体球面上。这种分布上的不规则,在噪声函数的结果中产生了不希望的高频

实现大理石材质

将柏林噪声运用在光线追踪中

  • 用一个随机生成的三维数组铺满整个空间,会得到重复的区块

  • 实现:

    1. 用哈希表完成贴瓷砖

      class noise_texture : public texture{    private:    perlin noise;    public:    noise_texture() = default;    virtual vec3 value(double u, double v, const vec3& p) const override    {        return vec3(1, 1, 1) * noise.noise(p);    }};//perlin noiseclass perlin{    private:    static const int point_count = 256;    double* ranfloat;//Gradients table    //Permutation Table    int* perm_x;    int* perm_y;    int* perm_z;    //make the value of permutation table out of order    static void permute(int* p, int n)    {        for (int i = n - 1; i > 0; --i)        {            int target = random_int(0, i);            int tmp = p[i];            p[i] = p[target];//得到随机值            p[target] = tmp;//防止排列表重复        }    }    //生成排列表对应的伪随机斜率    static int* perlin_generate_perm()    {        auto p = new int[point_count];        for (int i = 0; i < perlin::point_count; ++i)        {            p[i] = i;        }        permute(p, point_count);        return p;    }        public:    perlin()    {        ranfloat = new double[point_count];        for (int i = 0; i < point_count; ++i)        {            ranfloat[i] = random_double();        }        perm_x = perlin_generate_perm();        perm_y = perlin_generate_perm();        perm_z = perlin_generate_perm();    }    ~perlin()    {        delete[] ranfloat;        delete[] perm_x;        delete[] perm_y;        delete[] perm_z;    }    //算得位置p的伪随机值    double noise(const vec3& p) const    {        //小数部分        auto u = p.x() - floor(p.x());        auto v = p.y() - floor(p.y());        auto w = p.z() - floor(p.z());        //整数部分        //限制索引值在[0,255]        auto i = static_cast(4 * p.x()) & 255;        auto j = static_cast(4 * p.y()) & 255;        auto k = static_cast(4 * p.z()) & 255;        return ranfloat[perm_x[i] ^ perm_y[j] ^ perm_z[k]];    }};hittable_list two_perlin_spheres(){    hittable_list objects;    auto pertext = make_shared();    objects.add(make_shared(vec3(0, -1000, 0), 1000, make_shared(pertext)));    objects.add(make_shared(vec3(0, 2, 0), 2, make_shared(pertext)));    return static_cast(make_shared(objects, 0.0, 1.0));}
      1. 输出:
    2. 我们可以让其更加平滑去掉格子,使用线性插值

      三线性插值:

      //三线性插值inline double trilinear_interp(double c[2][2][2], double u, double v, double w){    auto lerp_x_y0_z0 = c[0][0][0] + (u - 0) / 1 * (c[1][0][0] - c[0][0][0]);    auto lerp_x_y1_z0 = c[0][1][0] + u * (c[1][1][0] - c[0][1][0]);    auto lerp_x_y0_z1 = c[0][0][1] + u * (c[1][0][1] - c[0][0][1]);    auto lerp_x_y1_z1 = c[0][1][1] + u * (c[1][1][1] - c[0][1][1]);    auto lerp_y_z0 = lerp_x_y0_z0 + v * (lerp_x_y1_z0 - lerp_x_y0_z0);    auto lerp_y_z1 = lerp_x_y0_z1 + v * (lerp_x_y1_z1 - lerp_x_y0_z1);    auto lerp_z = lerp_y_z0 + w * (lerp_y_z1 - lerp_y_z0);    return lerp_z;}class perlin{    //...    //Calculate the random value of position P    double noise(const vec3& p) const    {        //decimal part        auto u = p.x() - floor(p.x());        auto v = p.y() - floor(p.y());        auto w = p.z() - floor(p.z());        //integral part        int i = floor(p.x());        int j = floor(p.y());        int k = floor(p.z());        double c[2][2][2];        //xor        //Record the Hasche transform to get the gradient of the reference point        //The limiting gradient value is [0,255]        for (int di = 0; di < 2; ++di)        {            for (int dj = 0; dj < 2; ++dj)            {                for (int dk = 0; dk < 2; ++dk)                {                    //hash transform                    //Finds the index of the gradient table bythe index values of the permutation table                    //The gradient value is found by the index of the gradient table                    //Here, the operands between Perm have no effect, just the color of the result                    c[di][dj][dk] = ranfloat[perm_x[(i + di) & 255] & perm_y[(j + dj) & 255] & perm_z[(k + dk) & 255]];                }            }        }        return trilinear_interp(c, u, v, w);    }};
      1. 输出:
    3. 从上图可以看出还是有明显的格子,其中有一部分是mach band效应,这是由线性变化的颜色构成的视觉感知效果。解决这一问题我们用:hermite cube来平滑差值

      1. 什么是mach bands效应?

        简单来说,mach band效应是一种边缘对比效应。当观察两块亮度不同的区域时,边界亮度对比加强,使得轮廓更为明显

      2. mach bands效应的原因

        人类的视觉系统可以增强边缘对比度

      3. 什么是hermite cube?

        就是我们在前面提到柏林噪声时说的替换t,\(x = x1 + t(x2 – x1),t = 3t^2 – 2t^3\)

        double noise(const vec3& p) const{    auto u = p.x() - floor(p.x());    auto v = p.y() - floor(p.y());    auto w = p.z() - floor(p.z());    u = 3 * std::pow(u, 2) - 2 * std::pow(u, 3);    v = 3 * std::pow(v, 2) - 2 * std::pow(v, 3);    w = 3 * std::pow(w, 2) - 2 * std::pow(w, 3);    //...}
      4. 输出:

    4. 加入花纹

      加入一个变量使得点的采样频率更高,实际上是索引值&255也有周期

      class noise_texture : public texture{    private:    perlin noise;    double scale;    public:    noise_texture() = default;    noise_texture( double sc ) : scale(sc) {}    virtual vec3 value(double u, double v, const vec3& p) const override    {        return vec3(1, 1, 1) * noise.noise(p * scale);    }};hittable_list two_perlin_spheres(){    auto pertext = make_shared(10.0);}
      1. 输出:
    5. 上图现在看起来还是有格子,这是因为实现的max和min总是精确地到整数。因此我们将使用梯度向量,并以点乘的方法将min和max的值退离网格点(也就是变为小数)

      将参考点到小数部分(u,v,w)的向量和伪随机单位向量进行点乘,再进行线性插值

      //vec3的三线性插值inline double vec3_trilinear_interp(vec3 c[2][2][2], double u, double v, double w){    auto uu = 3 * std::pow(u, 2) - 2 * std::pow(u, 3);    auto vv = 3 * std::pow(v, 2) - 2 * std::pow(v, 3);    auto ww = 3 * std::pow(w, 2) - 2 * std::pow(w, 3);    auto lerp_x_y0_z0 = dot(vec3(u, v, w), c[0][0][0]) + (dot(vec3(u - 1, v, w), c[1][0][0]) - dot(vec3(u, v, w), c[0][0][0])) * uu;    auto lerp_x_y1_z0 = dot(vec3(u, v - 1, w), c[0][1][0]) + (dot(vec3(u - 1, v - 1, w), c[1][1][0]) - dot(vec3(u, v - 1, w), c[0][1][0])) * uu;    auto lerp_x_y0_z1 = dot(vec3(u, v, w - 1), c[0][0][1]) + (dot(vec3(u - 1, v, w - 1), c[1][0][1]) - dot(vec3(u, v, w - 1), c[0][0][1])) * uu;    auto lerp_x_y1_z1 = dot(vec3(u, v - 1, w - 1), c[0][1][1]) + (dot(vec3(u - 1, v - 1, w - 1), c[1][1][1]) - dot(vec3(u, v - 1, w - 1), c[0][1][1])) * uu;    auto lerp_y_z0 = lerp_x_y0_z0 + (lerp_x_y1_z0 - lerp_x_y0_z0) * vv;    auto lerp_y_z1 = lerp_x_y0_z1 + (lerp_x_y1_z1 - lerp_x_y0_z1) * vv;    auto lerp_z = lerp_y_z0 + (lerp_y_z1 - lerp_y_z0) * ww;    return lerp_z;}class perlin{    private:    vec3* ran_vec3;//Get a random vec3    //...    public:    perlin()    {        ran_vec3 = new vec3[point_count];        for (int i = 0; i < point_count; ++i)        {            ran_vec3[i] = unit_vector(vec3::random(-1,1));        }        //...    }    ~perlin()    {        delete[] ran_vec3;        //...    }    //Calculate the random value of position P    double noise(const vec3& p) const    {        //decimal part        auto u = p.x() - floor(p.x());        auto v = p.y() - floor(p.y());        auto w = p.z() - floor(p.z());        //integral part        int i = floor(p.x());        int j = floor(p.y());        int k = floor(p.z());        vec3 c[2][2][2];        //xor        //Record the Hasche transform to get the gradient of the reference point        //The limiting gradient value is [0,255]        for (int di = 0; di < 2; ++di)        {            for (int dj = 0; dj < 2; ++dj)            {                for (int dk = 0; dk < 2; ++dk)                {                    //hash transform                    //Finds the index of the gradient table by the index values of the permutation table                    //The gradient value is found by the index of the gradient table                    //Here, the operands between Perm have no effect, just the color of the result                    c[di][dj][dk] = ran_vec3[perm_x[(i + di) & 255] ^ perm_y[(j + dj) & 255] ^ perm_z[(k + dk) & 255]];                }            }        }        return vec3_trilinear_interp(c, u, v, w);    }};//The output of the Berlin interpolation may be negative,These negative numbers become NaN when Gama corrects them by sqrt()//We should map the output to between 0 and 1class noise_texture : public texture{    //...    public:    virtual vec3 value(double u, double v, const vec3& p) const override    {        return vec3(1, 1, 1) * ( 1.0 + noise.noise( p * scale) ) * 0.5;    }};
      1. 输出:
    6. 噪声叠加(扰动,turblence)

      使用多个频率相加得到复合噪声

      1. 关于倍频系数weight的算法:这个weight是会不断变化的(否则并不会那么自然,也就是起伏过于平缓)

        \[\LARGE frequency = 2^n\\\LARGE amplitude = persistence^i\\其中,frequency为频率,amplitude为振幅\]

      class perlin{    public:    //turbulence    double turb(const vec3& p, int depth = 7)const    {        auto accum = 0.0;        vec3 temp_p = p;        auto weight = 1.0;        for (int i = 0; i < depth; ++i)        {            accum += weight * noise(temp_p);            weight *= 0.5;//amplitude            temp_p *= 2;//frequence        }        return std::fabs(accum);    }}
      1. 输出:
    7. 大理石纹理

      1. 调整颜色:让颜色与sin函数的值成比例

      2. 条纹起伏波荡:并使用turblence函数去调整相位(平移sin函数的参数),使得带状条纹起伏波荡

      3. 可以看到上图的频率很高,若我们只是降低频率是否可以实现大理石纹理呢?

        class noise_texture : public texture{    //...    public:    virtual vec3 value(double u, double v, const vec3& p) const override    {        return vec3(1, 1, 1) * (1 + 5 * noise.turb(p)) * 0.5;    }};
        1. 输出:可以看到并没有大理石纹理,亮度特别高,颜色值不够平缓
      4. 加入sin函数

        class noise_texture : public texture{    //...    public:    virtual vec3 value(double u, double v, const vec3& p) const override    {        return vec3(1, 1, 1) * ( sin( 5 * noise.turb(p) ) ) * 0.5;    }};
        1. 输出:可以看到还差点亮度,应该更白
      5. 实现大理石纹理

        class noise_texture : public texture{    //...    public:    virtual vec3 value(double u, double v, const vec3& p) const override    {        return vec3(1, 1, 1) * ( 1 + sin( scale*p.z() + 10 * noise.turb(p) ) ) * 0.5;    }};

纹理贴图

  • 在这之前我们实现了程序式纹理,我们也能使用图像来进行改变物体外观。也就是,读取一张图片,将2D(u,v)坐标系映射在图片上

  • 纹理贴图可以通过一个通常的纹理管线来描述:

    如下对上图进行概述:

    1. 通过投影方程空间中的点转换为一组参数空间值(也就是u,v坐标),这个值和纹理有关
    2. 使用映射函数将参数空间值转换到纹理空间
    3. 使用纹理空间值从纹理中获取相应的
    4. 使用值变换函数对检索结果进行值变换,再使用得到的新值改变表面属性
    5. 来看个例子:
      1. 找到物体在空间中的位置(x,y,z)
      2. 使用投影函数转换为二元向量(u,v)
      3. 我们假设图像分辨率为(256*256),则分别乘以(u,v)上两个维度的坐标,得到纹理坐标
      4. 通过纹理坐标,在纹理贴图上查找到对应颜色值
      5. 假设找到的颜色值太暗,我们可以使用值变换函数,对每个向量乘以1.1
      6. 最后,我们就可以将作用于值变换函数的向量 用于着色方程
  • 映射函数的实现:

    • 使用(u,v)的直接想法:

      将u和v调整比例后取整,再将其对应到像素坐标(i,j)上

      可行,但比较麻烦。试想每次图片分辨率发生变化,我们都要修改代码

    • 更广泛的做法:

      采用纹理坐标系代替图像坐标系,范围是[0,1]

      如:对一般的图像来说,宽\(N_x\),高\(N_y\),他的一像素(i,j),在纹理坐标系中的坐标为:$$\large u = \frac{i}{N_x – 1}\
      \large v = \frac{j}{N_y – 1}$$

  • 球体uv坐标:

    1. 对于球体来说,uv的计算是基于球面坐标的,只需按比例转化即可得到uv坐标
      当我们有一个球面坐标:\(\large (\theta,\phi)\),其中\theta是朝下距离极轴的角度,\phi是绕极轴旋转的角度

      将球面坐标映射到[0,1]:\(\large u = \frac{\phi}{2 \pi} v = \frac{\theta}{\pi}\)

    2. 计算\(\large \theta 和 \phi\)

      将球面坐标系转化为直角坐标系:

      \[\large x = cos(\theta) cos(\theta) \\\large y = sin(\theta) cos(\phi) \\\large z = sin(\theta)\]

    3. 将第二步倒过来,根据将直角坐标系转化为球面坐标系

      利用库函数atan2(),给出任意一角度的sin和cos值,即可得出这个角的角度值

      //inverse functionstheta = atan2(y,x);//range [-pi, pi]phi = atan2(z);//range [-pi/2, pi/2]
    4. 实现工具函数

      void get_sphere_uv(const vec3& p, double& u, double& v){    auto phi = std::atan2(p.z(), p.x());    auto theta = std::asin(p.y());    u = 1 - (phi + pi) / (2 * pi);//Consider  phi + pi = 2pi    v = (theta + pi / 2) / pi;}

      更新sphere::hit():

      //spere::hitget_sphere_uv( (rec.p - center) / radius, rec.u, rec.v );
  • 存放图片

    我们使用stb_image.h——将图片信息读入一个unsigned char(8bit)的数组中

    //Texture mappingclass image_texture : public texture{    private:    unsigned char* data;    int nx, ny;    public:    image_texture() = default;    image_texture(unsigned char* pixels, int A, int B)        : data(pixels), nx(A), ny(B) {}    ~image_texture()    {        delete data;    }    virtual vec3 value(double u, double v, const vec3& p) const override    {        if (data == nullptr)        {            return vec3(0, 1, 1);        }        //texture-space locations        auto i = static_cast(u * nx);        auto j = static_cast((1 - v) * ny - 0.001); //because v doesn't take into account -2/Pi        if (i < 0)            i = 0;        if (j  nx - 1)            i = nx - 1;        if (j > ny - 1)            j = ny - 1;        auto r = static_cast(data[3 * i + 3 * nx * j + 0] / 255.0);        auto g = static_cast(data[3 * i + 3 * nx * j + 1] / 255.0);        auto b = static_cast(data[3 * i + 3 * nx * j + 2] / 255.0);        return vec3(r, g, b);    }};//main.cpp#define STB_IMAGE_IMPLEMENTATION#include "stb_image.h"hittable_list texture_mapping(){int nx, ny, nn;unsigned char* texture_data = stbi_load("Bronya.jpg", &nx, &ny, &nn, 0);auto bronya_surface = make_shared(make_shared(texture_data, nx, ny));auto bronya = make_shared(vec3(0, 0, 0), 3, bronya_surface);return hittable_list(bronya);}
  • 输出:
    这是我的素材

    这是我的输出结果

矩阵和光源发射光线的材质

  • 加入发射函数。这个材质只要考虑自己发射的光线的颜色,无需考虑反射折射
class material{    public:    //scatter    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const = 0;    //emit    virtual vec3 emitted(double u, double v, const vec3& p) const    {        return vec3(0, 0, 0);    }};class diffus_light : public material{    private:    shared_ptr emit;    public:    diffus_light(shared_ptr a) : emit(a) {}    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const override    {        return false;    }    virtual vec3 emitted(double u, double v, const vec3& p) const    {        return emit->value(u, v, p);    }};
  • 实现一个纯黑背景,并让所有光线都来自光源材质

    在ray_color()中加入背景色变量,并由emitted()递归产生新颜色值

    vec3 ray_color(const ray& r, const vec3& background,  const hittable& world, int depth){    hit_record rec;    if (depth emitted(rec.u, rec.v, rec.p);    if (rec.mat_ptr->scatter(r, rec, attenuation, scattered) == false)    {        return emitted;    }    return emitted + attenuation * ray_color(scattered, background, world, depth - 1);}void output_image(){    //...    const vec3 background(0, 0, 0);    //...    color += ray_color(r, background, world, max_depth);}

加入矩形光源

  1. 实现轴对齐的矩形

    1. 将矩形放在xy平面

      使用z值确定平面,xy确定平面范围。如:若z = k,一个轴对齐的矩形由\(\large x = x_0, x = x_1, y = y_0, y = y_1\)四条直线构成

    2. 判断光线是否于矩形相交

      \[\large 射线的z值为\; z(t) = a_z + t · \overrightarrow{b} \\\large 而z = k时 t的值为\; t = \frac{k – a_z}{\overrightarrow{b_z}} \\\large 即可求解x和y:\; \\\large x = a_x + t · \overrightarrow{b_x} \\\large y = a_y + t · \overrightarrow{b_y} \\\large 其中,若x在[x_0, x_1]且y在[y0,y1],射线就击中了矩阵\]

      //xyclass xy_rect : public hittable{    private:    double x0, x1, y0, y1, k;    shared_ptr mp;    public:    xy_rect() = default;    xy_rect( double _x0, double _x1, double _y0, double _y1, double _k, shared_ptr mat)        : x0(_x0), x1(_x1), y0(_y0), y1(_y1), k(_k), mp(mat) {}    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const override;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const override;};bool xy_rect::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const{    auto t = (k - r.get_origin().z()) / r.get_direction().z();    if (t  t_max)    {        return false;    }    auto x = r.get_origin().x() + t * r.get_direction().x();    auto y = r.get_origin().y() + t * r.get_direction().y();    if (x  x1 || y  y1)    {        return false;    }    rec.u = (x - x0) / (x1 - x0);    rec.v = (y - y0) / (y1 - y0);    rec.t = t;    vec3 outward_normal = vec3(0, 0, 1);    rec.set_face_normal(r, outward_normal);    rec.mat_ptr = mp;    rec.p = r.at(t);    return true;}bool xy_rect::bounding_box(double t0, double t1, aabb& output_box) const{    output_box = aabb(vec3(x0, y0, k - 0.0001), vec3(x1, y1, k + 0.0001));    return true;}//main.cpphittable_list simple_light(){    hittable_list objects;    auto pertext = make_shared(4);    objects.add(make_shared(vec3(0, -1000, 0), 1000, make_shared(pertext)));    objects.add(make_shared(vec3(0, 2, 0), 2, make_shared(pertext)));    auto diff_light = make_shared(make_shared(vec3(4, 4, 4)));    objects.add(make_shared(vec3(0, 7, 0), 2, diff_light));    objects.add(make_shared(3, 5, 1, 3, -2, diff_light));    return static_cast(make_shared(objects, 0.0, 1.0));}
  2. 输出:

    1. 球形光源:
    2. 矩形光源:

Cornell Box

  1. xz平面和yz平面

    //xzclass xz_rect : public hittable{    private:    double x0, x1, z0, z1, k;    shared_ptrmp;    public:    xz_rect() = default;    xz_rect(double _x0, double _x1, double _z0, double _z1, double _k, shared_ptr mat)        : x0(_x0), x1(_x1), z0(_z0), z1(_z1), k(_k), mp(mat) {}    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const override;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const override;};bool xz_rect::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const{    auto t = (k - r.get_origin().y()) / r.get_direction().y();    if (t  t_max)    {        return false;    }    auto x = r.get_origin().x() + t * r.get_direction().x();    auto z = r.get_origin().z() + t * r.get_direction().z();    if (x  x1 || z  z1)    {        return false;    }    rec.u = (x - x0) / (x1 - x0);    rec.v = (z - z0) / (z1 - z0);    rec.t = t;    vec3 outward_normal = vec3(0, 1, 0);    rec.set_face_normal(r, outward_normal);    rec.mat_ptr = mp;    rec.p = r.at(t);    return true;}bool xz_rect::bounding_box(double t0, double t1, aabb& output_box) const{    output_box = aabb(vec3(x0, k - 0.0001, z0), vec3(x1, k + 0.0001, z1));    return true;}//yzclass yz_rect : public hittable{    private:    double y0, y1, z0, z1, k;    shared_ptrmp;    public:    yz_rect() = default;    yz_rect(double _y0, double _y1, double _z0, double _z1, double _k, shared_ptr mat)        : y0(_y0), y1(_y1), z0(_z0), z1(_z1), k(_k), mp(mat) {}    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const override;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const override;};bool yz_rect::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const{    auto t = (k - r.get_origin().x()) / r.get_direction().x();    if (t  t_max)    {        return false;    }    auto y = r.get_origin().y() + t * r.get_direction().y();    auto z = r.get_origin().z() + t * r.get_direction().z();    if (y  y1 || z  z1)    {        return false;    }    rec.u = (y - y0) / (y1 - y0);    rec.v = (z - z0) / (z1 - z0);    rec.t = t;    vec3 outward_normal = vec3(1, 0, 0);    rec.set_face_normal(r, outward_normal);    rec.mat_ptr = mp;    rec.p = r.at(t);    return true;}bool yz_rect::bounding_box(double t0, double t1, aabb& output_box) const{    output_box = aabb(vec3(k - 0.0001, y0, z0), vec3(k + 0.0001, y1, z1));    return true;}
  2. 实现墙壁

    hittable_list cornell_box(){    hittable_list objects;    auto red = make_shared(make_shared(vec3(0.65, 0.05, 0.05)));    auto white = make_shared(make_shared(vec3(0.73, 0.73, 0.73)));    auto green = make_shared(make_shared(vec3(0.12, 0.45, 0.15)));    auto light = make_shared(make_shared(vec3(15, 15, 15)));    objects.add(make_shared(0, 555, 0, 555, 555, green));    objects.add(make_shared(0, 555, 0, 555, 0, red));    objects.add(make_shared(213, 343, 227, 332, 554, light));    objects.add(make_shared(0, 555, 0, 555, 0, white));    objects.add(make_shared(0, 555, 0, 555, 555, white));    objects.add(make_shared(0, 555, 0, 555, 555, white));    return static_cast(make_shared(objects, 0.0, 1.0));}void output_image(){    //...    auto world = cornell_box();    vec3 lookfrom(278, 278, -800);    vec3 lookat(278, 278, 0);    vec3 up_vector(0, 1, 0);    auto dist_to_focus = 10.0;    auto aperture = 0.0;    const auto aspect_ratio = double(image_width) / image_height;    camera cam(lookfrom, lookat, up_vector, 40.0, aspect_ratio, aperture, dist_to_focus, 0.0, 1.0);    //...}
  3. 输出:

  4. 加入两个长方体

    Cornell Box里一般含有两个相对墙面来说有一定角度的长方体

    1. 我们先实现轴对齐的长方体,向其加入六个面

      //Axis-aligned cuboidclass box : public hittable{    private:    vec3 box_min;    vec3 box_max;    hittable_list sides;    public:    box() = default;    box(const vec3& p0, const vec3& p1, shared_ptr ptr);    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const override    {        return sides.hit(r,t_min, t_max, rec);    }    virtual bool bounding_box(double t0, double t1, aabb& output_box) const override    {        output_box = aabb(box_min, box_max);        return true;    }};inline box::box(const vec3& p0, const vec3& p1, shared_ptr ptr){    box_min = p0;    box_max = p1;    sides.add(make_shared(p0.x(), p1.x(), p0.y(), p1.y(), p1.z(), ptr));    sides.add(make_shared(p0.x(), p1.x(), p0.y(), p1.y(), p0.z(), ptr));    sides.add(make_shared(p0.x(), p1.x(), p0.z(), p1.z(), p1.y(), ptr));    sides.add(make_shared(p0.x(), p1.x(), p0.z(), p1.z(), p0.y(), ptr));    sides.add(make_shared(p0.y(), p1.y(), p0.z(), p1.z(), p1.x(), ptr));    sides.add(make_shared(p0.y(), p1.y(), p0.z(), p1.z(), p0.x(), ptr));}
      1. 输出:
    2. 旋转&平移长方体

      在光追中,我们通常用instance来完成旋转、平移

      instance是一种经过旋转或平移等操作的几何图元

      1. 平移

        我们并不需要去移动图元,我们只需对ray进行操作。如下图中将粉色盒子平移至右边的盒子,将位于原点的粉色盒子所有的组成部分的x+2。或者在hit函数中,把射线的原点-2计算得到的是左边的粉色,最后再+2得到右边盒子

        class translate : public hittable{    private:    shared_ptrptr;    vec3 offset;    public:    translate(shared_ptr p, const vec3& displacement) : ptr(p), offset(displacement) {}    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const override;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const override;};bool translate::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const{    ray moved_r(r.get_origin() - offset, r.get_direction(), r.get_time());    if (ptr->hit(moved_r, t_min, t_max, rec) == false)    {        return false;    }    rec.p += offset;    rec.set_face_normal(moved_r, rec.normal);    return true;}bool translate::bounding_box(double t0, double t1, aabb& output_box) const{    if (ptr->bounding_box(t0, t1, output_box) == false)    {        return false;    }    output_box = aabb(output_box.get_min() + offset, output_box.get_max() + offset);    return true;}
      2. 旋转

        从零手写一个软渲染器day3 MVP变换 – 爱莉希雅 – 博客园 (cnblogs.com)一些推导过程

        1. 绕z轴旋转

          \[\large x’ = cos(\theta)·x – sin(\theta)·y \\\large y’ = sin(\theta)·x + cos(\theta)·y\]

        2. 绕y轴旋转

          \[\large x’ = cos(\theta)·x + sin(\theta) · z \\\large z’ = -sin(\theta) · x + cos(\theta)·z\]

        3. 绕x轴旋转

          \[\large y’ = cos(\theta) · y – sin(\theta) · z \\\large z’ = sin(\theta) · y + cos(\theta) · z\]

        4. 旋转时,图元表面的法向量也发生了改变,因此还需重新计算法向量

          class rotate_y : public hittable{    private:    shared_ptr ptr;    double sin_theta;    double cos_theta;    bool hasbox;    aabb bbox;    public:    rotate_y(shared_ptr p, double angle);    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const override;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const override    {        output_box = bbox;        return hasbox;    }};//Calculate the coordinates of the six rotating points of the axially aligned cubeinline rotate_y::rotate_y(shared_ptr p, double angle) : ptr(p){    auto radians = degrees_to_radians(angle);    sin_theta = sin(radians);    cos_theta = cos(radians);    hasbox = ptr->bounding_box(0, 1, bbox);    vec3 min(infinity, infinity, infinity);    vec3 max(-infinity, -infinity, -infinity);    for (int i = 0; i < 2; ++i)    {        for (int j = 0; j < 2; ++j)        {            for (int k = 0; k < 2; ++k)            {                auto x = i * bbox.get_max().x() + (1 - i) * bbox.get_min().x();                auto y = j * bbox.get_max().y() + (1 - j) * bbox.get_min().y();                auto z = k * bbox.get_max().z() + (1 - k) * bbox.get_min().z();                auto newx = cos_theta * x + sin_theta * z;                auto newz = -sin_theta * x + cos_theta * z;                vec3 tester(newx, y, newz);                for (int c = 0; c hit(rotated_r, t_min, t_max, rec) == false)    {        return false;    }    vec3 normal = rec.normal;    vec3 p = rec.p;    p[0] = cos_theta * rec.p[0] + sin_theta * rec.p[2];    p[1] = -sin_theta * rec.p[0] + cos_theta * rec.p[2];    normal[0] = cos_theta * rec.normal[0] + sin_theta * rec.normal[2];    normal[0] = -sin_theta * rec.normal[0] + cos_theta * rec.normal[2];    rec.p = p;    rec.set_face_normal(rotated_r, normal);    return true;}
        5. 输出:

Volume

  • volume rendering(体渲染):用于显示离散三维采样数据集的二维投影的技术。主要用于渲染体积云、体积光

  • volume表示为一个随机表面
    加入volume部分内容会导致代码结构混乱,因为volume和平面表面是完全不同的两种东西。但我们可以将volume表示为一个随机表面。如一团烟雾可以用在概率上不确定位置的平面来代替

  • 实现固定密度的volume

    1. ray可以在volume内部任意一点发生散射,也可以直接穿过volume,volume越薄,直接穿过的概率更大;volume越浓,越可能发生散射
    2. 在任意微小的距离差\(\large \Delta L\)发生散射的概率:
      1. \[\large probability = C· \Delta L \\其中,C是volume的光学密度比例常数\]
    3. 经过一系列不同的等式运算,将会随机得到一个光线发生散射的距离值。根据这个距离值,若散射点在volume外部,则认为没有相交
    4. 对于静态的体积体,我们只需它的密度C和边界
    //The boundary of a volumeclass constant_medium : public hittable{    private:    shared_ptr boundary;    shared_ptr phase_function;    double neg_inv_density;    constant_medium() = default;    public:    constant_medium(shared_ptr b, double d, shared_ptr a)        : boundary(b), neg_inv_density(-1 / d)        {            phase_function = make_shared(a);        }    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const override;    virtual bool bounding_box(double t0, double t1, aabb& output_box) const override    {        return boundary->bounding_box(t0, t1, output_box);    }};bool constant_medium::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const{    //Print  samples when debugging. To enable, set enableDebug true.    const bool enableDebug = false;    const bool debugging = enableDebug && random_double() hit(r, -infinity, infinity, rec1) == false)    {        return false;    }    if (boundary->hit(r, rec1.t + 0.0001, infinity, rec2) == false)    {        return false;    }    if (debugging == true)    {        std::cerr << std::endl << "t0 = " << rec1.t << ", t1 = " << rec2.t << std::endl;    }    if (rec1.t  t_max) rec2.t = t_max;    if (rec1.t >= rec2.t) return false;    if (rec1.t  distance_inside_boundary)    {        return false;    }    rec.t = rec1.t + hit_distance / ray_length;    rec.p = r.at(rec.t);    if (debugging) {        std::cerr << "hit_distance = " << hit_distance << '\n'            << "rec.t = " << rec.t << '\n'            << "rec.p = " << rec.p << '\n';    }    rec.normal = vec3(1, 0, 0);    rec.front_face = true;    rec.mat_ptr = phase_function;    return true;}class isotropic : public material{    private:    shared_ptr albedo;    public:    isotropic(shared_ptr a) : albedo(a) {}    virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const override    {        scattered = ray(rec.p, random_in_unit_sphere(), r_in.get_time());        attenuation = albedo->value(rec.u, rec.v, rec.p);        return true;    }};

最终成品

  • 将所有实现的内容渲染出来,并用体积雾盖住,再加入一个在电解质内部充满volume的蓝色球体
static hittable_list final_scene() {    hittable_list boxes1;    auto ground =        make_shared(make_shared(vec3(0.48, 0.83, 0.53)));    const int boxes_per_side = 20;    for (int i = 0; i < boxes_per_side; i++) {        for (int j = 0; j < boxes_per_side; j++) {            auto w = 100.0;            auto x0 = -1000.0 + i * w;            auto z0 = -1000.0 + j * w;            auto y0 = 0.0;            auto x1 = x0 + w;            auto y1 = random_double(1, 101);            auto z1 = z0 + w;            boxes1.add(make_shared(vec3(x0, y0, z0), vec3(x1, y1, z1), ground));        }    }    hittable_list objects;    objects.add(make_shared(boxes1, 0, 1));    auto light = make_shared(make_shared(vec3(7, 7, 7)));    objects.add(make_shared(123, 423, 147, 412, 554, light));    auto center1 = vec3(400, 400, 200);    auto center2 = center1 + vec3(30, 0, 0);    auto moving_sphere_material =        make_shared(make_shared(vec3(0.7, 0.3, 0.1)));    objects.add(make_shared(center1, center2, 0, 1, 50, moving_sphere_material));    objects.add(make_shared(vec3(260, 150, 45), 50, make_shared(1.5)));    objects.add(make_shared(        vec3(0, 150, 145), 50, make_shared(vec3(0.8, 0.8, 0.9), 10.0)    ));    auto boundary = make_shared(vec3(360, 150, 145), 70, make_shared(1.5));    objects.add(boundary);    objects.add(make_shared(        boundary, 0.2, make_shared(vec3(0.2, 0.4, 0.9))    ));    boundary = make_shared(vec3(0, 0, 0), 5000, make_shared(1.5));    objects.add(make_shared(        boundary, .0001, make_shared(vec3(1, 1, 1))));    int nx, ny, nn;    auto tex_data = stbi_load("Bronya.jpg", &nx, &ny, &nn, 0);    auto emat = make_shared(make_shared(tex_data, nx, ny));    objects.add(make_shared(vec3(400, 200, 400), 100, emat));    auto pertext = make_shared(0.1);    objects.add(make_shared(vec3(220, 280, 300), 80, make_shared(pertext)));    hittable_list boxes2;    auto white = make_shared(make_shared(vec3(0.73, 0.73, 0.73)));    int ns = 1000;    for (int j = 0; j < ns; j++) {        boxes2.add(make_shared(vec3::random(0, 165), 10, white));    }    objects.add(make_shared(        make_shared(            make_shared(boxes2, 0.0, 1.0), 15),        vec3(-100, 270, 395)    )               );    return static_cast(make_shared(objects, 0.0, 1.0));}
  • 输出:借用一下原教程的图,因为太慢了!

reference

Ray Tracing: The Next Week V3.0中文翻译(上) – 知乎 (zhihu.com)

QianMo/Real-Time-Rendering-3rd-CN-Summary-Ebook: 电子书 -《Real-Time Rendering 3rd》提炼总结 | 全书共9万7千余字。你可以把它看做中文通俗版的《Real-Time Rendering 3rd》,也可以把它看做《Real-Time Rendering 3rd》的解读版与配套学习伴侣,或者《Real-Time Rendering 4th》的前置阅读材料。 (github.com)

GAMES101:现代计算机图形学入门 – 计算机图形学与混合现实在线平台 (games-cn.org)

Wikipedia, the free encyclopedia (jinzhao.wiki)

[Nature of Code] 柏林噪声 – 知乎 (zhihu.com)

(26条消息) 什么是马赫带效应_李守聪的博客-CSDN博客_马赫带效应

3D Math Primer for Graphics and Game Development 2nd

gpu gems1