3D图形算法.doc

3D简介
　　我们首先从坐标系统开始。你也许知道在2D里我们经常使用Ren?笛卡儿坐标系统在平面上来识别点。我们使用二维(X,Y):X表示水平轴坐标,Y表示纵轴坐标。在3维坐标系，我们增加了Z，一般用它来表示深度。所以为表示三维坐标系的一个点，我们用三个参数(X,Y,Z)。这里有不同的笛卡儿三维系统可以使用。但是它们都是左手螺旋或右手螺旋的。右手螺旋是右手手指的卷曲方向指向Z轴正方向，而大拇指指向X轴正方向。左手螺旋是左手手指的卷曲方向指向Z轴负方向。实际上，我们可以在任何方向上旋转这些坐标系，而且它们仍然保持本身的特性。在计算机图形学，常用坐标系为左手坐标系，所以我们也使用它。：

X 正轴朝右
Y 正轴向上
Z 正轴指向屏幕里

矢量
什么是矢量？几句话，它是坐标集合。首先我们从二维矢量开始，(X,Y)：例如矢量P(4,5)（一般，我们用->表示矢量）。我们认为矢量P代表点(4,5)，它是从原点指向(4,5)的有方向和长度的箭头。我们谈论矢量的长度指从原点到该点的距离。二维距离计算公式是
| P | = sqrt( x^2 + y^2 )
这里有一个有趣的事实：在1D（点在单一的坐标轴上），平方根为它的绝对值。让我们讨论三维矢量：例如P(4, -5, 9)，它的长度为
| P | = sqrt( x^2 + y^2 + z^2 )
它代表在笛卡儿3D空间的一个点。或从原点到该点的一个箭头代表该矢量。在有关操作一节里，我们讨论更多的知识。

矩阵
开始，我们从简单的开始：我们使用二维矩阵4乘4矩阵，为什么是4乘4？因为我们在三维坐标系里而且我们需要附加的行和列来完成计算工作。在二维坐标系我们需要3乘3矩阵。着意味着我们在3D中有4个水平参数和4个垂直参数，一共16个。例如：
4x4单位矩阵
| 1 0 0 0 |
| 0 1 0 0 |
| 0 0 1 0 |
| 0 0 0 1 |
因为任何其它矩阵与之相乘都不改变，所以称之为单位阵。又例如有矩阵如下：
| 10 -7 22 45 |
| sin(a) cos(a) 34 32 |
| -35 28 17 6 |
| 45 -99 32 16 |

有关矢量和矩阵的操作
我们已经介绍了一些非常简单的基本概念，那么上面的知识与三维图形有什么关系呢？
本节我们介绍3D变换的基本知识和其它的一些概念。它仍然是数学知识。我们要讨论
有关矢量和矩阵操作。让我们从两个矢量和开始：

( x1 , y1 , z1 ) + ( x2 , y2 , z2 ) = ( x1 + x2 , y1 + y2 , z1 + z2 )

很简单，现在把矢量乘于系数：

k ?( x, y, z ) = ( kx, ky, kz )
可以把上面的公式称为点积，如下表示：
(x1 , y1 , z1 ) ?( x2 , y2 , z2 ) = x1x2 + y1y2 + z1z2
实际上，两个矢量的点积被它们的模的乘积除，等于两个矢量夹角的余弦。所以
cos (V ^ W) =V ?W / | V | | W |

注意^并不表示指数而是两个矢量的夹角。点积可以用来计算光线于平面的夹角，我们在计算阴影一节里会详细讨论。
现在讨论叉乘：

( x1 , y1 , z1 ) X ( x2 , y2 , z2 ) =
( y1z2 - z1y2 , z1x2 - x1z2 , x1y2 - y1x2 )

叉乘对于计算屏幕的法向量非常有用。

