从0开始的图形学大师（1~4章）

从0开始的图形学大师

标签（空格分隔）：图形学 peter Shirley <计算机图形学>

这篇blog是看《计算机图形学》（绿皮书）的笔记，博主个人保存。

第1章引言

1.1 图形学领域

建模。对形状和外观性质进行数学定义，而且能够存储在计算机中。例如：一只咖啡杯可以表示为有序的三维点集，这些三维点按一定的插值规则连接起来，还可以建立咖啡杯受光线照射的反射模型。
绘制。该术语来自艺术领域，根据三维计算机模型生成带阴影的图像。
动画。大家都知道的，其中用到了建模和绘制。

。。。

1.2 主要应用

仿真。
CAD。
游戏/动画。

1.3 图像学API

解释API是什么，简单介绍两种两种API模式

##1.4 三维几何模型

三维图像采集有时要借用距离扫描仪（淘宝上有卖），这些距离扫描仪得到的就是三维几何模型，最常用的模型由三维三角形组成，这些三角形共享顶点，常被称作三角形网格。美术设计人员用设计软件也可以生成这样的三维几何模型。

Ps：为何常见的是三维三角形呢？

三角形重心好算，建系可以根据三条边来建，以三条边作为基准向量。

##1.5 图形流水线（graphic pipline）

图形流水线是现在计算机都有的一个软硬件子系统。包含很多对矩阵和向量进行高效处理和组合运算的机制。

多数图形流水线的速度都与正在画的三角形的数量有关。因为交互能力要比视觉效果更重要，所以在生成模型的过程中需要将三角形数量减到最少。另外，如果模型出现在不远处，所需要的三角形数量就比模型在近处时少。这就意味着，在表示模型时要采用不同的细节等级（LOD）（原书摘抄）

1.6 数值问题

图像学程序实际是三维数值代码，所以表示数值很重要（这让我想起图书馆的《c语言数值算法》这本书，还没看过）。幸运的是，现代计算机都服从IEEE浮点标准。后面举了个例子说明IEEE的高效。

1.7 效率

效率往往是各种利弊的权衡，往往没有超级规则。

1.8 软件工程（重点）

图形学关键部分都有比较好的几何基本代码。哪些声明为内联函数，哪些用单精度哪些用双精度不清楚，哪些用const，需要包含很多assert宏？这一节是重点。

第2章数学知识

作为大学生，这些高中知识我们应该都知道了。这里就不讲了。

线性插值：告诉你几个点，用直线连起来。

关于三角形重心要求看懂

第3章光栅算法

3.1 光栅显像

光栅=像素矩阵

3.2 显示器亮度和γ值

我们显示器显示的信号和我们输入的信号往往不同，于是我们要进行伽马矫正，

3.3 RGB

24位系统 = 1个字节（0~255）* 8位 * 3基色

3.4 α通道

如图

3.5 直线绘制（重点）

