tinyrender games101笔记

curry dog

已于 2022-02-05 14:29:25 修改

阅读量311

点赞数

分类专栏：视觉 3D视觉文章标签： 3d 图形学

于 2021-10-02 15:46:12 首次发布

本文链接：https://blog.csdn.net/weixin_39849839/article/details/120549398

版权

视觉同时被 2 个专栏收录

14 篇文章 0 订阅

订阅专栏

3D视觉

5 篇文章 0 订阅

订阅专栏

draw line

给定(x0,y0) (x1,y1)两个点，画出一条直线。首先通过交换x,y轴使得变化较快的那个轴变为x轴，再交换x0,x1使得x1永远大于x0。令dx=x1-x0;dy=y1-y0;
从(x0,y0)开始，x轴每增加一个像素，因为y轴变化较慢所以增加的值dy/dx小于1，把增加的值累积起来用变量error记录，每当error>0.5就令y的值增加1，同时error的值减去1. 但是这样做的缺点是需要使用0.5和dy/dx两个浮点数，所以直接两者同时乘以2dx。

void line(int x0, int y0,int x1, int y1, TGAImage &image,TGAColor color){
    bool steep=false;
    if(abs(y1-y0)>abs(x1-x0)){
        steep=true;
        std::swap(x0,y0);
        std::swap(x1,y1);
    }
    if(x1<x0){
        std::swap(x0,x1);
        std::swap(y0,y1);
    }
    int dx=x1-x0;
    int dy=y1-y0;

    int derror=abs(dy)*2;
    int error=0;
    for(int x=x0, y=y0;x<=x1;x++){
        error+=derror;
        if(error>dx){
            y+=dy>0?1:-1;
            error-=2*dx;
        }
        if(steep){
            image.set(y,x, color);
        }else{
            image.set(x,y,color);
        }
    }
}

模型导入
.obj文件的每一行要么是点，要么表示一个面

v 0.296 -0.070 0.953 v -0.319 -0.065 0.946

f 1193/1240/1193 1180/1227/1180 1179/1226/1179

第一行表示一个点，使用归一化坐标，第二行表示由第1193,1180和1179个点组成的面。这些序号都是从1开始。

triangle

思路，求出三角形的bbox后，需要判断其中的每一个点是否在三角形内。这里用到了重心坐标的概念，即如果P在三角形ABC内，则存在三个分别放在ABC点上的质量，使得P成为三角形的重心。
即P点满足
$\overrightarrow{P A}+v \overrightarrow{P B}+(1-u-v) \overrightarrow{PC}=0$
现在给定P点需要求出u,v,(1-u-v),如果其中一个值小于0则说明P不在三角形内。
整理得
$\left\{\begin{array}{l} {\left[\begin{array}{lll} u & v & 1 \end{array}\right]} & {\left[\begin{array}{c} \overrightarrow{CA}_{x} \\ \overrightarrow{CB}_{x} \\ \overrightarrow{P C}_{x} \end{array}\right]=0} \\ {\left[\begin{array}{lll} u & v & 1 \end{array}\right]} & {\left[\begin{array}{c} \overrightarrow{CA}_{y} \\ \overrightarrow{CB}_{y} \\ \overrightarrow{PC}_{y} \end{array}\right]=0} \end{array}\right.$

可以看出(u,v,1)和两个向量(ABx,ACx,PAx) and (ABy,ACy,PAy)正交，所以可以先求出两个向量的叉乘，然后除以z轴的值，得到（u,v,1）,注意这里叉乘结果的z值可能为0，这种情况可以直接归为P不在三角形内。

Vec3f barycentric(Vec2i* pts,Vec2i P){
    Vec3i v1(pts[1].x-pts[0].x,pts[2].x-pts[0].x,pts[0].x-P.x);
    Vec3i v2(pts[1].y-pts[0].y,pts[2].y-pts[0].y,pts[0].y-P.y);
    auto u=v1.cross(v2);
    if(u.z==0){
        return Vec3f(-1,1,1);
    }else{
        Vec3f v=u;
        return Vec3f(v.x/v.z,v.y/v.z,1.0f-(v.x+v.y)/v.z);
    }
}

void triangle(Vec2i* pts, TGAImage& image,TGAColor color){
    Vec2i boxmax(0,0);
    Vec2i boxmin(image.get_width()-1,image.get_height()-1);
    for(int i=0;i<3;i++){
        for(int j=0;j<2;j++){
            if(pts[i].raw[j]>boxmax.raw[j]){
                boxmax.raw[j]=pts[i].raw[j];
            }
            if(pts[i].raw[j]<boxmin.raw[j]){
                boxmin.raw[j]=pts[i].raw[j];
            }
        }
    }
    for(int x=boxmin.x;x<=boxmax.x;x++){
        for(int y=boxmin.y;y<=boxmax.y;y++){
            Vec2i  P(x,y);
            auto b= barycentric(pts,P);
            if(b.x>=0 && b.y>=0&& b.z>=0){
                image.set(x,y,color);
            }
        }
    }
}

坑点：注意u,v和x,y,z的对应关系。

Zbuffer

上面的实现的triangle没有考虑深度的问题，即使三角形A的z值比三角形B大，但是渲染时还是可能被B给覆盖，解决方法是利用一个zbuffer[image_width*imag_height],来缓存渲染过程中，三角形上每个点的z值，这样渲染结束后，zbuffer中每一个z值都是所有对应点中最大的。三角形上某点z值的计算过程是以该点的重心坐标为权重，乘以三角形三个顶点的z坐标。

void triangle(Vec3i* pts, float* z_buffer, TGAImage& image,TGAColor color){
    Vec2i boxmax(0,0);
    Vec2i boxmin(image.get_width()-1,image.get_height()-1);
    for(int i=0;i<3;i++){
        for(int j=0;j<2;j++){
            if(pts[i].raw[j]>boxmax.raw[j]){
                boxmax.raw[j]=pts[i].raw[j];
            }
            if(pts[i].raw[j]<boxmin.raw[j]){
                boxmin.raw[j]=pts[i].raw[j];
            }
        }
    }
    for(int x=boxmin.x;x<=boxmax.x;x++){
        for(int y=boxmin.y;y<=boxmax.y;y++){
            Vec2i  P(x,y);
            auto b= barycentric(pts,P);
            if(b.x>=0 && b.y>=0&& b.z>=0){
                float z=0;
                for(int i=0;i<3;i++){
                    z+=b.raw[i]*pts[i].z;
                }
                int index=P.y*image.get_width()+P.x;
                if(z_buffer[index]<z){
                    z_buffer[index]=z;
                    image.set(x,y,color);
                }
            }
        }
    }
}

projection texture

v 0.296 -0.070 0.953 v -0.319 -0.065 0.946

f 1193/1240/1193 1180/1227/1180 1179/1226/1179

第二行除了表示由第1193,1180和1179个点组成的面之外，还表示三个顶点分别对应第1240,1227,1226个uv坐标点(vt)和第1193,1189,1179个法向量点(vn)。
在这里插入图片描述
和视觉里的相机模型不同，这里原点在像平面上，假设相机坐标为(0,0,c), 那么物体上某点在像平面上的投影坐标可以用下式计算：
$\left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & -1 / c & 1 \end{array}\right]\left[\begin{array}{c} x \\ y \\ z \\ 1 \end{array}\right]=\left[\begin{array}{c} x \\ y \\ z \\ 1-z / c \end{array}\right]$
计算投影矩阵的函数一般以1/c为参数，因为可以令1/c为0，这样就等于正射投影。
相机移动需要提供三个参数，一个是相机的位置eye，一个是目标点center,即新的坐标系原点，以及up向量。up向量在像平面的投影就是新的y轴。由center到eye的向量就是新的z轴方向。利用新z轴和up的叉乘可以得到新的x轴方向，最后用z和x的叉乘可以得到新y轴方向。
新坐标系下的坐标 $x^{\prime},y^{\prime},z^{\prime}$ 可以用下式计算
$\left[\begin{array}{l} x \\ y \\ z \end{array}\right]=\left[\begin{array}{l} O_{x}^{\prime} \\ O_{y}^{\prime} \\ O_{z}^{\prime} \end{array}\right]+M\left[\begin{array}{l} x^{\prime} \\ y^{\prime} \\ z^{\prime} \end{array}\right] \quad \Rightarrow \quad\left[\begin{array}{l} x^{\prime} \\ y^{\prime} \\ z^{\prime} \end{array}\right]=M^{-1}\left(\left[\begin{array}{l} x \\ y \\ z \end{array}\right]-\left[\begin{array}{c} O_{x}^{\prime} \\ O_{y}^{\prime} \\ O_{z}^{\prime} \end{array}\right]\right)$
其中M矩阵是新坐标系三个轴对应的向量组成标准正交基。 $O^{\prime}$ 是给定的center点的坐标。
在视角变换，投影变换之前，模型的坐标值在定义域 $[- 1, 1] * [- 1, 1] * [- 1, 1]$ 中，需要利用viewport矩阵将转换后的定义域变为 $[x, x + w] * [y, y + h] * [0, d]$ ，这里d一般设置为255，方便将zbuffer映射为一个图片:
$\left[\begin{array}{cccc} \frac{w}{2} & 0 & 0 & x+\frac{w}{2} \\ 0 & \frac{h}{2} & 0 & y+\frac{h}{2} \\ 0 & 0 & \frac{d}{2} & \frac{d}{2} \\ 0 & 0 & 0 & 1 \end{array}\right]$

