Games101 学习笔记

第一章

Rasterization 光栅化

将三维空间的几何形体显示在屏幕空间上

实时显示（高于30fps，低于则是离线）

课程主页：https://sites.cs.ucsb.edu/~lingqi/teaching/games101.html

计算机图形学与混合现实在线平台GAMES: http://games-cn.org

http://games-cn.org/forums/topic/allhw/，作业链接地址

虎书https://pan.baidu.com/s/1WSwD7HJv1QBEXBCo2ZMENg 提取：n2rx

第二章 线性代数

向量

单位向量：

向量相加：

笛卡尔坐标系：

向 量乘法：

点乘

点乘dot ,向量点乘得到一个数，点乘可以得到两个向量之间的夹角

作用：1.找到两个向量之间的余弦夹角 2.一个向量投影到另一个向量上的长度

分解向量，可以用来计算这个b向量是靠近两边哪一个向量

3.判断两个向量的前后或者距离

a和b点乘是正数，方向比较一致，a和c点乘是负数，方向差不多相反，如果依据中线对称，则点乘为0，则以此来判断向量的前后关系。

还可以用来判断两个向量的距离，重合为1，然后逐渐远离，由1到0，再远离直到完全相反，由0到-1。可以用于镜面反射，从出射光方向看过去会非常亮，入射光则不然，金属高光同理

叉乘

右手螺旋定则，叉乘方向为手指指向a旋转到b时大拇指的方向，不满足交换律,如果从b旋转到a，则叉乘结果朝下

如果一个坐标系，那就称这个坐标系是右手系，注意UE是左手系

一个向量叉乘他自己等于0的向量，注意叉乘一定得到的是向量

作用：1.判定左和右, 2.判定内与外

上图,得到的向量是正的，所以b在a的左侧，反之则在右侧

对于点P，AP在AB左侧，BP在BC左侧，CP在CA左侧，所以P点位于三边的同一侧，P点在三角形内，不全为同一侧则在三角形外

对于在边上concer case的情况，自己决定内外

标准正交基orthonormal

矩阵

m x n的矩阵是指 m行n列

矩阵乘法

两个矩阵如果想要相乘必须要满足：(m x p)(p x n) = (m x n)

矩阵a和b相乘时先得出矩阵c的大小m x n，第i 行j列元素为a的第i行点乘b的第j列

性质：

1.不满足交换律，AB和BA在一般情况下不相等、

2.具有分配律和结合律，

3.将一个向量看作是m x 1的列矩阵

矩阵乘法作用例子：2D相对y轴镜面反射

矩阵的转置：

转置公式：

单位矩阵和逆矩阵

用于向量

第三章 Transform

缩放

切变shear

旋转矩阵

在数学上，如果一个矩阵的逆等于它的转置，则将这个矩阵称为正交矩阵

所以

线性变换

上面三种变换均为线性变换，可以总结为

平移

齐次坐标

这里平移使用普通的线性变换已经不能满足，需要额外加上一个矩阵才能达到目的。就引入了齐次坐标以将平移写成一个矩阵乘以一个向量的形式。

2D 点 =

2D向量= 之所以2d向量的w分量是0，是因为向量代表方向，具有平移不变性，当这个2d向量进行平移以后，与原来的向量并没有区别，这个0是为保护这个性质(向量平移或者加上一个向量它还是向量，加一个点w值为1就变成了点)。

仿射变换Affine Transformation

仿射变换=线性变换+平移注意做变换一定是先做线性变换，再平移

使用齐次坐标就是：

仿射变换最后一行都是001
变换矩阵都是左乘

这里先旋转再平移

第四章 Transform Cont

3D旋转

对于旋转有：

称之为欧拉角

viewing transformation观测变换

变换顺序：简称MVP变换

• model transformation
• view（这里其实是camera） transformation
• projection transformation

针对一个相机，需要：1.位置 2.相机朝向 3.相机的向上方向 ；去标定相机的状态

通常情况下，将相机置于标准位置，坐标0，0，0，朝向-z，向上方向为y，右手系

projection transformation

orthographic projection

将视锥体（立方体）先平移到0，0，0原点，再缩放成[-1,1]的标准立方体

注意这里摄像机看向的是-z轴，所以其实f（远平面）是小于n（近平面）的（所以OpenGL投影坐标是左手系）

这里缩放矩阵中2/r-l是因为要把边长缩放到-1到1，乘以这个系数因子以缩放到长度为2。

perspective projection

透视投影

考虑到透视投影的矩阵直接写出有难度，可将平截头体frustum的远端向内压缩，形成红色部分的立方体，之后再进行正交投影即可，注意压缩前后n不变，同时中心点蓝色点也是不变的。

将这个压缩操作的矩阵称为

来看这个平截头体的侧视图：

寻找n平面上的点与f平面上的点的对应关系，依据相似三角形得出：

viewport transformation

假定l=-r，b=-t，红线夹角为fovY，宽高比aspect ratio = width/height

由上图，通过fovY和aspect即可得到相应的lrbt

屏幕Screen

什么是屏幕：

• 像素的数组
• 数组大小称为分辨率（1080p，p指的是逐行扫描）
• 是一种典型的光栅化显示
• 光栅化就是在屏幕上进行绘制的过程

补充H1

罗德里格旋转公式

贴一下H1的代码：

float get_degree(float angle)
{
    return (angle / 180) * std::acos(-1);
}

Eigen::Matrix4f get_model_matrix(float rotation_angle)
{
    Eigen::Matrix4f model = Eigen::Matrix4f::Identity();

    // TODO: Implement this function
    // Create the model matrix for rotating the triangle around the Z axis.
    // Then return it.
    float rotation_degree = get_degree(rotation_angle);
    float aa = std::cos(rotation_degree);
    float bb = std::sin(rotation_degree);
    model << aa, -bb, 0, 0, bb, aa, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1;
    return model;
}

Eigen::Matrix4f get_projection_matrix(float eye_fov, float aspect_ratio,
                                      float zNear, float zFar)
{
    // Students will implement this function

    Eigen::Matrix4f projection = Eigen::Matrix4f::Identity();

    float t = zNear * std::tan(get_degree(eye_fov / 2));
    float r = t * aspect_ratio;
    float b = -t;
    float l = -r;
    float halfz = (zNear + zFar) / 2;
    Eigen::Matrix4f TO;
    TO << 1, 0, 0, -(l + r), 
          0, 1, 0, -(t + b), 
          0, 0, 1, -halfz, 
          0, 0, 0, 1;
    Eigen::Matrix4f SO;
    SO << 2 / (r - l) ,0 ,0, 0, 
          0, 2 / (t - b), 0, 0, 
          0, 0, 2/(zNear - zFar),0, 
          0, 0, 0, 1;
    Eigen::Matrix4f Persp;
    Persp << zNear, 0, 0, 0, 
             0, -zNear, 0, 0, //这里为负是为解决左手系三角形显示是倒着的问题
             0, 0, zNear + zFar, -zNear * zFar, 
             0, 0, 1, 0;
    projection = SO * TO *Persp ;

    // TODO: Implement this function
    // Create the projection matrix for the given parameters.
    // Then return it.

    return projection;
}

Eigen::Matrix4f get_rotation(Vector3f axis, float angle)
{
    float rotation_degree = get_degree(angle);
    Eigen::Matrix4f rotation = Eigen::Matrix4f::Identity();
    //std::cout<< rotation <<std::endl;
    Eigen::Matrix4f times;
    float x = axis[0];
    float y = axis[1];
    float z = axis[2];
    times << 0, -z, y, 0,
             z, 0, -x, 0,
             -y, x, 0, 0,
             0, 0, 0,  1;
    Vector4f taxis = {x, y, z, 1};
    rotation = std::cos(rotation_degree) * rotation + (1 -std::cos(rotation_degree))*taxis*taxis.transpose() 
    + std::sin(rotation_degree) * times;
              
    return rotation;
}

第五章Rasterization(Triangle)光栅化-三角形

三角形

为什么是三角形：

• 最简单的多边形
• 可以将其他多边形拆解成三角形

性质：

• 一个三角形一定是一个平面
• 一个三角形可以利用叉乘很明确的定义内外
• 三角形可以依据三个顶点在三角形内部进行插值计算

光栅化

采样：采样就是将一个函数离散化的过程

通过将三角形作为限定条件，对像素中心进行采样，显示中心在三角形内的像素点，从而实现光栅化

（注意对中心点在三角形边上的情况，自定义处理即可）

注意这里是将屏幕上所有像素点都进行了采样，考虑到一个三角形只在一部分范围内存在，使用一个axis along bounding box （AABB）进行包围，并在其中进行采样，可以减少采样成本

第二种优化方法是考虑只处理从边界最左侧到最右侧内的像素点

第六章Rasterization(Antialiasing，z-buffering)反走样和深度缓冲

走样

在计算机图形学中的采样瑕疵sampling artifacts

• jaggies 锯齿-空间中的采样
• moire 摩尔纹-图片降采样
• wagon wheel effect-车轮效应-对时间的采样

本质是：相对来说信号变换过快（频率过高）而采样速度太慢

可以看见，随着信号频率的增加，两次采样之间的变化越多，已逐渐无法通过逆傅里叶变换用两次采样的数据恢复成原有的信号。如下图，恢复出来的信号（黑线）与实际的完全不同了。将黑色线与白色线看作两个函数，会发现采样得到的结果完全相同。

走样：相同的采样方法采样两种不同频率的函数得到完全相同的结果，这种情况就称为走样

如何反走样：在采样前进行一次模糊Blurring（Pre-Filtering滤波）

、

滤波

将左图做一个傅里叶变换就可以得到右图，图中红线方向频率由低变高，而白色部分代表信号（这个信号可以看成两个相邻像素点之间的差异，差异越明显，频率越高，反之，频率越低），可以看出来左图大部分是低频信号。

做一个高通滤波，去掉图中的低频信号，只留下高频信号（即相邻的像素变化差异较大，这种像素一般是边缘部分），做一个逆傅里叶变换，可以得到大概的轮廓图。

做一个低通滤波，去掉图中的高频信号，即去掉边缘细节，只保留低频信息，则整个图像变得模糊，因为失去了边缘和细节部分。

图像可以通过傅里叶变换进行时域到频域的转换，通过逆傅里叶变换进行频域到时域的转换。时域上的卷积等于频域上的乘积，可以看见上图，做一个3x3的卷积核卷积，最后结果类似于低通滤波，只留下了低频信号。

靠。。。这部分有点复杂。。。

采样可以理解为对原始频谱在频域上的重复

稀疏采样的情况下，两次频谱之间重复之间会有部分重合，从而导致走样。

简单理解下：

–对于一个给定的三角形，采样用的像素点越小，需要采样的像素就越多，采样密度就越大，越不容易走样，而采样的像素过大，采样密度就过小，导致了走样。

–所以对于同一张图片，在不同分辨率的屏幕上，高分辨率的屏幕采样密度更大，更不容易走样

反走样

两种减少走样的方法：

• 提高采样率：即增加频域上两个频谱之间的距离，这样在稀疏采样时频谱之间更不容易重叠。实际上就是提高画面分辨率，但开销很高
• 反走样：使用低通滤波去掉高频信号，相当于减少频域上频谱的宽度（见下图），这样在稀疏采样时也不容易重叠。所以要先滤波再采样，反之则不行

反走样的实际应用：使用一个1x1的卷积核对单个像素进行卷积（或者说滤波、平均）

（上图中，黑色部分为在三角形内部的部分，将这部分与这个像素剩下的部分进行平均）

Multi Sample Antialiasing （MSAA）

MSAA在一个像素点设置多个采样点，通过这多个采样点去近似计算一个像素内在三角形内部的比例。

上图是进行2x2MSAA的示意图，对每个在三角形内部和边缘的像素内进行2x2=4次采样，蓝框中的那个像素点，像素比例即为75%

所以MSAA并没有提高分辨率，实际上做的是采样前进行模糊的那一部分工作。

注意和SSAA区别，SSAA在每一个像素内部采样点维护了深度和色值，这个像素色值就是内部采样点的加权平均，等于是一个像素四次shading

而MSAA像素内部只维护深度值，并依据像素内采样点在三角形的比例对像素色值加权，只做一次shading

FXAA只是一个后处理技术，对原图绘制完成后，通过算法识别边缘，然后以像素级别进行混合；

https://blog.csdn.net/weixin_40273050/article/details/121886452

其他

Fast Approximate AA （FXAA）快速近似抗锯齿

与采样无关，在图像层面上进行的抗锯齿，是一种后处理

Temporal AA（TAA）

复用上一帧的信息

DLSS（Deep Learning Super Sampling）

Games101中提及的光栅化是三角形光栅化，直线光栅化没讲，DDA和中点Bresenham也不复杂，先挖个坑后续再看

Z-buffering

为每一个采样的像素点存储一个当前最小的（最近的)z值

