光线追踪渲染器

Created by miccall (转载请注明出处 miccall.tech)

目录

  • 1.基本的向量数学
  • 2.光线追踪入门
  • 3.开始渲染图片
  • 4.渲染一个球
  • 5.使用OpenGL承载
  • 6.法线的计算
  • 7.为了以后方便,我们重构一下代码架构
  • 8.使用摄像机类
  • 9.抗锯齿
  • 10.漫反射
  • 11.伽马 gamma 矫正
  • 12.材质 material
  • 13.镜面模糊
  • 14.透明折射
  • 15.相机角度
  • 16.openGL重现采样过程
  • 17.motion blur

一. 基本的向量数学基础

我们只 使用一个Vec3类表示一个三维向量 包括xyz位置矢量和rgb颜色矢量
代码还包括他的性质和一些代数运算规则

#pragma once
#include <iostream>
#include <math.h>
#include <stdlib.h>
class Vec3
{
public:
    float e[3];
public:
    Vec3() {};
    Vec3(float e0, float e1, float e2) { e[0] = e0; e[1] = e1; e[2] = e2; }
    inline float x() const { return e[0]; }
    inline float y() const { return e[1]; }
    inline float z() const { return e[2]; }
    inline float r() const { return e[0]; }
    inline float g() const { return e[1]; }
    inline float b() const { return e[2]; }

    inline const Vec3 & operator+() const { return *this; }
    inline Vec3 operator-() const { return Vec3(-e[0], -e[1], -e[2]); }
    inline float operator[] (int i) const { return e[i]; }
    inline float& operator[] (int i) { return e[i]; }
    inline Vec3& operator+=(const Vec3 & v2);
    inline Vec3& operator*=(const Vec3 & v2);
    inline Vec3& operator-=(const Vec3 & v2);
    inline Vec3& operator/=(const Vec3 & v2);
    inline Vec3& operator*=(const float t);
    inline Vec3& operator/=(const float t);

    inline float length() const {
        return sqrt(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);
    }

    inline float squared_length() const {
        return e[0] * e[0] + e[1] * e[1] + e[2] * e[2];
    }

    inline void make_unit_vector();
};

//输入流 输入格式
inline std::istream& operator>>(std::istream &is, Vec3 &t) {
    is >> t.e[0] >> t.e[1] >> t.e[2];
    return is;
}

//输出流 输出格式
inline std::ostream& operator<<(std::ostream &os, const Vec3 &t) {
    os << t.e[0] << " " << t.e[1] << " " << t.e[2];
    return os;
}

//生成单位向量 向量/向量的模
inline void Vec3::make_unit_vector() {
    float k = 1.0 / sqrt(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);
    e[0] *= k; e[1] *= k; e[2] *= k;
}

//两向量相加
inline Vec3 operator+(const Vec3 &v1, const Vec3 &v2) {
    return Vec3(v1.e[0] + v2.e[0], v1.e[1] + v2.e[1], v1.e[2] + v2.e[2]);
}

//两向量相减
inline Vec3 operator-(const Vec3 &v1, const Vec3 &v2) {
    return Vec3(v1.e[0] - v2.e[0], v1.e[1] - v2.e[1], v1.e[2] - v2.e[2]);
}

//两向量相乘(for colors)
inline Vec3 operator*(const Vec3 &v1, const Vec3 &v2) {
    return Vec3(v1.e[0] * v2.e[0], v1.e[1] * v2.e[1], v1.e[2] * v2.e[2]);
}

//两向量相除(for colors)
inline Vec3 operator/(const Vec3 &v1, const Vec3 &v2) {
    return Vec3(v1.e[0] / v2.e[0], v1.e[1] / v2.e[1], v1.e[2] / v2.e[2]);
}

//一个标量乘一个向量
inline Vec3 operator*(float t, const Vec3 &v) {
    return Vec3(t*v.e[0], t*v.e[1], t*v.e[2]);
}

//一个向量除以一个标量
inline Vec3 operator/(Vec3 v, float t) {
    return Vec3(v.e[0] / t, v.e[1] / t, v.e[2] / t);
}

//一个向量乘一个标量
inline Vec3 operator*(const Vec3 &v, float t) {
    return Vec3(t*v.e[0], t*v.e[1], t*v.e[2]);
}

//两向量点乘(求内积,得一个常数)(for locations)
inline float dot(const Vec3 &v1, const Vec3 &v2) {
    return v1.e[0] * v2.e[0] + v1.e[1] * v2.e[1] + v1.e[2] * v2.e[2];
}

//两向量叉乘(求外积,得一个向量)(for locations)
inline Vec3 cross(const Vec3 &v1, const Vec3 &v2) {
    return Vec3((v1.e[1] * v2.e[2] - v1.e[2] * v2.e[1]),
        (-(v1.e[0] * v2.e[2] - v1.e[2] * v2.e[0])),
        (v1.e[0] * v2.e[1] - v1.e[1] * v2.e[0]));
}

// += 一个向量
inline Vec3& Vec3::operator+=(const Vec3 &v) {
    e[0] += v.e[0];
    e[1] += v.e[1];
    e[2] += v.e[2];
    return *this;
}

// *= 一个向量
inline Vec3& Vec3::operator*=(const Vec3 &v) {
    e[0] *= v.e[0];
    e[1] *= v.e[1];
    e[2] *= v.e[2];
    return *this;
}

// /= 一个向量
inline Vec3& Vec3::operator/=(const Vec3 &v) {
    e[0] /= v.e[0];
    e[1] /= v.e[1];
    e[2] /= v.e[2];
    return *this;
}

// -= 一个向量
inline Vec3& Vec3::operator-=(const Vec3& v) {
    e[0] -= v.e[0];
    e[1] -= v.e[1];
    e[2] -= v.e[2];
    return *this;
}

// *= 一个数
inline Vec3& Vec3::operator*=(const float t) {
    e[0] *= t;
    e[1] *= t;
    e[2] *= t;
    return *this;
}