直接对法向量进行变换的结果不一定正确(当变换矩阵为正交矩阵是可以直接对法向量变换)。假设变换前法向量为[A,B,C,0]（在齐次坐标系中，向量的最后一维为0）, 它对应一个平面，对平面上任何点（x,y,z）都有
$\left[\begin{array}{llll} A & B & C & 0 \end{array}\right] \times\left[\begin{array}{l} x \\ y \\ z \\ 1 \end{array}\right]=0$
在上式中间插入单位矩阵 $M^{-1}M$ 得到
$\left(\left[\begin{array}{llll} A & B & C & 0 \end{array}\right] \times M^{-1}\right) \times\left(M \times\left[\begin{array}{l} x \\ y \\ z \\ 1 \end{array}\right]\right)=0$
整理得
$\left(\left(M^{\top}\right)^{-1} \times\left[\begin{array}{l} A \\ B \\ C \\ 0 \end{array}\right]\right)^{\top} \times\left(M \times\left[\begin{array}{l} x \\ y \\ z \\ 1 \end{array}\right]\right)=0$
可以发现右边正好是变换后的（x,y,z）点。所以应该对法向量进行的变换是 $\left(M^{\top}\right)^{-1}$

注意：我们真正关心的重心坐标是根据三维空间中的顶点得到的，而不是根据投影到屏幕上的点得到的，但是对三角形进行采样时必须使用屏幕坐标，可以通过修正算法将屏幕上的重心坐标变成三维空间中的重心坐标：
$\begin{gathered} z_{P}=\frac{1}{\frac{a}{z_{A}}+\frac{b}{z_{B}}+\frac{c}{z_{C}}} \\ I_{P}=\left(a \frac{I_{A}}{z_{A}}+b \frac{I_{B}}{z_{B}}+c \frac{I_{C}}{z_{C}}\right) / \frac{1}{z_{P}} \end{gathered}$

其中a,b,c是根据屏幕上的点得到的重心坐标,z是三维空间的深度，I是某个属性值
在这里插入图片描述
这个原理是，假设三维空间点(x,y,z)的投影点为 $P_{x},P_{y},0)$ ,根据投影矩阵可得
$x=(P_{x} C-P_{x} z)/C$ 同理可得y的表示方法，代入三维空间中的平面方程可以发现 $1 / z$ 可以用 $P_{x}$ ， $P_{y}$ 线性表示,所以在投影平面的插值也可以应用到 $1 / z$ 上