需要额外的buffer去存储这个深度值

–frame buffer存储颜色值

–depth buffer存储深度值

（注意我们这里规定，z值永远为正，且z越小，离摄像机越近，越大就越远，且初始化时z值为无限大）

伪代码：针对每一个被光栅化的三角形的其内的一个采样点，记录最小的z值

所以z-buffer永远维护的当前最近的像素的深度信息。

z-buffer计算的时间复杂度：我们只需要遍历场景中的每一个三角形，对其中的采样像素点记录最小z值，所以对于n个三角形，复杂度为O(n)

通过GPU硬件进行计算

注意在使用MASS时，z-buffer就不单纯是针对像素进行，而是针对其中的每个采样点进行带宽会急剧增加，这是个人觉得为啥延迟渲染不适用MSAA的原因

补充H2

H2 作业 光栅化

1. 创建三角形的 2 维 bounding box。
2. 遍历此 bounding box 内的所有像素（使用其整数索引）。然后，使用像素中
   心的屏幕空间坐标来检查中心点是否在三角形内。
3. 如果在内部，则将其位置处的插值深度值 (interpolated depth value) 与深度
   缓冲区 (depth buffer) 中的相应值进行比较。
4. 如果当前点更靠近相机，请设置像素颜色并更新深度缓冲区 (depth buffer)。

static bool insideTriangle(float x, float y, const Vector3f* _v)
{   
    // TODO : Implement this function to check if the point (x, y) is inside the triangle represented by _v[0], _v[1], _v[2]
    Vector3f P = {x, y, 1};
    Vector3f AB, BC, CA, AP, BP, CP;
    
    AB = _v[1] - _v[0];
    BC = _v[2] - _v[1];
    CA = _v[0] - _v[2];
    AP = P - _v[0];
    BP = P - _v[1];
    CP = P - _v[2];

    
    auto r0 = AB.cross(AP);
    auto r1 = BC.cross(BP);
    auto r2 = CA.cross(AP);
    
    return (r0[2] > 0 && r1[2] > 0 && r2[2] > 0) || 
        (r0[2] < 0 && r1[2] < 0 && r2[2] < 0);
}

void rst::rasterizer::rasterize_triangle(const Triangle& t) {
    auto v = t.toVector4();
    
    // TODO : Find out the bounding box of current triangle.
      Eigen::Vector4f v1[] = {
                Vector4f(t.v[0][0], t.v[0][1], t.v[0][2], 1.0),
                Vector4f(t.v[1][0], t.v[1][1], t.v[1][2], 1.0),
                Vector4f(t.v[2][0], t.v[2][1], t.v[2][2], 1.0)
        };
     float tmaxX, tminX, tmaxY, tminY;
    tmaxX = tminX = v1[0][0];
    tmaxY = tminY = v1[0][1];
    for (int i = 1; i < 3; i++)
    {
        if (v1[i][0] > tmaxX)
        {
            tmaxX = v1[i][0];
        }
        if (v1[i][0] < tminX)
        {
            tminX = v1[i][0];
        }
        if (v1[i][1] > tmaxY)
        {
            tmaxY = v1[i][1];
        }
        if (v1[i][1] < tminY)
        {
            tminY = v1[i][1];
        }
        
    }

    int maxX = std::ceil(tmaxX);
    int maxY = std::ceil(tmaxY);
    int minX = std::floor(tminX);
    int minY = std::floor(tminY);

    Eigen::Vector3f v2[] = {
        {v1[0][0], v1[0][1], 1},
         {v1[1][0], v1[1][1], 1},
          {v1[2][0], v1[2][1], 1}
    };
    
    for (int x = minX; x < maxX; x++)
    {
        for (int y = minY; y < maxY; y++)
        {
            if (insideTriangle(x, y, v2))
            {
                /* code */
                auto[alpha, beta, gamma] = computeBarycentric2D(x, y, t.v);
                float w_reciprocal = 1.0/(alpha / v[0].w() + beta / v[1].w() + gamma / v[2].w());
                float z_interpolated = alpha * v[0].z() / v[0].w() + beta * v[1].z() / v[1].w() + gamma * v[2].z() / v[2].w();
                z_interpolated *= w_reciprocal;
                if(depth_buf[get_index(x, y)] > z_interpolated)
                {
                    depth_buf[get_index(x, y)] = z_interpolated;
                    set_pixel({x, y, 1}, t.getColor());
                }
                //set_pixel({x, y, 1}, t.getColor());
            }
             //set_pixel({x, y, 1}, t.getColor());
        }
        
    }
    

    // iterate through the pixel and find if the current pixel is inside the triangle

    // If so, use the following code to get the interpolated z value.
    

    // TODO : set the current pixel (use the set_pixel function) to the color of the triangle (use getColor function) if it should be painted.
}

mass版本

void rst::rasterizer::rasterize_triangle(const Triangle& t) {
    auto v = t.toVector4();
    
    // TODO : Find out the bounding box of current triangle.
      Eigen::Vector4f v1[] = {
                Vector4f(t.v[0][0], t.v[0][1], t.v[0][2], 1.0),
                Vector4f(t.v[1][0], t.v[1][1], t.v[1][2], 1.0),
                Vector4f(t.v[2][0], t.v[2][1], t.v[2][2], 1.0)
        };
    int maxX, minX, maxY, minY;
    maxX = minX = v1[0][0];
    maxY = minY = v1[0][1];
    for (int i = 1; i < 3; i++)
    {
        if (v1[i][0] > maxX)
        {
            maxX = v1[i][0];
        }
        if (v1[i][0] < minX)
        {
            minX = v1[i][0];
        }
        if (v1[i][1] > maxY)
        {
            maxY = v1[i][1];
        }
        if (v1[i][1] < minY)
        {
            minY = v1[i][1];
        }
        
    }

    Eigen::Vector3f v2[] = {
        {v1[0][0], v1[0][1], 1},
         {v1[1][0], v1[1][1], 1},
          {v1[2][0], v1[2][1], 1}
    };
    
    for (int x = minX; x < maxX; x++)
    {
        for (int y = minY; y < maxY; y++)
        {
             Vector3f P[] = {
                {x + 0.25f , y + 0.25f, 1},
                {x + 0.75f , y + 0.75f, 1},
                {x + 0.75f , y + 0.25f, 1},
                {x + 0.25f , y + 0.75f, 1}
                };
            int val = 0;
            float minZ = std::numeric_limits<float>::infinity();
            for (int i = 0; i < 4; i++)
            {
                if (insideTriangle(P[i][0], P[i][1], v2))
                {
                /* code */
                    auto[alpha, beta, gamma] = computeBarycentric2D(P[i][0], P[i][1], t.v);
                    float w_reciprocal = 1.0/(alpha / v[0].w() + beta / v[1].w() + gamma / v[2].w());
                    float z_interpolated = alpha * v[0].z() / v[0].w() + beta * v[1].z() / v[1].w() + gamma * v[2].z() / v[2].w();
                    z_interpolated *= w_reciprocal;
                    if(minZ > z_interpolated)
                    {
                        minZ = z_interpolated;
                    }
                    ++val;
                }
            }
            if(val > 0 && depth_buf[get_index(x, y)] > minZ)
            {
                depth_buf[get_index(x, y)] = minZ;
                set_pixel({x, y, 1}, t.getColor() * val * 0.25f);
            }
        }
    }
}

可以看见锯齿好了很多，要注意的是如果出现黑线是因为深度检测有问题，应该用采样点的最低深度值来代替这个像素点的深度值。

第七、八、九章Shading1

Shading

定义：对一个物体应用一个材质的一个过程

(Blinn-Phong Reflectance Model)冯氏光照模型

Diffuse Lambertian Shading漫反射着色

先定义如下几个参数：

光线与物体表面的碰撞点 shading point

单位向量l，代表光线的光线

单位向量n，物体表面的法线

单位向量v，观测方向

Lamberts 余弦定理，单位面积上接受光线的多少与光线方向l和表面法线n的余弦值有关，所以当两个方向垂直时，接受的光照几乎为0

考虑在空间的光源，每一次发出的能量大小固定，依据能量守恒定律，在能量大小一致的情况下，得到光的强度与传播半径r的平方成反比

综合上面光照接收和光照强度衰减的定义，得到如下公式：

注意这里Kd指的的0-1的rgb值，代表表面吸收的颜色；而max(0, n dot 1)是考虑cos值为负时光线方向为物体表面的背面，在不考虑折射的情况下不具有物理意义。

同时，漫反射会在物体表面均匀反射，与观测方向无关，所以公式与向量v无关

Specular Term （Blinn-Phong）

物体表面看到高光，实际上是观测方向与光线反射的出射方向重合，记出射方向为r，考虑观测方向v与出射方向r的cos值，越大说明越接近，高光越强，即为冯氏高光

但出射光方向r不好计算，这时考虑入射方向l和观测方向v的半程向量h，h平分l和v形成的夹角，同时h十分好计算，通过向量加法l+v再归一化即可。

通过比较这个h向量与法线n的夹角cos值，可得到同样的结果，为布林冯高光，是一种优化。

这里高光颜色ks一般记为白色，同时指针指向的系数p代表了高光系数.

因为在经验中，高光只会在一个很小的角度内生成，而cos值在0-90度内都有效，所以加了一个高光系数做指数去收敛高光有效范围

可以看见，随着高光系数p的减少，高光范围也在逐步缩减

Ambient Term

环境光与光线方向l，观测方向v，表面法线n无关，可以说是一个常数

Reflectance Model

Shading Frequencies

下图分别是对同一个球体：逐面、逐顶点、逐像素进行shading

逐面（Flat shading）

依据每一个三角形面片的法线n去着色，所以这个三角形面片内部颜色不会产生变化，无法表现光滑的表面

逐顶点（Gouraud shading）

求出三角形面片三个顶点的法线，对这三个顶点着色，再依据这三个颜色对三角形进行插值计算

逐像素（Phong shading）

求出三角形面片三个顶点的法线，再依此插值计算出三角形面片内部每一个像素自己的法线，对每一个像素进行着色

每个着色频率与三角形面片数量是息息相关的，并不存在说某种频率效果一定更好，在高面下flat shading也可以有很好的效果，同时如果面数超过像素数，flat shading开销反而会比phong shading大

计算法线

通过求共享顶点v的平面的法线的平均（或依据面大小加权平均）即可求出该顶点的法线

如何计算每个像素的法线见下图

Graphics（Real-time Rendering） Pipeline

可编程渲染管线就是指可以控制顶点阶段和像素阶段怎么着色

一个用于实时预览shader的网站IQ

纹理映射 Texture Mapping

纹理映射定义任意一个点的属性（包括颜色，金属度等）

认为3D空间中的物体表面是一个2D的贴图

对于每一个三角形的顶点，都会被映射到纹理坐标上的某一个点（u，v）

这里uv范围0-1

四方连续纹理tileable texture，可以进行uv贴图的重复

三角形内部插值：重心坐标（Barycentric Coordinates）

对于一个三角形所在平面的点：有这样的性质

如果α β γ均大于0 ，则这个点在三角形内部

对于三角形内部的任一点，可以通过面积比的方式求出来

通过任意的重心坐标点xy，以及三角形三个顶点，可以求出对应的α β γ，

以此就可以把三个顶点的值线性组合起来表示三角形内部任一点的值，这个值可以是位置、纹理坐标、颜色、法线、深度、材质属性等

注意，3维的三角形投影到2维上重心坐标可能会发生变化，所以要先在3d空间中算出重心坐标再投影到2d

纹理映射

纹理坐标定义在三角形面片顶点上，再通过重心坐标插值，可以得到三角形内部任意点的uv，通过这个uv值从纹理中拿到对应的像素颜色，这个颜色就是漫反射系数kd

纹理放大

纹理上的一个像素点称为一个  纹理像素、纹素、texel

Nearset方法：放大后，着色时将一个像素以及周围一定范围的像素如3x3 5x5映射到同一个纹素上。如上图所示就是简单的将目标纹素等同于最近的标准纹素

Bilinear Interpolation 双线性插值：定义两个纹素之间的距离为1，在水平方向s上做两次插值，在竖直方向t上做一次插值，综合考量目标纹素点周围的像素值。

Bicubic Interpolation 双三次插值：双线性是取周围四个像素点，而双三次是取16个像素点，水平竖直都做三次插值

纹理过滤

纹理过大，要将其缩小，反而会有更严重的问题

下图红框部分被称为摩尔纹

实质上缩放后，一个像素点会包含多个纹素，但是采样点是逐像素进行的，所以范围内纹素的信息丢失了，产生了走样。

这时可以将点查询（采样）替换为范围查询

MipMap

允许使用：快速的、近似的、方形的，范围查询

在渲染前先对贴图做预处理，对边长为二次幂的正方形贴图，每次取其边长的一半，做新的分辨率更低的贴图,称为l层，就等于是已经提前做好了一个范围内的像素求值并存储在一张新的贴图内。