// /= 一个数
inline Vec3& Vec3::operator/=(const float t) {
    float k = 1.0 / t;

    e[0] *= k;
    e[1] *= k;
    e[2] *= k;
    return *this;
}

//归一化向量
inline Vec3 unit_vector(Vec3 v) {
    return v / v.length();
}

二. 光线追踪入门

光线使用一个数学模型 R = A + tB 表示 。其中 :

  • R : 射线到达的位置
  • A : origin 表示光线发射的原点。
  • B : Direction 表示光线发射的方向 。
  • t : 一个变量参数 表示光线发射的时间 。

#include "Vec3.h"
class Ray
{
    public:

        Ray();
        Ray(const Vec3 & a , const Vec3 & b );
        Vec3 Origin() const;
        Vec3 Direction() const; 
        Vec3 point_at_parameter(float t) const;

        Vec3 A;
        Vec3 B;
};

Ray::Ray()
{}

Ray::Ray(const Vec3 & a , const Vec3 & b )
{
    A = a;
    B = b;
}

Vec3 Ray::Origin() const 
{
    return A;
}

Vec3 Ray::Direction() const
{
    return B;
}

Vec3 Ray::point_at_parameter( float t ) const
{
    return  A + t * B ;
}

三. 开始渲染图片

我们从一个简单的图片着手,来认识一下光线追踪的基本方法和思路
首先我们选用了一个简单的图片格式 ppm格式的图片作为基本。


#include <fstream>
#include "Ray.h"

Vec3 color(const Ray& r)
{

    Vec3 unit_direction = unit_vector(r.Direction());

    // 因为y的值是 -1到1 ,我们这里使用了一个常见的saturate技巧,重新映射到 (0,1) 
    float t = 0.5 * (unit_direction.y() + 1.0) ;

    // 这是另一个常规技巧 当 t 等于 0 是 Vec3(0.0, 0.0, 0.0) 当t 等于 1 时是Vec3(1.0, 1.0, 1.0) ,  介于 0到1之间。就是两个颜色的混合 。
    return (1.0 - t) * Vec3(1.0, 1.0, 1.0) + t * Vec3(0.0, 0.0, 0.0);

}

int main()
{
    int nx = 200;
    int ny = 100;

    std::ofstream outf;
    outf.open("C:/Users/miccall/Desktop/abc.ppm");
    outf << "P3\n" << nx << " " << ny << "\n255\n";

    // 左下角的坐标 
    Vec3 lower_left_corner(-2.0 , -1.0 , - 1.0 );

    // 水平方向的大小
    Vec3 horizontal(4.0, 0.0, 0.0);

    // 垂直方向的大小 
    Vec3 vertical(0.0,2.0,0.0);

    // 原点坐标 
    Vec3 origin(0.0, 0.0, 0.0);

    // 我们用坐标在原点的光线,射向一个z在 -1 的平面 。 从第一行第一列开始,依次遍历每个像素
    for (int j = ny - 1  ; j  >= 0 ; j -- )
    {
        for (int i = 0; i < nx; i++)
        {

            // UV 计算 
            float u = float(i) / float(nx);
            float v = float(j) / float(ny);

            // 射线
            Ray r(origin, lower_left_corner + u* horizontal + v * vertical);

            // 追踪 
            Vec3 col = color(r);

            // RGB 
            int ir = int(255.99 * col[0]);
            int ig = int(255.99 * col[1]);
            int ib = int(255.99 * col[2]);

            // 数据写入文件 
            outf << ir << " " << ig << " " << ib << " " << "\n";

        }
    }

    outf.close();
    std::cout << "finished" << std::endl;
    getchar();
    return 0;
}

四. 渲染一个球

球的方程
R^2 = x^2 + y^2 + z^2
这是球心在000的情况,如果我们假设球心的坐标为 C = (cx,cy,cz) , 那么等式变为
R^2 = (x-cx)^2 + (y-cy)^2 + (z-cz)^2
在图形学中,我们可以描述为
R^2 = dot( p - C , p - C )
p又是关于t的函数 P = A + t * B ,这是光线追踪方程 。
所以带入可以得到

R^2 = dot( (A + t B - C) , (A + t B - C) )
最后,R,A , B , C都是我们指定的参数,我们要知道的是 t 的值. 只要解出这个方程 .
根据代数知识 ,可以推出

dot(B,B) t ^ 2 + 2 dot( B , A-C ) * t + dot(A-C,A-C) - R^2 = 0
可以看出,这是一个关于t的 一元二次方程 ( ax^2 + bx + c = 0 ) , 我们以前就有办法解出这个答案

bool ishit_sphere(const Vec3 & center ,float Radius , const Ray & r )
{
    Vec3 oc = r.Origin() - center;
    float a = dot(r.Direction(), r.Direction());
    float b = 2.0 * dot( oc , r.Direction());
    float c = dot( oc , oc ) - Radius * Radius ;
    float discriminant = b* b - 4 * a * c;
    return (discriminant > 0);

}

Vec3 color(const Ray& r)
{
    // 射线画一个圆 
    if (ishit_sphere(Vec3(0, 0, -1), 0.5, r))
        return Vec3(1.0,0.0,0.0);

    Vec3 unit_direction = unit_vector(r.Direction());
    float t = 0.5 * (unit_direction.y() + 1.0) ;
    return (1.0 - t) * Vec3(1.0, 1.0, 1.0) + t * Vec3(0.5, 0.7, 1.0);
}

五. 使用OpenGL承载

openGL - glut 基本框架