OK，我们已经讲完了矢量的基本概念。我们开始两个矩阵的和。它与矢量相加非常相似，这里就不讨论了。设I是矩阵的一行，J是矩阵的一列，(i,j)是矩阵的一个元素。我们讨论与3D变换有关的重要的矩阵操作原理。两个矩阵相乘，而且M x N <> N x M。例如：
A 4x4矩阵相乘公式
如果 A=(aij)4x4， B=(bij)4x4, 那么
A x B=
| S> a1jbj1 S> a1jbj2 S> a1jbj3 S> a1jbj4 |
| |
| S> a2jbj1 S> a2jbj2 S> a2jbj3 S> a2jbj4 |
| |
| S> a3jbj1 S> a3jbj2 S> a3jbj3 S> a3jbj4 |
| |
| S> a4jbj1 S> a4jbj2 S> a4jbj3 S> a4jbj4 |

其中 j=1,2,3,4

而且如果 AxB=(cik)4x4 那么我们可以在一行上写下：
cik = S>4, j=1 aijbjk
( a1, a2, a3 ) x B =
(Sum(aibi1) + b4,1, Sum(aibi2) + b4,2, Sum(aibi3) + b4,3 )

现在,我们可以试着把一些矩阵乘以单位阵来了解矩阵相乘的性质。我们把矩阵与矢量相乘结合在一起。下面有一个公式把3D矢量乘以一个4x4矩阵（得到另外一个三维矢量）如果B=(bij)4x4，那么：
( a1, a2, a3 ) x B = (S>aibi1 + b4,1, S>aibi2 + b4,2, S>aibi3 + b4,3 )>

这就是矢量和矩阵操作公式。从这里开始，代码与数学之间的联系开始清晰。

变换
我们已经见过象这样的公式：
t( tx, ty ): ( x, y ) ==> ( x + tx, y + ty )

这是在二维笛卡儿坐标系的平移等式。下面是缩放公式：

s( k ): ( x, y ) ==> ( kx, ky )

旋转等式:
r( q ): ( x, y ) ==> ( x cos(q) - y sin(q), x sin(q) + y cos(q) )
以上都是二维公式，在三维里公式的形式仍然很相近。
平移公式：
t( tx, ty, tz ): ( x, y, z ) ==> ( x + tx, y + ty, z + tz )
缩放公式：
s( k ): ( x, y, z ) ==> ( kx, ky, kz )
旋转公式（围绕Z轴）：
r( q ): ( x, y, z ) ==> ( x cos(q) - y sin(q), x sin(q) + y cos(q), z )
所以我们可以写出像在二维中同样的变换公式。我们通过乘以变换矩阵而得到新的矢量，新矢量将指向变换点。下面是所有三维变换矩阵：

平移(tx, ty, tz)的矩阵

| 1 0 0 0 |
| 0 1 0 0 |
| 0 0 1 0 |
| tx ty tz 1 |

缩放(sx, sy, sz)的矩阵
| sz 0 0 0 |
| 0 sy 0 0 |
| 0 0 sx 0 |
| 0 0 0 1 |

绕X轴旋转角q的矩阵
| 1 0 0 0 |
| 0 cos(q) sin(q) 0 |
| 0 -sin(q) cos(q) 0 |
| 0 0 0 1 |

绕Y轴旋转角q的矩阵:
| cos(q) 0 -sin(q) 0 |
| 0 1 0 0 |
| sin(q) 0 cos(q) 0 |
| 0 0 0 1 |

绕Z轴旋转角q的矩阵:
| cos(q) sin(q) 0 0 |
|-sin(q) cos(q) 0 0 |
| 0 0 1 0 |
| 0 0 0 1 |

所以我们已经可以结束关于变换的部分.通过这些矩阵我们可以对三维点进行任何变换.

平面和法向量
平面是平坦的，无限的，指向特定方向的表面可以定义平面如下：
Ax + By + Cz + D = 0