可以得出的是，这样做只会增大原有图片内存的三分之一：

那么如何通过mipmap进行范围求值呢

通常，先将屏幕像素点映射回到纹素点上，如下

然后求出映射到纹素上后，目标纹素与周围纹素的距离L，这个L可能有多个值，取最大值，然后将这个L作为边长，就近似得到了一个边长为L的正方形来表示一个区域，这个区域代表着该屏幕像素映射到纹素上所对应的范围（可能包含多个纹素）。

如果使用普通的方法，那这时就是点查询，会丢失这个范围内的其他数据，所以使用去计算出在第D层这个范围会变为一个纹素大小，此时这个像素就代表这片范围的数据平均。这时就是屏幕像素、纹素一一对应了。

但是这个D值会截断成整数，所以实际应用中会出现纹理的接缝断层

所以如果要求非整数层Dx的值，可以对D层和D+1层两层对应范围的纹素各做一个双线性插值，再将两层的纹素值根据层数Dx做一次插值，即可得到平滑的过渡值，这称为trilinear 三线性插值

各向异性过滤 Anisotropic Filtering

Mipmap虽然十分巧妙但也不是完美的，其范围计算是近似的正方形区域，对于纹理映射到斜着的或者不规则的情况下，效果不是很好，所以有了各向异性过滤。

在mipmap中，我们考虑正方形近似计算，所以两个边是等比压缩，但在实际情况中，很多情况下是不会等比的，所以将这部分贴图也预计算进去，mipmap的部分其实就是下图的左上到右下的对角线。

所以理论上各向异性过滤增加了纹理内存三倍的大小。

各向异性过滤同mipmap一样，可以压缩多层，一般为16层，游戏中选择最大第几层其实在显存足够的情况下并无太大区别，因为纹理内存都近似增大三倍。

各向异性过滤因为可以在单方向上压缩，所以对于处理这种纹理映射到一个矩形区域的情况十分有优势，可以计算更精确的矩形范围，相反mipmap是正方形区域，误差更多。

除此之外还有一个EWA filtering使用椭圆进行多次查询以求近似范围，多次查询会很慢，这里不做展开。

纹理应用

可以将纹理理解为 纹理=内存+范围查询（过滤）

环境贴图

假设环境内光照距离无限远，将环境图像映射到贴图上以显示反射。

在一个房间内放一个光滑的球，这个球的表面就可以完整的表现出房间内的各个方向的画面

将环境光映射到球面称为Sphere Map。

但是将SphereMap展开为2dmap时，球的上下两极会出现压缩。

所以可以将球分割为六个部分，映射到正方体的六个面上，称为cube map

bump/normal Map 凹凸/法线贴图

在计算光照时会使用几何表面的法线，我们可以通过将法线信息存储在一张贴图上，在计算光照时使用这额外的信息去扰动正常的法线，从而在不改变几何形体的情况下去改变物体表面的光照从而达到凹凸的效果

在正常情况下，法线为p，使用凹凸贴图后，法线高度被提高，变成了法线n，自然光照结果也会改变

通过导数得出新shading point的切线方向，通过旋转矩阵，旋转90°得到切线的法线。

这里C值为影响因子，是一个系数。

注意计算新法线时，是先将坐标系转换到法线方向为（0，0，1）的切线空间内，计算完新法线后，再切换回原有坐标系。

注意我们通过s移动一个单位后得到一个点，该点减去p点的向量的u和v坐标的值是和p无关的，因为这个点相对于p的v方向的高度就是dp/dv,u方向的高度就是dp/du

Displacement mapping位移贴图

相比于凹凸贴图，位移贴图是真实的改变几何形体，改变了三角形的顶点。凹凸贴图只是模拟光照的变化，对于阴影、物体阴影投影到自身、以及轮廓等效果不尽人意。位移贴图可以很好的做到这些，但要求几何体顶点数足够多以反应位移贴图的变化，DX中可以使用动态曲面细分

AO贴图

3D纹理

3d空间内的一组值

Shadow Mapping

shadow mapping是一种图像空间的做法，计算时其实并不知道场景中的几何信息。同时会产生走样的问题。

主要思想：一个不在阴影中的点，必定能同时被光源和eye看见

主要是两步：

1.将点光源视为一个camera，从光源的方向，将场景投影到一个光源的近平面上（图中红色箭头指指向的平面），相当于光栅化了一遍场景，不过不进行着色，仅存储光源看得见的点的深度值，这张深度图即为shadow map。

2.从eye处进行光栅化，将光栅化时某一个像素点P投影回3D场景P1，再将这个P1点投影回光源的深度图上的某个像素点P2，如果此时P2存储的深度值与P1投影回shadowmap的深度值相同，则这个点能被光源看见（下图橙色的线），则它不在阴影里，反之深度值不同，则这个点P1在阴影里（下图红色的线）

shadowmap相关的问题：

• 只能处理点光源，生成硬阴影（即一个点要么在阴影里，要么不在）
• 由于光源光栅化生成的深度图shadowmap也具有分辨率，所以当其分辨率低时，生成的阴影会产生走样
• 深度比较，由于是浮点值比较，较为困难。

对于有范围的光源，会产生软阴影。

从物理上来说，是本影umbra（完全无光源）和伴影penumbra（有部分光源）的现象。

补充H3

H3作业，做了好久，呜呜

法向量、纹理坐标插值

static Vector4f getInterpolatedVector(const std::tuple<float, float, float>& par, const Vector4f* verts)
{
    Vector4f newVec;
    for(int i = 0; i < 3; ++i)
    {
       newVec[i] = getInterpolatedValue(par, verts, i);
    }
    return newVec;
}

//Screen space rasterization
void rst::rasterizer::rasterize_triangle(const Triangle& t, const std::array<Eigen::Vector3f, 3>& view_pos) 
{
 auto v = t.toVector4();
    
    // TODO : Find out the bounding box of current triangle.
      Eigen::Vector4f v1[] = {
                Vector4f(t.v[0][0], t.v[0][1], t.v[0][2], 1.0),
                Vector4f(t.v[1][0], t.v[1][1], t.v[1][2], 1.0),
                Vector4f(t.v[2][0], t.v[2][1], t.v[2][2], 1.0)
        };
    float tmaxX, tminX, tmaxY, tminY;
    tmaxX = tminX = v1[0][0];
    tmaxY = tminY = v1[0][1];
    for (int i = 1; i < 3; i++)
    {
        if (v1[i][0] > tmaxX)
        {
            tmaxX = v1[i][0];
        }
        if (v1[i][0] < tminX)
        {
            tminX = v1[i][0];
        }
        if (v1[i][1] > tmaxY)
        {
            tmaxY = v1[i][1];
        }
        if (v1[i][1] < tminY)
        {
            tminY = v1[i][1];
        }
        
    }

    int maxX = std::ceil(tmaxX);
    int maxY = std::ceil(tmaxY);
    int minX = std::floor(tminX);
    int minY = std::floor(tminY);
 
    Eigen::Vector4f v2[] = {
        {v1[0][0], v1[0][1], 1, 1},
         {v1[1][0], v1[1][1], 1, 1},
          {v1[2][0], v1[2][1], 1, 1}
    };
    
    for (int x = minX; x < maxX; x++)
    {
        for (int y = minY; y < maxY; y++)
        {
             Vector3f P[] = {
                {x + 0.25f , y + 0.25f, 1},
                {x + 0.75f , y + 0.75f, 1},
                {x + 0.75f , y + 0.25f, 1},
                {x + 0.25f , y + 0.75f, 1}
                };
            int val = 0;
            float minZ = std::numeric_limits<float>::infinity();
            for (int i = 0; i < 4; i++)
            {
                if (insideTriangle(P[i][0], P[i][1], v2))
                {
                /* code */
                    auto[alpha, beta, gamma] = computeBarycentric2D(P[i][0], P[i][1], t.v);
                    float w_reciprocal = 1.0/(alpha / v[0].w() + beta / v[1].w() + gamma / v[2].w());
                    float z_interpolated = alpha * v[0].z() / v[0].w() + beta * v[1].z() / v[1].w() + gamma * v[2].z() / v[2].w();
                    z_interpolated *= w_reciprocal;
                    if(minZ > z_interpolated)
                    {
                        minZ = z_interpolated;
                    }
                    ++val;
                }
            }
            if(val > 0 && depth_buf[get_index(x, y)] > minZ)
            {
                auto par = computeBarycentric2D(x, y, t.v);
                Vector4f normals[] = {
                    {t.normal[0][0], t.normal[0][1], t.normal[0][2], 1},
                    {t.normal[1][0], t.normal[1][1], t.normal[1][2], 1},
                    {t.normal[2][0], t.normal[2][1], t.normal[2][2], 1},
                };
                Vector4f colors[] = {
                    {t.color[0][0], t.color[0][1], t.color[0][2], 1},
                    {t.color[1][0], t.color[1][1], t.color[1][2], 1},
                    {t.color[2][0], t.color[2][1], t.color[2][2], 1},
                };
                Vector4f txs[] = {
                    {t.tex_coords[0][0], t.tex_coords[0][1], 1, 1},
                    {t.tex_coords[1][0], t.tex_coords[1][1], 1, 1},
                    {t.tex_coords[2][0], t.tex_coords[2][1], 1, 1},
                };
                Vector4f normal = getInterpolatedVector(par, normals);
                Vector4f col4 = getInterpolatedVector(par, colors);
                Vector2f tx = {getInterpolatedValue(par, txs, 0), getInterpolatedValue(par, txs, 1)};
               
                fragment_shader_payload shader({col4[0],col4[1],col4[2]}, {normal[0],normal[1],normal[2]}, tx, &texture.value());
                Vector3f col = fragment_shader(shader);
                depth_buf[get_index(x, y)] = minZ;
                set_pixel({x, y}, col * val * 0.25f);
            }
        }
    }
}

blinn-phong反射

static Vector3f vectorScalarMul(const Vector3f& l, const Vector3f& r)
{
    return {l.x() * r.x(), l.y() * r.y(), l.z() * r.z()};
}
Eigen::Vector3f phong_fragment_shader(const fragment_shader_payload& payload)
{
    Eigen::Vector3f ka = Eigen::Vector3f(0.005, 0.005, 0.005);
    Eigen::Vector3f kd = payload.color;
    Eigen::Vector3f ks = Eigen::Vector3f(0.7937, 0.7937, 0.7937);

    auto l1 = light{{20, 20, 20}, {500, 500, 500}};
    auto l2 = light{{-20, 20, 0}, {500, 500, 500}};

    std::vector<light> lights = {l1, l2};
    Eigen::Vector3f amb_light_intensity{10, 10, 10};
    Eigen::Vector3f eye_pos{0, 0, 10};

    float p = 150;

    Eigen::Vector3f color = payload.color;
    Eigen::Vector3f point = payload.view_pos;
    Eigen::Vector3f normal = payload.normal;


    Eigen::Vector3f result_color = {0, 0, 0};
    for (auto& light : lights)
    {
         
        Eigen::Vector3f viewDir = eye_pos - point;
        Eigen::Vector3f lightDir = light.position - point;
        float dis = lightDir.dot(lightDir);
        Eigen::Vector3f hi = (viewDir.normalized() + lightDir.normalized()).normalized();
        Eigen::Vector3f lighta = vectorScalarMul(ka, amb_light_intensity); 

        Eigen::Vector3f lightd = vectorScalarMul(kd , (light.intensity / dis)) * 
        std::max(0.0f, normal.normalized().dot(lightDir.normalized()));

        Eigen::Vector3f lights = vectorScalarMul(ks , (light.intensity / dis)) * 
        std::pow(std::max(0.0f, normal.normalized().dot(hi)), p);
        
        result_color += (lighta + lightd + lights);
    }

    return result_color * 255.f;
}

texture