const GLint     ImageWidth = 800 ;
const GLint     ImageHeight = 400 ;
GLint     PixelLength;
GLbyte * PixelData;
void init(int argc, char** argv, unsigned int mode, int x_position, int y_position, int width, int heigth, const char * title)
{
    // GL 初始化 
    glutInit(&argc, argv);
    // 开启双缓冲机制 和 颜色缓冲 
    glutInitDisplayMode(mode);
    // 窗口位置 
    glutInitWindowPosition(x_position, y_position);
    // 窗口尺寸
    glutInitWindowSize(width, heigth);
    // 窗口名称 
    glutCreateWindow(title);

    // 清空颜色缓存 
    glClearColor(1.0, 0.0, 1.0, 1.0);
    // 矩阵模式 
    /*
    GL_MODELVIEW: 应用这个参数后,表示接下来的矩阵操作都是针对模型视景矩阵堆栈 。  直到下一次调用这个函数并更改参数为止。
    GL_PROJECTION:应用这个参数后,表示接下来的矩阵操作都是针对投影矩阵堆栈    。  直到下一次调用这个函数并更改参数为止。
    GL_TEXTURE :  应用这个参数后,表示接下来的矩阵操作都是针对纹理矩阵堆栈    。  直到下一次调用这个函数并更改参数为止。
    */
    glMatrixMode(GL_PROJECTION);
    //特殊的正射投影函数 left right bottom top 
    gluOrtho2D(0.0, 200.0, 0.0, 150.0);
}

void data()
{
    // 参见下文
}

void pixel()
{
    // 清除颜色缓存
    glClear(GL_COLOR_BUFFER_BIT);

    // 数据加载 
    glDrawPixels(ImageWidth, ImageHeight,
        GL_RGB, GL_UNSIGNED_BYTE, PixelData);

    // 渲染缓存 
    glutSwapBuffers();
}

int main(int argc, char** argv)
{
    data();
    // GL 初始化 自定义初始化内容
    init(argc, argv, GLUT_DOUBLE | GLUT_RGBA, 100, 100, ImageWidth, ImageHeight, "miccall");
    glutDisplayFunc(pixel);
    glutMainLoop();
    return 0;
    return 0;
}

数据加载 :

void data()
{
    // 计算像素数据长度
    PixelLength = ImageWidth * 3;
    while (PixelLength % 3 != 0)
        ++PixelLength;
    PixelLength *= ImageHeight;

    Vec3 lower_left_corner(-2.0, -1.0, -1.0);
    Vec3 horizontal(4.0, 0.0, 0.0);
    Vec3 vertical(0.0, 2.0, 0.0);
    Vec3 origin(0.0, 0.0, 0.0);

    const GLint rgbwidth = ImageWidth * 3;
    PixelData = new GLbyte [PixelLength];

    for (int j = ImageHeight - 1 ; j >= 0; j--)
    {
        for (int i = 0 , k = 0 ; i < rgbwidth; i+=3, k++ )
        {
            float u = float(k) / float(ImageWidth );
            float v = float(j) / float(ImageHeight );

            Ray r(origin, lower_left_corner + u* horizontal + v * vertical);
            Vec3 col = color(r);

            int ir = int(255.99 * col[0]);
            int ig = int(255.99 * col[1]);
            int ib = int(255.99 * col[2]);

            PixelData[ j *(rgbwidth)+(i + 0 )] = (GLbyte) ir;  // R
            PixelData[ j *(rgbwidth)+(i + 1 )] = (GLbyte) ig;  // G
            PixelData[ j *(rgbwidth)+(i + 2 )] = (GLbyte) ib;  // B

        }
    }
}

六.法线计算

我们适当的更改一下代码 , 之前是判断是否射中一个球,是的话,直接填充了一个红色,
那么现在,我们变更为,继续获取他射中的点,然后求出这个点的法线。
法向量的计算方法就是,这个点矢量减去球心矢量 。

实现方法也很简单,把之前那方法稍微改一下,我们直接返回找的的 t 的值 (其中一个根)。

float hit_sphere(const Vec3 & center, float Radius, const Ray & r)
{

    Vec3 oc = r.Origin() - center;
    float a = dot(r.Direction(), r.Direction());
    float b = 2.0 * dot(oc, r.Direction());
    float c = dot(oc, oc) - Radius * Radius;
    float discriminant = b* b - 4 * a * c;
    if (discriminant < 0)
        return -1;
    else
        return (-b - sqrt(discriminant)) / (2.0 * a);
}

随后,我们有了所要的 t 值 ,也就是 t 时 ,光线能够达到的位置 R(t) 。 为了我们方便直观看到法线是否正确,我们直接使用该点的法向量作为他的颜色矢量 。

Vec3 color(const Ray& r)
{
    Vec3 sphereorgin(0,0,-1);
    /* 射线画一个圆 
    if (ishit_sphere(Vec3(0, 0, -1), 0.5, r))
        return Vec3(1.0, 0.0, 0.0);
    */
    float t = hit_sphere(sphereorgin, 0.5, r);
    if (t > 0.0) 
    {
        // 计算法线 
        Vec3 N = unit_vector(r.point_at_parameter(t) - sphereorgin );
        // 同样把 -1 - 1 重映射到 0-1
        return 0.5 * Vec3(N.x()+1 , N.y()+1 ,N.z()+1 );
    }

    Vec3 unit_direction = unit_vector(r.Direction());
    t = 0.5 * (unit_direction.y() + 1.0);
    return (1.0 - t) * Vec3(1.0, 1.0, 1.0) + t * Vec3(0.5, 0.7, 1.0);
}

然后我们就得到了法线信息。

为什么这个图比之前那个图清晰了那么多? 是因为我

const GLint     ImageWidth = 800 ;
const GLint     ImageHeight = 400 ;

而之前那是 200 * 100 的 ,确实提高了很多 。

七. 为了以后方便,我们重构一下代码架构

这样我们就用射线射出了一个球,以后就要反复射其他物体,物体一旦多起来,这确实是一个代码届的苦力活,所以我们重构一下代码,使得复用性更强一些 。

抽象 -> hitable类
这里方便数据传递,做了一个结构体 。表示当光线击中物体时,我们所需要得到的信息 。