其中 A, B, C称为平面的法向量，D是平面到原点的距离。我们可以通过计算平面上的两个矢量的叉积得到平面的法向量。为得到这两个矢量，我们需要三个点。P1，P2，P3逆时针排列，可以得到：
矢量1 = P1 - P2

和

矢量2 = P3 - P2

计算法向量为：
法向量 = 矢量1 X 矢量2

把D移到等式的右边得到：
D = - (Ax + By + Cz)

或

D = - (A??1.x + B??2.y + C??3.z)>

或更简单:

D = - Normal ?P1>

但是为计算A,B,C分量。可以简化操作按如下等式：
A = y1 ( z2 - z3 ) + y2 ( z3 - z1 ) + y3 ( z1 - z2 )
B = z1 ( x2 - x3 ) + z2 ( x3 - x1 ) + z3 ( x1 - x2 )
C= x1 ( y2 - y3 ) + x2 ( y3 - y1 ) + x3 ( y1 - y2 )
D = - x1 ( y2z3 - y3z2 ) - x2 ( y3z1 - y1z3 ) - x3 ( y1z2 - y2z1 )

三维变换
存储坐标
实现矩阵系统
实现三角法系统
创建变换矩阵
如何创建透视
变换对象

存储坐标
首先可以编写星空模拟代码。那么我们基本的结构是什么样？每一个对象的描述是如何存储的？为解决这个问题，首先我们思考另一个问题：我们需要的是什么样的坐标系？最明显的答案是：
屏幕坐标系：相对于显示器的原点的2D坐标系
本地坐标系：相对于对象的原点的3D坐标系
但是我们不要忘记变换中间用到的坐标系，例如：
世界坐标系：相对于3D世界的原点三维坐标系
对齐(视点)坐标系：世界坐标系的变换，观察者的位置在世界坐标系的原点。

下面是坐标的基本结构：

// 二维坐标
typedef struct
{
short x, y;
}_2D;

//三维坐标
typedef struct
{
float x, y, z;
}_3D;

这里,我们定义了称为顶点的坐标结构。因为顶点一词指两个或两个以上菱形边的
交点。我们的顶点可以简单地认为是描述不同系统的矢量。

//不同的坐标系的坐标
typedef struct
{
_3D Local;
_3D World;
_3D Aligned;
}Vertex_t;

实现矩阵系统
我们需要存储我们的矩阵在4x4浮点数矩阵中。所以当我们需要做变换是我们定义如下矩阵：
float matrix[4][4];
然后我们定义一些函数来拷贝临时矩阵到全局矩阵:

void MAT_Copy(float source[4][4], float dest[4][4])
{
int i,j;
for(i=0; i<4; i++)
for(j=0; j<4; j++)
dest[i][j]=source[i][j];
}

很简单！现在我们来写两个矩阵相乘的函数。同时可以理解上面的一些有关矩阵相乘的公式代码如下：

void MAT_Mult(float mat1[4][4], float mat2[4][4], float dest[4][4])
{
int i,j;
for(i=0; i<4; i++)
for(j=0; j<4; j++)
dest[i][j]=mat1[i][0]*mat2[0][j]+mat1[i][1]*mat2[1][j]+mat1[i][2]*mat2[2][j]+mat1[i][3]*mat2[3][j];
}
//mat1----矩阵1
//mat2----矩阵2
//dest----相乘后的新矩阵

现在你明白了吗？现在我们设计矢量与矩阵相乘的公式。
void VEC_MultMatrix(_3D *Source,float mat[4][4],_3D *Dest)
{
Dest->x=Source->x*mat[0][0]+Source->y*mat[1][0]+Source->z*mat[2][0]+mat[3][0];
Dest->y=Source->x*mat[0][1]+Source->y*mat[1][1]+Source->z*mat[2][1]+mat[3][1];
Dest->z=Source->x*mat[0][2]+Source->y*mat[1][2]+Source->z*mat[2][2]+mat[3][2];
}
//Source-----源矢量(坐标)
//mat--------变换矩阵
//Dest-------目标矩阵(坐标)