Eigen::Vector3f texture_fragment_shader(const fragment_shader_payload& payload)
{
    Eigen::Vector3f return_color = {0, 0, 0};
    if (payload.texture)
    {
        // TODO: Get the texture value at the texture coordinates of the current fragment
        return_color = payload.texture->getColorBilinear(payload.tex_coords.x(), payload.tex_coords.y());
    }
    Eigen::Vector3f texture_color;
    texture_color << return_color.x(), return_color.y(), return_color.z();

    Eigen::Vector3f ka = Eigen::Vector3f(0.005, 0.005, 0.005);
    Eigen::Vector3f kd = texture_color / 255.f;
    Eigen::Vector3f ks = Eigen::Vector3f(0.7937, 0.7937, 0.7937);

    auto l1 = light{{20, 20, 20}, {500, 500, 500}};
    auto l2 = light{{-20, 20, 0}, {500, 500, 500}};

    std::vector<light> lights = {l1, l2};
    Eigen::Vector3f amb_light_intensity{10, 10, 10};
    Eigen::Vector3f eye_pos{0, 0, 10};

    float p = 150;

    Eigen::Vector3f color = texture_color;
    Eigen::Vector3f point = payload.view_pos;
    Eigen::Vector3f normal = payload.normal;

    Eigen::Vector3f result_color = {0, 0, 0};

    for (auto& light : lights)
    {
        // TODO: For each light source in the code, calculate what the *ambient*, *diffuse*, and *specular* 
        // components are. Then, accumulate that result on the *result_color* object.
        
        Eigen::Vector3f viewDir = eye_pos - point;
        Eigen::Vector3f lightDir = light.position - point;
        float dis = lightDir.dot(lightDir);
        Eigen::Vector3f hi = (viewDir.normalized() + lightDir.normalized()).normalized();
        Eigen::Vector3f lighta = vectorScalarMul(ka, amb_light_intensity); 

        Eigen::Vector3f lightd = vectorScalarMul(kd , (light.intensity / dis)) * 
        std::max(0.0f, normal.normalized().dot(lightDir.normalized()));

        Eigen::Vector3f lights = vectorScalarMul(ks , (light.intensity / dis)) * 
        std::pow(std::max(0.0f, normal.normalized().dot(hi)), p);
        
        result_color += (lighta + lightd + lights);
    }

    return result_color * 255.f;
}

bump

Eigen::Vector3f bump_fragment_shader(const fragment_shader_payload& payload)
{
    
    Eigen::Vector3f ka = Eigen::Vector3f(0.005, 0.005, 0.005);
    Eigen::Vector3f kd = payload.color;
    Eigen::Vector3f ks = Eigen::Vector3f(0.7937, 0.7937, 0.7937);

    auto l1 = light{{20, 20, 20}, {500, 500, 500}};
    auto l2 = light{{-20, 20, 0}, {500, 500, 500}};

    std::vector<light> lights = {l1, l2};
    Eigen::Vector3f amb_light_intensity{10, 10, 10};
    Eigen::Vector3f eye_pos{0, 0, 10};

    float p = 150;

    Eigen::Vector3f color = payload.color; 
    Eigen::Vector3f point = payload.view_pos;
    Eigen::Vector3f normal = payload.normal;


    float kh = 0.2, kn = 0.1;

    // TODO: Implement bump mapping here
    // Let n = normal = (x, y, z)
    // Vector t = (x*y/sqrt(x*x+z*z),sqrt(x*x+z*z),z*y/sqrt(x*x+z*z))
    // Vector b = n cross product t
    // Matrix TBN = [t b n]
    // dU = kh * kn * (h(u+1/w,v)-h(u,v))
    // dV = kh * kn * (h(u,v+1/h)-h(u,v))
    // Vector ln = (-dU, -dV, 1)
    // Normal n = normalize(TBN * ln)
    Eigen::Vector3f n = normal;
    Eigen::Vector3f t = {n.x() * n.y() / std::sqrt(n.x() * n.x() + n.z() * n.z()),
                        std::sqrt(n.x() * n.x() + n.z() * n.z()),
                        n.z() * n.y() / std::sqrt(n.x() * n.x() + n.z() * n.z())
    };
    Eigen::Vector3f b = n.cross(t);
    Eigen::Matrix3f TBN ;
    TBN << t , b, n;
    float w = payload.texture->width;
    float h = payload.texture->height;
    float u = payload.tex_coords.x();
    float v = payload.tex_coords.y();
    float du = kh * kn * (payload.texture->getColor(u + 1.0f/w, v).norm() -
    payload.texture->getColor(u, v).norm());
    float dv = kh * kn * (payload.texture->getColor(u, v + 1.0f/h).norm() -
    payload.texture->getColor(u, v).norm());
    Eigen::Vector3f ln = {-du, -dv, 1.0f};
    Eigen::Vector3f result_color = (TBN * ln).normalized();

    return result_color * 255.f;
}

displacement

Eigen::Vector3f displacement_fragment_shader(const fragment_shader_payload& payload)
{
    
    Eigen::Vector3f ka = Eigen::Vector3f(0.005, 0.005, 0.005);
    Eigen::Vector3f kd = payload.color;
    Eigen::Vector3f ks = Eigen::Vector3f(0.7937, 0.7937, 0.7937);

    auto l1 = light{{20, 20, 20}, {500, 500, 500}};
    auto l2 = light{{-20, 20, 0}, {500, 500, 500}};

    std::vector<light> lights = {l1, l2};
    Eigen::Vector3f amb_light_intensity{10, 10, 10};
    Eigen::Vector3f eye_pos{0, 0, 10};

    float p = 150;

    Eigen::Vector3f color = payload.color; 
    Eigen::Vector3f point = payload.view_pos;
    Eigen::Vector3f normal = payload.normal;

    float kh = 0.2, kn = 0.1;
    Eigen::Vector3f n = normal;
    Eigen::Vector3f t = {n.x() * n.y() / std::sqrt(n.x() * n.x() + n.z() * n.z()),
                        std::sqrt(n.x() * n.x() + n.z() * n.z()),
                        n.z() * n.y() / std::sqrt(n.x() * n.x() + n.z() * n.z())
    };
    Eigen::Vector3f b = n.cross(t);
    Eigen::Matrix3f TBN ;
    TBN << t , b, n;
    float w = payload.texture->width;
    float h = payload.texture->height;
    float u = payload.tex_coords.x();
    float v = payload.tex_coords.y();
    float du = kh * kn * (payload.texture->getColor(u + 1.0f/w, v).norm() -
    payload.texture->getColor(u, v).norm());
    float dv = kh * kn * (payload.texture->getColor(u, v + 1.0f/h).norm() -
    payload.texture->getColor(u, v).norm());
    Eigen::Vector3f ln = {-du, -dv, 1.0f};
    point = point + kn * n * payload.texture->getColor(u, v).norm();
    normal = (TBN * ln).normalized();
    // TODO: Implement displacement mapping here
    // Let n = normal = (x, y, z)
    // Vector t = (x*y/sqrt(x*x+z*z),sqrt(x*x+z*z),z*y/sqrt(x*x+z*z))
    // Vector b = n cross product t
    // Matrix TBN = [t b n]
    // dU = kh * kn * (h(u+1/w,v)-h(u,v))
    // dV = kh * kn * (h(u,v+1/h)-h(u,v))
    // Vector ln = (-dU, -dV, 1)
    // Position p = p + kn * n * h(u,v)
    // Normal n = normalize(TBN * ln)

    Eigen::Vector3f result_color = {0, 0, 0};

    for (auto& light : lights)
    {
        Eigen::Vector3f viewDir = eye_pos - point;
        Eigen::Vector3f lightDir = light.position - point;
        float dis = lightDir.dot(lightDir);
        Eigen::Vector3f hi = (viewDir.normalized() + lightDir.normalized()).normalized();
        Eigen::Vector3f lighta = vectorScalarMul(ka, amb_light_intensity); 

        Eigen::Vector3f lightd = vectorScalarMul(kd , (light.intensity / dis)) * 
        std::max(0.0f, normal.normalized().dot(lightDir.normalized()));

        Eigen::Vector3f lights = vectorScalarMul(ks , (light.intensity / dis)) * 
        std::pow(std::max(0.0f, normal.normalized().dot(hi)), p);
        
        result_color += (lighta + lightd + lights);
    }

    return result_color * 255.f;
}

纹理采样双线性插值

Eigen::Vector3f getColorBilinear(float u, float v)
    {
        u = std::max(0.0f, u);
        v = std::max(0.0f, v);
        u = std::min(1.0f, u);
        v = std::min(1.0f, v);
        auto lerp = [](float x, const Eigen::Vector3f& v0, const Eigen::Vector3f& v1)
        {
            return v0 + x * (v1 - v0);
        };
        auto getCol = [&](float x, float y)
        {
            auto color = image_data.at<cv::Vec3b>(y, x);
            return Eigen::Vector3f(color[0], color[1], color[2]);
        };

        float u_img = u * width;
        float v_img = (1 - v) * height;
        float lx = std::max(0.0f, std::floor(u_img));
        float rx = std::min(float(width), std::ceil(u_img));
        float by = std::max(0.0f, std::floor(v_img));
        float ty = std::min(float(height), std::ceil(v_img));
        //std::cout<<"lx: " << lx <<"  rx: " <<rx << " by: " << by <<"  ty:  "<<ty  <<" u : "<<u <<" v: " <<v << std::endl;
        Eigen::Vector3f u0 = lerp(u_img - lx, getCol(lx, ty), getCol(rx, ty));
        Eigen::Vector3f u1 = lerp(u_img - lx, getCol(lx, by), getCol(rx, by));

        return lerp(ty - v_img, u0, u1);
    }

第十、十一、十二章Geometry几何

隐式几何，给出一组特殊的点满足的关系式而不给出特定具体的点。

如下图给定一个圆球的方程f，这个f指定了一个圆球的几何形体。

缺点：很难知道整个几何面上的点，直观来说很难知道这个f指定的几何的形状

有点：可以快速的算出一个点是否在这个几何表面上

**显式几何，**给出所有点或者通过参数映射的方式指定，如下图给定参数u和v可以映射到几何形体表面对应的点xyz。

难以判断一个点是否在表面

**距离函数：**定义出了距离一个物体表面任意一点最近的点的距离（有正负号）

距离函数中，值为0的点组合在一起就形成了物体表面，通过对两个物体的距离函数（SDF，signed distance function）做blend，可以根据正负号和0值求出混合后距离函数，再反推出其表示的新的几何形体

点云（显式）：一个包含大量点xyz的列表

点数足够多以形成面，一般是3d扫描出来的形体用点云表示

多边形面（显示）：通常用三角形或者四边形表示

曲线

Bezier Curves 贝塞尔曲线

考虑空间中有三个点b0、b1、b2，有一条关于时间t[0, 1]的曲线，起始点一定在b0，结束点一定在b2，对于给定一个时间t，根据其比例求出线段b0b1上对应的点b01，b1b2线段上对应的点b11，此时得到一个新的线段b01b11，再在此线段上根据t求出对应点b02，这个点就是该贝塞尔曲线在时间为t时对应的值点。

算法以此类推，类似递归，可以处理n个点的情况。

在代数上有：

下图例子是3阶的伯恩斯坦公式，对于0-1的上任意点t，其四条伯恩斯坦曲线值之合一定为1

piecewise bezier curves 分段贝塞尔曲线

可以使用多条3次贝塞尔曲线连接起来表达复杂的曲线，通过四个控制点进行控制

c1连续：两条曲线连续且可导（切线方向和大小是一样的）

Spline 样条

样条：一条通过给定点集合的，连续且导数连续的曲线

曲面

贝塞尔曲面

对于1维的贝塞尔曲线，仅有一个控制参数t，而对于曲面，则有两个0-1范围的参数u、v来控制

以4x4的贝塞尔曲面为例，四条贝塞尔曲线在u0点可以得到4个确定的点，通过这四个点又可以组成一条新的关于时间v的贝塞尔曲线，即图中两条蓝线，以此类推得到所有关于uv的曲线，则确定一个贝塞尔曲面。

几何操作

• 曲面细分 subdivision
• 曲面简化 simplification
• 曲面正规化 regularization（减少非均匀长条状三角形）

曲面细分

两步：

• 创建更多的三角形
• 调整它们的位置

loop细分

先将每一个三角形按如下方式分成四个，然后根据权重调整新顶点的位置

对于新的顶点（白色的点），有四个关联的旧顶点ABCD，经验公式认为A和B对新顶点的影响贡献度更大，所以直接对ABCD点加权平均即得到新顶点位置

对于旧的顶点，根据旧点的度（相连的线段数，图中为6），由经验公式得到u（如果n=3，u=16，否则u=3/8n)，再对旧点和其周围关联的点加权平均。

loop细分只能细分三角形网格

Catmull-Clark 细分

catmull-clark细分可以用于常规的几何体，定义非四边形面（下图中三角形所在的面），定义 奇异点：点的度不为4。

做一次Catmull细分，取面中任意一点，将其连接到面的每一条边的中点上，在一次细分后，所有非四边形面都会消失，同时增加非四边形数个奇异点。

增加新的顶点后按照如下公式加权平均计算新的位置即可

细分结果

曲面简化

使用一种叫边坍缩edge collapsing的方法

如何求出新的顶点位置呢，如果单纯使用平均，效果并不是很好，所以这里引入了Quadric Error Metric 二次度量误差的概念，求出一个新点，它距离相关联的四个平面的距离的平分和有最小值。

求出一条边中点的二次度量误差值，这个值就是这条边的二次度量误差值。在实际细分时，求出模型所有边的二次度量误差，使用优先队列、堆来存储，每次坍缩二次度量误差值最小的边，然后重新计算相关被影响的边的值，再次坍缩当前值最小的边，以此循环。这是一种局部最优的做法。

补充H4

