3D图像算法

http://dev.gameres.com/Program/Visual/3D/3Darit.htm
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[ABS((int)x&255)]
#define COS(x) CosTable[ABS((int)x&255)]

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

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 )
的相反值。
        Lx = p1d.x + ((i - p1d.y) * xd1) / yd1;
        Rx = p1d.x + ((i - p1d.y) * xd2) / yd2;
如果不在同一个点，绘制线段，按次序传递这两个点：

        if(Lx!=Rx)
            VID_HLine(MIN(Lx,Rx),MAX(Lx,Rx),i,c);
    } 
现在我们重新计算第一个增量，而且计算第二条边
    xd1=p3d.x-p2d.x;
    yd1=p3d.y-p2d.y;

    if(yd1)
    for(i = p2d.y; i <= p3d.y; i++)
    {
        Lx = p1d.x + ((i - p1d.y) * xd2) / yd2;
        Rx = p2d.x + ((i - p2d.y) * xd1) / yd1;
        if(Lx!=Rx)
            VID_HLine(MIN(Lx,Rx),MAX(Lx,Rx),i,c);
    }
}

以上我们已经得到多边形填充公式，对于平面填充更加简单：
void VID_HLine(int x1, int x2, int y, char c)
{
    int x;
    for(x=x1; x<=x2; x++)
        putpixel(x, y, c);
}

Sutherland-Hodgman剪贴
概述
Z-剪贴
屏幕剪贴 

概述
一般地，我们更愿意剪贴我们的多边形。必须靠着屏幕的边缘剪贴，但也必须在观察的前方（我们不需要绘制观察者后面的事物，当z左边非常小时）。当我们剪贴一个多边形，并不考虑是否每一个点在限制以内，而我们更愿意增加必须的顶点，所以我们需要一个第三个多边形结构：
typedef struct
{
    int Count;
    _3D Vertex[20];
}CPolygon_t;

由于我们有附加的顶点来投影，我们不再投影顶点，而是投影剪贴的3D多边形到
2D多边形。
void M3D_Project(CPolygon_t *Polygon,Polygon2D_t *Clipped,int focaldistance)
{
    int v;
    for(v=0; v<Polygon->Count; v++)
    {
        if(!Polygon->Vertex[v].z)Polygon->Vertex[v].z++;
        Clipped->Points[v].x=Polygon->Vertex[v].x*focaldistance/Polygon->Vertex[v].z+Origin.x;
        Clipped->Points[v].y=Polygon->Vertex[v].y*focaldistance/Polygon->Vertex[v].z+Origin.y;
    }
    Clipped->PointsCount=Polygon->Count;
}

Z-剪贴
首先我们定义计算在第一个点和第二个点之间以及在第一个点和最小z值的z增量的宏。
然后，我们计算比例，注意不要被零除。
WORD ZMin=20;
#define INIT_ZDELTAS dold=V2.z-V1.z; dnew=ZMin-V1.z;
#define INIT_ZCLIP INIT_ZDELTAS if(dold) m=dnew/dold;

我们建立一个函数，它主要剪贴多边形指针的参数（它将记下作为结果的剪贴的顶点），第一个顶点（我们剪贴的边的开始）和第二个顶点（最后）：
void CLIP_Front(CPolygon_t *Polygon,_3D V1,_3D V2)
{
    float dold,dnew, m=1;
    INIT_ZCLIP

现在我们必须检测边缘是否完全地在视口里，离开或进入视口。如果边缘没有完全地
在视口里，我们计算视口与边缘的交线，用m值表示，用INIT_ZCLIP计算。

如果边缘在视口里： 
if ( (V1.z>=ZMin) && (V2.z>=ZMin) )
    Polygon->Vertex[Polygon->Count++]=V2;

如果边缘正离开视口：
if ( (V1.z>=ZMin) && (V2.z<ZMin) )
{
    Polygon->Vertex[Polygon->Count ].x=V1.x + (V2.x-V1.x)*m;
    Polygon->Vertex[Polygon->Count ].y=V1.y + (V2.y-V1.y)*m;
    Polygon->Vertex[Polygon->Count++ ].z=ZMin;
}

如果边缘正进入视口：
if ( (V1.z<ZMin) && (V2.z>=ZMin) )
{
    Polygon->Vertex[Polygon->Count ].x=V1.x + (V2.x-V1.x)*m;
    Polygon->Vertex[Polygon->Count ].y=V1.y + (V2.y-V1.y)*m;
    Polygon->Vertex[Polygon->Count++ ].z=ZMin;
    Polygon->Vertex[Polygon->Count++ ]=V2;
}
这就是边缘Z-剪贴函数
}