我们已经得到了矩阵变换函数，不错吧！！
//注意,这里的矩阵变换与我们学过的矩阵变换不同
//一般的,Y=TX,T为变换矩阵，这里为Y = XT，
//由于矩阵T为4x4矩阵

实现三角法系统
几乎每一个C编译器都带有有三角函数的数学库，但是我们需要简单的三角函数时，不是每次都使用它们。正弦和余弦的计算是阶乘和除法的大量运算。为提高计算速度，我们建立自己的三角函数表。首先决定你需要的角度的个数，然后在这些地方用下面的值代替：
float SinTable[256], CosTable[256];
然后使用宏定义，它会把每一个角度变成正值，并对于大于360度的角度进行周期变换，然后返回需要的值。如果需要的角度数是2的幂次，那么我们可以使用&代替%，它使程序运行更快。例如256。所以在程序中尽量选取2的幂次。

三角法系统：
#define SIN(x) SinTable
#define COS(x) CosTable

一旦我们已经定义了需要的东西，建立初始化函数，并且在程序中调用宏。

void M3D_Init()
{
int d;
for(d=0; d<256; d++)
{
SinTable[d]=sin(d*PI/128.0);
CosTable[d]=cos(d*PI/128.0);
}
}

建立变换矩阵
下面使用C编写的变换矩阵代码

float mat1[4][4], mat2[4][4];

void MAT_Identity(float mat[4][4])
{
mat[0][0]=1; mat[0][1]=0; mat[0][2]=0; mat[0][3]=0;
mat[1][0]=0; mat[1][1]=1; mat[1][2]=0; mat[1][3]=0;
mat[2][0]=0; mat[2][1]=0; mat[2][2]=1; mat[2][3]=0;
mat[3][0]=0; mat[3][1]=0; mat[3][2]=0; mat[3][3]=1;
}
//定义单位阵

void TR_Translate(float matrix[4][4],float tx,float ty,float tz)
{
float tmat[4][4];
tmat[0][0]=1; tmat[0][1]=0; tmat[0][2]=0; tmat[0][3]=0;
tmat[1][0]=0; tmat[1][1]=1; tmat[1][2]=0; tmat[1][3]=0;
tmat[2][0]=0; tmat[2][1]=0; tmat[2][2]=1; tmat[2][3]=0;
tmat[3][0]=tx; tmat[3][1]=ty; tmat[3][2]=tz; tmat[3][3]=1;
MAT_Mult(matrix,tmat,mat1);
MAT_Copy(mat1,matrix);
}
//tx,ty.tz------平移参数
//matrix--------源矩阵和目标矩阵
//矩阵平移函数

void TR_Scale(float matrix[4][4],float sx,float sy, float sz)
{
float smat[4][4];
smat[0][0]=sx; smat[0][1]=0; smat[0][2]=0; smat[0][3]=0;
smat[1][0]=0; smat[1][1]=sy; smat[1][2]=0; smat[1][3]=0;
smat[2][0]=0; smat[2][1]=0; smat[2][2]=sz; smat[2][3]=0;
smat[3][0]=0; smat[3][1]=0; smat[3][2]=0; smat[3][3]=1;
MAT_Mult(matrix,smat,mat1);
MAT_Copy(mat1,matrix);
}
//矩阵缩放