这里的光照模型比较简单，由漫反射diffuse和全反射specular组成，前者用法向量n和入射光线l的内积表示，而后者由反射光线r的z值表示,这里的z轴指的就是由相机指向原点的轴，所以要求以上所有的方向向量全部都是经过了变换之后的：
在这里插入图片描述
这里入射光线是远离物体的。

alpha 是当前待绘制颜色和缓冲区颜色的混合系数，alpha越接近1越不透明。
不翻转的话，屏幕的y轴向下，物体坐标系的y轴往物体上方纹理图本身的y轴向下
不管模型坐标系中,x指向左边还是右边，y轴指向上面还是下面，只要把该轴作为up向量，那么它在屏幕中都指向上方。深度测试优先显示z值大即离相机近的点。

实现
着色器分为两个部分，vertex和fragment，前者负责计算三角形顶点的坐标，以及获取对应的uv坐标，每个面对应的所有uv坐标缓存在shader中，方便fragment使用。 fragment负责根据当前内部点的重心坐标获取颜色信息。两者每次都只负责处理一个点。

struct MyShader : public IShader{

    explicit MyShader(Model*m){
        M=ViewPort*Projection*ModelView;
        invertTransposeM=M.invert_transpose();
        model=m;
    }

    Vec4f vertex(int fid,int nthVertex) override{
        Vec4f vert=proj<4>(model->vert(fid,nthVertex),1.0f);
        originPts.set_col(nthVertex,vert);
        vert=ViewPort*Projection*ModelView*vert;
        varying_uv.set_col(nthVertex,model->uv(fid,nthVertex));
        return vert;
    }

    bool fragment(Vec3f bar,TGAColor& color) override{
        Vec2f uv=varying_uv*bar;
        Vec3f n= proj<3,float>(invertTransposeM*proj<4,float>(model->normal(uv)),0).normalize();
        Vec3f l=proj<3,float>(M*proj<4,float>(light_dir),0).normalize();
        float diff=std::max(0.f, n*l);


        Vec3f r= (n*(n*l*2.f) - l).normalize().normalize();
        float spec=pow(std::max(0.f,r.z),model->specular(uv));
        auto c=model->diffuse(uv);
        for(int i=0;i<3;i++){
            color.raw[i]=std::min<float>(5+c.raw[i]*(diff+spec),255);
        }
        return true;
    }
};

求矩阵的逆，需要分解为求矩阵的伴随矩阵adjugate和求行列式，而求伴随矩阵相当于求每个元素的余子式cofactor，而求余子式又需要求出低阶的行列式。所以整个过程可以看成是cofactor函数和det函数的递归调用，这里用模板的特化来终止递归。

template<size_t DimCols,size_t DimRows,typename T> class mat;
template<size_t DIM,typename T> struct dt {
    static T det(const mat<DIM,DIM,T>& src) {
        T ret=0;
        for (size_t i=DIM; i--; ret += src[0][i]*src.cofactor(0,i));
        return ret;
    }
};

template<typename T> struct dt<1,T> {
    static T det(const mat<1,1,T>& src) {
        return src[0][0];
    }
};
mat<DimRow-1,DimCol-1,T> get_minor(int row, int col) const{
        mat<DimRow-1,DimCol-1,T> ret;
        for(int i=0;i<DimRow-1;i++){
            for(int j=0;j<DimCol-1;j++){
                ret[i][j]= rows_[i>=row?i+1:i][j>=col?j+1:j];
            }
        }
        return ret;
    }

    T cofactor(int row,int col) const{
        return (get_minor(row,col).det())*((row+col)%2?-1:1);
    }
    T det() const{
        return dt<DimCol,T>::det(*this);
    }

    mat<DimRow,DimCol,T> adjugate(){
        mat<DimRow,DimCol,T> ret;
        for(int i=0;i<DimRow;i++){
            for(int j=0;j<DimCol;j++){
                ret[i][j]= cofactor(i,j);
            }
        }
        return ret;
    }

坑点：法向量贴图在读取时要注意变换成[-1,1], 在做深度测试时要对z值进行裁剪z= std::max(0, std::min(255, int(z+.5)));

在这里插入图片描述

tangent space

之前的法向量贴图是建立在对象空间上的，如果想把纹理用到其它模型上就不行了，而建立在切线空间上的法向量贴图可以用在不同的模型上。

有了切线坐标系后，从法向量纹理图中读取的就是在切线坐标系中的坐标值。在切线坐标系中，z的值永远大于0并且在[0.1]之间，其它两轴在[-1,1]之间，用(normal + 1) / 2*255 转换成RGB格式保存到文件中，可以看出由于z大于0所以轴对应的blue通道总是大于255/2,所以切线纹理图像中出现大片蓝色。

三角形上每一个点都有自己的切向量空间，它的x轴(T)就是当前点u变化最快的方向，y轴(B)是纹理坐标v变化最快的方向，z轴（N）可以是插值得到的法向量方向。
以求T为例，假设u是顶点坐标的线性函数Ax+By+Cz+D那么其变化最快的方向就是梯度(A,B,C)，它的解法为：