现在我们可以写下完整的多边形Z-剪贴程序。为了有代表性，定义一个宏用来
在一个对象结构中寻找适当的多边形顶点。
#define Vert(x) Object->Vertex[Polygon->Vertex[x]]

下面是它的函数：
void CLIP_Z(Polygon_t *Polygon,Object_t *Object,CPolygon_t *ZClipped)
{
    int d,v;
    ZClipped->Count=0;
    for (v=0; v<Polygon->Count; v++)
    {
        d=v+1;
        if(d==Polygon->Count)d=0;
            CLIP_Front(ZClipped, Vert(v).Aligned,Vert(d).Aligned);
    }
}

这个函数相当简单：它仅仅调用FrontClip函数来做顶点交换。

剪贴屏幕 
剪贴屏幕的边缘同Z-剪贴相同，但是我们有四个边缘而不是一个。所以我们需要四个
不同的函数。但是它们需要同样的增量初始化：
#define INIT_DELTAS dx=V2.x-V1.x; dy=V2.y-V1.y;
#define INIT_CLIP INIT_DELTAS if(dx)m=dy/dx;

边缘是：
_2D TopLeft, DownRight;

为了进一步简化_2D和 _3D结构的使用，我们定义两个有用的函数：
_2D P2D(short x, short y)
{
    _2D Temp;
    Temp.x=x;
    Temp.y=y;
    return Temp;
}
_3D P3D(float x,float y,float z)
{
    _3D Temp;
    Temp.x=x;
    Temp.y=y;
    Temp.z=z;
    return Temp;
}

然后使用这两个函数来指定视口：
TopLeft=P2D(0, 0);
DownRight=P2D(319, 199);

下面是四个边缘剪贴函数：
/*
=======================
剪贴左边缘
=======================
*/
void CLIP_Left(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
    float dx,dy, m=1;
    INIT_CLIP

    // ************OK************
    if ( (V1.x>=TopLeft.x) && (V2.x>=TopLeft.x) )
        Polygon->Points[Polygon->PointsCount++]=V2;
    // *********LEAVING**********
    if ( (V1.x>=TopLeft.x) && (V2.x<TopLeft.x) )
    {
        Polygon->Points[Polygon->PointsCount].x=TopLeft.x;
        Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(TopLeft.x-V1.x);
    }
    // ********ENTERING*********
    if ( (V1.x<TopLeft.x) && (V2.x>=TopLeft.x) )
    {
        Polygon->Points[Polygon->PointsCount].x=TopLeft.x;
        Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(TopLeft.x-V1.x);
        Polygon->Points[Polygon->PointsCount++]=V2;
    }
}
/*
=======================
剪贴右边缘
=======================
*/

void CLIP_Right(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
    float dx,dy, m=1;
    INIT_CLIP
    // ************OK************
    if ( (V1.x<=DownRight.x) && (V2.x<=DownRight.x) )
        Polygon->Points[Polygon->PointsCount++]=V2;
    // *********LEAVING**********
    if ( (V1.x<=DownRight.x) && (V2.x>DownRight.x) )
    {
        Polygon->Points[Polygon->PointsCount].x=DownRight.x;
        Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(DownRight.x-V1.x);
    }
    // ********ENTERING*********
    if ( (V1.x>DownRight.x) && (V2.x<=DownRight.x) )
    {
        Polygon->Points[Polygon->PointsCount].x=DownRight.x;
        Polygon->Points[Polygon->PointsCount++].y=V1.y+m*(DownRight.x-V1.x);
        Polygon->Points[Polygon->PointsCount++]=V2;
    }
}
/*
=======================
剪贴上边缘
=======================
*/
void CLIP_Top(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
    float dx,dy, m=1;
    INIT_CLIP
    // ************OK************
    if ( (V1.y>=TopLeft.y) && (V2.y>=TopLeft.y) )
        Polygon->Points[Polygon->PointsCount++]=V2;
    // *********LEAVING**********
    if ( (V1.y>=TopLeft.y) && (V2.y<TopLeft.y) )
    {
        if(dx)
            Polygon->Points[Polygon->PointsCount].x=V1.x+(TopLeft.y-V1.y)/m;
        else
            Polygon->Points[Polygon->PointsCount].x=V1.x;
        Polygon->Points[Polygon->PointsCount++].y=TopLeft.y;
    }
    // ********ENTERING*********
    if ( (V1.y<TopLeft.y) && (V2.y>=TopLeft.y) )
    {
        if(dx)
            Polygon->Points[Polygon->PointsCount].x=V1.x+(TopLeft.y-V1.y)/m;
        else
            Polygon->Points[Polygon->PointsCount].x=V1.x;
        Polygon->Points[Polygon->PointsCount++].y=TopLeft.y;
        Polygon->Points[Polygon->PointsCount++]=V2;
    }
}