Hw4是做贝塞尔曲线以及反走样，还是比较简单的，使用了递归和代数公式两种办法：

//递归的做法
cv::Point2f recursive_bezier(const std::vector<cv::Point2f> &control_points, float t) 
{
    // TODO: Implement de Casteljau's algorithm
    if (control_points.size() == 1)
    {
       return control_points[0];
    }
    else if (control_points.size() > 1)
    {
        std::vector<cv::Point2f> newPoints;
        for (size_t i = 0; i < control_points.size() - 1; ++i)
        {
           cv::Point2f newPoint = (1 - t) * control_points[i] + t * control_points[i + 1];
           newPoints.push_back(std::move(newPoint));
        }
        return recursive_bezier(std::move(newPoints), t);
    }
    
    return cv::Point2f();
}

//代数的做法
cv::Point2f recursive_bezierMath(const std::vector<cv::Point2f> &control_points, float t) 
{
    auto factorial = [](int x)
    {
        int res = 1;
        for (int i = 2; i <= x; ++i)
            res *= i;
        return res;
    };
    cv::Point2f newPoint;
    int n = control_points.size() - 1;
    for (int i = 0; i < control_points.size(); ++i)
    {
        newPoint += (factorial(n) / (factorial(i) * factorial(n - i))) *
        std::pow(t, i) * std::pow(1 - t, n - i) * control_points[i];
    }
    
    return newPoint;
}

画出曲线以及反走样：我对曲线像素点周围5x5的范围做了颜色插值，缩小看起来效果还是比较明显的

void bezier(const std::vector<cv::Point2f> &control_points, cv::Mat &window) 
{
    //for (double t = 0.0; t <= 1.0; t += 0.0001) 
    {
       // auto point = recursive_bezierMath(control_points, t);

       // window.at<cv::Vec3b>(point.y, point.x)[1] = 255;
    }

    static int max = 8;
    static int pad = 2;
    for (double t = 0.0; t <= 1.0; t += 0.0001) 
    {
        auto point = recursive_bezierMath(control_points, t);
        window.at<cv::Vec3b>(point.y, point.x)[1] = 255;

        for (int i = -pad; i <= pad; i++)
        {
            for (int j = -pad; j <= pad; j++)
            {
                cv::Point2f newPoint(point.x + i, point.y + j);
                if (newPoint.x < 0 || newPoint.x > window.cols || 
                    newPoint.y < 0 || newPoint.y > window.rows||
                    (window.at<cv::Vec3b>(newPoint.y, newPoint.x)[1]!= 0))
                continue;
                float dis = std::pow(newPoint.x - point.x, 2) + 
                std::pow(newPoint.y - point.y, 2);
                window.at<cv::Vec3b>(newPoint.y, newPoint.x)[1] = 
                std::max(255 * (1 - dis / max), 
                (float)window.at<cv::Vec3b>(newPoint.y, newPoint.x)[1]);
            }
        }
    }
}

第十三、十四、十五、十六章 Ray Tracing 光追

光栅化的局限：软阴影，glossy反射，Indirect illumin间接光照

光栅化很快，是实时的，但是光照质量很低，是一种近似的结果

光追很慢，一般是离线的，但光照很准确

做如下假设-光线：

• 光线沿直线传播
• 各条光线之间不会产生碰撞
• 光线具有可逆性，当某条光线从光源出发沿着某条光路到达人眼，人眼也可以从这条光路逆向找到光源

光线投射

光线投射，从眼睛处看见的像素点发射视线 eye ray，碰撞到场景中物体最近的一个点（自然就解决了深度测试的问题），然后求这个点对于光源的方向，以及这个点的法线，然后就可以求得这个像素的颜色。

Recursive(Whitted-Style) Ray Tracing

一条eye ray会生成多条光路，将多条光路的值加权之和即为这个像素的颜色

primary ray：第一次反射或者折射前的eye ray

secondary ray：第二次以及多次反射折射的eye ray

shadow ray：连接光源用于判定可见性的ray

Ray-Surface Intersection

光线追踪的第一个问题就是如何检测eye ray与物体表面的碰撞点呢

一条射线由它的起点和方向进行定义

ray = [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2bey640U-1651142592611)(https://cdn.nlark.com/yuque/__latex/a7a9188995b4ceebdd1d2f577cb15dfa.svg)] [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ObpHgnPG-1651142592611)(https://cdn.nlark.com/yuque/__latex/b29e76e863aec456181fa7b5b1eca702.svg)]

对于一个球体，有[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ybiwsCxx-1651142592611)(https://cdn.nlark.com/yuque/__latex/0af2899450d32db88a3145da8142666e.svg)]

将球的方程与射线方程联立求解（r(t)带入p点），得出的实数非负数根即可求出交点。

所以对于隐式表面，可以很简单的将ray方程带入求解：

如何求解三角形面与ray的交点呢：对每个三角形面做碰撞求解，求最小的t值

每个三角形都在一个确定的平面上，所以转而是求解ray和一个平面的交点。

一个平面可以用一个点和一条法线来确定

同样联立求解，有：

Möller-Trumbore算法

更简单的方法是根据三角形的重心坐标直接求解：

1.求得的t 必须是一个大于0的数。

2.求得的b1 必须是一非负的数且值小于等于1。

3.求得的b2 必须是一非负的数且值小于等于1。

4.求得的b1 + b2 必须是一个小于等于1的数。

补充H5

H5比较简单。主要要将场景像素映射到[-1,1]的屏幕空间（注意扫描从屏幕左上角开始，所以是1-jx），然后用Trumbore求解即可

void Renderer::Render(const Scene& scene)
{
    std::vector<Vector3f> framebuffer(scene.width * scene.height);

    float scale = std::tan(deg2rad(scene.fov * 0.5f));
    float imageAspectRatio = scene.width / (float)scene.height;

    // Use this variable as the eye position to start your rays.
    Vector3f eye_pos(0);
    int m = 0;
    for (int j = 0; j < scene.height; ++j)
    {
        for (int i = 0; i < scene.width; ++i)
        {
            // generate primary ray direction
            float x = (2 * (i+0.5) / scene.width - 1)  * scale * imageAspectRatio;
            float y = (1 - 2 * (j+0.5) / scene.height) * scale;
            // TODO: Find the x and y positions of the current pixel to get the direction
            // vector that passes through it.
            // Also, don't forget to multiply both of them with the variable *scale*, and
            // x (horizontal) variable with the *imageAspectRatio*            

            Vector3f dir = normalize(Vector3f(x, y, -1)); // Don't forget to normalize this direction!
            framebuffer[m++] = castRay(eye_pos, dir, scene, 0);
        }
        UpdateProgress(j / (float)scene.height);
    }

    // save framebuffer to file
    FILE* fp = fopen("binary.ppm", "wb");
    (void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);
    for (auto i = 0; i < scene.height * scene.width; ++i) {
        static unsigned char color[3];
        color[0] = (char)(255 * clamp(0, 1, framebuffer[i].x));
        color[1] = (char)(255 * clamp(0, 1, framebuffer[i].y));
        color[2] = (char)(255 * clamp(0, 1, framebuffer[i].z));
        fwrite(color, 1, 3, fp);
    }
    fclose(fp);    
}

Trumbore

bool rayTriangleIntersect(const Vector3f& v0, const Vector3f& v1, const Vector3f& v2, const Vector3f& orig,
                          const Vector3f& dir, float& tnear, float& u, float& v)
{
    // TODO: Implement this function that tests whether the triangle
    // that's specified bt v0, v1 and v2 intersects with the ray (whose
    // origin is *orig* and direction is *dir*)
    // Also don't forget to update tnear, u and v.
    Vector3f E1 = v1 - v0;
    Vector3f E2 = v2 - v0;
    Vector3f S = orig - v0;
    Vector3f S1 = crossProduct(dir, E2);
    Vector3f S2 = crossProduct(S, E1);
    Vector3f temp(dotProduct(S2, E2), dotProduct(S1, S), dotProduct(S2, dir));
    Vector3f res = (1 / dotProduct(S1, E1)) * temp;
    tnear = res.x;
    u = res.y;
    v = res.z;

    return tnear > 0 && u >= 0 && u <= 1 && v >= 0 && v <= 1 && (u + v) <= 1;
}

Accelerating Ray-Surface Intersection

AABB

对于光追的加速，任意条光线对所有几何体都进行判定是很慢的，可以通过包围盒进行碰撞检测的加速。

我们可以使用一些简单的形体将几何体包围起来，称为包围盒，对于一条光线，可以先检测它是否与包围盒碰撞，如果有碰撞，再检测其与几何体的碰撞。

一般包围盒是一个长方体，我们有如下定义：

一个box是空间中三组对面的交集

如下图所示，box六个面可以组成3组相对平行面，这三组面在3d空间中的交集定义了这个box 的大小范围。

如果这个box任意一条边都与x、y或者z轴平行，则这个box称为 Axis-Aligned Bounding Box（AABB）

那么如何求AABB与射线的交点呢，以2维情况类推

2维的AABB有两组对面，可以求出射线穿过对面的两组 时间tmin和tmax，一组tmin和tmax可以定义一条线段，对两组t求出的线段求交集，这个交集的线段即为射线在这个box内的距离。

在类比到3维的情况

关键是：

• 只有当射线进入所有对面的时候才算进入了box
• 当射线离开任一对面的时候就算离开了box

对于轴对齐面，我们就可以很简单的通过轴分量以及斜率求出t

Uniform grids

有了AABB怎么加速光线，可以使用均匀网格划分的方法。

• 先找出场景的包围盒
• 将其划分成统一的网格
• 储存网格与物体表面相交的信息

假定光线与网格求交远比与物体求交快，当光线进入包围盒后，先求光线与网格的交点并判定这个网格内是否有物体表面，有的话再进行光线与物体表面的求交。

这里判断光线与哪些网格相交的算法类似于直线的光栅化。

对于如何划分一个包围盒，这里有经验值 网格数cells=c * objNums c约为27。

对于场景物体分布不均的情况，uniform grids不是很适合，可能会出现很多额外的网格相交的计算。

Spatial Partitions

大多数场景物体分布不均，更适合于Spatial Partitions 空间划分，一般有三种方法划分：

对于八叉树和BSP树，存在高纬度运算量大或者不是轴对齐的问题，一般采用KD tree

将盒子沿轴方向依次循环划分，直到划分的盒子内物体为空或者其内物体数量小于规定值就停止。

所以所有的物体信息都存储在叶子节点上，中间节点不包含信息。

然后计算光线与叶子节点相交结果。上图方块1应该继续划分的，这里只是展示作用。

但KD tree仍有两个问题：

• 同一个物体可能存在于多个叶子节点（如图中hit的圆存在于3，4，5节点）
• 三角形面与划分盒子的相交十分难计算

Object Partitions & Bounding Volume Hierarchy(BVH)

KD树是考虑将空间进行划分，而BVH是将物体分组进行划分。

考虑有一组物体，先求出他们的包围盒，然后将这组物体分为两部分，再重新求每组的包围盒，这样进行递归计算，直到包围盒满足条件（为空或者物体数量少）。将物体信息存储在叶子节点上。

那么如何划分呢：

• 每次选择一个维度进行划分
• 启发式1：总是选择最长的轴进行划分
• 启发式2：将中间的物体作为分界线进行划分（将三角形重心坐标作为值，用O(n)快速选择算法选择出中位数，进行划分）

BVH数据结构：

• 中间节点：存储包围盒以及左右子节点指针
• 叶子节点：存储包围盒以及物体列表
• 所有的物体信息都存储在叶子节点中

BVH检测伪代码，先求出BVH树。先检测出光线是否与叶子节点相交，如果这个节点不是叶子节点，则检测光线是否与当前节点的两个左右子节点相交，递归计算。

空间划分（KD tree）：

• 划分区域不会重叠
• 一个物体可能属于多个区域

物体划分（BVH）：

• 划分区域可能重叠
• 一个物体只可能属于一个区域

补充H6（SAH）

=-=太菜了，这次作业也磨蹭好久

inline Intersection Triangle::getIntersection(Ray ray)
{
    Intersection inter;

    if (dotProduct(ray.direction, normal) > 0)
        return inter;
    double u, v, t_tmp = 0;
    Vector3f pvec = crossProduct(ray.direction, e2);
    double det = dotProduct(e1, pvec);
    if (fabs(det) < EPSILON)
        return inter;

    double det_inv = 1. / det;
    Vector3f tvec = ray.origin - v0;
    u = dotProduct(tvec, pvec) * det_inv;
    if (u < 0 || u > 1)
        return inter;
    Vector3f qvec = crossProduct(tvec, e1);
    v = dotProduct(ray.direction, qvec) * det_inv;
    if (v < 0 || u + v > 1)
        return inter;
    t_tmp = dotProduct(e2, qvec) * det_inv;

    // TODO find ray triangle intersection
    if (t_tmp < 0)
    {
        return inter;
    }
    
    inter.happened = true;
    inter.coords = ray(t_tmp);
    inter.distance = t_tmp;
    inter.normal = normal;
    inter.m = m;
    inter.obj = this;
    return inter;
}
void Renderer::Render(const Scene& scene)
{
    std::vector<Vector3f> framebuffer(scene.width * scene.height);

    float scale = tan(deg2rad(scene.fov * 0.5));
    float imageAspectRatio = scene.width / (float)scene.height;
    Vector3f eye_pos(-1, 5, 10);
    int m = 0;
    for (uint32_t j = 0; j < scene.height; ++j) {
        for (uint32_t i = 0; i < scene.width; ++i) {
            // generate primary ray direction
            float x = (2 * (i + 0.5) / (float)scene.width - 1) *
                      imageAspectRatio * scale;
            float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
            // TODO: Find the x and y positions of the current pixel to get the
            // direction
            //  vector that passes through it.
            // Also, don't forget to multiply both of them with the variable
            // *scale*, and x (horizontal) variable with the *imageAspectRatio*

            // Don't forget to normalize this direction!
            Vector3f dir = normalize(Vector3f(x, y, -1)); // Don't forget to normalize this direction!
            framebuffer[m++] = scene.castRay(Ray(eye_pos, dir), 0);
        }
        UpdateProgress(j / (float)scene.height);
    }
    UpdateProgress(1.f);

    // save framebuffer to file
    FILE* fp = fopen("binary.ppm", "wb");
    (void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);
    for (auto i = 0; i < scene.height * scene.width; ++i) {
        static unsigned char color[3];
        color[0] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].x));
        color[1] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].y));
        color[2] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].z));
        fwrite(color, 1, 3, fp);
    }
    fclose(fp);    
}