$\left\{\begin{array}{rl} \overrightarrow{p_{0} p_{1}} & \cdot(A & B & C) & =f_{1}-f_{0} \\ \overrightarrow{p_{0} p_{2}} & \cdot(A & B & C) & =f_{2}-f_{0} \\ \overrightarrow{n} & \cdot (A & B & C) & =0 \end{array}\right.$

实际上这样求出来的三个轴不一定正交，所以可以先求出u的梯度作为T,然后和v的梯度叉乘得到N轴，最后NxT得到B轴.

沿着任意方向变化时的导数，可以看出梯度方向时，变化最快
$D_{u} f(x, y)=f_{x}(x, y) \cos \theta+f_{y}(x, y) \sin \theta$

    bool fragment(Vec3f bar,TGAColor& color) override{
        Vec2f uv=varying_uv*bar;

        Vec3f n=(tri_normals*bar).normalize();

        mat<3,3,float> A;
        A[0]=tri_vertexes.get_col(1)-tri_vertexes.get_col(0);
        A[1]=tri_vertexes.get_col(2)-tri_vertexes.get_col(0);
        A[2]=n;
        A=A.invert();
        Vec3f bu(varying_uv[0][1]-varying_uv[0][0],varying_uv[0][2]-varying_uv[0][0],0);
        Vec3f bv(varying_uv[1][1]-varying_uv[1][0],varying_uv[1][2]-varying_uv[1][0],0);
        auto T=(A*bu).normalize();
        auto B=(A*bv).normalize();

        mat<3,3,float> cord;
        cord.set_col(0,T);
        cord.set_col(1,B);
        cord.set_col(2,n);

        auto read_n=model->normal(uv);

        n=(cord*read_n).normalize();


        float diff=std::max(0.f, n*light_dir);


        Vec3f r= (n*(n*light_dir*2.f) - light_dir).normalize();
        float spec=pow(std::max(0.f,r.z),model->specular(uv));
        auto c=model->diffuse(uv);
        for(int i=0;i<3;i++){
            color.raw[i]=std::min<float>(c.raw[i]*(diff+spec),255);
//            color.raw[i]=std::min<float>((diff+spec)*255,255);
        }
        return true;
    }
};

上面所有的向量和顶点都统一只经过了projection和modelview。
在这里插入图片描述

shadow

这里引入了视锥的概念，并且相机点为原点。原因是我们想在三维空间内找到一个感兴趣的空间，这个空间内的点可以尽可能的填充屏幕。

视锥有近平面n和远平面f，近平面相当于像平面，所以近平面上所有点的齐次坐标在经过透视投影后不变，远平面的中心点(0,0,f,1)经过透视投影后不变。这两个性质保证这是一个过原点的视锥。
所以可以得到投影矩阵
$\begin{array}{cccc} n & 0 & 0 & 0 \\ 0 & n & 0 & 0 \\ 0 & 0 & n+f & -n f \\ 0 & 0 & 1 & 0 \end{array}$

然后视锥被压缩成一个长方体[l,r][b,t][f,n]，然后可以用正交投影装换成正方体[-1,1][-1,1][-1,1]

$M_{\text {ortho }}=\left[\begin{array}{cccc} \frac{2}{r-l} & 0 & 0 & 0 \\ 0 & \frac{2}{t-b} & 0 & 0 \\ 0 & 0 & \frac{2}{n-f} & 0 \\ 0 & 0 & 0 & 1 \end{array}\right]\left[\begin{array}{cccc} 1 & 0 & 0 & -\frac{r+l}{2} \\ 0 & 1 & 0 & -\frac{t+b}{2} \\ 0 & 0 & 1 & -\frac{n+f}{2} \\ 0 & 0 & 0 & 1 \end{array}\right]$

模拟阴影的思路是将相机移动到光源方向(C1坐标系)，渲染然后计算shadow_buffer, 在正常渲染时（C2坐标系），把当前点变换到C1坐标系，然后比较它的z值和shadow_buffer中的z值，如果小于说明该点在光源方向被遮挡。已知屏幕上的点(带z值)，因为齐次坐标的尺度不变性，可以通过乘以变换矩阵的逆矩阵得到空间点的准确坐标。

精度问题
经过变换后的空间称为NDC(除以第四维w),变换前为viewspace,如果NDC中近平面为0，远平面为1，用浮点数来保存NDC中的深度值时，1附近保存了大量的远点，而浮点数靠近1时的精度反而不如靠近0的精度，所以一般令近平面为1，远平面为0.
另一方面，在像素平面上，同一块面积，能包含的远点比近点要多得多，所以像素点移动一格可能对应两个隔得非常远的远点。因此远点的贴图一般会用高斯卷积核进行模糊。在把C2中的坐标转到光源坐标并利用深度测试来判断是否遮挡时，可以把当前点的z值加上一个常数，以减少z-fighting的情况。
透视修正
前面只修复了z值，并没有修复重心坐标。前面得到z的导数是线性插值，又
因为其它变换前的属性和z值成线性关系：
$a m^{}+b z^{}=c$
所以属性值除以z也是线性插值
$\frac{a m^{}}{c z^{}}+\frac{b}{c}=\frac{1}{z^{}}$
因此已知三个顶点投影前的z值，以及属性值，可以求得中间某点经过修正后的属性值：
$v_{n}=Z_{n}\left(\frac{1-u-v}{Z_{1}} u v_{1}+\frac{u}{Z_{2}} u v_{2}+\frac{v}{Z_{3}} u v_{3}\right)$
其中 $Z_{n}=\frac{1}{\frac{1-u-v}{Z_{1}}+\frac{u}{Z_{2}}+\frac{v}{Z_{3}}}$

struct MyShader : public IShader{
    Mat4f otherM;
    Mat4f shadowM;
    float * shadowBuffer;
    explicit MyShader(Model*m,float * buffer,Mat4f otherTransM):IShader(m),otherM(otherTransM)
    ,shadowBuffer(buffer){
        M=Projection*ModelView;
        invertTransposeM=M.invert_transpose();
        shadowM=otherM*(M.invert());
    }

    Vec4f vertex(int fid,int nthVertex) override{
        Vec4f vert=proj<4>(model->vert(fid,nthVertex),1.0f);
        originPts.set_col(nthVertex,vert);
        vert=M*vert;
        varying_uv.set_col(nthVertex,model->uv(fid,nthVertex));
        tri_vertexes.set_col(nthVertex, vert/vert[3]);
        tri_normals.set_col(nthVertex, proj<3,float>(invertTransposeM* proj<4,float>(model->normal(fid,nthVertex),0.0f)));
        return ViewPort*vert;
    }