/*
=======================
剪贴下边缘
=======================
*/

void CLIP_Bottom(Polygon2D_t *Polygon,_2D V1,_2D V2)
{
    float dx,dy, m=1;
    INIT_CLIP
    // ************OK************
    if ( (V1.y<=DownRight.y) && (V2.y<=DownRight.y) )
        Polygon->Points[Polygon->PointsCount++]=V2;
    // *********LEAVING**********
    if ( (V1.y<=DownRight.y) && (V2.y>DownRight.y) )
    {
        if(dx)
            Polygon->Points[Polygon->PointsCount].x=V1.x+(DownRight.y-V1.y)/m;
        else
            Polygon->Points[Polygon->PointsCount].x=V1.x;
        Polygon->Points[Polygon->PointsCount++].y=DownRight.y;
    }
    // ********ENTERING*********
    if ( (V1.y>DownRight.y) && (V2.y<=DownRight.y) )
    {
        if(dx)
            Polygon->Points[Polygon->PointsCount].x=V1.x+(DownRight.y-V1.y)/m;
        else
            Polygon->Points[Polygon->PointsCount].x=V1.x;
        Polygon->Points[Polygon->PointsCount++].y=DownRight.y;
        Polygon->Points[Polygon->PointsCount++]=V2;
    }
}

为了得到完整的多边形剪贴函数，我们需要定义一个附加的全局变量：
polygon2D_t TmpPoly;

void CLIP_Polygon(Polygon2D_t *Polygon,Polygon2D_t *Clipped)
{
    int v,d;

    Clipped->PointsCount=0;
    TmpPoly.PointsCount=0;

    for (v=0; v<Polygon->PointsCount; v++)
    {
        d=v+1;
        if(d==Polygon->PointsCount)d=0;
        CLIP_Left(&TmpPoly, Polygon->Points[v],Polygon->Points[d]);
    }
    for (v=0; v<TmpPoly.PointsCount; v++)
    {
        d=v+1;
        if(d==TmpPoly.PointsCount)d=0;
        CLIP_Right(Clipped, TmpPoly.Points[v],TmpPoly.Points[d]);
    }
    TmpPoly.PointsCount=0;
    for (v=0; v<Clipped->PointsCount; v++)
    {
        d=v+1;
        if(d==Clipped->PointsCount)d=0;
        CLIP_Top(&TmpPoly, Clipped->Points[v],Clipped->Points[d]);
    }
    Clipped->PointsCount=0;
    for (v=0; v<TmpPoly.PointsCount; v++)
    {
        d=v+1;
        if(d==TmpPoly.PointsCount)d=0;
        CLIP_Bottom(Clipped, TmpPoly.Points[v],TmpPoly.Points[d]);
    }
}

程序原理同Z-剪贴一样，所以我们可以轻松地领会它。


隐面消除
Dilemna
底面消除
Z-缓冲 
The Dilemna
三维引擎的核心是它的HSR系统，所以我们必须考虑选择那一种。一般来说，最流行
的几种算法是：
画笔算法
需要的时间增长更快
难以实现（尤其重叠测试）
不能准确地排序复杂的场景
字节空间分区树
特别快
难以实现
仅仅能排序静态多边形
需要存储树
Z-缓存
需要的时间随多边形的数目线性地增加
在多边形大于5000后速度比画笔算法快
能够完美地渲染任何场景，即使逻辑上不正确
非常容易实现
简单
需要大量的内存
很慢
所以我们的选择是Z-缓存。当然也可以选择其他算法。