画一个点集

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 画一个点集*/
#include<iostream>
#include<cstdio>using namespace std;class Vec
{public:int x, y, z;
};int f(int x, int y)
{return (2*x-3*y);
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];for(int y = 0; y < h; y++){fprintf(stderr, "\rRendering %5.2f%%", 100.0*y/(h-1));for(int x = 0; x < w; x++){int i = (h - y - 1) * w + x;if(f(x, y) == 0){c[i].x = 255;c[i].y = 255;c[i].z = 255;}else{c[i].x = 0;c[i].y = 0;c[i].z = 0;}}}FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for (int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

3.5.1 基于隐式方程绘制直线

中点算法（Pitteway，1967）

重要假设：我们能绘出没有间隔的最细的直线。两对角像素之间的连接被认为不产生间隔。

我们先讨论 0<m(斜率)<1 的情况。其他3种情况（-1~0 和 1~正无穷和负无穷~-1）需要另外讨论。

我们知道，把任意一个点带入隐式直线方程，根据正负可以判断点在直线上方还是下方。在我们讨论的情况基础上，我们假设p1（x0,y0）点在p2（x1,y1）点的左下方。

我们考虑x = x0+1的这个点的位置，由于直线像素点相邻，我们有两种选择，为（x0+1,y0）和（x0+1,y0+1）即（x0,y0）向右或向右上方加一个像素，选择哪个呢？我们把中点（x0+1,y0+0.5）代入隐式直线方程，如果大于0就往右走，小于等于0往右上走。

即：往左遍历，根据条件向上走。

得出

 y = y0d = f(x0+1, y0+0.5)for x = x0 to x1 dodraw(x,y)if d < 0 theny = y+1;d = d+(x1-x0)+(y0-y1)elsed = d+(y0-y1)

这里我们用了增量算法，不必在循环内每次都计算d的值，我们用上一次的d的值通过加法计算d的值，以代替原来算d时的乘法运算。

因为我们可以推导出下面两个公式:

f ( x + 1 , y ) = f ( x , y ) + ( y 0 − y 1 ) f(x+1,y) = f(x,y)+(y0-y1) f(x+1,y)=f(x,y)+(y0−y1)

f ( x + 1 , y + 1 ) = f ( x , y ) + ( y 0 − y 1 ) + ( x 1 − x 0 ) f(x+1,y+1) = f(x,y) + (y0-y1) + (x1-x0) f(x+1,y+1)=f(x,y)+(y0−y1)+(x1−x0)

完整代码为

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 中点法画直线*/
#include<iostream>
#include<cstdio>using namespace std;class Vec
{public:int x, y, z;};class Point2
{public:int x, y;Point2(int x_, int y_){x = x_;y = y_;}
};double f(double x, double y)
{return (2*x-3*y+100);
}void drawline(Point2 p1, Point2 p2, Vec * c, int h, int w)
{int y = p1.y;double d = f(p1.x+1, p1.y+0.5);for(int x = p1.x; x <= p2.x; x++){int i = (h - y - 1) * w + x;if(d < 0){y = y + 1;d = d + (p2.x - p1.x) + (p1.y - p2.y);}elsed = d + (p1.y - p2.y);c[i].x = 255;c[i].y = 255;c[i].z = 255;}
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];for(int y = 0; y < h; y++){fprintf(stderr, "\rRendering %5.2f%%", 100.0*y/(h-1));for(int x = 0; x < w; x++){int i = (h - y - 1) * w + x;if(f(x, y) == .0f){c[i].x = 255;c[i].y = 255;c[i].z = 255;}else{c[i].x = 0;c[i].y = 0;c[i].z = 0;}}}Point2 p1(0, 0),p2(1000, 200);drawline(p1, p2, c, h, w);FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for (int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

我们用两点式，这样不用写f()函数：

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 中点法画直线*/
#include<iostream>
#include<cstdio>using namespace std;class Vec
{public:int x, y, z;};class Point2
{public:int x, y;Point2(int x_, int y_){x = x_;y = y_;}
};void drawline(Point2 p1, Point2 p2, Vec * c, int h, int w)
{int y0 = p1.y, x0 = p1.x;int y1 = p2.y, x1 = p2.x;int y = y0;double d = (y0-y1)*(x0+1)+(x1-x0)*(y0+0.5)+x0*y1-x1*y0;for(int x = x0; x <= x1; x++){int i = (h - y - 1) * w + x;if(d < .0f){y = y + 1;d = d + (x1 - x0) + (y0 - y1);}elsed = d + (y0 - y1);c[i].x = 255;c[i].y = 255;c[i].z = 255;}
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];for(int y = 0; y < h; y++){fprintf(stderr, "\rRendering %5.2f%%", 100.0*y/(h-1));for(int x = 0; x < w; x++){int i = (h - y - 1) * w + x;if(2*x-3*y == 0){c[i].x = 255;c[i].y = 255;c[i].z = 255;}else{c[i].x = 0;c[i].y = 0;c[i].z = 0;}}}Point2 p1(0, 0),p2(1000, 200);drawline(p1, p2, c, h, w);FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for (int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

增量算法会使代码变快，但结果会使误差得以积累。

有时，如果只有整数进行运算，算法会执行得更快。但是我们初始化时用到了 d = f(x0+1, y0+0.5)于是我们用2f(x,y)代替f(x,y)。

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 中点法画直线*/
#include<iostream>
#include<cstdio>using namespace std;class Vec
{public:int x, y, z;};class Point2
{public:int x, y;Point2(int x_, int y_){x = x_;y = y_;}
};void drawline(Point2 p1, Point2 p2, Vec * c, int h, int w)
{int y0 = p1.y, x0 = p1.x;int y1 = p2.y, x1 = p2.x;int y = y0;int d = 2*(y0-y1)*(x0+1)+(x1-x0)*(2*y0+1)+2*x0*y1-2*x1*y0;for(int x = x0; x <= x1; x++){int i = (h - y - 1) * w + x;if(d < 0){y = y + 1;d = d + 2*(x1 - x0) + 2*(y0 - y1);}elsed = d + 2*(y0 - y1);c[i].x = 255;c[i].y = 255;c[i].z = 255;}
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];for(int y = 0; y < h; y++){fprintf(stderr, "\rRendering %5.2f%%", 100.0*y/(h-1));for(int x = 0; x < w; x++){int i = (h - y - 1) * w + x;if(2*x-3*y == 0){c[i].x = 255;c[i].y = 255;c[i].z = 255;}else{c[i].x = 0;c[i].y = 0;c[i].z = 0;}}}Point2 p1(0, 0),p2(1000, 200);drawline(p1, p2, c, h, w);FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for (int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

还有就是我们在运行这段代码时，要检查0<m<1这一前提条件，原来我们的假设是遍历x选择一个y，如果条件为1<m<正无穷，直线在一个x上会有多个像素点，就不是扩展一个点的问题了，每次只多选一的算法会使直线变短。（你可以画一画验证一下），我们判断0<m<1一般用下面的条件：

( ( y 1 > = y 0 ) a n d ( x 1 − x 0 > y 1 − y 0 ) ) = ( m ∈ [ 0 , 1 ) ) ((y1>=y0)and(x1-x0>y1-y0)) = (m \in [0,1)) ((y1>=y0)and(x1−x0>y1−y0))=(m∈[0,1))

3.5.2 基于参数方程绘制直线

我们考虑-1<m<1的情况，那么for循环可以沿x轴进行。

参数方程画法就是遍历x，用x求t，然后用t求y。

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 参数方程画直线*/
#include<iostream>
#include<cstdio>
#include<cmath>using namespace std;class Vec
{public:int x, y, z;
};class Point2
{public:int x, y;Vec color;Point2(int x_, int y_){x = x_;y = y_;}
};//double round(double r)
//{//    return (r > 0.0)?floor(r + 0.5):ceil(r - 0.5);
//}void drawline(Point2 p1, Point2 p2, Vec * c, int h, int w)
{int y0 = p1.y, x0 = p1.x;int y1 = p2.y, x1 = p2.x;double y = y0;double yy = 1.0*(y1-y0)/(x1-x0);for(int x = x0; x <= x1; x++){int i = (h - round(y) - 1) * w + x;c[i].x = 255;c[i].y = 255;c[i].z = 255;y = y + yy;}
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];for(int y = 0; y < h; y++){fprintf(stderr, "\rRendering %5.2f%%", 100.0*y/(h-1));for(int x = 0; x < w; x++){int i = (h - y - 1) * w + x;if(2*x-3*y == 0){c[i].x = 255;c[i].y = 255;c[i].z = 255;}else{c[i].x = 0;c[i].y = 0;c[i].z = 0;}}}Point2 p1(0, 0),p2(1000, 200);drawline(p1, p2, c, h, w);FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for (int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

我们还可以画一条彩色的渐变线

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 参数方程画直线*/
#include<iostream>
#include<cstdio>
#include<cmath>using namespace std;class Vec
{public:int x, y, z;Vec(int x_ = 0, int y_ = 0, int z_= 0):x(x_),y(y_),z(z_){}
};class Point2
{public:int x, y;Vec color;Point2(int x_, int y_, int r, int g, int b):x(x_),y(y_),color(r,g,b){}
};//double round(double r)
//{//    return (r > 0.0)?floor(r + 0.5):ceil(r - 0.5);
//}void drawline(Point2 p1, Point2 p2, Vec * c, int h, int w)
{int y0 = p1.y, x0 = p1.x;int y1 = p2.y, x1 = p2.x;int r1 = p2.color.x, g1 = p2.color.y, b1 = p2.color.z;int r0 = p1.color.x, g0 = p1.color.y, b0 = p1.color.z;double y = y0, r = r0, g = g0, b = b0;double yy = 1.0*(y1-y0)/(x1-x0);double rr = 1.0*(r1-r0)/(x1-x0);double gg = 1.0*(g1-g0)/(x1-x0);double bb = 1.0*(b1-b0)/(x1-x0);for(int x = x0; x <= x1; x++){int i = (h - round(y) - 1) * w + x;r = r + rr;g = g + gg;b = b + bb;c[i].x = round(r);c[i].y = round(g);c[i].z = round(b);y = y + yy;}
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];for(int y = 0; y < h; y++){fprintf(stderr, "\rRendering %5.2f%%", 100.0*y/(h-1));for(int x = 0; x < w; x++){int i = (h - y - 1) * w + x;if(2*x-3*y == 0){c[i].x = 255;c[i].y = 255;c[i].z = 255;}else{c[i].x = 0;c[i].y = 0;c[i].z = 0;}}}Point2 p1(0, 0, 255, 0, 0),p2(1000, 200, 0, 0, 255);drawline(p1, p2, c, h, w);FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for (int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

3.6 三角形光栅化(重点)

我们用重心坐标系可以很清楚的画出一个三角形，而且如果我们画一个最简单的彩色三角形，也可以用顶点颜色（假设为c0,c1,c2）表示：c = αc0+βc1+γc2。这种颜色插值方法为Gouraud插值法（作者Gouraud，1971）。

暴力光栅化伪代码

for all x dofor all y docompute(α,β,γ)for(x,y)if(0<=α<=1 and 0<=β<=1 and 0<=γ<=1)thenc = αc0+βc1+γc2drawpixel(x,y) with color c

我们不用遍历所有像素点，只需要遍历三角形最小包围矩形。

x_min = floor(x0,x1,x2)//floor指求最小
x_max = ceiling(x0,x1,x2)//ceiling指求最大
y_min = floor(y0,y1,y2)
y_max = ceiling(y0,y1,y2)
for y=y_min to y_max do//一般都是从y开始遍历，一行一行来for x=x_min to x_maxi doα = f12(x,y)/f12(x0,y0)β = f20(x,y)/f20(x1,y1)γ = f01(x,y)/f01(x2,y2)if(α>0 and β>0 and γ>0)thenc = αc0+βc1+γc2drawpixel(x,y) with color c其中fij表示编号为i与j的直线，α>0表示点在(x1,y1)与(x2,y2)形成的这条直线靠近x0的五边形区域。
f01(x,y) = (y0-y1)x+(x1-x0)y+x0y1-x1y0xb
f01(x,y) = (y1-y2)x+(x2-x1)y+x1y2-x2y1
f01(x,y) = (y2-y0)x+(x0-x2)y+x2y0-x0y2

完整代码为：

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 三角形光栅化（画一个彩色三角形）*/#include<iostream>
#include<cstdio>
#include<cmath>using namespace std;
int max(int a, int b, int c)
{if(a >= b && a >= c)return a;else if(b >= a && b >= c)return b;elsereturn c;
}
int min(int a, int b, int c)
{if(a <= b && a <= c)return a;else if(b <= a && b <= c)return b;elsereturn c;
}class Vec
{public:int x, y, z;Vec(int x_ = 0, int y_ = 0, int z_= 0):x(x_),y(y_),z(z_){}Vec operator+(const Vec &b) const{ return Vec(x + b.x, y + b.y, z + b.z); }
};class Point2
{public:int x, y;Vec color;Point2(int x_, int y_, int r, int g, int b):x(x_),y(y_),color(r,g,b){}
};//double round(double r)
//{//    return (r > 0.0)?floor(r + 0.5):ceil(r - 0.5);
//}int fij(int x, int y, int x0, int y0, int x1, int y1)
{return (y0-y1)*x+(x1-x0)*y+x0*y1-x1*y0;
}void drawtriangle(Point2 p0, Point2 p1, Point2 p2, Vec *c, int w, int h)
{int y0 = p0.y, x0 = p0.x;int y1 = p1.y, x1 = p1.x;int y2 = p2.y, x2 = p2.x;int x_min = min(x0, x1, x2);int x_max = max(x0, x1, x2);int y_min = min(y0, y1, y2);int y_max = max(y0, y1, y2);double a,b,r;for(int y = y_min; y <= y_max; y++){for(int x = x_min; x <= x_max; x++){int i = (h - y - 1) * w + x;a = 1.0*fij(x,y,x1,y1,x2,y2)/fij(x0,y0,x1,y1,x2,y2);b = 1.0*fij(x,y,x2,y2,x0,y0)/fij(x1,y1,x2,y2,x0,y0);r = 1.0*fij(x,y,x0,y0,x1,y1)/fij(x2,y2,x0,y0,x1,y1);if(a > 0 && b > 0 && r > 0){c[i].x = round(a*p0.color.x + b*p1.color.x + r*p2.color.x);c[i].y = round(a*p0.color.y + b*p1.color.y + r*p2.color.y);c[i].z = round(a*p0.color.z + b*p1.color.z + r*p2.color.z);}
//            else
//            {//                c[i].x = 0;
//                c[i].y = 0;
//                c[i].z = 0;
//            }}}
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];Point2 p0(30, 30, 255, 0, 0),p1(30, 300, 0, 255, 0),p2(300, 30, 0, 0, 255);drawtriangle(p0, p1, p2, c, w, h);FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for(int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

处理三角形边上的像素

我们发现fij(x,y)>0（x,y随便定）和α>0在两个相邻三角形中往往只有一个满足。

于是我们规定在边上的那些像素点必须满足fij(-1,-1)>0才画，伪代码如下：

x_min = floor(x0,x1,x2)//floor指求最小
x_max = ceiling(x0,x1,x2)//ceiling指求最大
y_min = floor(y0,y1,y2)
y_max = ceiling(y0,y1,y2)
for y=y_min to y_max do//一般都是从y开始遍历，一行一行来for x=x_min to x_maxi doα = f12(x,y)/f12(x0,y0)β = f20(x,y)/f20(x1,y1)γ = f01(x,y)/f01(x2,y2)if(α>0 and β>0 and γ>0)thenif(α>0 or f12(-1,-1)>0) and (β>0 or f20(-1,-1)>0) and (γ>0 or f01(-1,-1)>0) thenc = αc0+βc1+γc2drawpixel(x,y) with color cf01(x,y) = (y0-y1)x+(x1-x0)y+x0y1-x1y0xb
f01(x,y) = (y1-y2)x+(x2-x1)y+x1y2-x2y1
f01(x,y) = (y2-y0)x+(x0-x2)y+x2y0-x0y2

作者指出，这些以上的这些代码还有很多需要改进的，比如如果α是负的，就不用计算β和γ了，对于输入为一条直线的三角形，可能出现除数为0的情况，应采用浮点误差条件。

3.7 简单反走样技术

用盒状滤波即可

/*** @title 从0开始的图形学大师* @author hezenggeng* @code 中点法画直线+不大恰当的均值滤波反走样*/
#include<iostream>
#include<cstdio>using namespace std;class Vec
{public:int x, y, z;Vec(int x_ = 255, int y_ = 255, int z_ = 255):x(x_),y(y_),z(z_){}
};class Point2
{public:int x, y;Point2(int x_, int y_){x = x_;y = y_;}
};void drawline(Point2 p1, Point2 p2, Vec * c, int h, int w)
{int y0 = p1.y, x0 = p1.x;int y1 = p2.y, x1 = p2.x;int y = y0;int d = 2*(y0-y1)*(x0+1)+(x1-x0)*(2*y0+1)+2*x0*y1-2*x1*y0;Vec *temp = new Vec[w*h];for(int x = x0; x <= x1; x++){int i = (h - y - 1) * w + x;if(d < 0){y = y + 1;d = d + 2*(x1 - x0) + 2*(y0 - y1);}elsed = d + 2*(y0 - y1);temp[i].x = 0;temp[i].y = 0;temp[i].z = 0;}//反走样的不大恰当的演示,我用九宫格的均值代替中间那个值（均值滤波）for(y = y0; y <= y1; y++){for(int x = x0; x <= x1; x++){int i00 = (h - y - 2) * w + x-1;int i01 = (h - y - 2) * w + x;int i02 = (h - y - 2) * w + x+1;int i10 = (h - y - 1) * w + x-1;int i11 = (h - y - 1) * w + x;int i12 = (h - y - 1) * w + x+1;int i20 = (h - y) * w + x-1;int i21 = (h - y) * w + x;int i22 = (h - y) * w + x+1;c[i11].x = (temp[i00].x + temp[i01].x + temp[i02].x +temp[i10].x + temp[i11].x + temp[i12].x +temp[i20].x + temp[i21].x + temp[i22].x) / 9;c[i11].y = c[i11].x;c[i11].z = c[i11].x;}}free(temp);
}int main(int argc, char *argv[])
{int w = 1024, h = 768;Vec *c = new Vec[w*h];Point2 p1(10, 10),p2(1000, 200);drawline(p1, p2, c, h, w);FILE *f = fopen("hh.ppm", "w");fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);for (int i = 0; i < h*w; i++){fprintf(f, "%d %d %d ", c[i].x, c[i].y, c[i].z);}free(c);return 0;
}

3.8 图像捕捉与储存

3.8.1 扫描仪和数码摄像机

看看就好

3.8.2 图像储存

gif。有损压缩，只能表示256种颜色。适用于自然图像。
jpeg。有损压缩，以人类视觉系统为阈值。适用于自然图像。
tiff。无损压缩，每个像素通常占24位（3*8）。
ppm。无损压缩，每个像素通常占24位（3*8）。
png。无损格式集合，具有很好的开放源码管理工具。

由于压缩和变形，对于非商业应用，如果没有读/写库可用的话，最好使用原始的ppm。

第4章信号处理

我们在上一章用相邻的像素点表示一条直线，这叫采样表示法。

这部分需要看看书p49页。讲的不错。

20世纪20年代，通信领域已经用到了采样与重构技术。到20世纪40年代，关于采样与重构理论的论述就和现在所用的形势一样了（摘录）

##4.1 数字音频

要接触这类知识，记得日本有一位科普漫画家画一个解释傅里叶变化的漫画，当时女猪脚的设定就是乐队的成员。

##4.2 卷积

符号 f ∗ g f*g f∗g

4.2.1 滑动平均

我们想要知道函数某一段的平均值，需要对该段积分再除以区间长度。

h ( x ) = 1 2 r ∫ x − r x + r g ( t ) d t h(x) = \frac{1}{2r} \int ^{x+r} _{x-r} g(t)dt h(x)=2r1∫x−rx+rg(t)dt
滑动平均可以做滤波，叫滑动平均滤波

4.2.2 离散卷积

( a ∗ b ) [ i ] = ∑ j a [ j ] ∗ b [ i − j ] (a*b)[i] = \sum \limits _{j} a[j]*b[i-j] (a∗b)[i]=j∑a[j]∗b[i−j]

在图形学中，假设a的支撑集有限，（使|j|>r时，都有a[j]=0）此时，上面的求和公式可写为：

( a ∗ b ) [ i ] = ∑ j = − r r a [ j ] ∗ b [ i − j ] (a*b)[i] = \sum \limits ^r _{j = -r} a[j]*b[i-j] (a∗b)[i]=j=−r∑ra[j]∗b[i−j]
理解1：（加权平均）就是b中以i为中心，取b[i+r,i-r]范围（我故意写反的）与a所有值a[-r,r]相乘，得到的变为新的(a*b)[i]。

一般来说，a是滤波器，提供权值，b是信号。

书上伪代码

function convolve(sequence a, sequence b, int r, int i)s = 0for j = -r to rs = s + a[j]b[i-j]
return

####卷积滤波器

当非零区间上存在一个常数，我们称之为盒式滤波（不是何氏滤波）。滑动平均公式也可以由盒式滤波得到：

KaTeX parse error: No such environment: equation at position 16: a[j] = \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ \begin{cases}…

代入上式得滑动平均公式。

卷积性质

卷积的一大优点是：滤波器和信号可以互换。我们用i-k替换j推导得到（推导得到或者通过图形想象）：

( a ∗ b ) [ i ] = ∑ j a [ j ] ∗ b [ i − j ] (a*b)[i] = \sum \limits _{j} a[j]*b[i-j] (a∗b)[i]=j∑a[j]∗b[i−j]
( a ∗ b ) [ i ] = ∑ i − k a [ i − k ] ∗ b [ i − ( i − k ) ] (a*b)[i] = \sum \limits _{i-k} a[i-k]*b[i-(i-k)] (a∗b)[i]=i−k∑a[i−k]∗b[i−(i−k)]
( a ∗ b ) [ i ] = ∑ k a [ i − k ] ∗ b [ k ] (a*b)[i] = \sum \limits _{k} a[i-k]*b[k] (a∗b)[i]=k∑a[i−k]∗b[k]（i-k和k和j都是全集R，于是可以替换）

因此，对于任意序列a和b， ( a ∗ b ) [ i ] = ( b ∗ a ) [ i ] (a*b)[i] = (b*a)[i] (a∗b)[i]=(b∗a)[i]，卷积满足交换律

卷积还满足结合律、对加法满足分配率。

一种简单的滤波器：离散脉冲d[i] = …,0,0,0,1,0,0,0…
( d ∗ b ) [ i ] = ∑ j = 0 j = 0 d [ j ] ∗ b [ i − j ] = b [ i ] (d*b)[i] = \sum \limits ^{j = 0} _{j = 0} d[j] * b[i-j] = b[i] (d∗b)[i]=j=0∑j=0d[j]∗b[i−j]=b[i]
所以d就相当于单位向量或者单位矩阵或者单位元，运算时可以任意加。

$c = b-ab = db-a*b = (d-a)*b $

4.2.3 把卷积看做移位滤波器之和

标题说明这里又产生了卷积的一种新解释：卷积就是b的位移序列之和。

首先定义位移序列： b → j = b [ i − j ] b_{\rightarrow j} = b[i-j] b→j=b[i−j]，这只是一种表示而已，把原来卷积中的 b [ i − j ] b[i-j] b[i−j]换成 b → j b_{\rightarrow j} b→j即可。得到新的公式：
( a ∗ b ) = ∑ j a [ j ] ∗ b → j (a*b) = \sum \limits _j a[j]*b_{\rightarrow j} (a∗b)=j∑a[j]∗b→j
其实按照理解，我们可以写一些不符合数学规范的公式，来更形象的说明公式：
( a ∗ b ) = ∑ j a 序列 ∗ b → j ( b 向右移 j 后的序列 ) (a*b) = \sum \limits _j a序列*b_{\rightarrow j}(b向右移j后的序列) (a∗b)=j∑a序列∗b→j(b向右移j后的序列)
理解2：之前我们是一段a[i+r,i-r]和一段b[-r,r]的乘积得到一个(a*b)[i]。现在我们的解释为，先把每一个a[j]与b[i]相乘（两序列相乘）得到第一个(a*b)，然后再平移b，得到b[i-1]，再相乘每一个a[j]与b[i-1]（b平移得到的序列）得到又一个(a*b)，重复j次平移相乘，得到j个(a*b)，相加所有(a*b)为最终结果。

（ps：）这样我们就有两种算卷积的算法了。

4.2.4 与连续函数的卷积

连续函数卷积用积分代替求和运算。
( f ∗ g ) ( x ) = ∫ − ∞ + ∞ f ( t ) g ( x − t ) d t (f*g)(x) = \int ^{+\infty} _{-\infty} f(t)g(x-t)dt (f∗g)(x)=∫−∞+∞f(t)g(x−t)dt
理解：和卷积的理解1相同，(f*g)(x)就是移动f使f(0)与g(x)对应后，两函数之积形成的曲线下面的面积。

连续函数的卷积也满足交换律、结合律、对加法满足分配率。同样，连续卷积也可以看做是对移位滤波器求和（理解2）。不同之处在于，在连续的情况下，有无限多的位移滤波器拷贝。

之后书中有一个简单的例题P57可以看一看。

狄拉克 δ \delta δ函数(狄拉克脉冲函数)
∫ + ∞ + ∞ δ ( x ) f ( x ) d x = f ( 0 ) \int ^{+\infty}_{+\infty} \delta(x)f(x)dx = f(0) ∫+∞+∞δ(x)f(x)dx=f(0)

( δ ∗ f ) ( x ) = ∫ + ∞ + ∞ δ ( t ) f ( x − t ) d t = f ( x ) 即 δ ∗ f = f (\delta*f)(x) = \int ^{+\infty}_{+\infty} \delta(t)f(x-t)dt = f(x)即\delta*f = f (δ∗f)(x)=∫+∞+∞δ(t)f(x−t)dt=f(x)即δ∗f=f

δ \delta δ函数在0处没有确定值，对非零x都有 δ ( x ) = 0 \delta(x) = 0 δ(x)=0

4.2.5 离散-连续卷积

连续序列变为离散序列：采样即可

a [ i ] = f ( i ) a[i] = f(i) a[i]=f(i)

离散序列变为连续序列：用连续滤波器对离散序列进行滤波，其中x不是离散的，而i是离散的。
( a ∗ f ) ( x ) = ∑ i a [ i ] f ( x − i ) (a*f)(x) = \sum \limits _i a[i]f(x-i) (a∗f)(x)=i∑a[i]f(x−i)
我们用伪代码表示滤波器半径为r的离散-连续卷积(好久没敲了-。-)
( a ∗ f ) ( x ) = ∑ i = ∣ x − r ∣ x + r a [ i ] f ( x − i ) (a*f)(x) = \sum \limits ^{x+r} _{i = |x-r|} a[i]f(x-i) (a∗f)(x)=i=∣x−r∣∑x+ra[i]f(x−i)

function reconstruct(sequence a, filter f, real x)s = 0r = f.radiusfor i=|x-r| to |x+r| dos = s + a[i]f(x-i)
return s

4.2.6 多维卷积

( a ∗ b ) [ i , j ] = ∑ i ′ ∑ j ′ a [ i ′ , j ′ ] ∗ b [ i − i ′ , j − j ′ ] (a*b)[i,j] = \sum \limits _{i'} \sum \limits _{j'} a[i',j']*b[i-i',j-j'] (a∗b)[i,j]=i′∑j′∑a[i′,j′]∗b[i−i′,j−j′]

( a ∗ b ) [ i , j ] = ∑ i ′ = − r i ′ = r ∑ j ′ = − r j ′ = r a [ i ′ , j ′ ] ∗ b [ i − i ′ , j − j ′ ] (a*b)[i,j] = \sum \limits ^{i'=r}_{i'=-r} \sum \limits ^{j'=r}_{j'=-r} a[i',j']*b[i-i',j-j'] (a∗b)[i,j]=i′=−r∑i′=rj′=−r∑j′=ra[i′,j′]∗b[i−i′,j−j′]

用伪代码表示

function convolve2d(filter2d a, filter2d b, int i, int j)s = 0;r = a.radiusfor i' = -r to r dofor j' = -r to r dos = s + a[i'][j']b[i-i'][j-j']
return s

解释：输出样本等于输入面积的加权平均值，其中将二维滤波器作为掩膜，以确定每个样本的权值。

继续推广，可推出二维的连续-连续卷积和离散-连续卷积
( f ∗ g ) ( x , y ) = ∫ ∫ f ( x ′ , y ′ ) ∗ g ( x − x ′ , y − y ′ ) d x ′ d y ′ (f*g)(x,y) = \int\int f(x',y')*g(x-x',y-y')dx'dy' (f∗g)(x,y)=∫∫f(x′,y′)∗g(x−x′,y−y′)dx′dy′

( a ∗ f ) ( x , y ) = ∑ i ∑ j a [ i , j ] ∗ g ( x − i , y − j ) (a*f)(x,y) = \sum \limits _{i} \sum \limits _{j} a[i,j]*g(x-i,y-j) (a∗f)(x,y)=i∑j∑a[i,j]∗g(x−i,y−j)

4.3 卷积滤波器

下面是图形学常用的特殊滤波器

4.3.1 各种卷积滤波器

盒式滤波器 f b o x f _{box} fbox，（存在跳变，边界要引起重视）
帐篷式滤波器 f t e n t f_{tent} ftent，（不存在跳变）
高斯滤波器 f g ( x ) f_g(x) fg(x)，（后面还会学，常用来采样）
三次B样条滤波器 f B ( x ) f_B(x) fB(x)（导数连续，感觉也很平滑，常用来重构）
。。。

为什么高斯滤波很平滑？

4.3.2 滤波器的性质

具有负值的滤波器存在振铃
重构函数继承滤波器连续性
构建可分滤波器最简单就是 f 2 ( x , y ) = f 1 ( x ) f 1 ( y ) f_2(x,y) = f_1(x)f_1(y) f2(x,y)=f1(x)f1(y)

比如高斯滤波器（具有可分性和径向对称性）

可分滤波器实现效率高
KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ \begin{aligned…
伪代码：

function filterImage(image I, filter f)r = f.radiusnx = I.widthny = I.heightallocate storage array S[0,...,nx-1]allocate image Iout[r,...,nx-r-1][r,...,ny-r-1]initialize S and Iout to all zerofor y=r to ny-r-1 dofor x=0 to nx-1 doS[x]=0for i=-r to r doS[x]=S[x]+f[i]I[x][y-i]for x=r to nx-r-1 dofor i=-r to r doIout[x][y]=Iout[x][y]+f[i]S[x-i]
return Iout

4.4 图像信号处理

4.4.1 离散图像滤波

图像模糊化和锐化（给了公式）

产生阴影（给了公式）
KaTeX parse error: No such environment: equation at position 23: …}(i,j) = \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ \begin{cases}…

KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲ \begin{aligned…

opencv3+python3代码实现阴影

import cv2
import numpy as npimg1 = cv2.imread("../images/favorite/luoxiaohei.jpg")
img2 = cv2.imread("../images/favorite/luoxiaohei.jpg")# kernel = np.array([(0, 0, 0, 0, 0),
#                    (0, 0, 0, 0, 0),
#                    (0, 0, 0, 0, 0),
#                    (0, 0, 0, 1, 0),
#                    (0, 0, 0, 0, 0)])kernel = np.zeros((205,205), np.float32)
kernel[100][85] = 1img1 = cv2.filter2D(img1, -1, kernel)
img1 = cv2.GaussianBlur(img1, (205,205), 0)# cv2.imshow("blur", img1)rows,cols,channels = img2.shape
roi = img1[0:rows, 0:cols ]img2gray = cv2.cvtColor(img2,cv2.COLOR_BGR2GRAY)
ret, mask = cv2.threshold(img2gray, 175, 255, cv2.THRESH_BINARY)
# cv2.imshow("mask",mask)mask_inv = cv2.bitwise_not(mask)
# cv2.imshow("mask_inv", mask_inv)img1_bg = cv2.bitwise_and(roi,roi,mask = mask)
# cv2.imshow("img1_bg", img1_bg)img2_fg = cv2.bitwise_and(img2,img2,mask = mask_inv)
# cv2.imshow("img2_fg", img2_fg)dst = cv2.add(img1_bg,img2_fg)
img1[0:rows, 0:cols ] = dstcv2.imshow("img1", img1)cv2.waitKey(20000)cv2.destroyAllWindows()

等以后再写一下c++的

4.4.2 图像采样中的反走样技术

我们在采样中会遇到莫尔图案。因此要用平滑滤波进行处理。

4.4.3 重构与重采样

重采样=重构+采样

我们改变图像大小时，可以用重构来实现重采样，就类似电视或者bilibili改变屏幕大小却不改变内容。

无论放大还是缩小，我们把它看做是像素密度的变化，而非尺寸的变化。（就像化学中相同容积，气体的分子数减小一样）这样我们就可以想出一个方法来求得输出图像的像素值：取图像中与样本点最近的值作为输出值。而这个方法与用一个像素宽的盒式滤波器（假设函数为f）重构图像然后点采样一样。（重构-重采样过程）

我们考虑一维的情况。

function resample(sequence a, float x0, float xx, int n, filter f)create sequence b of length nfor i=0 to n-1 dob[i]=reconstruct(a,f,x0+i*xx)
return b

我们上一节讲过，采样前要进行平滑滤波，假设函数为g。

于是我们把结果表示为a*f*g，我们把f和g先运算了，（f*g）叫做重采样滤波器

在图像的边缘处，这种简单的处理会超过序列的范围，有三种方案：

在序列末尾停止循环，这等于在图像所有边缘上补0
修改序列末尾可以访问，如a[-1] = a[0]
当接近边缘时修改滤波器，使之不会超过序列范围

第一种方法会产生模糊边缘，第二种方法容易实现，第三种方法最好。在图像边缘处，修改滤波器最好的方法是对滤波器重新规范化。

4.5 采样理论

4.5.1 傅里叶变换

看不懂，呜呜呜