void TR_Rotate(float matrix[4][4],int ax,int ay,int az)
{
float xmat[4][4], ymat[4][4], zmat[4][4];
xmat[0][0]=1; xmat[0][1]=0; xmat[0][2]=0;
xmat[0][3]=0;

xmat[1][0]=0; xmat[1][1]=COS(ax); xmat[1][2]=SIN(ax);
xmat[1][3]=0;

xmat[2][0]=0; xmat[2][1]=-SIN(ax); xmat[2][2]=COS(ax); xmat[2][3]=0;
xmat[3][0]=0; xmat[3][1]=0; xmat[3][2]=0; xmat[3][3]=1;

ymat[0][0]=COS(ay); ymat[0][1]=0; ymat[0][2]=-SIN(ay); ymat[0][3]=0;
ymat[1][0]=0; ymat[1][1]=1; ymat[1][2]=0; ymat[1][3]=0;
ymat[2][0]=SIN(ay); ymat[2][1]=0; ymat[2][2]=COS(ay); ymat[2][3]=0;
ymat[3][0]=0; ymat[3][1]=0; ymat[3][2]=0; ymat[3][3]=1;

zmat[0][0]=COS(az); zmat[0][1]=SIN(az); zmat[0][2]=0; zmat[0][3]=0;
zmat[1][0]=-SIN(az); zmat[1][1]=COS(az); zmat[1][2]=0; zmat[1][3]=0;
zmat[2][0]=0; zmat[2][1]=0; zmat[2][2]=1; zmat[2][3]=0;
zmat[3][0]=0; zmat[3][1]=0; zmat[3][2]=0; zmat[3][3]=1;

MAT_Mult(matrix,ymat,mat1);
MAT_Mult(mat1,xmat,mat2);
MAT_Mult(mat2,zmat,matrix);
}
//ax------绕X轴旋转的角度
//ay------绕Y轴旋转的角度
//az------绕Z轴旋转的角度
//矩阵旋转

如何建立透视
如何建立对象的立体视觉,即显示器上的一些事物看起来离我们很近,而另外一些事物离我们很远。透视问题一直是困绕我们的一个问题。有许多方法被使用。我们使用的3D世界到2D屏幕的投影公式：
P( f ):(x, y, z)==>( f*x / z + XOrigin, f*y / z + YOrigin )

其中f是焦点距离，它表示从观察者到屏幕的距离，一般在80到200厘米之间。XOrigin和YOrigin是屏幕中心的坐标，(x,y,z)在对齐坐标系上。那么投影函数应该是什么样？

#define FOCAL_DISTANCE 200
//定义焦点距离
void Project(vertex_t * Vertex)
{
if(!Vertex->Aligned.z)
Vertex->Aligned.z=1;
Vertex->Screen.x =FOCAL_DISTANCE * Vertex->Aligned.x / Vertex->Aligned.z + XOrigin;
Vertex->Screen.y =FOCAL_DISTANCE * Vertex->Aligned.y / Vertex->Aligned.z + YOrigin;
}
//得到屏幕上的投影坐标
因为0不能做除数，所以对z进行判断。

变换对象
既然我们已经掌握了所有的变换顶点的工具，就应该了解需要执行的主要步骤。
一、初始化每一个顶点的本地坐标
二、设置全局矩阵为单位阵
三、根据对象的尺寸缩放全局矩阵
四、根据对象的角度来旋转全局矩阵
五、根据对象的位置移动全局矩阵
六、把本地坐标乘以全局矩阵来得到世界坐标系
七、设置全局矩阵为单位阵
八、用观测者的位置的负值平移全局矩阵
九、用观测者的角度的负值旋转全局矩阵
十、把世界坐标系与全局矩阵相乘得到对齐坐标系
十一、投影对齐坐标系来得到屏幕坐标
即：本地坐标系-->世界坐标系-->对齐坐标系-->屏幕坐标系

多边形填充
多边形结构
发现三角形
绘制三角形

多边形结构
我们如何存储我们的多边形？首先，我们必须知道再这种状态下多边形是二维多边形，而且由于初始多边形是三维的，我们仅需要一个临时的二维多边形，所以我们能够设置二维顶点的最大数为一个常量，而没有浪费内存：

2D结构：
typedef struct
{
_2D Points[20];
int PointsCount;
int Texture;
}Polygon2D_t;