底面消除
除了这些方法，我们可以很容易地消除多边形的背面来节省大量的计算时间。首先
我们定义一些有用的函数来计算平面和法向量以及填充。然后，我们给这个函数增加
纹理和阴影计算。这些变量为全局变量：
float A,B,C,D;
BOOL backface;
下面是我们的引擎函数，每一个坐标都是浮点变量：
void ENG3D_SetPlane(Polygon_t *Polygon,Object_t *Object)
{
    float x1=Vert(0).Aligned.x;
    float x2=Vert(1).Aligned.x;
    float x3=Vert(2).Aligned.x;
    float y1=Vert(0).Aligned.y;
    float y2=Vert(1).Aligned.y;
    float y3=Vert(2).Aligned.y;
    float z1=Vert(0).Aligned.z;
    float z2=Vert(1).Aligned.z;
    float z3=Vert(2).Aligned.z;


然后我们计算平面等式的每一个成员：
    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*(y2*z3-y3*z2)-x2*(y3*z1-y1*z3)-x3*(y1*z2-y2*z1);

再检查是否它面朝我们或背朝：
    backface=D<0;
}


Z-缓存
Z-缓存是把显示在屏幕上的每一个点的z坐标保持在一个巨大的数组中，并且当我们
我们检查是否它靠近观察者或是否在观察者后面。我们仅仅在第一种情况下绘制它。所以我们不得不计算每一个点的z值。但是首先，我们定义全局树组和为他分配空间。
（内存等于追至方向与水平方向的乘积）：
typedef long ZBUFTYPE;
ZBUFTYPE *zbuffer;
zbuffer=(ZBUFTYPE *)malloc(sizeof(ZBUFTYPE)*MEMORYSIZE);
我们使用长整形作为z-缓存类型，因为我们要使用定点数。我们必须记住设置每一个z坐标来尽可能得到更快的速度：
int c;
for(c=0; c<MEMORYSIZE; c++)
    zbuffer[c]=-32767;

下面是数学公式。如何才能发现z坐标？我们仅仅已经定义的顶点，而不是多边形的
每一个点。实际上，我们所需要做的是投影的反变换，投影公式是：
u = f ?x / z

和

v = f ?y / z


其中u是屏幕上x的坐标，最小值为XOrigin，v是屏幕上的y的坐标，最小值YOrigin。
平面公式是：

Ax + By + Cz + D = 0

一旦我们已经得到分离的x和y，有：
x = uz / f

和

y = vz / f


如果我们在平面等式中替代变量，公式变为：

A(uz / f) + B(vz / f) + Cz = -D 

我们可以提取z分量：

z(A(u / f) + B(v / f) + C) = -D

所以我们得到z：
z = -D / (A(u / f) + B(v / f) + C)

但是由于对于每一个像素我们需要执行以上的除法，而计算1/z将提高程序的速度：
1 / z = -(A(u / f) + B(v / f) +C) / D

1 / z = -(A / (fD))u - (B / (fD))v - C / D

所以在一次像素程序运行的开始：

1 / z = -(A / (fD))u1 - (B / (fD))v - C / D

对于每一个像素，增量为：
-(A / (fD))>


下面是程序：
#define FIXMUL (1<<20)

int offset=y*MODEINFO.XResolution+x1;
int i=x1-Origin.x, j=y-Origin.y;
float z_,dz_;
ZBUFTYPE z,dz;

//初始化 1/z 值 (z: 1/z)
dz_=((A/(float)Focal_Distance)/-D);
z_=((dz_*i)+( (B*j/(float)Focal_Distance) + C) /-D);
dz=dz_*FIXMUL;
z=z_*FIXMUL;

然后，对于每一个像素，我们简单的计算：
if(z>ZBuffer[offset])
{
    zbuffer[offset]=z;
    SCREENBUFFER[offset]=color;
}
z+=dz;


3D纹理映射
概述
魔幻数字
纹理映射的透视修正

概述
在做纹理映射时首先考虑的是建立纹理数组和初始化3D纹理坐标。纹理将存储在：
#define MAXTEXTURES 16
bitmap_t Textures[MAXTEXTURES];
我们从PCX文件分配和加载纹理。这里假设纹理大小为64x64。我们使用polygon_t
结构的纹理坐标：
vertex_t P,M,N;

