投影矩阵推导过程

3d程序经常要用到矩阵算法，
比较常见的如旋转矩阵，平移矩阵，以及投影矩阵
opengl与d3d均有对应的api进行相应的操作。

本文主要介绍一下投影矩阵,
(在阅读irricht与ogre代码时碰到了一些问题,发现视截体是根据投影矩阵计算出来的,

其实也可以根据视角与摄影机的位置与朝向，计算出视截体的,quake就是这么做的,

另外处理纹理阴影时,也需要对投影矩阵有一定程度的了解)

程序内存中保存的是x y z这样的三维数据,最终的显示结果,却只有x y这两个维度。

将三维空间的点转换成屏幕(显存)的二维数据,就是通过投影矩阵来完成的.
投影矩阵其实很简单。

想像一下，近平面与远平面，与原点组成了一个金字塔,

所有的视截体内在三维点,全部投影到近平面上,得到一个二维点
(z变成深度值,对点投影在二维平面上没有影响,这点很重要,在后继的阐述中将会发现,z的值是非线性的,
越接近近平面,值被放大,接近远平面的值则放大系数变小)

再来理一下三维点转换成二维点的过程：
1 世界坐标转换成观察坐标
2 投影变换
3 从坐标中去除齐次坐标的成份(表达的不清楚,实际上是除以点的z坐标值,后面会阐述)

经过旋转平移，我们可以直接假定摄影机就位于0,0,0 对着 0,0,-1(即z轴负向)
而视截体则为一个缺顶的金字塔,由6个平面构成(上下左右,远平面z=n,近平面z=f)
我们假设视角为a(在计算中,这个好像没什么用), 近平面宽w 高 h (远平面不用管,只要知道远平面的z=f即可)
则位于视截体内部的一个点 (x, y, z) 它在近平面上的投影 (x1, y1)
可以表示如下(根据相似三角形)
x1 = x * (n/z)
y1 = y * (n/z)

然后x1 y1 和 z(没有z1) 再转换成标准视空间(一个从 -1,-1,0 延伸到 1,1,1的轴对齐矩形)
即最终的数据(x2, y2, z2)　x2 y2取值一定是从-1,1之间 z2的取值则为0,1
按比例变换即可
x2 = x1 / w
y2 = y1 / h

即
x2 = (x * (n/z))/w
y2 = (y * (n/z))/h
所以投影矩阵其实是一个缩放矩阵,
这里分母包含z,也就是上文提到的齐次坐标的问题,投影矩阵里都是常数(w,h给定情况下)
改写上面的的等式
x2 * z = (x*n)/w
y2 * z = (y*n)/w

然后是z2
因为引入了齐次坐标,自然也希望来个 z2 * z = p*z + q,以便投影矩阵就是一个缩放矩阵
由于 z2取值是 0,1
即 z=n 时 z2等于 0
即 z=f 时 z2等于 1
我们假设这样的一个等式
z2 * z = p * z + q (p q是常数)
分别将z2 = 0, z = n  与 z2 = 1, z = f 代入,即可解出这个方程p 与 q 的值
可以得到
p = f / (f - n)
q = - (f * n) / (f - n)

上面提到过 z2是非线性的,即，接近近平面，被放大,接近远平面，放大倍数变小,在远平面上的点无放大(z2=1)
如果按严格的比例变换,z2的真实值应该是
z2_real = (z - n)/(f - n)  (z2_real严格处在0,1之间)
再来看我们假设求得的z2
p q代入
z2 = (f/z) *( (z-n)/ (f-n) )  ---- 这里笔者做了一些变换,以便说明问题

