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];
我們得到了三維明暗。