#include "Ray.h"

struct hit_record
{
    float t;
    Vec3 p;
    Vec3 normal;
};

class hitable {
public:
    virtual bool hit(const Ray & r, float t_min, float t_max, hit_record & rec) const = 0;
};

然后我们抽象一个球类,这样封装更符合C++特性些


class sphere : public hitable
{
public:
    Vec3 center;
    float Radius;
public:
    sphere() {};
    sphere(Vec3 cen, float r) : center(cen), Radius(r) {}
    virtual bool hit(const Ray & r , float tmin , float tmax , hit_record & rec ) const ;
};

当然这样,我们重写父类的方法,就可以得到球的一个射中情况了。
写之前,几个小方面要注意下。我们在公式的几个部分做了割舍 ,首先计算 b 的时候 ,我们舍弃了前面的系数 2 . 并且在 计算 discriminant的时候,同样舍弃了系数 4 .
随后,在计算 temp的时候 ,同样舍弃了系数 2 . 这样舍弃的目的就是化简计算,结果是没有变化的.

inline bool sphere::hit(const Ray & r, float tmin, float tmax, hit_record & rec) const
{
    Vec3 oc = r.Origin() - center;
    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(b*b - a*c)) / (a) ;
        if (temp < tmax && temp > tmin)
        {
            rec.t = temp;
            rec.p = r.point_at_parameter(rec.t);
            rec.normal = (rec.p - center) / Radius;
            return true;
        }
        temp = (-b + sqrt(b*b - a*c)) / ( a);
        if (temp < tmax && temp > tmin)
        {
            rec.t = temp;
            rec.p = r.point_at_parameter(rec.t);
            rec.normal = (rec.p - center) / Radius;
            return true;
        }
    }
    return false;
}

可以看到,我们不仅仅做了前向的 temp1 而且1在背面做了 temp2 即使我们看不到背面.

可以说到这里就应该可以了,但是为了方便管理,我们多做一步,把射线的所有物体,管理到一个list


#include "hitable.h"

class hitable_list:public hitable
{
public:
    hitable **list;
    int list_size;
public: 
    hitable_list(){}
    hitable_list(hitable ** l, int n) 
    {
        list = l;
        list_size = n;
    }

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

为什么要这样做呢 , 因为物体有远有近 , 我们要按照顺序渲染物体,才不会出现穿插现象 .

bool hitable_list::hit(const Ray & r, float tmin, float tmax, hit_record & rec) const
{
    hit_record temp_rec;
    bool hit_anything = false;
    double closest_so_far = tmax;
    for (int i = 0; i < list_size ; i++)
    {
        if (list[i]->hit(r, tmin, closest_so_far, temp_rec))
        {
            hit_anything = true;
            closest_so_far = temp_rec.t;
            rec = temp_rec; 
        }
    }
    return hit_anything;
}

经历了这么多的变化,我们终于迎来熟悉的步骤,就是重新规划一下Color方法
首先,我们既然用list管理渲染的物体,Color肯定要获取这个数据,那么,我们就得在Color中添加一个参数 .

Vec3 color(const Ray& r, hitable * world)
{
    hit_record rec;
    if (world->hit(r, 0.0, FLT_MAX, rec))
        return 0.5 * Vec3(rec.normal.x() + 1, rec.normal.y() + 1, rec.normal.z() + 1);
    else
    {
        Vec3 unit_direction = unit_vector(r.Direction());
        float t = 0.5 * (unit_direction.y() + 1.0);
        return (1.0 - t) * Vec3(1.0, 1.0, 1.0) + t * Vec3(0.5, 0.7, 1.0);
    }
}

好了,看过代码之后就很清楚了,我们使用world-list里面的物体对射线进行判断,射中的话,就进行法线计算,把法线信息作为他的颜色值.射不中的话,还是一个渐变的背景色.我们已经看过很多遍了 .

好了,一切准备就绪,我么可以开始渲染了 . 第一步,开始加载数据. 这里我们除了渲染之前那个球之外,我们多渲染一个大的球,当作我们的表面 earth


    hitable *list[2];
    list[0] = new sphere(Vec3(0, 0, -1), 0.5);
    list[1] = new sphere(Vec3(0, -100.5, -1), 100);
    hitable *world = new hitable_list(list, 2);

当然,只要在主方法当中调用我们的新方法就ok啦 .

Vec3 col = color(r,world);

八 .使用摄像机类

现在呢 ,同样的为了减少主函数的代码整洁。也为了以后方便,我们新增一个Camera类

#include "Ray.h" 
class Camera
{
public:
    Vec3 origin;
    Vec3 lower_left_corner;
    Vec3 horizontal;
    Vec3 vertical;
public:
    Camera();
    Ray getRay(float u,float v);
};

Camera::Camera()
{
    lower_left_corner = Vec3(-2.0, -1.0, -1.0);
    horizontal = Vec3(4.0,0.0,0.0);
    vertical = Vec3(0.0,2.0,0.0);
    origin = Vec3(0.0, 0.0, 0.0);
}

inline Ray Camera::getRay(float u, float v)
{
    return Ray(origin , lower_left_corner + u*horizontal + v * vertical - origin );
}

gerRay 方法要再细说一下 ,对比之前的方法 ,我们的光线仅仅获取到了 lower_left_corner + uhorizontal + v vertical 方向 ,这里我们在减去一个 orgin 方向 。
当然,之前的理解也当然没问题了,只不过我们现在的摄像机类使用了一个坐标,光线的发射就要以摄像机的位置考虑。

好了,我们就可以在主方法里面使用了 。


    // 声明
    Camera cam;