我们在函数中初始化纹理，该函数在建立多边形后被调用。P是纹理的原点，M是
纹理的水平线末端，N是垂直线的末端。
void TEX_Setup(Polygon_t * Polygon, Object_t *Object)
{
    Polygon->P.Local=P3D(Vert(1).Local.x,Vert(1).Local.y,Vert(1).Local.z);
    Polygon->M.Local=P3D(Vert(0).Local.x,Vert(0).Local.y,Vert(0).Local.z);
    Polygon->N.Local=P3D(Vert(2).Local.x,Vert(2).Local.y,Vert(2).Local.z);
}

我们需要象任何其他对象的顶点一样变换纹理坐标，所以我们需要建立世界变换和
一个对齐变换函数：
void TR_Object(Object_t *Object, float matrix[4][4])
{
    int v,p;
    for(v=0; v<Object->VertexCount; v++)
        VEC_MultMatrix(&Object->Vertex[v].Local,matrix,&Object->Vertex[v].World);
    for(p=0; p<Object->PolygonCount; p++)
    {
        VEC_MultMatrix(&Object->Polygon[p].P.Local,matrix,&Object->Polygon[p].P.World);
        VEC_MultMatrix(&Object->Polygon[p].M.Local,matrix,&Object->Polygon[p].M.World);
        VEC_MultMatrix(&Object->Polygon[p].N.Local,matrix,&Object->Polygon[p].N.World);
    }
}

void TR_AlignObject(Object_t *Object, float matrix[4][4])
{
    int v,p;
    for(v=0; v<Object->VertexCount; v++)
        VEC_MultMatrix(&Object->Vertex[v].World,matrix,&Object->Vertex[v].Aligned);
    for(p=0; p<Object->PolygonCount; p++)
    {
        VEC_MultMatrix(&Object->Polygon[p].P.World,matrix,&Object->Polygon[p].P.Aligned);
        VEC_MultMatrix(&Object->Polygon[p].M.World,matrix,&Object->Polygon[p].M.Aligned);
        VEC_MultMatrix(&Object->Polygon[p].N.World,matrix,&Object->Polygon[p].N.Aligned);
    }
}

魔幻数

既然我们已经得到了变幻的纹理坐标，我们的目标是发现在纹理位图上的像素
水平和垂直的坐标在屏幕上如何绘制。纹理坐标称为u和v。下面的公式给出坐标：
u = a * TEXTURE_SIZE / c

和 
v = b * TEXTURE_SIZE / c


a，b，c满足下面的等式：
a = Ox + Vx j + Hx i

b = Oy + Vy j + Hy i

c = Oz + Vz j + Hz i


其中O,H,V数是魔幻数。它根据下面的公式由纹理坐标计算得到：
Ox = NxPy - NyPx
Hx = NyPz - NzPy
Vx = NzPx - NxPz
Oy = MxPy - MyPx
Hy = MyPz - MzPy
Vy = MzPx - MxPz
Oz = MyNx - MxNy
Hz = MzNy - MyNz
Vz = MxNz - MzNx

这里，我们不解释魔幻数的原因。它看起来像奇怪的叉积。


纹理映射透视修正
O,H,V数的计算需要一些修正，所以我们增加下面到ENG3D_SetPlane:
//用于修正当数字变得太大时候的错误
#define FIX_FACTOR 1/640

//初始化纹理矢量
P=Polygon->P.Aligned;
M=VEC_Difference(Polygon->M.Aligned,Polygon->P.Aligned);
N=VEC_Difference(Polygon->N.Aligned,Polygon->P.Aligned);

P.x*=Focal_Distance;
P.y*=Focal_Distance;

M.x*=Focal_Distance;
M.y*=Focal_Distance;

N.x*=Focal_Distance;
N.y*=Focal_Distance;

下面是VEC_Difference的实现：
_3D VEC_Difference(_3D Vector1, _3D Vector2)
{
return P3D(Vector1.x-Vector2.x,Vector1.y-Vector2.y,Vector1.z-Vector2.z);
}

然后计算魔幻数。
_3D O, H, V;
ENG3D_SetPlane: 
H.x=(N.y*P.z-N.z*P.y)*FIX_FACTOR;
V.x=(N.z*P.x-N.x*P.z)*FIX_FACTOR;
O.x=(N.x*P.y-N.y*P.x)*FIX_FACTOR;

H.z=(M.z*N.y-M.y*N.z)*FIX_FACTOR;
V.z=(M.x*N.z-M.z*N.x)*FIX_FACTOR;
O.z=(M.y*N.x-M.x*N.y)*FIX_FACTOR;