f/z就是放大的倍数(就是上文反复提到的) 而且 z2也一定处于(0,1）之间,并且单调递增
因为z2并不影响投影结果，它只表征深度，所以这样的变换满足3d绘图需求的。

于是有

x2 * z = (x*n)/w
y2 * z = (y*n)/h
z2 * z = pz + q  (p,q已经在上文求出)　---- <标记1>

则点(x,y,z,1) 这里加上w,因为要解决齐次坐标问题,经过这样一个投影矩阵,变成

(x2*z,  y2*z, z2*z,  z),  去除齐次坐标后 (x2, y2, z2, 1) 就变换到标准视空间了(也就是最终绘图的依据)

现在给出投影矩阵

  n/w, 0,   0, 0    --- 此处的w是宽除2
  0,   n/h, 0, 0     ----  此处h同上
  0,   0,   p, q
  0,   0,   1, 0
 
前面的三行都很好理解,根据<标记1>处的三个等式得来的,
最后一行的1,就是为了引入齐次坐标z.
以便能opengl能正确地去除z

这样的一个矩阵,最后计算结果x2 y2受z影响,就能表现出透视的效果,z2也正确地在0,1之间单调递增(非线性)。

ps:
d3d的标准视空间是 -1,-1,0 1,1,1
opengl的标准的视空间好像是 -1,-1,-1 到 1,1,1 (z的取值好像是-1,1)但这不影响相应的公式推导

w, h的值可以用a与n组合得出(还有纵横比),这也是opengl与d3d那个计算投影矩阵函数要传的参数 a, 纵横比, n, f

附上代码片断与注释

void matrix4::gen_pers_proj_matrix(ant_f32 fovy, ant_f32 aspect, ant_f32 z_near, ant_f32 z_far, matrix4& out_mtx)
{
	/*
	* 假设视截体内的点(x, y, z) 投影到标准视空间里的(x_p, y_p, z_p)
	* 需要把视截体里的点,映身到标准视空间里(-1, -1, 0) (1, 1, 1)
	* 即投影之后 x_p y_p取值为(-1, 1) z_p的取值为(0, 1)
	* 根据相似三角形,投影形成的 x_p y_p肯定与z有关,所以矩阵中需要包含一个齐次空间因子,以便做矩阵乘法后w取值为z
	*
	* 2n/w, 0,    0, 0
	* 0,    2n/h, 0, 0   --- x y只是单纯的缩放(根据相似三角形)
	* 0,    0,    p, q   --- 这两个参数纯粹为了矩阵运算是造出来的,在几何上不存在意义,它们能保证变换后的z_p在(0,1)区间单调递增,并且z会与比例值放大(f/z)倍,这也是很多3d编程书籍中提到的放大倍数
	* 0,    0,    1, 0   --- 1添加齐次空间因子,最终(x, y, z)经过矩阵运算会变成(x_p * z, y_p * z, z_p * z, z),去除齐次坐标因子后得到(x_p, y_p, z_p, 1)
	*
	* p = f / (f - n)
	* q = - (f * n) / (f - n)
	*
	* 以上是d3d的情况,opengl的标准视空间是 (-1, -1, -1) (1, 1, 1),即 z_p的取值范围是(-1, 1)
	* 并且z是负值,因为opengl默认是位于(0, 0, 0)并且看向z负方向,相应的缩放系数会变成
	*  x_p * (-z) = (2n/w) * x
	*  y_p * (-z) = (2n/h) * y
	*  z_p * (-z) = p_o * z + q_o   --- 这里要做一个变通 z虽然是负的,但解的时候,要假设 z在(n,f) 然后z_p对应(-1,1)
	*
	* 解得 p_o = (f+n)/(f-n)  q_o=(-2fn)(f-n),  然后考虑z实际是处于(-f, -n)这个区间,则 p_o要取 -(f+n)/(f-n)
	* 对应的矩阵
	*
	*  2n/w, 0,    0,   0
	*  0,    2n/h, 0,   0
	*  0,    0,    p_o, -1  --- 这个取-1,是因为齐次坐标因子此时是 -z;
	*  0,    0,    q_o, 0   
	*/

	/* 2n/h 实际上就是tan(fovy/2) */
	ant_f32 h = tan(fovy * ANT_AR_CONVERT * 0.5f);
	ant_f32 w = h / aspect;//aspect不可为零
	ant_f32 p = (z_far + z_near) / (z_near - z_far);//分母不可为零
	ant_f32 q = (2 * z_far * z_near) / (z_near - z_far);//分母不可为零

	out_mtx.m_[0] = w;
	out_mtx.m_[1] = 0;
	out_mtx.m_[2] = 0;
	out_mtx.m_[3] = 0;

	out_mtx.m_[4] = 0;
	out_mtx.m_[5] = h;
	out_mtx.m_[6] = 0;
	out_mtx.m_[7] = 0;

	/*  */
	out_mtx.m_[8] = 0;
	out_mtx.m_[9] = 0;
	out_mtx.m_[10] = p;
	out_mtx.m_[11] = -1;

	out_mtx.m_[12] = 0;
	out_mtx.m_[13] = 0;
	out_mtx.m_[14] = q;
	out_mtx.m_[15] = 0;
}