    // 渲染循环 
    Ray r = cam.getRay(u, v);

九 . 抗锯齿

如果我们放大我们的渲染图片,其实还是可以看到一个明显的锯齿的。
因为我们的渲染器分别渲染了不同的物体,以至于他不知道他们之间的交接处应该是什么,直接把最前面的物体覆盖了上去,就显得和突兀。

那么我们怎么来处理他呢,首先我们想到的就是增加一个采样率,使他获得更多的信息 。

采样率的算法很简单,对于每一个像素,随机产生一个数,利用这个数作为我们新的uv值进行重计算,这样得到的结果就会融合各个角度得到的信息 。

随机函数

static unsigned long long seed = 1;

float drand48()
{
    const long long  m = 0x100000000LL , ra  = 0x5DEECE66DLL ;
    int rc = 0xB16;
    seed = (ra * seed + rc) & 0xFFFFFFFFFFFFLL;
    unsigned int x = seed >> 16;
    return (float) ((double)x / (double)m);
}

增加采样率

    const GLint SamplingRate = 100;


    for (int j = ImageHeight - 1 ; j >= 0; j--)
    {
        for (int i = 0 , k = 0 ; i < rgbwidth; i+=3, k++ )
        {
            Vec3 col(0, 0, 0);
            for (int s = 0; s < SamplingRate ; s++)
            {
                float u = float(k + drand48()) / float(ImageWidth);
                float v = float(j + drand48()) / float(ImageHeight);

                Ray r = cam.getRay(u, v);
                Vec3 p = r.point_at_parameter(2.0);

                col += color(r, world);
            }
            col /= float(SamplingRate);

            // 像素数据 
            PixelData[ j *(rgbwidth)+(i + 0 )] = (GLbyte) int(255.99 * col[0]);  // R
            PixelData[ j *(rgbwidth)+(i + 1 )] = (GLbyte) int(255.99 * col[1]);  // G
            PixelData[ j *(rgbwidth)+(i + 2 )] = (GLbyte) int(255.99 * col[2]);  // B

        }
    }


可以看到,我们使用一个drand48()函数对uv做了更改。并且增加了一层循环 。
这样就会导致运算量得指数增加 。原来几秒算完的结果,此时就变得以分钟衡量 。

放大来看,似乎还是有锯齿,但是我们已经混合了背景色,看出来了吗 ?

真实得照片很难有锯齿,就是因为他的边缘混杂得是背景的颜色,缩放之后就很难发掘。

十. 漫反射

好了,到了我们真正大展身手得时候了,前面铺垫太多,都很难做下去了。今天鼓足了力气,把这部分搞完了。

之前得光线我们只要一次发射,然后就得到颜色,而现在,我们有两个物体,他们之间得空隙,肯定是有光线来回反射的,然而,光线的反射往往不是固定的,就像之前我们所做的那样,来一些随机数是必不可少的。

为了研究方便,我们假设光线在碰到物体之后,有50%的光被吸收,而剩下的光线被反射。我们现在考虑的就是这个反射光线的方向。

那么这个方向到底怎么算呢 ? 查了查网上通用的算法,总结了一个比较简单的方案。
我们假设有一个单位球和一个单位立方体,我们随机生成一个点,要满足这个点在球之外,还必须在立方体之内,否则的话,我们就放弃这个点 。


Vec3 Random_in_unitt_cube()
{
    Vec3 p;
    do
    {
        p = 2.0 * Vec3(drand48(), drand48(), drand48()) - Vec3(1,1,1);

    } while (p.squared_length() >= 1.0 );
    return p;
}

这个点还不能直接用,我们还得考虑法线位置和原始点的位置 。

我们先这样来

    Vec3 target = rec.p + rec.normal + Random_in_unitt_cube();

好了, 目标点有了 ,我们的光线就可以按这个点拐弯了 。


    if (world->hit(r, 0.001, FLT_MAX, rec))
    {
        Vec3 target = rec.p + rec.normal + Random_in_unitt_cube();
        return 0.5 * color(Ray(rec.p, target - rec.p), world);
    }

emmmm , 应该不错的吧 ,跑一个试试 。

emmmm , 跑了不知道多久,反着等的我不耐烦了 。 我们先来写一个方法记录一下过程,不然死等还真烦 。


void showProgress(int num , int sum )
{
    system("cls");
    std::cout << (sum - num) * 100 / sum << "%" << std::endl;
}

    for (int j = ImageHeight - 1 ; j >= 0; j--)
    {
        showProgress(j,ImageHeight);
        for (int i = 0 , k = 0 ; i < rgbwidth; i+=3, k++ )
        {
            // code 
        }
    }

好了,再来一次 。emmm 不多,我的电脑也就跑了两分半钟吧 。。。。

十一 . 伽马 gamma 矫正

由于我们使用的RGB显示方式,当然界外的值我们就不好表示了,画面就会偏黑或者偏白.
这里我们使用gamma矫正的方法来使得画面更加美观一些


Vec3 Gamma_Correct(Vec3 col)
{
    float gammafactor = float(SamplingRate);
    col /= gammafactor ;
    return Vec3(sqrt(col[0]), sqrt(col[1]), sqrt(col[2]));
}


col = Gamma_Correct(col); 

可以对比一些两幅图 ~ 今天的内容很少,我先休息了 ….

十二. 材质 material

终于迎来了我们第一幅真实的画面,想想还是很激动的。但是我们想要更细节表现的话,还要从更多的方向入手,那么,我们就来模拟真实的材质吧。

今天,我们要引入的是一半的漫反射材质和金属材质 ,漫反射材质也就是我们之前写的那样。光的随机方向反射,但是金属就不一样了。他相对磨砂材质,光的反射相对的比较集中,一个极端情况就是完全向预定的反射方向反射,那样就会显得十分的光滑,就跟镜子一样,嗯,你可以想象的到。

好了,我们一步一步来。先来写一个材质类

#include "ray.h"
#include "hitable.h"

class material {
public:
    virtual bool scatter(const Ray& r_in, const hit_record& rec, Vec3& attenuation, Ray& scattered) const = 0;
};

显而易见,这是个抽象方法。我们要实现的是一个具体的材质类,那么就要来继承这个方法,只要调用这个抽象方法,就显得很方便了。
从简单的入手,我们把之前写的漫反射反射封装一下。来看看基本架构

class lambertian : public material
{
public:
    Vec3 albedo;
public:
    lambertian(const Vec3 & a) : albedo(a) {}
    virtual bool scatter(const Ray & r_in , const hit_record& rec , Vec3 & attenuation , Ray & scatterd ) const {
        Vec3 target = rec.p + rec.normal + Random_in_unitt_cube() ;
        scatterd = Ray(rec.p, target - rec.p);
        attenuation = albedo;
        return true;
    }
};

可以看到,其实核心并没有什么变化,变化的是。我们引入了一些额外的变量 。

albedo 反照率
attenuation 光照衰减

这些我们先不用细究太深,我们以后对他进行细化的时候,再来考虑他 。

然后,我们就来模拟金属反射。金属反射可以看作是:

反射的向量是红色所示,我们知道了入射光v和法线n ,根据向量的概念,可以推出反射的光
R = v - 2 B
B = dot( v , n ) * n