H.y=(M.y*P.z-M.z*P.y)*FIX_FACTOR;
V.y=(M.z*P.x-M.x*P.z)*FIX_FACTOR;
O.y=(M.x*P.y-M.y*P.x)*FIX_FACTOR; 
下面为TEX_HLine改变VID_HLine以便使用纹理映射（难以理解的部分），首先我们
必须初始化魔幻数坐标：
a=-((long)O.x+((long)V.x*(long)j)+((long)H.x*(long)i))*64L;
b= ((long)O.y+((long)V.y*(long)j)+((long)H.y*(long)i))*64L;
c= ((long)O.z+((long)V.z*(long)j)+((long)H.z*(long)i));
long Hx,Hy,Hz;
int u,v;
BYTE color=0;
BYTE *mapping=Textures[texture].Picture;
然后把H.x 和 H.y乘以64，这样我们不需要为每一个像素做计算。我们用长整形数代替
浮点数：
Hx=H.x*-64;
Hy=H.y*64;
Hz=H.z;

对于每一个像素，改变最后一个参数并且替代绘制原来的参数：
if(c)
{
    u=a/c;
    v=b/c;
    color=mapping[((v&63)<<6)+(u&63)];
    if(color)
    {
        zbuffer[offset]=z;
        SCREENBUFFER[offset]=LightTable[light][color];
    }
}
a+=Hx;
b+=Hy;
c+=Hz;

现在我们得到自己的纹理映射！




三维明暗
计算法向量
计算叉乘
使用光线表

计算法向量
在3D数学教程里我们已经深入讨论了矢量和法向量，这里我们给出一些实现：
float VEC_DotProduct(_3D Vector1, _3D Vector2)
{
    return (Vector1.x*Vector2.x+Vector1.y*Vector2.y+Vector1.z*Vector2.z);
}

_3D VEC_CrossProduct(_3D Vector1, _3D Vector2)
{
    return P3D(Vector1.y*Vector2.z-Vector1.z*Vector2.y,Vector1.z*Vector2.x-Vector1.x*Vector2.z,Vector1.x*Vector2.y-Vector1.y*Vector2.x);
}
void VEC_Normalize(_3D * Vector)
{
    float m=sqrt(Vector->x*Vector->x+Vector->y*Vector->y+Vector->z*Vector->z);
    Vector->x/=m;
    Vector->y/=m;
    Vector->z/=m;
}

对于3D明暗，需要向ENG3D_SetPlane增加：
//计算平面的法向量
PlaneNormal=VEC_CrossProduct(P3D(x2-x1,y2-y1,z2-z1),P3D(x3-x1,y3-y1,z3-z1));
VEC_Normalize(&PlaneNormal);

计算叉积
正如在数学部分所说的，两个矢量间的夹角等于他们的点积（当他们归一化后）。
为发现到达平面的光线，我们简单地把环境光和归一化的光源的世界坐标系的叉积
与最大的光强度的乘积相加。下面是代码：

全局变量定义：
WORD AmbientLight=20;
#define MAXLIGHT 32
static Vertex_t LightSource;
WORD light;

在SetPlane函数里：
//计算法向量和光源的点积的强度
light=MAXLIGHT*VEC_DotProduct(PlaneNormal,LightSource.World)+AmbientLight;
if(light>MAXLIGHT)light=MAXLIGHT;
if(light<1)light=1;


明显地，我们需要初始化光源，正如计算法向量顶点。
//初始化光源
LightSource.Local=P3D(0,0,0);
MAT_Identity(matrix);
TR_Translate(matrix, 10,10,10);
TR_Rotate(matrix, 0,128-32,128);
VEC_MultMatrix(&LightSource.Local,matrix,&LightSource.World);
VEC_Normalize(&LightSource.World);

使用光线表

光线表是在基于调色板技术上使用光线强度的一种方法。在每一个强度上可以发现
最好的颜色匹配。Christopher Lampton在他的幻想花园一书中设计了大量的光源表
发生器。不幸的是，他使用的是自己的格式，所以我们可以选择自己的或他人的光线表。一旦我们已经有了光线表，一旦我们有了广西表，就可以在全局数组中加载它。
BYTE LightTable[MAXLIGHT+1][256];
一旦已经加载，我们在TEX_HLine函数中改变如下：
screen[offset]=LightTable[light][color];
我们得到了三维明暗。