与面相交用前文的射线与点法式联立即可，但要主义这个dirisneg变量。刚开始没有用这个变量，结果出来有问题，后面发现是因为要同时考虑光线源点在pmin左侧，或者pmax右侧的情况。对于第二种情况，射线是反向的则有值，同时tmin和tmax应该交换。

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
                                const std::array<int, 3>& dirIsNeg) const
{
    // invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division
    // dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic
    // TODO test if ray bound intersects
    float tenter = -std::numeric_limits<float>::max();
    float texit = std::numeric_limits<float>::max();
    std::vector<Vector3f> normals = {Vector3f(0, 0, 1), Vector3f(0, 1, 0), Vector3f(1, 0, 0)};
    for (int i = 0; i < 3; ++i)
    {
        /* code */
        float tmax = (float)(pMax[i] - ray.origin[i]) * (float)invDir[i];
        float tmin = (float)(pMin[i] - ray.origin[i]) * (float)invDir[i];
        if (!dirIsNeg[i])
        {
            std::swap(tmax, tmin);
        }
        
        tenter = std::max(tenter, tmin);
        texit = std::min(texit, tmax);
        //std::cout<<"enter : " <<tenter << "  texit: " <<texit <<std::endl;
    }
    return texit > 0 && texit > tenter;
}

递归求解参照伪代码即可

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    // TODO Traverse the BVH to find intersection
    Intersection isect;
    if (node == nullptr)
    {
       return isect;
    }
    
    std::array<int, 3> dirIsNeg = {ray.direction[0] > 0, ray.direction[1] > 0, ray.direction[2] > 0};
    if (!node->bounds.IntersectP(ray, ray.direction_inv, dirIsNeg))
    {
        return isect;
    }
    if (node->object != nullptr)
    {
        return node->object->getIntersection(ray); 
    }
    
    Intersection left = getIntersection(node->left, ray);
    Intersection right = getIntersection(node->right, ray);

    if (left.distance < right.distance && left.happened)
    {
       return left;
    }
    return right;
}

重点是最后的BVH的SAH启发式加速算法，但实现也不是很难。

对于BVN使用轴分割的方法，可能在物体分布不均时造成较差的性能，因为一边密集一边稀疏。

而SAH会综合考虑场景中每一个物体以及分布。

有一个包围盒N，其内有两个包围盒A和B，一条射线进入N，击中A的概率为 N表面积/A表面积 ，B同理，再乘上包围盒A内的物体数NA，B同理，相加即得到击中物体数的期望值。假设一次遍历开销为Ctrav（一般认为是碰撞检测的1/8），一次碰撞检测的开销为Cisect，则总的开销为上图所示。

SAH就是遍历物体数组，每次找到当前的最小开销值进行物体划分，是一种贪心算法。

为了方便计算，我们会先将物体按照某个轴的顺序分组，即分成一个一个的桶，一般小于32个桶（我代码中为12个桶），将这些桶作为划分的最小的子元素进行计算，比如计算b0 以及 b1-b7 ，b0b1 和b2-b7 …b0-b6和b7

的开销，找出最小的开销，即为最优的划分组合。

实现按照上图伪代码即可。

这篇文章解析SVH很好https://blog.csdn.net/qq_39300235/article/details/106999514

SVH只是一个划分的优化方法，其他递归构建BVH的地方没有差异，实现较为简单。

BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects)
{
    BVHBuildNode* node = new BVHBuildNode();

    // Compute bounds of all primitives in BVH node
    Bounds3 bounds;
    for (int i = 0; i < objects.size(); ++i)
        bounds = Union(bounds, objects[i]->getBounds());
    if (objects.size() == 1) {
        // Create leaf _BVHBuildNode_
        node->bounds = objects[0]->getBounds();
        node->object = objects[0];
        node->left = nullptr;
        node->right = nullptr;
        return node;
    }
    else if (objects.size() == 2) {
        node->left = recursiveBuild(std::vector{objects[0]});
        node->right = recursiveBuild(std::vector{objects[1]});

        node->bounds = Union(node->left->bounds, node->right->bounds);
        return node;
    }
    else {
        Bounds3 centroidBounds;
        for (int i = 0; i < objects.size(); ++i)
            centroidBounds =
                Union(centroidBounds, objects[i]->getBounds().Centroid());
        std::vector<Object*> leftshapes;
        std::vector<Object*> rightshapes;
        if (splitMethod == SplitMethod::NAIVE)
        {
            int dim = centroidBounds.maxExtent();
        switch (dim) {
        case 0:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().x <
                       f2->getBounds().Centroid().x;
            });
            break;
        case 1:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().y <
                       f2->getBounds().Centroid().y;
            });
            break;
        case 2:
            std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                return f1->getBounds().Centroid().z <
                       f2->getBounds().Centroid().z;
            });
            break;
        }

        auto beginning = objects.begin();
        auto middling = objects.begin() + (objects.size() / 2);
        auto ending = objects.end();

        leftshapes = std::vector<Object*>(beginning, middling);
        rightshapes = std::vector<Object*>(middling, ending);
        }
        else if (splitMethod == SplitMethod::SAH)
        {
            float minCost = std::numeric_limits<float>::infinity();
            int minAxis = 0;
            int MaxBucketCount = 12;
            const float costIter = 0.125f;
            struct Bucket
            {
                Bounds3 bound;
                int objNums = 0;
            };

            struct BucketAxis
            {
                int i = 0; //012 ---xyz
                std::vector<Bucket> buckets;
                std::map<int, int> objToBuckets;
                int minDivideIndex = 0;
                BucketAxis(int axis, int bucketCount):i(axis)
                {
                    buckets.resize(bucketCount);
                    Bucket newBucket;
                    std::fill(buckets.begin(), buckets.end(), newBucket);
                }
            };

            std::vector<BucketAxis> bucketAxises;

            for (int i = 0; i < 3; ++i)
            {
                BucketAxis newBucketAxis(i, MaxBucketCount);
                for (int j = 0; j < primitives.size(); ++j)
                {
                    int objIndex = centroidBounds.Offset(primitives[i]->getBounds().Centroid())[i] * MaxBucketCount;
                    if (objIndex == MaxBucketCount)
                    {
                        objIndex = MaxBucketCount - 1;
                    }
                    Bounds3 bound = newBucketAxis.buckets[objIndex].bound;
                    newBucketAxis.buckets[objIndex].bound = Union(bound, primitives[j]->getBounds());
                    newBucketAxis.buckets[objIndex].objNums++;
                    newBucketAxis.objToBuckets.insert(std::make_pair(j, objIndex));
                }
                
                for (int k = 1; k < MaxBucketCount; ++k) //divide the buckets
                {
                    float costA = 0, costB = 0;
                    Bounds3 boundA, boundB;
                    int sizeA = 0, sizeB = 0;
                    for (int m = 0; m < k; ++m)
                    {
                        boundA = Union(boundA, newBucketAxis.buckets[m].bound);
                        sizeA += newBucketAxis.buckets[m].objNums;
                    }
                    for (int m = k; m < MaxBucketCount; ++m)
                    {
                        boundB = Union(boundB, newBucketAxis.buckets[m].bound);
                        sizeB += newBucketAxis.buckets[m].objNums;
                    }
                    
                    float cost = costIter + (boundA.SurfaceArea() * sizeA + 
                        boundB.SurfaceArea() * sizeB) / centroidBounds.SurfaceArea();
                    if (cost < minCost)
                    {
                        minAxis = i;
                        newBucketAxis.minDivideIndex = k;
                    }
                }
                bucketAxises.push_back(std::move(newBucketAxis));
            }
            for (int i = 0; i < primitives.size(); ++i)
            {
                bucketAxises[minAxis].objToBuckets[i] < bucketAxises[minAxis].minDivideIndex ? 
                leftshapes.push_back(primitives[i]) :  rightshapes.push_back(primitives[i]);
            }
        }
        assert(objects.size() == (leftshapes.size() + rightshapes.size()));

        node->left = recursiveBuild(leftshapes);
        node->right = recursiveBuild(rightshapes);

        node->bounds = Union(node->left->bounds, node->right->bounds);          
    }

    return node;
}

附上邦尼的照片~

Basic radiometry 辐射度量学

闫老师的学习论~

Radiant Energy and Flux(Power)

定义：辐射的能量就是电磁辐射的能量，使用单位 焦来定义：

定义：Radiant flux(power)，是单位时间内发射、反射、传递或者接收的能量，定义为

相关的物理量：

Radiant Intensity

定义：radiant（luminous）intensity是一个点光源在单位立体角（solid power）内所发出的能量。

二维空间内的角是弧长除以半径，推广到三维就是面积除以半径平方。

单位立体角对整个球面做积分，结果为4pi

所以对于一个点光源，其intensity为Radiant flux(power) / 4pi

例子：对于一个815lumens的点光源有

Irradiance

定义：irradiance是单位面积上垂直（非垂直要投影到垂直方向）方向的能量。

所以phong光照模型中，说光照的intensity衰减其实是不对的，因为随着能量传播半径增加，其单位立体角大小并没有变换，所以intensity没有变，但是单位面积A变大了，所以衰减的其实是irradiance。

Radiance

定义：radiance是表面每单位面积上每单位立体角所发射、反射、传递或者接收的能量

简单理解就是具有方向性的Irradiance，在单位面积上对radiance做积分，得到的就是irradiance。

E是irradiance，L是radiance

Bidirectional Reflectance Distribution Function（BRDF）双向反射分布函数

BRDF描述的是从每个入射的光会反射多少光到出射方向w上，简单理解就是描述物体表面如何进行反射。

对所有入射方向的光照的radiance应用BRDF（fr）做球面积分，则可以得到出射的总的radiance（Lr）。

The Rendering Equation 渲染公式

考虑BRDF以及物体自发光的情况，做球面积分最后得到总的radiance，即渲染公式。

针对面光源，将其视为无数个点光源的集合做积分，而其他物体表面反射过来的光，也可以将那个物体的表面作为面光源进行计算。

中间推导好复杂。。。略过