    Vec3 reflect(const Vec3& v, const Vec3& n) {
        return v - 2 * dot(v, n)*n;
    }

有了反射光,我们就可以让光线按照反射光的方向继续追踪了。那么金属材料就可以用这个方法写了 。

class metal : public material
{
public:
    Vec3 albedo;
public:
    metal(const Vec3 & a , float f ) : albedo(a),fuzz(f) {}
    virtual bool scatter(const Ray & r_in, const hit_record& rec, Vec3 & attenuation, Ray & scatterd) const {
        Vec3 reflected = reflect(unit_vector(r_in.Direction()), rec.normal);
        scatterd = Ray(rec.p, reflected );
        attenuation = albedo;
        return (dot(scatterd.Direction(), rec.normal) > 0);
    }
};

上面的内容不是很难理解,那我们应该怎么使用这个材质呢 ??

还记得我们定义了一个结构体来保存射中点的信息吗 ?
struct hit_record;
我们为他新添加一个成员,来保存材质的信息 。

struct hit_record
{
    float t;
    Vec3 p;
    Vec3 normal;
    material * mat_ptr;
};

这个你看我是保存了一个材质指针,因为材质是通用的,不用频繁的更改信息,也节省了时间空间。
既然这样,当我们射中这个点时,就会为他赋值不是? 那这个值应该就是物体本身具有的属性了,只有这样,当射中时,我们才能在其之中找到它本身的材质特性。
我们修改球的构造方法 :

sphere(Vec3 cen, float r , material * m ) : center(cen), Radius(r),mat_ptr(m) {}

还有当射中时 :

    if (discriminant > 0)
    {
        float temp ;
        if (temp < tmax && temp > tmin)
        {
            // ... code 
            rec.mat_ptr = mat_ptr;
        }

        if (temp < tmax && temp > tmin)
        {
            // ... code
            rec.mat_ptr = mat_ptr;
            return true;
        }
    }

省略了很多,希望可以看的懂 ~~
修改了球的构造方法,那么我们就要重新声明一个球了 。。


    list[0] = new sphere(Vec3(0, 0, -1), 0.5 , new lambertian(Vec3(0.8,0.3,0.3)));
    list[1] = new sphere(Vec3(0, -100.5, -1), 100 , new lambertian(Vec3(0.8,0.8,0.0)));

有了球,我们还要他的颜色信息,来重写 Color 方法

Vec3 color(const Ray& r, hitable * world,int depth)
{
    hit_record rec;
    if (world->hit(r, 0.001, FLT_MAX, rec))
    {
        Ray scattered;
        Vec3 attenuation; 
        if (depth < 50 && rec.mat_ptr->scatter(r, rec, attenuation, scattered))
        {
            return attenuation * color(scattered, world, depth + 1);
        }
        else
        {
            return Vec3(0, 0, 0);
        }

    }
    else
    {
        Vec3 unit_direction = unit_vector(r.Direction());
        float t = 0.5 * (unit_direction.y() + 1.0);
        return (1.0 - t) * Vec3(1.0, 1.0, 1.0) + t * Vec3(0.5, 0.7, 1.0);
    }
}

为了我们比较材质,和更清晰的认识光线追踪,我们再来添加两个球 。


    hitable *list[4];
    list[0] = new sphere(Vec3(0, 0, -1), 0.5 , new lambertian(Vec3(0.8,0.3,0.3)));
    list[1] = new sphere(Vec3(0, -100.5, -1), 100 , new lambertian(Vec3(0.8,0.8,0.0)));
    list[2] = new sphere(Vec3(1, 0, -1), 0.5, new metal(Vec3(0.8,0.6,0.2),0.3));
    list[3] = new sphere(Vec3(-1, 0, -1), 0.5, new metal(Vec3(0.8, 0.8, 0.8),0.9));
    hitable *world = new hitable_list(list, 4);

好了,一切就绪,又是一个大工程,如果不出意外的话 ,六分钟以后,你会得到一副很漂亮的画面 ;

两个金属球亮的就跟镜面反射一样,我们可以降低一下他的反射效果。

十三. 镜面模糊

我们之前的反射是直接拿公式算出来的,这样未免有一些的规律化,所以会看出,他反射的倒影都是平直的,我们对他进行细微的模糊化。

就像这样,在一定范围内对他随机处理 。

添加方法也很简单,我们只要传递一个模糊参数进来,就在反射光的部分加上一个随机数 。


class metal : public material
{
public:
    Vec3 albedo;
    float fuzz;
public:
    metal(const Vec3 & a , float f ) : albedo(a),fuzz(f) {}
    virtual bool scatter(const Ray & r_in, const hit_record& rec, Vec3 & attenuation, Ray & scatterd) const {
        Vec3 reflected = reflect(unit_vector(r_in.Direction()), rec.normal);
        scatterd = Ray(rec.p, reflected + fuzz * Random_in_unitt_cube());
        attenuation = albedo;
        return (dot(scatterd.Direction(), rec.normal) > 0);
    }
};

我们传一个稍大的和一个稍微小的数试试 。。