3D 结构:
typedef struct
{
int Count;
int * Vertex;
int Texture;

Vertex_t P,M,N;
}Polygon_t;

为什么顶点数组包含整数值呢?仔细思考一下,例如在立方体内,三个多边形公用同一个
顶点,所以在三个多边形里存储和变换同一个顶点会浪费内存和时间。我们更愿意存储
它们在一个对象结构里，而且在多边形结构里，我们会放置相应顶点的索引。请看
下面的结构：
typedef struct
{
int VertexCount;
int PolygonCount;
Vertex_t * Vertex;
Polygon_t * Polygon;
_3D Scaling;
_3D Position;
_3D Angle;
int NeedUpdate;
}Object_t;

发现三角形
因为绘制一个三角形比绘制任意的多边形要简单，所以我们从把多边形分割成
三顶点的形状。这种方法非常简单和直接：
void POLY_Draw(Polygon2D_t *Polygon)
{
_2D P1,P2,P3;
int i;

P1 = Polygon->Points[0];
for(i=1; i < Polygon->PointsCount-1; i++)
{
P2=Polygon->Points[i];
P3=Polygon->Points[i+1];
POLY_Triangle(P1,P2,P3,Polygon->Texture);
}
}
//上面的算法，对于凹多边形就不太适用
_____
|\ |
| \ |
|____\|

绘制三角形
现在怎样得到三角形函数?我们怎样才能画出每一条有关的直线,并且如何发现
每一行的起始和结实的x坐标。我们通过定义两个简单有用的宏定义开始来区别
垂直地两个点和两个数：

#define MIN(a,b) ((a<b)?(a):(b))
#define MAX(a,b) ((a>b)?(a):(b))
#define MaxPoint(a,b) ((a.y > b.y) ? a : b)
#define MinPoint(a,b) ((b.y > a.y) ? a : b)

然后我们定义三个宏来区别三个点：

#define MaxPoint3(a,b,c) MaxPoint(MaxPoint(a,b),MaxPoint(b,c))
#define MidPoint3(a,b,c) MaxPoint(MinPoint(a,b),MinPoint(a,c))
#define MinPoint3(a,b,c) MinPoint(MinPoint(a,b),MinPoint(b,c))

你也许注意到MidPoint3宏不总是正常地工作，取决于三个点排列的顺序，
例如,a<b & a<c 那么由MidPoint3得到的是a,但它不是中间点。
我们用if语句来修正这个缺点，下面为函数的代码：

void POLY_Triangle(_2D p1,_2D p2,_2D p3,char c)
{
_2D p1d,p2d,p3d;
int xd1,yd1,xd2,yd2,i;
int Lx,Rx;

首先我们把三个点进行排序：
p1d = MinPoint3(p1,p2,p3);
p2d = MidPoint3(p2,p3,p1);
p3d = MaxPoint3(p3,p1,p2);

当调用这些宏的时候为什么会有点的顺序的改变？（作者也不清楚）可能这些点被逆时针传递。试图改变这些宏你的屏幕显示的是垃圾！现在我们并不确定中间的点，所以我们做一些检查，
而且在这种状态下，得到的中间点有似乎是错误的，所以我们修正：

if(p2.y < p1.y)
{
p1d=MinPoint3(p2,p1,p3);
p2d=MidPoint3(p1,p3,p2);
}
这些点的排列顺序看起来很奇怪，但是试图改变他们那么所有的东西就乱套了。只有理解或
接受这些结论。现在我们计算增量

xd1=p2d.x-p1d.x;
yd1=p2d.y-p1d.y;
xd2=p3d.x-p1d.x;
yd2=p3d.y-p1d.y;

OK，第一步已经完成，如果有增量 y：
if(yd1)
for(i=p1d.y; i<=p2d.y; i++)
{

我们用x的起始坐标计算x值，在当前点和起始点之间加上增量 y，乘以斜率( x / y )
的相反值。