最后可以将公式简化为[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-JHx4l4fD-1651142592623)(https://cdn.nlark.com/yuque/__latex/a0d6fb5d821f5638275ff3420e819cbb.svg)]的形式

展开计算

可以看出渲染公式考虑了n次弹射的情况。当光线不弹射，即[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Jpz7JD0F-1651142592623)(https://cdn.nlark.com/yuque/__latex/8e7b54d418e531d3484258633d4d270e.svg)]时，只有自发光。当光源只反射一次，即[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XpI1D5xF-1651142592624)(https://cdn.nlark.com/yuque/__latex/0616d99692b476f05c57597c6aaaff54.svg)],则只有直接光照，而在这之上的反射之和，称之为间接光照。

全局光照 = 直接光照 + 间接光照

只有直接光照，只反射一次，所以光照无法照射到的地方一定是黑的。

这里额外弹射一次，使用全局光照，P点就可见了。注意图片上面的灯是黑的，这是因为其材质是透明的，光线会进入其中，需要多次反射才会出来，现在额外一次弹射（两次反射一共）只会照亮灯的外表面，还需要一次反射进入，再需要一次反射射出。可以看见四次反射时灯的物体就被照亮了。

随着反射次数的增多，场景会越来越亮，最终收敛到一个恒定的亮度。

概率论回顾：对离散的随机变量，期望等于加权平均，对于概率密度函数（Probability Density Functions， PDF），期望等于积分

Monte Carlo Path Tracy 蒙特卡洛路径追踪

Monte Carlo Integration蒙特卡洛积分

对函数值进行随机采样从而得到近似的积分结果

通过采样某个点（上图红点、绿点），然后将其近似为矩形，重复n次，对矩形面积和求平均，即可得到积分值。

又因为概率函数p是固定的 1 / （b - a），所以可以写成：

蒙特卡洛要求：在X上采样，就在X上积分

Path Tracing

whitted-style ray tracing有两个问题：

• 只进行镜面反射、界面折射
• 遇见漫反射表面会停止

whitted不能实现图二这种glossy反射

whitted也无法做到全局光照，没办法实现光渗现象（Light Bleeding）

所以使用渲染方程去实现：

但渲染方程也有两个待解决问题：

• 求半球积分
• 递归执行

现在假设从像素投影到场景中的某个点，求直接光照，即这个点点某条光线直接打到光源

前提：不考虑自发光，只考虑反射光

这个点会在w0到wi方向之间发出均匀radiance的光线，把这样一条光线当作蒙特卡洛积分的一次采样，可以看作是在半球上随机取一个点，pdf = 1/2pi，将渲染方程当作fx带入蒙特卡洛积分求解，这样就得到了直接光照的正确算法

对这个方程写出伪代码，对每个方向wi进行一次采样，然后求和得到这个点的radiance

上面的代码只考虑了直接光照

如果求间接光照呢，这个问题可以转换为如果P点测试光线没有打到光源而打到另一个物体O呢？则可以将o看作是一个光源，其发出强度为Li的光线。更进一步，可以看成是从P观察O点，求O点的直接光照。

这样就可以写成一个递归，求出间接光照

但这样还是有问题，一个点对光线进行N次采样，假如一条光线打到另一个问题，又会进行N次采样，如此重复会造成指数爆炸

只有当N=1时，指数爆炸不会有影响。这个N是蒙特卡洛积分的采样次数，1也是合理值，但会有很大的噪声。带入N=1时的伪代码

N=1，这时就叫path tracing 路径追踪，反之就是分布式光线追踪。

如何解决N=1带来的噪声呢，当N=1时一个点只会采样一条光线，那么我们可以对一个像素采样多个点求平均即可

对一个像素点均匀采样N次，进行N次路径追踪得到结果最后求平均

但现在还有一个问题，计算式是递归的，却没有停止条件。这里使用了一种俄罗斯轮盘赌（RR）的方法，简化来说就是，一条光线每次弹射，有P的概率会进行下一次弹射，有1-P的概率会停止返回结果0。

这样做可以保证最终结果的期望是不变的。

加上RR后得到如下伪代码：

现在我们得到了正确的光照算法，但是效率很低。下图对每个像素的采样次数关系着图像质量，SPP越高质量越好

如何提高效率呢。计算直接光照时，如果光源很小，则需要更多的像素点采样才能投射到光源上

但实际上有很多的path都是无用的，被浪费掉了，可以考虑反向直接从光源处采样投射到半球采样点上。

将光源视作一个面光源，其面积为A，则采样概率P为1 / A，其积分域是dA。而渲染函数是对立体角dw上进行积分，需要将积分域转换一下，直接套用下图的结果（考虑是先将光源的面转向采样点，然后等比缩放距离的平方）

变换完成带入渲染方程

至此，可以将光照分为两个部分：直接光照（从光源采样，无RR）和间接光照（从半球采样点采样，RR）

伪代码：

如何知道光源的采样点和半球采样点之间有没有被阻挡呢，可以在两点之间发出一条射线，检测其有没有碰撞到物体，得到最终的路径追踪代码：

path tracing可以达到照片级的光照效果

最后一些还没有解决的问题：

H7补充

有点难。。。做了一天

有些坑：

• 伪代码有个地方有错误，采样光照投影到半球上的点乘应该为 dot(-ws, NN)

shade(p, wo)
sampleLight(inter , pdf_light)
Get x, ws, NN, emit from inter
Shoot a ray from p to x
If the ray is not blocked in the middle
L_dir = emit * eval(wo, ws, N) * dot(ws, N) * dot(-ws, NN) / |x-p|^2 / pdf_light
L_indir = 0.0
Test Russian Roulette with probability RussianRoulette
wi = sample(wo, N)
Trace a ray r(p, wi)
If ray r hit a non -emitting object at q
L_indir = shade(q, wi) * eval(wo, wi, N) * dot(wi, N)
/ pdf(wo, wi, N) / RussianRoulette

Return L_dir + L_indir

• 画面出现黑色条纹可能是判断光线相交的精度问题，我把EPSILON 设置成了0.01
• 直接光照判断阻挡物如果条件写得有问题的话可能会出现天花板是黑色的，猜测可能是因为天花板和灯处于同一平面的缘故
• 微表面和多线程没做，留个坑

// Implementation of Path Tracing
Vector3f Scene::castRay(const Ray &ray, int depth) const
{
    // TO DO Implement Path Tracing Algorithm here

    Intersection intersection = Scene::intersect(ray);

    
    if(intersection.happened) {
        Material *m = intersection.m;
        if (m->hasEmission())
        {
            return m->getEmission();
        }
    
        Vector3f hitPoint = intersection.coords;
        Vector3f N = intersection.normal.normalized(); // normal
        Vector3f L_dir(0);

        Intersection lightInter;
        float lightPdf = 0;
        sampleLight(lightInter, lightPdf);
        Vector3f lightDir = normalize(lightInter.coords - hitPoint);
        
        Ray ray_p2x(hitPoint , lightDir);
        float dis = (lightInter.coords - hitPoint).norm();
        Intersection dLightRay = Scene::intersect(ray_p2x);
        if (dLightRay.happened &&  (dLightRay.coords - lightInter.coords).norm() <= EPSILON)
        {
            Vector3f wo = ray_p2x.direction;
            L_dir = lightInter.emit * m->eval(ray.direction, wo, N) * dotProduct(wo, N)
                * dotProduct(-wo, lightInter.normal)/ std::pow(dis, 2) / lightPdf;
        }
        Vector3f L_indir(0);
        float curP = get_random_float();
        if (curP > RussianRoulette)
        {
            return L_indir + L_dir;
        }

        Vector3f wo = m->sample(ray.direction, N).normalized();
        Ray ray_obj(hitPoint, wo);
        Intersection objRayInter = intersect(ray_obj);
        if (objRayInter.happened && !objRayInter.m->hasEmission())
        {
            L_indir = castRay(ray_obj, depth+1) * m->eval(ray.direction, wo, N) * dotProduct(wo, N)
                / m->pdf(ray.direction, wo, N) / RussianRoulette;
        }

        return L_indir + L_dir;
    }
    else
        return {};
}

spp为50，23min

第十七章 Materials and Appearances

Material == BRDF

对于一个漫反射表面，假设其入射光在各个立体角上的radiance是相等，同时因为是漫反射，所以出射光在各个立体角上的radiance也是相等的，根据能量守恒定律，在表面不吸收能量的情况下，两个radiance相等，所以Lo和Li应该相等。

BRDF和Li都是常数，改写后BRDF fr = 1/pi ，所以BRDF是一个[0, 1/pi]的值，再额外定义一个反射率-albedo参数，以此控制表面的漫反射行为。

反射

完美反射，入射角等于出射角，根据平行四边形法则可以求出出射角

折射

对于折射，有如下性质

介质折射率的不同定义了光的折射行为

求这个出射角的余弦值，可以得出一个带根号的式子，唯有 其值大于0才有物理意义，即

当比例大于1时，代表入射的介质密度大于出射的介质，这时就不会发生折射，而是产生全反射现象（光纤的原理）。

对于折射应该是BTDF

BSDF = BRDF + BTDF

Fresnel Reflection / Term 菲涅尔项

可以发现，随着观察角度的不同，书的反射也不同，和桌面法线呈90度时，反射十分明显。

反射比定义了光线在表面折射与反射的比例

图中有3条线是因为 光具有波动性，蓝线和绿线是光极化的结果，一般考虑两值的均值

正常情况下，反射比Reff十分难计算，可以将其视为一条值域在0-1之间的曲线，通过Schlick表达式去拟合。R0称为基准反射比

Microfacet Material 微表面材质

对于一个表面，我们认为微观下是凹凸不平的，称为微表面。宏观上物体表面是平滑和粗糙的，如上图macrosurface指向的虚线；微观上物体表面是凹凸并且光滑的，如上图microsurface指向的部分。

每个微表面像镜子一样光滑并且具有其自己的法线。

宏观上微表面是材质，微观上微表面是几何

微表面BRDF其实描述的是微表面法线的贡献程度，如上图所示，微表面法线的不同影响了材质整体的效果。

• D(h)用到了半程向量，这里D其实描述的是 微表面法线方向为h方向的分布贡献值（即法线为h方向，反射方向就为wo方向）
• F为菲涅尔项，反应了表面的反射比
• G为几何阴影遮蔽项shadowing-masking：对于微表面，当光线从接近平行表面的方向摄入时，可能会出现部分微表面被另一部分微表面遮挡从而失效，所以对于一个圆球，它的边缘看起来应该比其他部分更暗一点，这个G就是描述这一效果的系数。

Isotropic/Anisotropic Materials(BRDFs) 各向同性/各向异性材质

基础表面的方向性不同，各向同性其法线从各个方向上看均是一样的，各向异性则不然

从BRDF的概念来说就是**各向异性材质的BRDF会随着方位角Φ变化而变化，**所以各向同性材质的BRDF其实是3维的，与方位角无关，与两个方位角之差有关

BRDF的性质

• 非负性

• 线性

• 可逆性

• 能量守恒

• 各向同性BRDF是三维的，两个方位角的正负没有关系

BRDF测量

对于实际的材质，物理上测出的菲涅尔项可能与理论曲线拟合合合的有很大出入

使用下图的测量方式，考虑所有的入射角和出射角，以及方位角，测量出射radiance，从而计算菲涅尔项

但4维枚举计算太复杂，可以这样加速：对于各向同性材质，四维变成三维；利用brdf可逆性，减少一半的方位角测量；使用采样计算，进行猜测。

第十八章 Advance Topics In Rendering

Advanced Light Transport

无偏差的蒙特卡洛估计：期望值总是正确的，与采样数无关

其他情况称为有偏差蒙特卡洛估计：一种特殊的情况，采样数无穷大时期望变为正确值，这种情况称为一致的。

Bidirectional Path Tracing（BDPT）双向路径追踪

从摄像机和光源处都进行路径追踪，并连接追踪路径的末尾点。无偏

传统路径追踪对于非均匀的高能量的区（光）需要要多次数采样才能有好的效果，但双向追踪也会在光源处投射追踪，增加了camera处path采样到光源的概率和效果。

但BDPT实现十分复杂，速度也非常慢。

Metropolis Light Transport（MLT）

在采样时使用了叫Markov Chain马尔可夫链 的概率性工具，可以提供相似的采样结果，生成更好的pdf。无偏的。

非常适用于复杂的光照场景，因为只要找到一条光照路径，剩下的路径可以通过马尔科夫链进行轻微扰动生成更多正确的光路。

缺点：

• 很难估算path计算收敛的时间
• 无法保证每个像素的收敛时间是一致的
• 得到的图像看起来很脏
• 无法保证速率也就无法应用于动画

Photon Mapping 光子映射

是一种有偏估计的方法，分为两个阶段。

非常适合于Specular-Diffuse-Specular（SDS）路径以及caustic现象

第一步：光子追踪，从光源处投射出一定数量的光子，在场景中进行相应的折射和反射，打到diffuse面的时候停止，并在diffuse面上记录这些光子的位置信息。

第二步：光子收集，从摄像机投射subpath，进行折射和反射，当打到diffuse面上时停止，并计算碰撞点周围一定面积的光子密度，密度越高，就越亮。

对于每个着色点，找到最邻近的N个光子，求出它们所占的表面积，再算出密度。

为什么是有偏的，是因为密度计算是估算的，，也就是说光子数量越多，dN/dA 和 ΔN / ΔA 的值越接近，结果也就越准确越清晰，所以光子映射是有偏的，但是也是一致的。

有偏会带来模糊，但保证一致性可以实现在样本足够多时结果清晰正确。

为什么不在diffuse面上固定一个面积A求密度呢，因为固定面积后随着光子数量的增加，ΔN会变大，但是ΔA不会变，最后与dN/dA 还是不等，这样做不仅是有偏的，同时是不一致的。

Vertex Connection and Merging （VCM）

同时应用了BDPT和光子映射的思想，对X1 X2 应用BDPT，对于较近的光路 X2*和X2 应用光子映射。

Instant Radiosity

将被光照亮的表面当作光源（虚拟）

从光源投射subpath，并在subpath 的末端生成虚拟点光源VPL，使用这些VPL去渲染场景。

效果一般很好，并且速度很快

缺点：

• 在接缝处会产生不正常光源，这是因为从光源采样是要除以光源和diffuse采样点的距离，距离很近时，除以一个极小的数会导致问题。
• 无法处理glossy材质

Advanced Appearance Modeling 表面模型

Non-Surface Models

Fog

当光线进入散射介质时，可能会出现被散射或者吸收的情况

类似于BRDF定义的diffuse的反射，相位方程phase function定义了在散射介质中散射的分布

雾的渲染过程与光追类似，见上图。

Hair

一般将头发视为一个圆柱体，kajiya-kay模型认为光打到毛发上会发生散射和反射，但实际效果更类似于布林冯模型

Marschner模型认为光打到毛发上，会产生散射R，两次穿透TT，穿透散射再穿透TRT三种类型。

认为毛发更类似于玻璃柱，光会穿透。

效果很好，但是因为会处理成百上千的头发，光线在期间散射，计算量非常大。

Fur

对于动物毛发则有些区别，动物毛发中的髓质比人毛发大得多，会产生额外的散射

闫老师提出的Double Cylinder Model额外考虑了动物毛发髓质的散射。

Granular Material 颗粒材质

Surface Models

Subsurface Scattering 次表面散射

BRDF假设光线射入和射出是在同一个点上，但是实际上有些材质并非如此。

BSSRDF在BRDF在方向上积分的基础上进行了面积的积分，考虑了出射点的位置信息。

对于一个半透明材质，光线射入后的结果就像是在表面的里外有两个对应的光源在照射一样，使用这种模型去模拟次表面散射

Cloth

cloth一般认为是：纤维fiber组成股ply，股组成线yarn

有几种渲染方式：

将织物看作表面使用BRDF渲染将织物表面作为一定体积的雾进行渲染

还原每一根fiber，ply进行渲染

Detial

真实世界的物体表面会有许多划痕。

修改微表面模型法线分布来模拟

对于十分微小的表面则要应用波动光学

Procedural Appearance 程序化生成表面

类似于3D纹理：有一个noise函数，给定3d空间的一个坐标值xyz，noise返回这个纹理的相应值

第十九章 Cameras，Lenses and Light Fields 相机，透镜和光场

Camera

相机成像分为：小孔成像、透镜成像

相机传感器记录的是光的irradiance，没有区分方向性，所以直接用传感器记录信息会使其接收各个方向的光，得到的东西是模糊的。

使用小孔成像生成的图像是没有景深效果的，path tracing使用的就是小孔成像的模型。

FOV

fov视场角与焦距和传感器大小有关，公式如下。

一般认为多少等效焦距是指在35mm大小传感器情况下等效的焦距大小，所以厚度不足10mm的手机也可以有28mm的等效焦距，因为手机的传感器更小。

传感器大小与画面大小也有关系

Exposure

H = T x E 曝光=时间 x irradiance

曝光主要是由三个要素控制：光圈大小，快门速度，ISO感光度

各要素效果预览：

ISO：是一种线性的增益，相当于直接将光信息乘以一定的倍数，是线性的，200iso效果一定是100iso 的两倍。同时由于是直接相乘增大，所以图片的噪声也会被增大，高iso会造成画面有很多噪点。

光圈：写作FN，N为光圈数，是光圈直径分之一。

快门：指的是快门打开关闭的时间，快门速度越快时间越短，慢的快门会导致画面因为物体高速移动而模糊，称为运动模糊motion blur，或者是高速物体扭曲（由于快门是缩放打开，画面其实类似于扫描成像，可能一部分画面的时间与另一部分的时间不等，而这期间物体运动了，从而产生扭曲）

以下快门和光圈的组合可以达到相近的曝光度：

同时 运动模糊（慢快门） 与 景深（小光圈）不可兼得。

Lens

在理论模型中我们假设透镜是没有厚度的薄透镜。

• 所有平行光经过透镜都会汇聚到焦点上
• 所有经过焦点的光穿过透镜都会变成平行光
• 透镜焦距可变（工程中通过几块透镜的组合可以实现）

为什么会产生景深呢？一个点的光线透过透镜汇聚到焦点，再继续延申打到成像平面上，这个点的光线会变成一个球即所谓的Cirlce of confusion 弥散圈，许多点的coc相互叠加从而导致模糊。

从下图公式可以看出coc与光圈大小有关。

所以光圈N/F的F数实际定义其实是 焦距 / 光圈直径。

Ray Tracing Ideal Thin Lenses

path tracing默认使用的模型是针孔摄像机模型，可以改为透镜模型。

• 决定传感器大小，焦距和光圈大小
• 决定物距
• 推算出像距zi

• 对于sensor上的每一个像素
• 得到光路信息，记录radiance在sensor上

Depth of Field景深

DOF指的就是像素空间中的光线穿过透镜形成的coc大小十分小，认为此时画面是清晰的。

可以看见随着光圈缩小，DOF也在增加，所以小孔时认为无限远，即没有景深

Light Fields/Lumigraph

Plenoptic Function

现实世界可以用一个七维的全光函数进行描述，θ和Φ描述一条光线的位置方向，λ描述光的波长即颜色，t描述时间，VxVyVz描述空间中的任意位置。

光场则可以认为是记录了任何一个点所接收的任何一个方向的irradiance。

即可以用一个四维函数去描述，两维描述物体表面上一个点的位置，两维描述这个点接收的光线的方向。

同时，两点确定一条直线，所以这个函数又可以通过两个平行平面改写为两个位置的函数

从uv平面看向st平面，得到的是这个物体从各个方向看的相应的图像，而从st看向uv，则得到的是这个物体在这个方向的irradiance的集合即radiance。

Light Field Camera

光场摄像机使用微透镜替换像素，通过微透镜拆分一个像素的irradiance并存储起来，从而可以实现相片的后期动态聚焦以及镜头的虚拟移动。

通过改变图像一个像素取的透镜中某一条光线的方向，可以模拟摄像机的移动，或者计算某些特定方向的光线，从而动态聚焦。

光场摄像机也有几个问题：

• 低分辨率，因为在相同的胶片上同时存储空间和方向信息。
• 高开销，微透镜十分复杂，同时信息量很大，很难存储。

第二十章 Color and Perception 颜色和感知

可见光在不同波长有不同颜色

光的Spectral Power Distribution （SPD）函数描述了它在不同波长上的强度大小

SPD是线性的

人的视网膜上分布着两种细胞：感知光的强弱的视杆细胞，感知光线颜色的视锥细胞。视锥细胞又分为三种S，M，L去感知不同波长的光照。下图是对应波长的感知强度曲线。

通过在波长上对感知曲线积分得到三个值SML，发送到大脑组合形成颜色。

Metamerism同色异谱

可以通过对感知曲线进行调整，达到不同的曲线积分得到相同的结果从而使SML值颜色完全相同。称为同色异谱。

Color Matching Function颜色匹配函数

由于spd是线性的，所以可以定义三个强度值RGB，去乘以三个标定的颜色红绿蓝，以此混合出不同的颜色。

实验得到下图的RGB颜色匹配函数，给出了任一波长的颜色三个系数强度值为多少。

Color Space色彩空间

与RGB同理，通过人为定义的颜色匹配函数，得到一个XYZ系统，其中y曲线分布均匀，所以将其视为颜色的亮度。

通过将Y值固定，对X和Z进行归一化，可以在二维空间下显示xyz系统所能表示的所有颜色。

因为xyz颜色系统是加色系统，所以白色是最不纯的颜色，而边缘的颜色是纯色。

色彩空间就是色彩系统所能表示的所有颜色，可以看见sRGB所能表示的色彩其实很局限。

CMYK是减色系统，随着混合，颜色会变为黑色。要额外一个黑色是为了考虑印刷的成本。

第二十章 Animation

Keyframe Animation

关键帧动画，原画师创造关键帧，动画师补出其中的过渡

关键帧动画插值不一定是线性的

Physical Simulation

通过力学方程去模拟物体真实的运动

Mass Spring System质点弹簧系统

质点弹簧系统的最小单元是两个有质量的点与一条弹簧，弹簧放松情况下长度为l，通过胡克定律可以得到两个点之间的力：

但是这个系统还没有考虑摩擦力，否则根据能量守恒定律在施加一个力后它会一直震荡下去。

在物理模拟系统中如下表示速度和加速度

增加一个摩擦力称为dumping，红框内是求出ab之间的相对速度并投影到ab之间的方向上去。

dump公式如下：

考虑如下一个布料结构的质点弹簧系统：

其无法抵挡斜切力以及弯折力

可以改进为如下结构：

蓝线为强力，红线为弱力

Particle Systems

可以通过计算一个个粒子之间的相互作用力（摩擦力，引力，排斥力。。。）去得到一个粒子集合的运动状态，再渲染出来。

一个简单的粒子系统分为以下步骤：

• 有需要就创建新的粒子
• 计算每个粒子之间的力
• 更新粒子的位置和速度
• 移除失效的粒子
• 渲染

但是粒子系统中的粒子一词并不是单纯的指粒子，也可以是一个类似个体组成的群体中的一个，如下图通过粒子系统定义一只鸟在鸟群中的运动规律从而模拟鸟群。

Forward Kinematics 正向运动学

主要是计算骨骼间的位置关系

十分容易计算，先确定各个骨骼的位置，再计算终点p位置

但是fk结果可能与物理不一致，同时不太适合于美工使用。

Inverse Kinematics 逆向运动学

IK先确定骨骼末端的位置p再计算确定中间骨骼的姿态

但有可能出现多解和无解的情况。

Rigging 绑定

Rigging是一组可以自由控制角色姿态骨骼的控制点。

Blend Shapes（BS）与Rig类似，不过控制的是模型的表面。

The Production Pipeline

第二十一章 Animation（cont）

Single Particle Simulation

计算单个粒子的运动，我们假设有一个速度场函数，可以通过粒子的位置与时间获得粒子的速度。

Euler’s Method

Explicit Euler’s Method 显式欧拉方法

最简单的方式是使用显式（前向）欧拉方法，始终使用上一个Δt的位置和加速度去估算粒子当前的位置和速度。

使用非常简单和广泛，但是会有问题：有误差，不稳定。

误差可以通过缩小Δt来减少

但是不稳定性很难解决，以此为基础的系统最终结果都会与正确结果有分歧。

显式欧拉方法不稳定，可以考虑使用后续几种方法提高稳定性。

Midpoint Method

中点法，先计算出下一个点a，求出当前点与a的中点b，再使用b的速度a的位置去计算下一个点c。

相比于显式欧拉方法，多了一个二次的项。显式欧拉可以认为是局部线性的方法，而中点法则模拟了抛物线，更稳定。

Adaptive Step Size 自适应步长

自适应步长方法每个点之间计算所使用的Δt是不同的。过程如下：

• 计算下一个点Xt，步长为t
• 使用t/2 计算另一个点Xt/2
• 计算两个点误差，如果误差大于阈值则缩减步长T并重新计算

Implicit methods 隐式欧拉方法

显式欧拉方法使用上一步的速度和位置去计算当前速度，而隐式方法使用下一步的位置速度计算当前的速度，更加准确，但是很难解。

Runge-Kutta Families

如何定义函数的稳定性呢，我们定义关于步长h的表达式：O(h^n)，n越高，则越稳定，可以理解为n=1，当步长h缩小一半，误差也缩小一半，当n=2，则误差缩小到四分之一。

隐式欧拉具有一阶的稳定性。

可以使用RK4方法去解常微分方程，这个方法是四阶的。

Position-Based / Verlet Integration

非物理的模拟

Rigid Body Simulation

刚体模拟，实际上是对单粒子模拟的补充。刚体不会发生形变，所以刚体内的所有粒子都具有相同的运动状态。

Fluid Simulation

流体模拟，基于position-based的方法：

• 假设水都是有一个个的刚体小球组成
• 水不可被压缩（所以任意位置密度应该一致）
• 当有地方密度不一致时，移动小球去修正位置
• 密度的梯度与小球位置有关
• 进行梯度下降

对于大规模的物体，我们有两种模型方法：

质点法（拉格朗日方法）：考虑模拟规模群体中的每一个个体，并考虑其间的作用力

网格法（欧拉方法）：将空间划分为一个一个的小格，模拟计算每一个网格内的运动

将两者结合的方法称为Material Point Method(MPM)材质点法。

在网格中算出运动信息并记录，然后写回到粒子中。

结语

H8的环境一直配不好，就先留个坑不做了，也不难。

零零散散学了一个多月，终于学完了，从Games101中真的学到了很多，感谢闫老师🙏

有机会二刷外加补坑，还待202和104去学习。