    bool fragment(Vec3f bar,TGAColor& color) override{
        Vec2f uv=varying_uv*bar;

        Vec3f n=(tri_normals*bar).normalize();


        mat<3,3,float> A;
        A[0]= proj<3,float>(tri_vertexes.get_col(1)-tri_vertexes.get_col(0));
        A[1]=proj<3,float>(tri_vertexes.get_col(2)-tri_vertexes.get_col(0));
        A[2]=n;
        A=A.invert();
        Vec3f bu(varying_uv[0][1]-varying_uv[0][0],varying_uv[0][2]-varying_uv[0][0],0);
        Vec3f bv(varying_uv[1][1]-varying_uv[1][0],varying_uv[1][2]-varying_uv[1][0],0);
        auto T=(A*bu).normalize();
        auto B=(A*bv).normalize();

        mat<3,3,float> cord;
        cord.set_col(0,T);
        cord.set_col(1,B);
        cord.set_col(2,n);

        auto read_n=model->normal(uv);

        n=(cord*read_n).normalize();


        float diff=std::max(0.f, n*light_dir);


        Vec3f r= (n*(n*light_dir*2.f) - light_dir).normalize();
        float spec=pow(std::max(0.f,r.z),model->specular(uv)+0.0001);
        auto c=model->diffuse(uv);

        auto p=tri_vertexes*bar;


        auto s_p=shadowM* p;
        s_p= s_p/s_p[3];
        float shadow=1.0f;
        int idx = int(s_p[0]) + int(s_p[1])*width;
        if(idx>=0&&idx<width*height&&shadowBuffer[idx]>s_p[2]+0.01){
            shadow-=0.8;
        }


        for(int i=0;i<3;i++){
            color.raw[i]=std::min<float>(1.3*c.raw[i]*shadow*(diff+spec),255);

        }
        return true;
    }
};

在这里插入图片描述