    list[2] = new sphere(Vec3(1, 0, -1), 0.5, new metal(Vec3(0.8,0.6,0.2),0.3));
    list[3] = new sphere(Vec3(-1, 0, -1), 0.5, new metal(Vec3(0.8, 0.8, 0.8),0.9));

好了,到目前位置,光线追踪也就算是基本入门了 。我们以后要做的,就是不断对这些内容进行细化和优化。

十四.透明折射

我们之前介绍了光的反射,在自然界中,光不仅仅只要反射,对于那些透明的物体,比如说电介质类,会对光进行透射,在其内部对光的方向做了折射 .

相信我们都知道,根据snell的折射定律,有 n1 * sin(theta1 ) = n2 sin( theta 2 )
其中 n 是在不同介质中的折射率 theta是 光线与法线的夹角 .

我们就是要根据这些信息,来计算出光的折射方向 ,这说简单也简单,说复杂还是有一些复杂的,我花了一些时间写了推导过程 :

总结一下 , 我们如何判断,对于一个半透明物体 , 光线什么时候发生折射现象呢 ?

我们考虑折射光和平面的夹角 theta2 . 如果他的 cos 值大于0 ,他在我们的物体内部,就会发生折射现象 .
如果小于等于0呢? 也就是说 , 如果入射光夹角增大, 大于某个临界角时 , 此时的折射光线会蹦出物体,发射到物体外面,我们称这个现象为全反射 ….

总结一下公社 , 我们用 w 代替 1/e ,就会得到

cos( theta2 ) = sqrt( 1 - pow (w,2) * (1 - pow(cos(theta1),2)));

我们只要判断他的大小就可以判断是否是折射现象
那么 , 如果折射现象,他的反向也可以推导一下 :
w * l + N (w * cos(theta1) - cos(theta2))


// 折射 
bool refract(const Vec3& v, const Vec3& n, float ni_over_nt, Vec3& refracted) {
    Vec3 uv = unit_vector(v);
    float dt = dot(uv, n);
    float discriminant = 1.0 - ni_over_nt*ni_over_nt*(1 - dt*dt);
    if (discriminant > 0) {
        refracted = ni_over_nt*(uv - n*dt) - n*sqrt(discriminant);
        return true;
    }
    else
        return false;
}

然而 ,我们知道了折射,但是,真实中,玻璃的反射率还是不同,因为他的反射率是随着观察角度不同而变化的。在图形学中,有一个简单的公式去计算它 , 叫做菲涅尔-托里克 多项式逼近函数 。

float schlick(float cosine, float ref_idx) {
    float r0 = (1 - ref_idx) / (1 + ref_idx);
    r0 = r0*r0;
    return r0 + (1 - r0)*pow((1 - cosine), 5);
}

然后,我们就可以开始考虑写电解质的材质了 。
首先我们知道,折射这个东西很奇怪,他能从一个面穿入,又能从一个面穿出,所以我们法线的计算就会变得稍微复杂一点,一个球的法线,正反面是相反的,而光线穿过有时候是和法线一致的,有时候又是相反的。所以我们还得根据他的投影来重新计算法线 。

说完了法线,我们还要讨论反射,之前说了反射的不确定性,所以我们得用一种反射范围来限定它,这也是菲尼尔方程的目的所在 。

class dielectric : public material {
public:
    float ref_idx;
public:
    dielectric(float ri) : ref_idx(ri) {}
    virtual bool scatter(const Ray& r_in, const hit_record& rec, Vec3& attenuation, Ray& scattered) const {
        Vec3 outward_normal;
        Vec3 reflected = reflect(r_in.Direction(), rec.normal);
        float ni_over_nt;
        attenuation = Vec3(1.0, 1.0, 1.0);
        Vec3 refracted;
        float reflect_prob;
        float cosine;
        if (dot(r_in.Direction(), rec.normal) > 0) {
            outward_normal = -rec.normal;
            ni_over_nt = ref_idx;
            cosine = dot(r_in.Direction(), rec.normal) / r_in.Direction().length();
            cosine = sqrt(1 - ref_idx*ref_idx*(1 - cosine*cosine));
        }
        else {
            outward_normal = rec.normal;
            ni_over_nt = 1.0 / ref_idx;
            cosine = -dot(r_in.Direction(), rec.normal) / r_in.Direction().length();
        }
        if (refract(r_in.Direction(), outward_normal, ni_over_nt, refracted))
            reflect_prob = schlick(cosine, ref_idx);
        else
            reflect_prob = 1.0;
        if (drand48() < reflect_prob)
            scattered = Ray(rec.p, reflected);
        else
            scattered = Ray(rec.p, refracted);
        return true;
    }
};

然后我们就可以使用这个材质创建新的球了 。

十五.相机角度

我们一般提起摄像机,就想到的就是相机的视野,只有在视野内部的图像,才会被相机所捕获,那我们如何去计算视野呢?


Camera::Camera(float vfov, float aspect)
{
    // 角度换弧度 
    float theta = vfov * M_PI / 180 ;
    float half_height = tan( theta / 2);
    float half_widgh = aspect * half_height;

    // 重新调整视野 
    lower_left_corner = Vec3(-half_widgh, -half_height, -1.0);
    horizontal = Vec3(2 * half_widgh, 0.0, 0.0);
    vertical = Vec3(0.0, 2 * half_height, 0.0);
    origin = Vec3(0.0, 0.0, 0.0);
}

然后我们调整摄像机的位置,就会得到下面的图像 。

但是我们这样不是很好控制摄像机的位置,我们用另一种思路实现一下,就是用一个点表示摄像机的位置,另一个点表示摄像机的观察位置,那么,他的方向就是两个值相减。

Camera::Camera(Vec3 lookfrom, Vec3 lookat, Vec3 vup, float vfov, float aspect)
{
    Vec3 u, v, w;
    float theta = vfov * M_PI / 180;
    float half_height = tan(theta / 2);
    float half_widgh = aspect * half_height;
    origin = lookfrom;
    w = unit_vector(lookfrom - lookat);
    u = unit_vector(cross(vup, w));
    v = cross(w, u);

    lower_left_corner = Vec3(-half_widgh, -half_height, -1.0);
    lower_left_corner = origin - half_widgh * u - half_height * v - w;
    horizontal = 2 * half_widgh * u ;
    vertical = 2 * half_height * v ;
}

然后设置摄像机的位置

Camera cam(Vec3(-2, 2, 1), Vec3(0.5, 0, -1), Vec3(0, 1, 0), 60, float(ImageWidth) / float(ImageHeight));

试着改变参数,你会得到下面的图像。。

十六.openGL重现采样过程

我们实现的是,当一开始运行程序就开始做opengl渲染,而不是全部计算结束之后再对图片进行加载。
我们之前的写法是 ,对于每一个像素,都进行一个SamplingRate的循环进行采样,然后把采样的颜色值进行平均。最后对他做了一个伽马矫正。那么我们现在要在第一时间得到全图,那么我们就先进行一个循环,把颜色值存储下来。以后的循环,在逐一遍历数组,对他进行采样和取平均。

那么我们的方法就改为:


// 主方法 
int main(int argc, char** argv)
{
    data();
    renderdata();
    std::thread t1(func1);
    // GL 初始化 自定义初始化内容
    init(argc, argv, GLUT_DOUBLE | GLUT_RGBA, 100, 100, ImageWidth, ImageHeight, "miccall");
    glutDisplayFunc(pixel);
    glutTimerFunc(10, timerProc, 1);
    glutMainLoop();
    t1.join();
    return 0;
}

主方法先取调用一次 data(); 然后调用 renderdata();
随后就对 pixel 做主循环 ,进行实时渲染 。
然后我们还有一个timerProc的时间方法来通知opengl重新加载数据。


// 数据计算加载 
void data( )
{
        // 计算像素数据长度
        PixelLength = ImageWidth * 3;
        while (PixelLength % 3 != 0)
            ++PixelLength;
        PixelLength *= ImageHeight;

        list[0] = new sphere(Vec3(0, 0, -1), 0.5, new lambertian(Vec3(0.1, 0.2, 0.5)));
        list[1] = new sphere(Vec3(0, -100.5, -1), 100, new lambertian(Vec3(0.8, 0.8, 0.0)));
        list[2] = new sphere(Vec3(1, 0, -1), 0.5, new metal(Vec3(0.8, 0.6, 0.2), 0.3));

        world = new hitable_list(list, 3);
        PixelData = new GLbyte[PixelLength];
        col = new Vec3[ImageHeight*ImageWidth];
}
//渲染数据 预处理
void renderdata()
{
    for ( int j = ImageHeight - 1; j >= 0; j-- )
    {
        for (int i = 0, k = 0; i < rgbwidth; i += 3, k++)
        {

            float u = float(k + drand48()) / float(ImageWidth);
            float v = float(j + drand48()) / float(ImageHeight);

            Ray r = cam.getRay(u, v);
            int index = j*ImageWidth + k;

            col[index] = Vec3(0, 0, 0);

            PixelData[j *(rgbwidth)+(i + 0)] = (GLbyte) int(255.99 * col[index][0]);  // R
            PixelData[j *(rgbwidth)+(i + 1)] = (GLbyte) int(255.99 * col[index][1]);  // G
            PixelData[j *(rgbwidth)+(i + 2)] = (GLbyte) int(255.99 * col[index][2]);  // B

        }
    }
}

这时我们已经有了一个粗略的图像,没有经过采样的 ,然后我们避免主线程阻塞,所以开一个子线程来加载数据 :


void func1()
{
    for (int i = 0; i < SamplingRate && count <= SamplingRate; i++)
    {
        setpixdata();
        count++;
    }
}


//采样 
void setpixdata()
{
    showProgress(SamplingRate - count,SamplingRate);
    std::ofstream outf;

    if (count == SamplingRate)
    {
        outf.open("C:/Users/miccall/Desktop/abc.ppm");
        outf << "P3\n" << ImageWidth << " " << ImageHeight << "\n255\n";
    }

    for (int j = ImageHeight - 1; j >= 0; j--)
    {
        for (int i = 0, k = 0; i < rgbwidth; i += 3, k++)
        {
            float u = float(k + drand48()) / float(ImageWidth);
            float v = float(j + drand48()) / float(ImageHeight);

            Ray r = cam.getRay(u, v);
            int index = j*ImageWidth + k;
            col[index] += color(r, 0);

            colorvec = col[index] / float(count);
            colorvec = Gamma_Correct(colorvec);

            int R = int(255.99 * colorvec[0]);
            int G = int(255.99 * colorvec[1]);
            int B = int(255.99 * colorvec[2]);

            if (count == SamplingRate)
            {
                outf << R << " " << G << " " << B << " \n";
            }

            // 像素数据 
            PixelData[j *(rgbwidth)+(i + 0)] = (GLbyte) R;  // R
            PixelData[j *(rgbwidth)+(i + 1)] = (GLbyte) G;  // G
            PixelData[j *(rgbwidth)+(i + 2)] = (GLbyte) B;  // B

        }
    }

    if (count == SamplingRate)
    {
        outf.close();
        system("cls");
        std::cout << "finished" << std::endl;
    }

}


//循环采样
void timerProc(int id)
{
    /*
    if (count <= SamplingRate)
    {
        setpixdata();
        count++;
    }
    */
    glutPostRedisplay();
    glutTimerFunc(10, timerProc, 1);//需要在函数中再调用一次,才能保证循环  
}

十六.Motion Blur