环境光
绘制环境光的方法是在一个半球面上均匀地采样点，对于每个光源点，用一张与纹理图大小相等的图来保存纹理图上点的是否可见的信息。最后把所有的图求平均。这里要注意的是，在球面上均匀采样是指，单位面积采样点数量相同，但是球面上单位面元为 $P_{A}=p(v) d A=\frac{1}{4 \pi} d A$ 和theta角相关，所以如果直接对theta和phi角采样会得到错误采样。积分坐标变换后微元要乘以一个积分后坐标对积分前坐标的雅克比式，用来表示单位面积。
在球面上采样的正确方法
$\left\{\begin{array}{l}\theta=\arccos \left(1-2 \xi_{x}\right) \\ \phi=2 \pi \xi_{y}\end{array}\right.$
其中x,y都是0,1上均匀分布的变量

Geometry

Implicit
给出点之间的关系 $f (x, y, z) = 0$ ，并不给出实际的点。隐式的好处是方便判断一个点是否在表面内。

algebraic surface
就用不同的曲面方程来表示
constructive solid geometry
通过一系列基本几何以及他们之间的布尔运算来定义几何。
signed distance function
空间中任何一个点的值定义为它到几何表面的距离，如果在表面外就是正数，如果在表面内就是负数。那既然已经可以用正负号来决定表面了为什么还需要距离呢，答案是方便进行插值，比如给定形状A和运动后的形状B，可以通过插值得到中间状态的形状。
TSDF
用于RGBD图的三维重建

把场景分成n个体素，对于体素的每一个顶点x，按照当前相机位姿计算深度，并和它在图片中的投影点p的深度相减，就得到了x和表面的距离sdf(x)，截断后变成了TSDF.
对多张图片计算tsdf并融合。对于每一个体素，知道了8个顶点的tsdf值后，找到值为0的等值面作为最后输出的mesh。
level sets
$\psi(\mathrm{x}, t)$ 是一个函数，x是N维向量。 $\gamma(t)$ 是固定t，并且令 $\psi(\mathrm{x}, t)=0$ 时得到的N-1维流型，也叫做level set。考虑level set上的一个点x(t)，有
$\psi(\mathbf{x}(t), t)=C$
The particle speed $\partial \mathbf{x} / \partial t$ in the direction $\mathbf{n}$ normal to $\gamma(t)$ is given by the speed function $F$ . Thus,
$\frac{\partial \mathbf{x}}{\partial t} \cdot \mathbf{n}=F$
where the normal vector $\mathbf{n}$ is given by $\mathbf{n}=\nabla \psi /|\nabla \psi| .$ By the chain rule we get,
$\psi_{t}+\frac{\partial \mathbf{x}}{\partial t} \cdot \nabla \psi=0$
and substitution yields
$\psi_{t}+F|\nabla \psi|=0$

上述构成一个hamilton jacobi equation，其中F是每一个点的变化速度，可以用曲率表示（用来平滑曲面），也可以通过学习得到。This equation can be solved using the stable, entropy-satisfying nite dierence schemes, borrowed from the literature on hyperbolic conservation laws
Explicit
所有的点都给出，或者通过 parameter mapping
比如uv空间中定义一些点，这些点可以通过函数映射到空间中。显示表示的好处是方便采样，但是判断一个点是否在表面内外就变得困难

point cloud
polygon mesh
subdivision, NURBS
bezier曲线的思想就是对初始集合中每一对前后点按照比例t进行线性插值，然后得到下一个集合，以此类推，直到最后只剩下唯一的点，得到的曲线是t的函数

如果直接用公式来计算bezier曲线的话，那么如果有n+1个控制点就得到n阶多项式
$\mathbf{b}^{n}(t)=\sum_{j=0}^{n} \mathbf{b}_{j} B_{j}^{n}(t)$
其中 $B_{i}^{n}(t)=\left(\begin{array}{l} n \\ i \end{array}\right) t^{i}(1-t)^{n-i}$ 可以看出多项式在t=0的值和b0相等，在t=1时的值和bn相等

从上图可以看出，组合数C(2,3)等于于b2到达 $b_{0}^{3}$ 的路径数。

bezier曲线还有一个性质是得到的曲线在控制点组成的凸包内，所以如果控制点在一条直线上，那么曲线也必须在同一条直线上。

对于分段bezier，一般去四个控制点为一段，所以每一段是3次多项式，前一段的最后一个控制点等于下一段第0个控制点。可以证明b0和b1点的切线(方向和大小)为 $\mathbf{b}^{\prime}(0)=3\left(\mathbf{b}_{1}-\mathbf{b}_{0}\right) ; \quad \mathbf{b}^{\prime}(1)=3\left(\mathbf{b}_{3}-\mathbf{b}_{2}\right)$ , 所以如果想让段与段之间的交接点的切线连续的话必须满足 $\mathbf{a}_{n}=\mathbf{b}_{0}=\frac{1}{2}\left(\mathbf{a}_{n-1}+\mathbf{b}_{1}\right)$

对于曲面来说，一个曲面由4x4个控制点表示，首先在水平方向对t采样得到4条样条曲线，然后固定不同的t值，利用四条曲线上对应的四个点对v进行采样从而得到垂直方向的曲线。

B样条的改进之处是某个控制点的改变只会影响某一范围内的样条曲线。

曲面细分和简化

细分
loop subvision
以每个三角形边的中点作为新顶点，并且用周围老顶点的加权平均更新新节点的坐标，
对于老节点，用自身和邻居老节点更新更新坐标。

catmull-clark subvision
能应用于任何mesh，比如四边形mesh
在每一面的中点和1/2边处增加新顶点。
在这里插入图片描述

上图本来有两个三角形，但是经过一次迭代后变成了6个四边形。

简化
每次选择一条边，然后把这条边坍缩成一个点。
坍缩后点移动到离坍缩之前相关面最近的位置

USTC ps1

坑点
必须setMouseTracking(true); 才能在不按下鼠标的情况下追踪鼠标移动事件。

USTC ps2

image warping 任务是在图像上画n条直线，寻找一个函数能把所有直线的原点转换到直线的终点。在Inverse distanc-weighted interpolation算法中，对于每一对控制点 $p i$ , $q i$ ，都有一个局部插值函数 $f_{i}(p)=q_{i}+D_{i}\left(p-p_{i}\right)$ ，整体的插值函数是各个局部插值函数的加权和，权重用 $w_{i}(p)=\frac{\sigma_{i}(p)}{\sum_{j=1}^{n} \sigma_{j}(p)}$ 表示，其中 $\sigma_{i}(\mathbf{x})=\frac{1}{\left(d_{i}(\mathbf{x})\right)^{\mu}}$ 。 $d_{i}(\mathbf{x})$ 表示当前待预测点和 $p i$ 的距离。
每个局部插值函数的损失函数为 $E_{i}(D)=\sum_{j=1, j \neq i}^{n} w_{i}\left(p_{j}\right)\left\|q_{i}+\left(\begin{array}{ll}d_{11} & d_{12} \\ d_{21} & d_{2}\end{array}\right)\left(p_{j}-p_{i}\right)-q_{j}\right\|^{2}$ ，这可以理解位每个局部插值函数应该尽量拟合其它控制点对，且离得越近的控制点对应该拟合地越好。应为这里 $w i (p j)$ 的分母可能为0，所以这里让 $w_{i}(p_j)=\sigma_{i}\left({p}_{j}\right)$ 。我们把w和p以及q的矩阵形式写在一起，就得到 $E_{i}\left(D_{i}\right)=\left\|\tilde{P}_{i} D_{i}^{T}-\tilde{Q}_{i}\right\|_{F}^{2}$ ，其中 $\tilde{P}_{i}=W_{i}\left(P-P_{i}\right), D=\left(\begin{array}{ll}d_{11} & d_{12} \\ d_{21} & d_{22}\end{array}\right), \tilde{Q_{i}}=W_{i}\left(Q-Q_{i}\right)$ ， $W_{i}=\left(\begin{array}{cccc}\sqrt{w_{i 1}} & & & 0 \\ & \sqrt{w_{i 2}} & & \\ & & \ddots & \\ 0 & & & \sqrt{w_{i n}}\end{array}\right)$ 。虽然这里 $D^T$ 有两列，但是仍然可以用最小二乘法解出 $D_{i}^{T}=\left(\tilde{P_{i}^{T}} \tilde{P_{i}}\right)^{-1} P_{i}^{T} \tilde{Q}_{i}$ .

//
// Created by van on 2021/12/5.
//
#include "math.h"
#include "IDWWarper.h"
#define DIST_2(X1,Y1,X2,Y2) ((X1-X2)*(X1-X2)+(Y1-Y2)*(Y1-Y2))

IDWWarper IDWWarper::inst;
void IDWWarper::initialize(vector<pair<QPoint, QPoint>> pairs, QPoint origin) {
    int n=pairs.size();
    P_points= MatrixXf::Zero(n,2);
    Q_points= MatrixXf::Zero(n,2);
    D_arrays=MatrixXf::Zero(2,2*n);

    for(int i=0;i<n;i++){
        P_points(i,0)=pairs[i].first.x()-origin.x();
        P_points(i,1)=pairs[i].first.y()-origin.y();
        Q_points(i,0)=pairs[i].second.x()-origin.x();
        Q_points(i,1)=pairs[i].second.y()-origin.y();
    }

    for(int i=0;i<n;i++){
        DiagonalMatrix<float, Eigen::Dynamic> Wi(n);
        float pix=P_points(i,0);
        float piy=P_points(i,1);
        float qix=Q_points(i,0);
        float qiy=Q_points(i,1);
        for(int j=0;j<n;j++){
            float pjx=P_points(j,0);
            float pjy=P_points(j,1);
            Wi.diagonal()[j]= sqrt(1.0 / DIST_2(pix,piy,pjx,pjy));
        }
        Wi.diagonal()[i]=0.0;
        Vector2f pi(pix, piy);
        Vector2f qi(qix, qiy);
        MatrixXf Pi=P_points;
        MatrixXf Qi=Q_points;
        Pi.rowwise()-=pi.transpose();
        Qi.rowwise()-=qi.transpose();

        Pi=Wi*Pi;
        Qi=Wi*Qi;
        MatrixXf DTi(2,2);
        DTi=((Pi.transpose()*Pi).inverse())*(Pi.transpose())*Qi;
        D_arrays.block<2,2>(0,2*i)=DTi;
    }

}

QPoint IDWWarper::trans(QPoint p,QPoint origin) {
    double w_sum=0;
    double px=p.x()-origin.x();
    double py=p.y()-origin.y();
    int n=P_points.rows();
    for(int i=0;i<n;i++){
        w_sum+= 1.0 /DIST_2(px,py,P_points(i,0),P_points(i,1));
    }
    DiagonalMatrix<float, Eigen::Dynamic> W(n);
    for(int j=0;j<n;j++){
        float pjx=P_points(j,0);
        float pjy=P_points(j,1);
        W.diagonal()[j]= (1.0 / DIST_2(px,py,pjx,pjy))/w_sum;
    }
    Vector2f p_vector(px, py);
    MatrixXf P=-P_points;
    P.rowwise()+=p_vector.transpose();
    MatrixXf PDT(n,2);
    for(int i=0;i<n;i++){
        PDT.block<1,2>(i,0)=(P.block<1,2>(i,0))*(D_arrays.block<2,2>(0,2*i));
    }
    PDT=W*(PDT+Q_points);
    Vector2f result=PDT.colwise().sum();
    return QPoint(result.x(),result.y());
}

对于变换之后产生的空洞，用ann库搜索出的最近邻有效点进行填充：

void ImageWidget::IDWWarp() {
    vector<QPoint> coods;
    int w = ptr_image_->width();
    int h = ptr_image_->height();
    IDWWarper::inst.initialize(pairs_,origin);
    QImage new_image(w,h,ptr_image_->format());
    new_image.fill(QColor(0,0,0));
    for(int i=0;i<w;i++){
        for(int j=0;j<h;j++){
            QColor color=ptr_image_->pixelColor(i,j);
            QPoint p(i,j);
            QPoint q;
            q=IDWWarper::inst.trans(p,QPoint(0,0));
            if(q.x()<0||q.x()>=w||q.y()<0||q.y()>=h){
//                qDebug("out range, at (%d,%d)",q.x(),q.y());
                continue;
            }
            coods.push_back(q);
            new_image.setPixelColor(q,color);
        }
    }
    ANNpointArray ptsArr = annAllocPts(coods.size(), 2);
    for(int i=0;i<coods.size();i++){
        ptsArr[i][0]=coods[i].x();
        ptsArr[i][1]=coods[i].y();
    }
    ANNbd_tree tree(ptsArr, coods.size(), 2);

    int K=1;
    for(int i=0;i<w;i++){
        for(int j=0;j<h;j++){
            if(new_image.pixelColor(i,j)!=QColor(0,0,0)){
                continue;
            }
            ANNpoint queryPt = annAllocPt(2);
            ANNidx idxArr[K];
            ANNdist distArr[K];
            queryPt[0] = i;
            queryPt[1] = j;
            tree.annkSearch(queryPt, K, idxArr, distArr);
            QPoint near(ptsArr[idxArr[0]][0],ptsArr[idxArr[0]][1]);
            QColor color=new_image.pixelColor(near);
            new_image.setPixelColor(i,j,color);
        }
    }
    *(ptr_image_) = new_image;
    update();
}

坑点，在建立ANNbd_tree时，应确保点集合中没有重复点，而当出现折叠情况时，会将源图像多个点映射到目标图像的同一个点，造成点集重复。

USTC ps3

poisson image editing 需要从source 图片中选出一个区域然后paste到target中的指定位置上。在mainwindow中的QMdiArea可以容纳多个childwindow子窗口，用来显示多张图片。当在主窗口按下rectChoose后,当前获得焦点的childwindow会被保存到主窗口的child_source_属性中，然后子窗口自己负责保存用户选择的矩形或多边形。当在主窗口按下paste后，会把child_source传给当前获得焦点的childwindow(target 图片)。poisson image editing 中从source图片复制到target图片中的区域 $\Omega$ 应该满足边界 $\partial \Omega$ 和目标图片相等，而内部的梯度尽量和source图片保持一致，即
$\min _{f} \iint_{\Omega}|\nabla f-\mathbf{v}|^{2}$ with $\left.f\right|_{\partial \Omega}=\left.f^{*}\right|_{\partial \Omega}$
它的解是带有 Dirichlet boundary 条件的泊松方程
$\Delta f=\operatorname{divv} \text { over } \Omega, \text { with }\left.f\right|_{\partial \Omega}=\left.f^{*}\right|_{\partial \Omega},$
where $\operatorname{div} \mathbf{v}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$ is the divergence of $\mathbf{v}=(u, v)$ .

它离散形式为
for all $\in \Omega, \quad\left|N_{p}\right| f_{p}-\sum_{q \in N_{p} \cap \Omega} f_{q}=\sum_{q \in N_{p} \cap \partial \Omega} f_{q}^{*}+\sum_{q \in N_{p}} v_{p q}$
其中 $N_{p}$ 是p点上下左右四个点组成的点集。for all $\langle p, q\rangle, v_{p q}=g_{p}-g_{q}$ ，g代表source图片。

剩下的工作室如何把用多边形表示的边界转换成一个mask.
这里用到多边形扫描转换算法。
首先建立一个边表(map数据结构)，key对应边下端点的y值，value是这个y值相等的所有边组成的集合。每条边保存以下数据。
在这里插入图片描述

假设key最小值为1，从y=1开始向上扫描，首先把key=1对应的边集作为初始活动边表。对活动边表按x值进行排序，每两个x值所夹的区域就是mask矩阵中应该设置为1的区域。y值加1时，应该考虑y是否超过了活动边表中某些边的ymax，如果是的话就去除这些边。另一方面，如果新的y值正好对应边表中的某个key，那么就应该把该key值对应的边集，加入活动边表中。那些本来就在边表中的边，按照x=x+ $\Delta x$ 更新x的值。

const Eigen::Matrix<int, -1, -1> &Polygon::getMaskMatrix() {
    int n=control_points.size();
    int h=max_y-min_y+1;
    int w=max_x-min_x+1;
    map<int,set<Triple>> edge_table;
    mask_mat = Eigen::Matrix<int, -1, -1>::Zero(h, w);
    for(int i=0;i<n;i++){
        QPoint s=control_points[i];
        QPoint e=control_points[(i+1)%n];
        QPoint low;
        QPoint up;
        if(s.y()==e.y()){
            continue;;
        }
        if(s.y()<e.y()){
            low=s;
            up=e;
        }else{
            low=e;
            up=s;
        }
        float step=float(up.x()-low.x())/float(up.y()-low.y());
        if(edge_table.find(low.y())!=edge_table.end()){
            edge_table[low.y()].insert(Triple(low.x(),step,up.y()));
        }else{
            edge_table[low.y()]=set<Triple>{Triple(low.x(),step,up.y())};
        }
    }

    list<Triple> active_edges;
    for(auto edge:edge_table[min_y]){
        active_edges.push_back(edge);
    }
    assert(active_edges.size()%2==0);
    for(int j=min_y+1;j<=max_y;j++){
        for(auto edge_it=active_edges.begin();edge_it!=active_edges.end();){
            if(j>=edge_it->ymax){
                active_edges.erase(edge_it++);
            }else{
                edge_it->x+=edge_it->step;
                edge_it++;
            }
        }
        if(edge_table.find(j)!=edge_table.end()){
            for(auto edge:edge_table[j]){
                active_edges.push_back(edge);
            }
        }
        assert(active_edges.size()%2==0);
        active_edges.sort();
        if(active_edges.size()==0){
            break;
        }

        for(auto edge_it=active_edges.begin();edge_it!=active_edges.end();){
            int left=edge_it->x;
            edge_it++;
            int right=edge_it->x;
            for(int i=left;i<=right;i++){
                if(i>min_x && i<max_x){
                    mask_mat(j-min_y,i-min_x)=1;
                }
            }
            edge_it++;
        }

    }
    return mask_mat;
}

坑点：应该确保活动边表中任意时刻都有偶数条边。

USTC PS4

$K=\lim _{\Delta s \rightarrow 0}\left|\frac{\Delta \alpha}{\Delta s}\right|=\left|\frac{d \alpha}{d s}\right|$ 曲率定义为切线转角和弧长之比极限的绝对值。
把M点的切线角度用一阶导数表示 $\alpha=arctany^{\prime}$
$\begin{aligned} \frac{d \alpha}{d x}=& \frac{y^{\prime \prime}}{1+\tan ^{2} \alpha}=\frac{y \alpha}{1+y^{\prime 2}} \\ & d \alpha=\frac{y^{\prime \prime}}{1+\left(y^{\prime 2}\right)} d x \end{aligned}$
又 $s=\sqrt{1+y^{\prime 2}} d x$ , 故曲线 $L$ 在 $M$ 点处的曲率为
$K=\frac{\left|y^{\prime \prime}\right|}{\left(1+y^{\prime 2}\right)^{\frac{3}{2}}}$

在平面上给定一些采样点，如何用三角mesh将这些点连接起来的问题叫做三角剖分。
在这里插入图片描述 voronoi图由相邻两点的垂直平分线组成的多边形构成。

delaunay是它的对偶图。 delaunay是比较好的三角剖分，可以证明delaunay网格中的最小角在所有三角剖分中最大，每个三角的外接圆内不包含任何其他点(empty sphere 性质)。并且delaunay图是唯一的。
可以从一般的三角剖分通过迭代的方法获得delaunay图
在这里插入图片描述
如果四边形中两个三角不满足empty sphere 性质，就改变对角线。

如果采样点本身分布很不均匀，那么即便使用delaunay也得不到很好的mesh。所以应该使采样点尽量和voronoi图中多边形的重心重合。
可以通优化能量函数来寻找最优的采样点
$F(X)=\sum_{i=1}^{N} \int_{V_{i}} \rho(\mathbf{x})\left\|\mathbf{x}-\mathbf{x}_{i}\right\|^{2} d \mathbf{x}$
在这里插入图片描述
欧拉公式：对于简单多面体来说，顶点数V，面数F，边数E，面上的洞数H，贯穿物体的洞数G满足
V+F-E=2+H-2G

二维流型：曲面上任意顶点一个无穷小的领域，拓扑同胚于平面上的一个圆盘

boundary edge：只有一个邻面的边
regular edge: 有两个邻面的边。
如果一个mesh中一条边有两个以上邻面，那么它就不属于流型。
下面这种情况也不属于流型，即在某个顶点上又长出多个三角形
在这里插入图片描述

curry dog

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
tinyrender games101笔记

draw line给定(x0,y0) (x1,y1)两个点，画出一条直线。首先通过交换x,y轴使得变化较快的那个轴变为x轴，再交换x0,x1使得x1永远大于x0。令dx=x1-x0;dy=y1-y0;从(x0,y0)开始，x轴每增加一个像素，因为y轴变化较慢所以增加的值dy/dx小于1，把增加的值累积起来用变量error记录，每当error>0.5就令y的值增加1，同时error的值减去1. 但是这样做的缺点是需要使用0.5和dy/dx两个浮点数，所以直接两者同时乘以2dx。void lin
复制链接

扫一扫