Graphics - 演算法筆記

2D Graphics

Graphics

Graphics是繪畫、製圖的意思。左圖是2D繪圖、右圖是3D繪圖，相信大家一眼就能看出差異：

Image的演算法著重於讓既有圖片產生變化，Graphics的演算法專注於從無到有產生圖片。兩者相輔相成、密不可分。

2D Graphics

先前已經介紹過「像素Pixel」。2D繪圖就是逐步設定每個像素的RGB值，呈現我們想要的圖案。

以人工方式，一點一點設定每個像素的RGB值，稱作Pixel Art，常常出現於電玩遊戲及牆壁磁磚。

以演算法，一點一點設定每個像素的RGB值，則是這裡要談的重點。

2D Graphics相關資源

C與C++沒有2D繪圖的函式庫。更高階的程式語言都有內建2D繪圖函式庫，例如C#的Graphics、Java的Graphics、Qt的QPaint。作業系統也有內建2D繪圖函式庫，例如Windows的GDI+。瀏覽器也有內建2D繪圖函式庫，例如HTML5 Canvas。

知名的2D繪圖軟體，諸如Microsoft Paint、AutoCAD、Microsoft Visio、Adobe Illustrator等等。

Rendering

Rendering是抹牆、展現的意思。計算學當中則是「填充像素產生圖片」的意思。中文翻譯千奇百怪，繁中「彩現」「算繪」，簡中「渲染」「绘制」。

Shape Rendering

前言

想要畫一張複雜的圖，無論是藝術家還是數學家，都是從簡單的幾何形狀：直線、方、圓開始。本篇文章將介紹如何繪製簡單的幾何形狀。

Shape Rendering的演算法已經在顯示卡上燒成電路、並且包裝成函式庫。前人種樹後人乘涼，我們大可不必痛下苦功編寫程式碼。然而本篇文章不談函式庫。

Shape Rendering擁有相關的檔案格式SVG。然而本篇文章不談資料結構。

操控鍵盤滑鼠即時繪製圖形，我們需要了解監聽和處理事件的機制。然而本篇文章不談事件處理。

繪製直線、曲線、方形、圓形等等圖形，管理鍵盤、滑鼠、視窗等等介面，我們需要規劃編排大量程式碼。然而本篇文章不談物件導向。

Shape Rendering的軟體已經發展成熟。然而本篇文章不談如何製作軟體，也不談如何使用軟體。

像素

直線

給定線段的兩個端點，如何找到線段對應的像素呢？

腦筋靈活的讀者肯定馬上聯想到一次內插。不過像素看起來零零散散的，根本沒有連成一線。

DDA Algorithm則是計算兩個端點在X軸、Y軸上的差距，取較大者作為像素數量，使得像素互相碰觸。另一方面，從逐點內插，改為預先計算斜率、逐步累加，大幅降低了運算量。

Bresenham's Algorithm搞定了浮點數誤差。

Wu's Algorithm加入了抗鋸齒效果。

曲線

大家習慣使用Bézier Curve。

遞推：窮舉各種參數值。相鄰控制點，反覆取加權平均數，直到得出一點。

分治法：貝茲曲線從中切開，分頭遞迴下去，直到兩端點位於相同像素。相鄰控制點，反覆取中點，以得到切開地點。

外插：以泰勒近似得到外插公式。計算速度極快，誤差稍高。

解方程式：窮舉水平線，求得水平線與曲線的交點。

還有加粗的貝茲曲線演算法、抗鋸齒的貝茲曲線演算法。

手繪線

其實就是一大堆的滑鼠座標位置，通常會斷斷續續不相連。

可視作Polyline的頂點。然後套用減少頂點的演算法，甚至以Bézier Curve作為模型。一方面平滑線條、一方面節省記憶體。

方

畫四條直線。或者視作簡單多邊形！

圓

Bresenham's Algorithm。可以推廣為橢圓和拋物線。

簡單多邊形（窮舉法）

分為空心和實心。空心很簡單，多邊形每條邊分別畫直線。

實心複雜多了。填充一個多邊形裡面所有的像素，讓整塊多邊形擁有顏色。多邊形以頂點座標的形式儲存，頂點座標是浮點數。像素座標則是整數。

最簡單的方式是試誤法：窮舉畫面上所有像素，運用「判斷點在多邊形內部」的演算法，若在內部則填色。

時間複雜度O(XYN)。X與Y是畫面的長與寬，N是簡單多邊形的頂點數目。

簡單多邊形（Flood Fill Algorithm）

還有沒有更好的方法呢？我們可以把問題分成兩個步驟：一、先鉤勒多邊形的邊界；二、再填充多邊形的內部。一旦有了邊界，就不必費心判斷像素是否在多邊形內。

鉤勒的部分，運用畫直線演算法即可。填墨的部分，讀者應該馬上就聯想到「Flood Fill Algorithm」，遞迴往四方向填墨。

鉤勒需時O(N+B)，填墨需時O(B+4A)，總時間複雜度O(N+B+A)。N是多邊形的頂點數目，B是多邊形邊界的像素數目，A是多邊形內部的像素數目。

凸多邊形（Scanline Fill Algorithm）

一、求出凸多邊形的最低像素和最高像素，把凸多邊形的邊，分為左邊界和右邊界。二、建立兩條陣列，陣列大小等同於凸多邊形的垂直高低差，用來儲存垂直方向各個像素的左邊界與右邊界。三、依序窮舉凸多邊形的邊，把每條邊碰到的像素位置算出來，作為邊界值，存到陣列。四、水平方向，一條一條填滿多邊形。

鉤勒需時O(N+2B)，填墨需時O(2B+A)，總時間複雜度O(N+B+A)。

雖然表面上Flood Fill Algorithm與Scanline Fill Algorithm的時間複雜度一模一樣，但是實際上Scanline Fill Algorithm的速度比較快。

const int X = 1024, Y = 768;
float xmin[Y], xmax[Y];

void fillEdge(Point p1, Point p2)
{
	int ymini = max(0,   (int)ceil (min(p1.y, p2.y)));
	int ymaxi = min(Y-1, (int)floor(max(p1.y, p2.y)));
	for (int y = ymini; y <= ymaxi; ++y)
	{
		float x = (p2.y == p1.y ? p1.x : (y - p1.y) / (p2.y - p1.y) * (x - p1.x) + p1.x);
		xmin[y] = min(xmin[y], x);	// 更新左邊界
		xmax[y] = max(xmax[y], x);	// 更新右邊界
	}
}

void fillPolygon(Polygon p)
{
	float ymin = +1e9, ymax = -1e9;
	for (int i=0; i<p.size(); ++i)
	{
		ymin = min(ymin, p[i].y);
		ymax = max(ymax, p[i].y);
	}

int ymini = max(0,   (int)ceil (ymin));
	int ymaxi = min(Y-1, (int)floor(ymax));
	for (int y = ymini; y <= ymaxi; ++y)
	{
		xmin[y] = +1e9;
		xmax[y] = -1e9;
	}

for (int i=0; i<p.size(); ++i)
		fillEdge(p[i], p[i+1]);
	fillEdge(p.back(), p[0]);

initImage();
	for (int y = ymini; y <= ymaxi; ++y)
	{
		int xmini = max(0,   (int)ceil (xmin[y]));
		int xmaxi = min(X-1, (int)floor(xmax[y]));
		for (int x = xmini; x <= xmaxi; ++x)
			drawPixel(x, y);
	}
}

UVa 143 356

簡單多邊形（Scanline Fill Algorithm）

Font Rendering

Character

Character就是電腦符號。譯作「字元」或「字符」。

電腦符號都有編號。世界公定的編號方式是Unicode。例如第30399號是「算」，第9775號是「☯」。電腦符號包括了世界各國文字、數字、圖示、表情符號等等。

電腦符號必須編碼，節省儲存空間。Unicode有著UTF-8和UTF-16兩種主流的編碼方式。

字元編碼，是編碼理論的範疇。字元操作，是字串理論的範疇。字元顯示，是計算機繪圖的範疇，是我們現在要談的內容。

Font

Font就是文字的造型。譯作「字體」或「字型」。

從書寫的觀點：中文字體有楷書、行書、草書、隸書、篆書五大類，有名的楷書如顏真卿、柳公權。字體是書法家發明的。

從電腦的觀點：繁體中文字型有標楷體、細明體、微軟正黑體。英文字型有Times New Roman、Helvetica、Palatino。字型是美術設計師和程式設計師發明的。

只要有文字的地方，通常都能更換字型，例如文字編輯器、瀏覽器、作業系統、命令提示字元。電腦可供使用的字型，主要是看你的電腦安裝了哪些字型。有專門製作字型、販售字型的公司，例如專精繁體中文字型的JustFont、文鼎、華康。亦有人號召製作字型，例如文泉驛。

電腦顯示字元，必須透過字型。每一套字型都只替一部分的字元設計造型，例如Times New Roman只設計了英文字母、拉丁字母、少量的特殊符號，例如微軟正黑體只設計了英文字母、中文字、少量的特殊符號。一套字型當中，沒有設計造型的字元，就無法顯示，或者是顯示問號、空白方格等等莫名其妙的符號。儘管現今沒有任何一套字型涵蓋Unicode規定的所有字元，不過其實也不需要一套萬能的字型。對於沒有設計造型的字元，作業系統將另尋其他字型。

字型由作業系統控管。想要讓視窗程式顯示特定字型，土法煉鋼的方式是Windows API，輕鬆寫意的方式是程式語言內建的函式庫。想要讓網頁顯示特定字型，利用CSS的font-family即可。

Glyph

接著到了本篇文章的主角「字形」。首先玩個小遊戲吧！拖曳控制點，按Enter評分。美感和眼力必須很傑出。設計師的日常。

「字形」就是單獨一個字元的形狀。由演算法自動調整形狀，或由文字設計師手動調整形狀。

Glyph的資料結構

三種：點陣字dot-matrix、筆劃字stroke、輪廓字outline。

點陣字：用黑色方格（像素）拼成字形。

電腦已經淘汰點陣字型了，例如Fixedsys。然而有些字型為了支援極小尺寸，仍會搭配點陣字型，例如文泉驛點陣宋。

點陣字使用最多的地方不是電腦，而是大眾運輸的LCD跑馬燈。

筆劃字：直覺的方式。特別適合漢字。

記錄起點座標、終點座標、寫法（永字八法：側勒努趯策掠啄磔）（英文字母筆劃剖析）、曲率參數、變換矩陣。最後以演算法加粗筆劃，生成字型。優點是容易創造新字型，缺點是美感不佳。累死工程師、悶死設計師。

可惜的是，我們無法獲得筆劃資料。中文筆劃資料，屬於商業機密，更有專利保護。大家不得其門而入。例如王漢宗字型侵權案，個人猜測是從神秘管道得到筆劃資料，再自創加粗筆劃的演算法。

可惜的是，我們無法得知筆劃加粗的演算法。即便是開源字型，也從未公開演算法。必須跟字型公司簽技術合作。

這多少阻礙了組字、造字、文字資料庫、筆劃檢索、手寫辨識的進展。

輪廓字：公定的方式。現今字型皆是輪廓字。

輪廓由許多線條組成，線條一律是「直線線段：兩個端點」和「Bézier Curve：兩個端點與兩個控制點」。最後以演算法填充輪廓內部，生成字型。

幸運的是，我們可以獲得輪廓資料。安裝免費開源的字形編輯軟體FontForge，下載免費開源的字型思源黑體、文泉驛正黑，將字形匯出成.eps或.svg檔案即得。上面的「算」就是這樣來的。

注意到，字型都有著作權，未經授權不可使用字型檔案裡面的數據，即使免費。讀者做實驗時，一定要記得選擇免費且開源的字型，以免觸法。

輪廓字的檔案格式主要有兩種：囊括所有格式的OpenType（.otf）、早期的TrueType（.ttc .ttf）。OpenType是當今主流，網路上有不少資訊。如果你有志鑽研字型渲染，可以自己寫程式存取OpenType檔案，也可以參考開源的字型渲染引擎FreeType。

em square：字形的編輯區域

活字印刷當中，em是指鉛字寬度。電腦字型當中，em square是指字形的編輯區域。就這樣。

outlining：描繪輪廓

運用上個章節的繪製直線與繪製曲線演算法即可。

描繪輪廓目前沒有公定標準。即便是相同字型相同字形，在Firefox、Chrome、Windows、OS X的顯示結果通常略有不同。不同的字型渲染引擎，呈現不同結果。你也可以自己寫一個。

為了讓輪廓座標變成像素座標，這裡介紹一下pt和dpi。

電腦字型大小，單位是pt (PostScript point)，是長度單位。72pt = 1英吋，pt除以72得到英吋。電腦字型通常預設12pt。

長度變成像素數量，單位是dpi (dot per inch)，每英吋多少點（像素）。電腦螢幕通常是96dpi，實體印刷至少是300dpi。

比方說，12pt字型，96dpi解析度，一個字形的邊長含有12 ÷ 72 × 96 = 16個像素。

假設em square的邊長是1024。輪廓座標(x,y)，等比例縮放得到像素座標(round(x÷1024×16), round(y÷1024×16))。

filling：填充輪廓

運用上個章節的Scanline Fill Algorithm。輪廓不能交叉重疊，以利判斷內外。

如果讓輪廓具有方向性，順時針為正，逆時針為負，則得以設計聯集與差集。

overlap removing：重疊的輪廓，簡化成一個輪廓

為了順利填充輪廓，凡是輪廓交叉重疊，就必須重新整合。

有兩種狀況：一個輪廓自身交叉、兩個輪廓互相重疊。

需要處理的事情：求交點、求聯集或差集、重新安排頂點順序。需要的演算法：線段交點、Bézier Curve交點、多邊形交集。不太容易實作，讀者可以挑戰看看。

字形編輯軟體有此功能，節省設計師的作業時間。

hinting / grid fitting：微調字形位置

像素座標只能是整數。座標四捨五入變成整數，可能往左往下跑、可能往右往上跑，筆劃粗細有所改變。尤其是字形縮放大小，例如從12pt到14pt，問題更加嚴重。

像素大小固定。當字形縮小至非常小，例如6pt，所有筆劃通通黏結在一起。此時整個字形都得重新設計，甚至改成點陣字。

目前已經有演算法，可以自動微調字形位置。不過凡事總有例外，仍需設計師人工檢視、手動微調。各種常見縮放尺寸，都得一一檢視。中文數萬字，一字一字做，一日一日做，愚公移山、精衛填海！

anti-aliasing：微調字形邊緣的像素數值

在字形邊緣設定漸層顏色，使得字形邊緣柔順、清晰。

像素整齊排列。對於一條斜線，黑白顯示器，無可避免地，一定有鋸齒；灰階顯示器、彩色顯示器，則可以微調鄰近像素數值，達到平滑效果。

三種演算法：面積、距離、分割像素。

面積：像素視作正方形。如果字形佔據40%的面積，就將數值設定為40%，也就是255 - (255 × 40%) = 153。

距離：像素到輪廓的最短距離。貝茲曲線轉換成多項式，再解方程式（求根）求得最短距離。也可以用分治法。

分割像素：彩色顯示器，一個像素由RGB三個元件組成，調整RGB亮度比重，強調邊緣。缺點是邊緣五顏六色，相當刺眼，不是每個人都喜歡。知名演算法是微軟ClearType。

ICPC 3702

metrics：設定字形間距

字形的部分已經介紹完畢了，接下來是排版的部分。

FreeType的說明文件已經有精美圖片。這裡就不重畫了。

kerning：微調字形間距

比方來說，AV合在一起時，視覺感受是留白太多，需要更緊密一點。AB、AD不需要更緊密，保持原本排版即可。

解法是事先設定特殊組合。OpenType支援此設定。屬於文字處理的範疇。

另外還有一種情況，出現於阿拉伯文、泰文等等。這些語言根據語義，銜接或斷開字形。

解法是剖析語意義。OpenType沒有此設定，由排版引擎自行實作。屬於自然語言處理的範疇。

ligature：連合字形

比方來說，fi合在一起時，書寫感受是筆劃太多，需要更簡潔一點，於是連合成新字形。解法同前。

專門用於程式設計的字體，也有連合功能。例如Cascadia Code、Fira Code，將!=改成≠，將->改成→。

未涉及的主題

因為我沒有跟字型公司簽技術合作，所以知道的事情有限。

一、輪廓字的編輯介面。

二、輪廓字、筆劃字、點陣字之間的轉換。華康有個專利，演算法是迴歸。

三、拆字組字。簡化字型製作流程。

四、調字。根據環境、根據版面，自動調整字形。

五、手寫字形。工程師（詳細演算法）和設計師（詳細流程）有著不同的策略。

六、壓縮。減少字型檔案大小。例如Compact Table Format和MicroType Express。

七、列印文字。印表機韌體設計。

Font Design

這是設計師的專業。已有人寫書拍片，我就不班門弄斧了。

演算法生成的字形太過呆板（例如墨字的四個點如何點）。目前的作業流程是：先由工程師產生大致輪廓，再由設計師調整細節。

繁體中文的筆劃寫法，台灣政府有訂立標準。不過政府標準不見得是好看的，也不代表是正統的。

Typeface Design / Typography Design

也許可以看看spiro curve、文字雲、筆劃裝飾。

Math Rendering

畫數學式子與畫字形，原理相同，只是還要再解決一個問題：數學式子的資料結構為何？目前最流行的方式，是用純文字記錄一道算式，所採用的文法是TeX的語法。

讀者可以參考MathJax。

Web Rendering

3D Graphics

當今科技，沒有任何設備可以清楚呈現3D物體。

退而求其次，以2D圖片呈現3D物體。目前公認的做法：將3D模型投射到2D圖片上面，逐一設定每個像素的RGB值。

3D物體投射到2D圖片，兩種方式。

Orthogonal Projection：平行投射至圖片，夾角是垂直90度。這個方式非常簡單，卻有一個嚴重的問題：物體不論遠近都一樣大，缺乏真實感。

Perspective Projection：聚焦於圖片後方一點：遠的物體小、近的物體大，稍微符合人類視覺觀感。它是主流的方式，此處只介紹它。

聰明的讀者肯定馬上聯想到物理課的「針孔成像」及「攝影機」。雖然Perspective Projection與「針孔成像」都有聚焦的概念，但是兩者不盡相同。

抑或聯想到美術課的「透視圖」。雖然Perspective Projection與「一點透視圖」的外觀相同、原理相似，但是兩者的製作方式截然不同。

Perspective Projection細分兩種實作方式：

相機到物體：窮舉每一個像素，從焦點朝像素射出射線，找到對應的物體。速度極慢，用於電影動畫、擬真圖片。

物體到相機：窮舉每一個物體，從物體射出射線至焦點，找到對應的像素。速度較快，用於電玩遊戲、虛擬實境。

「Camera」與「Model」的觀念，源自其他領域。

ICPC 5100

3D Graphics相關資源

3D繪圖廣泛應用於電腦遊戲、電影特效、電視廣告、模型設計、建築設計、醫學影像等等。像是李安引進台灣的Rhythm & Hues Studios公司、比如CG Taiwaner教學網站，都屬於3D繪圖領域。

知名的3D繪圖函式庫，例如OpenGL與GLEW與GLFW、Vulkan、Direct3D。瀏覽器的3D繪圖函式庫，例如WebGL、three.js。程式語言通常沒有內建3D繪圖函式庫，必須另外安裝。

知名的3D模型資料，例如Infinite-Realities。

知名的3D渲染軟體，例如Maxwell Render、LuxRender、Octane Render。

Mesh Rendering: Camera to Model

Mesh Rendering

這裡再提供一個更清楚的動畫。

Model

將物體表面簡化成大量的平坦的多邊形，再運用「多邊形三角剖分」切割成大量的三角形，得到Mesh。

替現實生活的物體建立Mesh，是一套複雜的學問。所幸我們已經有現成的模型可以使用。例如此處的模型是Utah Teapot。

Mesh由許多三角形組成。一個三角形擁有三個頂點座標（暨法向量）、正面顏色、背面顏色。

三角形的頂點順序，決定了三角形的正面：正視三角形的正面，我們習慣讓三角形頂點呈逆時針順序。依照頂點順序計算「叉積」，得到三角形的正面的法向量。這是大家約定俗成、心照不宣的規矩。

Mesh的資料結構通常是：一個陣列記錄每個點（暨法向量）、一個陣列記錄每個三角形的三個點的編號。由於三角形經常共用頂點，因此兩層式的資料結構，得以節省記憶體空間。

不過為了方便起見，以下採用笨蛋方法。

struct Point {float x, y, z;};

Point operator+(Point p1, Point p2)
{
	return {p1.x + p2.x, p1.y + p2.y, p1.z + p2.z};
}

Point operator-(Point p1, Point p2)
{
	return {p1.x - p2.x, p1.y - p2.y, p1.z - p2.z};
}

Point operator*(Point p, float s)
{
	return {p.x * s, p.y * s, p.z * s};
}

Point operator/(Point p, float s)
{
	return {p.x / s, p.y / s, p.z / s};
}

istream& operator>>(istream& in, Point p)
{
	return in >> p.x >> p.y >> p.z;
}

ostream& operator<<(ostream& out, Point p)
{
	return out << '(' << p.x  << ',' << p.y << ',' << p.z << ')';
}

Point min(Point p1, Point p2)
{
	return {min(p1.x, p2.x), min(p1.y, p2.y), min(p1.z, p2.z)};
}

Point max(Point p1, Point p2)
{
	return {max(p1.x, p2.x), max(p1.y, p2.y), max(p1.z, p2.z)};
}

float dot(Point p1, Point p2)
{
	return p1.x * p2.x + p1.y * p2.y + p1.z * p2.z;
}

Point cross(Point a, Point b)
{
	return {a.y * b.z - b.y * a.z,
			a.z * b.x - b.z * a.x,
			a.x * b.y - b.x * a.y};
}
/*
float length(Point p)
{
	return sqrt(dot(p, p));
}
*/
float InvSqrt(float x)
{
	float xhalf = 0.5f * x;
	int i = *(int*)&x;
	i = 0x5f3759df - (i >> 1);
	x = *(float*)&i;
	x = x * (1.5f - xhalf * x * x);
	return x;
}

Point normalize(Point p)
{
//	return p / length(p);
	return p * InvSqrt(dot(p, p));
}

UVa 10711

Camera

替相機設定三個單位向量：長的方向、寬的方向、面對方向。

求得Mesh中心點，倒退一萬步，作為焦點的位置。

Ray Casting

建立焦點往像素的射線，找出第一個射中的三角形。

一、以相機的三個單位向量，建立焦點往像素的射線。

二、窮舉三角形。針對每一個三角形，以「平面與直線交點」求得交點和距離，隨時記錄最短距離。為了避免緩慢的sqrt運算，我們不求直線距離，而是求dirz方向的距離、也就是深度。

最後一行程式碼可以簡化。* dt / dv可以提到dot外面。view是由dirx、diry、dirz組成，這三個向量互相垂直；view投影到dirz上，不需要看dirx和diry，只需要看dirz方向有多少投影量，點積結果顯然是depth。

depth是定值。有些人為了加速，甚至不乘上depth，只求倍率。倍率1.0剛好位於圖片上。

以倍率輕鬆求得交點。

三、以「判斷點在凸多邊形內部」，判斷交點是否在三角形內部。三維版本，叉積改為三重積，三角形面積改為四面體體積。

// 擊中的三角形、交點、深度、擊中正面或反面。
struct Hit {Triangle t; Point hit; float z; bool side;};

Hit Camera::cast(int x, int y)
{
	// 預設沒有擊中三角形
	Hit ans = {-1, eye, 1e9, true};

// 焦點往像素的射線
	Point view = (dirx * (x - X/2))
			   + (diry * (y - Y/2))
			   + (dirz * depth);

// 焦點與三角形距離，投影到dirz上的長度。即深度。
	// 目標是找到距離最短的。（已簡化成倍率）
	float zvalue = 1e9;

// 窮舉所有三角形，判斷三角形與射線相交。
	// 找到第一個射中的三角形。
	for (int i=0; i<(int)triangleList.size(); ++i)
	{
		Triangle& t = triangleList[i];

// plane-ray intersection
		float dv = dot(view, t.normal);
		float dt = dot(t.vertex[0] - eye, t.normal);
		if (dv == 0) continue;
		float z = dt / dv;
		if (!(z > 0 && z < zvalue)) continue;
		Point hit = eye + (view * z);

// point in triangle test
		Point v0 = t.vertex[0] - hit;
		Point v1 = t.vertex[1] - hit;
		Point v2 = t.vertex[2] - hit;
		float c2 = dot(cross(v0, v1), t.normal);
		float c0 = dot(cross(v1, v2), t.normal);
		float c1 = dot(cross(v2, v0), t.normal);
		if (c0 < 0) c0 = -c0, c1 = -c1, c2 = -c2;
		if (!(c0 > 0 && c1 > 0 && c2 > 0)) continue;

// 確定擊中三角形，更新深度。（已簡化成倍率）
		zvalue = z;

// 儲存結果
		ans = (Hit){t, hit, z, dv < 0};
	}
	return ans;
}

Ray Casting: Bounding Volume Hierarchy

若想加速，將所有三角形存入「Bounding Volume Hierarchy」資料結構，依照區域分類所有三角形。

若想繼續加速，每個像素平行計算。實務上是運用顯示卡的GPU，撰寫CUDA程式碼，大幅加速。

struct Box {Point a, b;};
Box operator+(Box b1, Box b2)
{
	return (Box){min(b1.a, b2.a), max(b1.b, b2.b)};
}

// 判斷射線是否穿過Bounding Box。
// 由於浮點數誤差，我們只能用複雜的方式，
// 判斷射線是否打中Bounding Box的六個表面。
#define check(x, y, z)						\
	if (view.x != 0)						\
	{										\
		t = (box.a.x - eye.x) / view.x;		\
		y = eye.y + view.y * t;				\
		z = eye.z + view.z * t;				\
		if (box.a.y <= y && y <= box.b.y	\
		 && box.a.z <= z && z <= box.b.z)	\
			return true;					\
											\
		t = (box.b.x - eye.x) / view.x;		\
		y = eye.y + view.y * t;				\
		z = eye.z + view.z * t;				\
		if (box.a.y <= y && y <= box.b.y	\
		 && box.a.z <= z && z <= box.b.z)	\
			return true;					\
	}

bool RayBoxIntersect(Point eye, Point view, Box box)
{
	float x, y, z, t;
	check(x, y, z);
	check(y, z, x);
	check(z, x, y);
	return false;
}

// 本來應該是放在struct BVH裡面的變數，
// 由於C++限制太多，只好改成全域變數，方便使用。
vector<Box> box;
vector<Point> center;
bool cmpx(int i, int j) {return center[i].x < center[j].x;}
bool cmpy(int i, int j) {return center[i].y < center[j].y;}
bool cmpz(int i, int j) {return center[i].z < center[j].z;}
const int bvh_limit = 1;
vector<int> idx;
vector<Box> bvhbox;

struct Hit {Triangle t; Point hit; float z; bool side;} ans;
struct BVH {Point eye, view;} bvh;

void BVH::init()
{
	int N = triangleList.size();
	idx.resize(N);
	box.resize(N);
	center.resize(N);
	bvhbox.resize(N * 4);

for (int i=0; i<(int)triangleList.size(); ++i)
	{
		Triangle& t = triangleList[i];
		idx[i] = i;
		box[i] = (Box){
			min(t.vertex[0], t.vertex[1], t.vertex[2]),
			max(t.vertex[0], t.vertex[1], t.vertex[2])
		};
		center[i] = (box[i].a + box[i].b) / 2;
	}

build(1, 0, N-1);
}

void BVH::build(int o, int L, int R)
{
	Box b = box[idx[L]];
	for (int i=L+1; i<=R; i++)
		b = b + box[idx[i]];
	bvhbox[o] = b;

if (R-L+1 <= bvh_limit) return;

Point d = b.b - b.a;
	if (d.x >= d.y && d.x >= d.z)
		sort(idx.begin()+L, idx.begin()+R+1, cmpx);
	else if (d.y >= d.z)
		sort(idx.begin()+L, idx.begin()+R+1, cmpy);
	else
		sort(idx.begin()+L, idx.begin()+R+1, cmpz);

int M = (L + R) / 2;
	build(o*2, L, M);
	build(o*2+1, M+1, R);
}

void BVH::query(int o, int L, int R)
{
	if (!RayBoxIntersect(eye, view, bvhbox[o]))
		return;

if (R-L+1 > bvh_limit)
	{
		int M = (L + R) / 2;
		query(o*2, L, M);
		query(o*2+1, M+1, R);
		return;
	}

for (int i=L; i<=R; i++)
	{
		Triangle& t = triangleList[idx[i]];

// plane-ray intersection
		float dv = dot(view, t.normal);
		float dt = dot(t.vertex[0] - eye, t.normal);
//		if (dv >= 0) continue;	// backface culling
		if (dv == 0) continue;
		float z = dt / dv;
		if (!(z > 0 && z < ans.z)) continue;
		Point hit = eye + (view * z);

// point in triangle test
		Point v0 = t.vertex[0] - hit;
		Point v1 = t.vertex[1] - hit;
		Point v2 = t.vertex[2] - hit;
		float c2 = dot(cross(v0, v1), t.normal);
		float c0 = dot(cross(v1, v2), t.normal);
		float c1 = dot(cross(v2, v0), t.normal);
		if (c0 < 0) c0 = -c0, c1 = -c1, c2 = -c2;
		if (!(c0 > 0 && c1 > 0 && c2 > 0)) continue;

ans = (Hit){t, hit, z, dv < 0};
	}
}

Hit BVH::cast(Point eye, Point view)
{
	this->eye = eye;
	this->view = view;
	ans = (Hit){-1, eye, 1e9, true};
	query(1, 0, triangleList.size()-1);
	return ans;
}

UVa 12312

繪製三角形！

此處的模型是Cornell Box。運用OpenGL建立視窗，將模型畫在圖片上。

窮舉每個像素，建立焦點往像素的射線；窮舉每個三角形，找出第一個射中的三角形，取得三角形顏色。

void display()
{
	// 清空所有像素
	glClear(GL_COLOR_BUFFER_BIT);

// 窮舉所有像素
	for (int x = 0; x < X; ++x)
		for (int y = 0; y < Y; ++y)
		{
			// 取得像素顏色
			Hit h = camera.cast(x, y);
			// 判斷焦點面對三角形的正面還是背面
			Point color = h.side ? h.t.foreColor : h.t.backColor;

// 一個一個像素慢慢畫，不使用進階的OpenGL函式。
			glColor3f(color.x, color.y, color.z);
			glBegin(GL_POINTS);
			glVertex2f(x, y);
			glEnd();
		}

// double buffering
	glutSwapBuffers();
}

// 設定視窗座標系統
void reshape(int w, int h)
{
	glViewport(0, 0, (GLsizei)w, (GLsizei)h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluOrtho2D(0, w, 0, h);
}

// 讀檔、設定視窗、設定繪圖細節
int main(int argc, char **argv)
{
	const char fileName[] = "csie.tri";
	bool load = LoadFile(fileName);
	if (!load) cout << "Cannot open: " << fileName << endl;

camera.init();

glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
	glutInitWindowSize(400, 300);
	glutInitWindowPosition(100, 100);
	glutCreateWindow("Mesh Rendering");

glClearColor(0., 0., 0., 0.);
	glShadeModel(GL_SMOOTH);

glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutMainLoop();
	return 0;
}

相機、矩陣、圖片函式庫，三種不同的座標系統，必須小心。

Mesh Rendering: Illumination

Illumination

直接取用三角形的原始顏色，外觀呆板，缺乏層次。真實世界，顏色繽紛，是因為光線照射。現在我們要引進一套光線系統，讓三角形華麗動人、閃閃發光！

一、Light。光線的來源。

二、Optics。光線射中物體之後，物體影響光線的機制。

三、Material。光線射中物體之後，光線亮度與彩度變化。

四、Ray Tracing。追蹤光線軌跡。運算量最大的地方。

五、Sampling。設定光線數量、方向，使圖片外觀柔順。

這套光線系統源自物理，源自人類對於世界現象的觀察。如果你熟悉物理理論、擅長物理實驗，可以創造更細膩的方式。

Light

point light (omni light)：點光源。例如燭火、燈泡。
cone light (spotlight)：圓椎光源。例如聚光燈。
directional light (collimated light)：有向光源。例如遙遠的太陽。

area light：區域光源。光源們組成一大片面積。
model：物體當作光源，表面發光。
environment light：環境光源。光源們組成球面，包覆在場景外圍。

實際範例請見此。數量、位置、方向、亮度、彩度，五個要素。光源還可進一步組成面，添增變化。

妥善設置光源，精心營造氣氛。甚至可以實地拍攝天空的照片作為光源，諸如朝霧夕霞、天光雲彩，營造真實感。

Optics

emission：發光。例如螢光棒。
absorption：吸收。日常生活常見物體都是。

reflect：反射。例如鏡子、水中倒影。
refract：折射。例如魚缸、晶鑽。
scatter：散射。例如葉隙光、雲隙光、頭髮。
diffract：繞射。很少見。

物體影響光線。或者簡單想成：光線碰撞物體，產生變化。

光波有兩個要素：頻率（彩度）和振幅（亮度）。

一道光通常由大量頻率組成，各種頻率有各種振幅。物體影響光線，不會改變頻率，而會改變各種頻率的振幅大小、行進方向。

發光吸收是振幅大小變化，反射折射散射繞射是行進方向變化。

Optics的數學模型

一、考慮反射角（等於入射角）、折射角（根據介質係數）。入射折射反射都在同一平面。適合當作業。

二、考慮各種反射方向、折射方向。此模型稱作BSDF = BRDF + BTDF。雖然光線遵守物理定律、光線隨從明確軌跡，但是混雜了大量複雜結果，才是人類視覺觀感：光線不是一條而是一束，照射處不是一點而是一塊，介面不是一個平面而是很多曲折，射出不僅一個方向而是很多方向。例如虹彩。

三、考慮地點差異。此模型稱作BSSRDF。光線在介質內不斷散射，從別處散逸而出。光線的散射、粒子的擴散，兩者有著異曲同工之妙。距離交點越遠，射出機率越小，近似常態分布。亦可想成表面套用Gauss filter。例如濁液、皮肉、頭髮。

四、考慮時間差異。光線被介質吸收，緩緩射出。例如磷光、螢光。

五、考慮各種頻率。白光反射時，部分頻率消失（光的吸收）；白光折射時，各種頻率錯開（光的色散）。接近真實了。例如dichroism、pleochroism。

考慮越多，運算量、記憶量越大。四和五，當今硬體設備難以即時實現，除非取巧。

Material

人類視覺透過複雜機制，揉合振幅和頻率，得到顏色。

工程師便宜行事，將光的要素，從頻率和振幅，換成顏色和亮暗程度。再瞎掰出：顏色×亮暗程度=呈現顏色。儘管不切實際，但是容易計算。

工程師便宜行事，反射光、折射光，顏色不變動，亮暗程度才變動。一種入射方向，擁有多種反射方向、折射方向，每種方向各有亮暗程度。方向以經緯度表示，亮暗程度以函數表示。函數有形容詞。

  diffuse / matte：漫射。霧面、粗糙表面。
                   反射光分散，朝不同方向。例如路面、牆壁、紙張。
specular / glossy：鏡射。鏡面、光滑表面。
                   反射光集中，朝同一方向。例如金屬、塑膠、磁磚。

給定一種材質，有物理儀器可以測量函數。採用實體的測量數據，便可模仿真實世界。已經有人收集測量數據、製作觀看軟體。

除了採用實體的測量數據，亦可採用虛擬的數學模型。現今已有多種數學模型，各有不同視覺效果。

Ray Tracing

一、光線接觸介面才有變化，產生反射折射。適用固體，例如居家擺設。演算法原理類似狀態空間搜尋。

   ray tracing: 焦點出發。
                光學：僅考慮反射角、折射角。
  path tracing: 焦點出發（亦可焦點和光源雙向出發）。
                光學：考慮各種反射方向、折射方向。
     radiosity: 光源出發，走K步，每步皆留下一些光子。
                焦點出發，走一步，以附近光子數量，決定像素亮度。
               （實作：每步結束後，另走一步到焦點，累加像素亮度。）
               （I^1 + I^2 + I^3 + ... + I^k）
photon mapping: 光源出發，走一步，留下一個光子。
                焦點出發，走K步，以附近光子數量，決定該步的光線亮度。

一、從焦點、光源出發，持續反射、折射，追蹤光線軌跡。
　　每次接觸介面當作一步。
二、從焦點出發時，每一步窮舉所有可能的入射方向。
　　從光源出發時，每一步窮舉所有可能的出射方向。
　　由於不可能完全窮舉，所以採用取樣，後面章節會介紹。
三、焦點、點光源僅有一點，通常瞄不準、追蹤不到。
　　因此最後一步必須無視物理定律，刻意走向焦點、點光源。(global illumination)
　　區域光源則無此問題。
四、為了減低計算量，通常只走K步。
　　又為了減少bias，以固定機率決定要不要走下一步。(Russian roulette)
　　又為了增加實感，累計一步、兩步、……、K步的結果。

二、光線接觸介質即有變化，於介質中不斷反射折射散射。適用氣體液體，例如雲霧、牛奶。依照濃度不同，有著各種演算法，但是皆昧於物理事實。

稀薄：視作濾心。穿過越多介質，亮暗程度變化越大。
　　　以介質厚度當作亮暗變化程度。
中等：視作粒子。氣體處處都是完美漫射。光線為直線，只散射（轉折）一次。
　　　窮舉視線的每一個位置，嘗試轉折至光源。
濃稠：視作固體。光學模型：考慮入射、出射地點差異。
　　　測量BSSRDF，套用固體的演算法。

Ray Tracing的常見伎倆

    ambient occlusion: 遮光。擋住光線，產生陰影。
                       例如鼻翼的陰影、門縫。
subsurface scattering: 透光。光線於半透明介質內散射，深層與淺層的光線重疊。
                       例如白裡透紅的肌膚、透光的窗簾、濁液、葉片、玉。

soft light/shadow: 柔光、柔影。
從交點發出大量射線，判斷打中多少光源，決定亮度。
volumetric light/shadow: 光芒、陰暗。例如百葉窗散落的光、積雨雲散落的陰影。
從光源發出大量光線。再從焦點發出大量射線，判斷打中多少光線，決定亮度。

environment mapping: 物體周圍擺放圖片，映射在物體上。
shadow mapping: 預先計算物體陰影範圍，貼上陰影圖片。

電腦所能繪製的現象十分稀少。至今仍有許多光學現象，繪製尚未成功，同志仍須努力。

UVa 12313

Sampling

一、反射折射：追蹤光線軌跡，列舉數種可能的方向。

二、發射：增加焦點射向像素的光線數量，射中的顏色們取平均數，產生平滑、模糊、抗鋸齒效果。

anti-alias: 增加亂數數量，取平均數。
motion blur: 根據物體速率，決定亂數偏移量。
depth of field: 根據物體遠近，決定亂數偏移量。

射線數量是一個自訂常數，射線方向的偏移量是一群二維亂數。光線反射折射，採用球形；光線發射，採用正方形。

繪製三角形！

光源：一個點光源。光學：僅反射（入射角等於反射角）、無折射。材質：完美漫射。光線追蹤演算法：Ray Tracing。無取樣。

圖片很陽春。如果想看更逼真的圖片，可以上網搜尋，或者自己動手畫畫看吧！

Mesh Rendering: Texture

Texture

以上介紹了光線照明系統，以下介紹表面紋路系統。

物體表面通常不是均一顏色、不是均一高低。工程師的解決之道：自訂表面紋理。紋理包括了顏色、高度、法向量等類型。

Texture Generation

製造紋理的方法有：請攝影師拍攝實際照片，請設計師繪製美工圖片，請工程師研發演算法。

製造紋理屬於圖片處理的範疇。原理是運用雜訊，產生紋理的顏色、高度、法向量等等，營造逼真畫面。要是不想自行製作紋理，亦可下載現成的紋理（但是只有顏色）。

Texture Mapping

將紋理貼上物體表面。物體表面、紋理，建立座標對應關係。請參考這篇文章。

簡易的方式是：物體的每個三角形（四邊形），各使用一片正方形紋理。採用三角形座標（四邊形座標），得到座標對應關係。請參考Image Warping，屬於圖片處理的範疇。

Texture Mapping: Texture Coordinates

將紋理布置於物體周圍。重設紋理的座標系統。

正方體：使用六片紋理，作為正方體的六個面。
圓柱體：使用一片紋理，左右銜接。（座標對應：長寬對應高度和角度）
半球體：使用一片紋理，中央凸出。（座標對應：長寬對應經緯）

Texture Mapping: Mesh Parameterization

將紋理映射至物體表面。重設物體表面的座標系統。

攤平：物體從某處剪開、攤平，得到對應關係。
　　　周界可以盡量調整成紋理形狀、內界可以盡量調整成整齊網格。
中心：源自物體中心的放射線，擊中表面和紋理，得到對應關係。
光學：源自焦點的射線，經由反射擊中紋理，得到對應關係。

攤平細分為兩種映射方式。

UV mapping：物體表面切割成棋盤方格，橫向座標稱作U，縱向座標稱作V。
UVW mapping：承上，處處追加深度數值。深度座標稱作W。

Texture Filtering

bilinear filtering：只有一張紋理。採用image warping改變紋理形狀。
trlinear filtering：考慮相機遠近，預先製作各種大小的紋理。scale transform。
anisotropic filtering：考慮相機角度，預先製作各種視角的紋理。shear transform。

二維紋理的映射當中，紋理與表面的大小形狀，經常不相符，必須調整紋理的大小形狀。

改變紋理的大小形狀，紋理可能出現鋸齒、可能失真。仿效內插的觀念，讓紋理柔順平滑。

Texture Data Structure / Texture Atlas

   mipmap: 紋理預先縮放成各種比率，以應付各種遠近距離的呈現觀感。
heightmap: 附帶高度資訊。

紋理的資料結構。我沒有研究，請參考這篇文章、Wavefront的obj檔案格式。

Texture Embossing

重新設定表面法向量、表面高度，令表面凹凸不平、有立體感、有陰影，宛如壓印。請參考這篇文章、這篇文章。

bump mapping: 表面是平面，重新設定每一處的高低。
normal mapping: 表面是平面，重新設定每一處的法向量。
displacement mapping: 三角形頂點沿法向量移動。
parallax mapping: http://exibeo.net/docs/parallax_mapping.pdf
relief mapping: http://graphics.cs.brown.edu/games/SteepParallax/index.html

Normal Interpolation

重新設定表面法向量，使其圓滑。

想得到三角形上某一點的法向量：

一、求得三角形的正面的法向量。三個頂點叉積，即得。

二、求得三角形的三個頂點的法向量。一個頂點有許多個相鄰的三角形。相鄰三角形的法向量相加之後，當作該頂點的法向量。然而模型本身多半已經提供精準的頂點的法向量，可以省略這兩步驟。

三、求得三角形的任意一點的法向量。想要求得內插比重，直覺的方式是：將一次內插重新表示成線性變換矩陣，然後求反矩陣──然而當三角形與原點共平面，反矩陣就不存在了。推薦的方式是「Barycentric Interpolation」，三塊小三角形的面積，就是內插比重。

法向量們方向漸層改變。打光後，三角形表面彷彿圓滑凸起。

形狀看似改變了，但是實際位置沒變。物體邊緣會產生破綻。

Normal Encoding

將一個法向量重新表示成一個浮點數，節省儲存空間。請參考這篇文章。

Mesh Rendering: Model to Camera

Perspective Projection

如何計算三維空間的一個點，在圖片上對應的XY座標呢？

先以「點積」求得三個方向的投影量、即是真實距離x y z。

再以「相似三角形邊長成比例」的原理，藉由z和depth的比例，求得投影距離x' y'。x z構成的三角形，縮小為depth/z倍，得到x'；y z構成的三角形，縮小為depth/z倍，得到y'。

注意到，真實距離和投影距離都是以圖片中心為基準，距離可以是負值。最後我們調整投影距離成為圖片座標。

繪製空心三角形！

此處的模型是台大資訊系。

窮舉每個三角形，把三角形三個頂點投射到圖片上，求得圖片座標。運用OpenGL畫出每個三角形的外框，最後形成了「線框（wireframe）圖」。

void display()
{
	glClear(GL_COLOR_BUFFER_BIT);
	for (int i=0; i<(int)triangleList.size(); ++i)
	{
		Triangle& t = triangleList[i];
		Point p0 = camera.project(t.vertex[0]);
		Point p1 = camera.project(t.vertex[1]);
		Point p2 = camera.project(t.vertex[2]);

glColor3f(1., 1., 1.);	// 白色
		glBegin(GL_LINES);
		glVertex2f(p0.x, p0.y); glVertex2f(p1.x, p1.y);
		glVertex2f(p1.x, p1.y); glVertex2f(p2.x, p2.y);
		glVertex2f(p2.x, p2.y); glVertex2f(p0.x, p0.y);
		glEnd();
	}
	glutSwapBuffers();
}

void reshape(int w, int h)
{
	glViewport(0, 0, (GLsizei)w, (GLsizei)h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluOrtho2D(0, w, 0, h);
}

int main(int argc, char **argv)
{
	const char fileName[] = "csie.tri";
	bool load = LoadFile(fileName);
	if (!load) cout << "Cannot open: " << fileName << endl;

camera.init();

glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
	glutInitWindowSize(400, 300);
	glutInitWindowPosition(100, 100);
	glutCreateWindow("Wireframe");

glClearColor(0., 0., 0., 0.);
	glShadeModel(GL_SMOOTH);

glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutMainLoop();
	return 0;
}

Polygon Filling

大家應該會畫空心三角形了。現在來畫實心三角形吧！

把三角形三個頂點投射到圖片上，求得圖片座標。套用填充多邊形演算法，把三角形畫到圖片上。

鉤勒線段的部分，採用了「一次內插」，求得準確的三角形邊界；然後運用ceil和floor函數，讓相鄰三角形的接縫密合。

Triangle* triangle[X][Y];

void FillEdge(Point p1, Point p2)
{
	int xmini = max(0,   (int)ceil (min(p1.x, p2.x)));
	int xmaxi = min(X-1, (int)floor(max(p1.x, p2.x)));
	for (int x = xmini; x <= xmaxi; ++x)
	{
		// 一次內插。用X座標求得Y座標。
		float y = (p2.x == p1.x ? p1.y : (p2.y - p1.y) / (p2.x - p1.x) * (x - p1.x) + p1.y);
		ymin[x] = min(ymin[x], y);
		ymax[x] = max(ymax[x], y);
	}
}

void FillTriangle(Triangle& t)
{
	// 三角形三個頂點，圖片座標
	Point p[3];
	for (int i=0; i<3; ++i)
		p[i] = camera.project(t.vertex[i]);

// 計算X座標的極小值、極大值
	float xmin = +1e9, xmax = -1e9;
	for (int i=0; i<3; ++i)
	{
		xmin = min(xmin, p[i].x);
		xmax = max(xmax, p[i].x);
	}

// 初始化Y座標的極小值、極大值
	int xmini = max(0,   (int)ceil (xmin));
	int xmaxi = min(X-1, (int)floor(xmax));
	for (int x = xmini; x <= xmaxi; ++x)
	{
		ymin[x] = +1e9;
		ymax[x] = -1e9;
	}

// 計算Y座標的極小值、極大值
	for (int i=0; i<3; ++i)
		FillEdge(p[i], p[(i+1)%3]);

// 填充
	for (int x = xmini; x <= xmaxi; ++x)
	{
		int ymini = max(0,   (int)ceil (ymin[x]));
		int ymaxi = min(Y-1, (int)floor(ymax[x]));
		for (int y = ymini; y <= ymaxi; ++y)
			triangle[x][y] = &t;
	}
}

Triangle* triangle[X][Y];

void FillLine(float x1, float x2, int y, Triangle& t)
{
	if (!(y >= 0 && y < Y)) return;
	int xmini = max(0,   (int)ceil (x1));
	int xmaxi = min(X-1, (int)floor(x2));
	for (int x = xmini; x <= xmaxi; x++)
		triangle[x][y] = &t;
}

void FillTriangle(Triangle& t)
{
	// 三角形三個頂點，圖片座標
	Point A = camera.project(t.vertex[0]);
	Point B = camera.project(t.vertex[1]);
	Point C = camera.project(t.vertex[2]);
	// 依照y座標排序
	if (A.y > B.y) swap(A, B);
	if (B.y > C.y) swap(B, C);
	if (A.y > B.y) swap(A, B);
	// 計算斜率
	float dx1 = (B.y-A.y > 0) ? (B.x-A.x)/(B.y-A.y) : 0;
	float dx2 = (C.y-A.y > 0) ? (C.x-A.x)/(C.y-A.y) : 0;
	float dx3 = (C.y-B.y > 0) ? (C.x-B.x)/(C.y-B.y) : 0;
	// 三角形的邊界
	float x1 = A.x;
	float x2 = A.x;
	int y  = (int)ceil(A.y);	// 導致x1 x2不夠精確，可能漏填
	int By = (int)floor(B.y);
	int Cy = (int)floor(C.y);
	// 填充
	if (dx1 > dx2) {
		for (; y<=By; y++, x1+=dx2, x2+=dx1)
			FillLine(x1, x2, y, t);
		x2 = B.x;
		for (; y<=Cy; y++, x1+=dx2, x2+=dx3)
			FillLine(x1, x2, y, t);
	} else {
		for (; y<=By; y++, x1+=dx1, x2+=dx2)
			FillLine(x1, x2, y, t);
		x1 = B.x;
		for (; y<=Cy; y++, x1+=dx3, x2+=dx2)
			FillLine(x1, x2, y, t);
	}
}

Visible Surface Determination

近的三角形，必須擋住遠的三角形。求得三角形與焦點的距離，才能判斷誰近誰遠。

從焦點往圖片像素的射線，射向三角形，以「三角形與直線交點」求得距離。為了避免緩慢的sqrt運算，我們不求直線距離，而是求dirz方向的距離、也就是深度。

最後一行程式碼可以簡化。* dt / dv可以提到dot外面。view是由dirx、diry、dirz組成，這三個向量互相垂直；view投影到dirz上，只需要看dirz方向有多少投影量，點積結果顯然是depth。

depth是定值。有些人為了加速，甚至不乘上depth，只求個縮放倍率。倍率1.0剛好位於圖片上。

一邊填充三角形，一邊判斷深度。三角形的交接處，可以優先選擇靠近焦點的三角形，或者為了爭取速度而不做判斷。

古人曾經使用另一種演算法：運用「BSP Tree」預先排序所有三角形，依照先後順序繪製，就不必判斷遠近。但是無法處理三角形循環互疊、無法處理相機移動，速度也極慢，現今沒人用。

Illumination

此處介紹懶惰的打光方式。

光線與表面垂直，感覺最亮；光線與表面平行，感覺最暗。光線與表面的夾角大小，決定了亮暗。更精準來說，光線的垂直分量大小，決定了亮暗。

光線向量與法向量的點積，即是亮暗程度。設定為單位向量，亮暗程度就是0.0到1.0之間的數字。至於法向量可以使用上個章節提到的演算法求得。

ambient：物體發光，亮度為定值。
diffuse：光線向量與法向量的點積，作為亮暗程度。
specular：視線向量與反射向量的點積，作為亮暗程度。
          取k次方，夾角越小、亮度驟升。
modified specular：「視線向量與光線向量的中垂線」與法向量的點積，作為亮暗程度。
                   取k次方，夾角越小、亮度驟升。

Lambert model = c1⋅ambient + c2⋅diffuse
Phong model = c1⋅ambient + c2⋅diffuse + c3⋅specular
Blinn-Phong model = c1⋅ambient + c2⋅diffuse + c3⋅modified specular

Point light = {1,0,0};	// 平行光源方向（單位向量）

float Camera::illuminate(Triangle t, Point hit,
						 Point view, Point light)
{
	view  = normalize(view);
	light = normalize(light);

// 交點的法向量
	Point weight = interpolate(t, hit);
	Point normal = (t.vnormal[0] * weight.x)
				 + (t.vnormal[1] * weight.y)
				 + (t.vnormal[2] * weight.z);
	// 反射向量
	Point reflect = normalize(normal * (dot(light, normal) * -2.0) + light);
	// 視線向量與光線向量的中垂線
	Point bisector = normalize(view + light);

// 全域基本亮度（0.0 ~ 1.0）
	float ambient  = 1.0;
	// 垂直分量亮度（0.0 ~ 1.0）
	float diffuse  = fabs(dot(light, normal));
	// 亮者越亮，產生亮點（0.0 ~ 1.0）
	float specular = pow(fabs(dot(view, reflect)), 5.0);
	// 亮者越亮，產生亮點（0.0 ~ 1.0）
	float specular2 = pow(fabs(dot(normal, bisector)), 20.0);
	// 光源照射的、焦點看到的，如果是不同面，就無亮度。
	// （缺陷：即將翻面的三角形整片漏失）
	if (dot(view, t.normal) * dot(light, t.normal) <= 0)
		diffuse = specular = specular2 = 0;
	// 光源照射的、焦點看到的，如果是不同面，就無亮度。
	// （缺陷：即將翻面的三角形邊緣滲光）
//	if (dot(view, normal) * dot(light, normal) <= 0)
//		diffuse = specular = specular2 = 0;
	// 自行調配比重，揉合三種亮度。
	return (ambient * 0.2)
		 + (diffuse * 0.4)
		 + (specular * 0.4);
}

Point Camera::interpolate(Triangle t, Point p)
{
	Point v0 = t.vertex[0] - p;
	Point v1 = t.vertex[1] - p;
	Point v2 = t.vertex[2] - p;
	float c2 = dot(cross(v0, v1), t.normal);
	float c0 = dot(cross(v1, v2), t.normal);
	float c1 = dot(cross(v2, v0), t.normal);
	float c = c0 + c1 + c2;
	return (Point){c0, c1, c2} / c;
}

Shading

素描的shading，著重線條樣式；電腦繪圖的shading，著重填充方式。

flat shading：求三角形任一點的顏色（例如重心），其餘通通一樣。
Gouraud shading：求三角形三頂點的顏色，其餘一次內插。
Phong shading: 求三角形每一點的顏色，法向量內插。

繪製實心三角形！

Culling

backface culling：預先剔除物體背後的三角形。
occlusion culling：預先剔除被其他物體遮住的三角形。
visibility culling：預先剔除視野範圍（field of view）之外的三角形。

culling有三種，這裡只介紹backface culling。

如果物體有著密封外殼，背向焦點的三角形一定被遮住；如果物體有著破洞外殼、物體是一張薄紙，背向焦點的三角形才有機會露臉。

如果物體有著密封外殼，我們可以剔除背向焦點的三角形，避免無謂的填充，加快程式速度。

UVa 12628

Clipping

裁切有兩種。第一種是設定視野範圍（Field of View）。

我們可以設定左、右、上、下、遠、近邊界，設定六個平面作為邊界、裁切物體。值得一提的是，以近平面裁切物體，讓我們可以瞧見物體內部（此時就不應該實施backface culling了）。我們習慣讓近平面等於圖片位置。

當然也可以設定不規則形狀的邊界。

裁切有兩種。第二種是阻止三角形伸展至焦點背後。

三角形頂點位於焦點背後，投射到圖片上就會歪斜！先前我們倒退一萬步，巧妙避開這個議題；現在我們想瞧瞧物體內部、讓圖片與焦點靠近物體、裁切物體，就必須解決這個議題。

解決方式是：運用「多邊形裁切」，以近平面裁切三角形，去除太近的部分，之後才投射、填充。裁切之後剩下的多邊形，通常會再剖分成三角形，方便套用原本的投射、填充程式碼。

第一種裁切，對象是圖片像素；第二種裁切，對象是三角形。

3D Transformation

物體位移、縮放、旋轉；相機位移、縮放、旋轉。一共六種。

物體旋轉。物體自轉，換個觀點來看，就是相機繞物體公轉。與其搬動大量三角形，不如只搬動一個相機。旋轉dirx、diry、dirz三個向量，重新倒退一萬步，就完成物體旋轉。

相機旋轉。相機自轉，換個觀點來看，就是物體繞相機公轉。其實不必想太多，直接旋轉相機就好了。

物體旋轉通常用於3D繪圖軟體，相機旋轉通常用於第一人稱射擊遊戲，焦點旋轉則不常使用。由於clipping的關係，焦點旋轉時會看不見近身景色。第一人稱射擊遊戲的角色是站在相機中心、採用相機旋轉，而不是站在焦點、採用焦點旋轉。

物體縮放。以物體中心為原點，縮放三角形頂點座標。

相機縮放。以相機焦點為原點，縮放三角形頂點投射之後的圖片座標。

物體縮放是先縮放再投影，相機縮放是先投影再縮放，兩者的視覺效果並不相同。物體縮放有鏡頭拉近拉遠的效果，相機縮放則是整張畫面等比例放大縮小。

物體位移、相機位移。懶的講了，兩者視覺效果相同。

繪製動感三角形！

// 用鍵盤的w與s控制遠近
void keyboard(unsigned char key, int x, int y)
{
	if		(key == 'w') camera.radius -= 50;
	else if (key == 's') camera.radius += 50;
	camera.rotate();	// 記得更新焦點位置
	glutPostRedisplay();
}

// 用滑鼠旋轉物體
int press_x, press_y;
Point press_angle;
void mouse(int button, int state, int x, int y)
{
	if (state == GLUT_DOWN)
	{
		press_x = x; press_y = y;
		press_angle = camera.angle;
	}
}

// 用滑鼠旋轉物體
void motion(int x, int y)
{
	camera.angle = press_angle + (Point){y - press_y, 0, x - press_x} / 40;
	camera.rotate();	// 記得更新焦點位置
	glutPostRedisplay();
}

int main(int argc, char **argv)
{
	const char fileName[] = "csie.tri";
	bool load = LoadFile(fileName);
	if (!load) cout << "Cannot open: " << fileName << endl;

camera.init();

glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
	glutInitWindowSize(400, 300);
	glutInitWindowPosition(100, 100);
	glutCreateWindow("Computer Graphics");

glClearColor(0., 0., 0., 0.);
	glShadeModel(GL_SMOOTH);

// 每當需要重新繪圖，自動呼叫display函式。
	glutDisplayFunc(display);
	// 每當改變視窗大小，自動呼叫reshape函式。
	glutReshapeFunc(reshape);
	// 每當按下鍵盤按鍵，自動呼叫keyboard函式。
	glutKeyboardFunc(keyboard);
	// 每當按下滑鼠按鍵，自動呼叫mouse函式。
	glutMouseFunc(mouse);
	// 每當移動滑鼠游標，自動呼叫motion函式。
	glutMotionFunc(motion);
	// 開始繪圖
	glutMainLoop();
	return 0;
}

Volume Rendering

Volume Rendering（Voxel Rendering）

Voxel

「像素Pixel」和「體素Voxel」兩者概念相仿，二維與三維的差別而已。

像素數值，通常代表RGB顏色；體素數值，通常代表物質密度。物質密度越高，體素數值越高。

替現實生活的物體建立體素，是一套複雜的學問。所幸網路上已經有電腦斷層掃描的資料。此處使用的CTHead是8bit TIFF圖片檔案，圖片一張張疊起來，像素成了體素。

這份資料中，皮肉密度低、體素數值低，骨骼密度高、體素數值高。我不清楚物質密度與體素數值是否恰好成正比。

Normal Vector

以「梯度gradient」的單位向量，當作法向量。梯度的計算過程很簡單：右邊減左邊、遠處減近處、上面減下面，除以兩個間距，分別得到XYZ三個方向的變化程度。

Ray Marching（Ray Casting）

體素視作固體solid、液體liquid、氣體gas，有著不同的光線追蹤演算法。此處先介紹固體的演算法：

一、找到射線從容積的哪裡射入、哪裡射出。

二、套用鉤勒直線演算法，求得射線碰到的體素們。

三、設定臨界值，求得射線無法穿越的那個體素，即為所求。

Point Camera::color(int x, int y)
{
	// 焦點往圖片像素的射線
	Point view = (dirx * (x - X/2))
			   + (diry * (y - Y/2))
			   + (dirz * depth);
	Point normalizeview = normalize(view);

// 求得射線打中容積的哪一個表面、求得z-value。
	float zmin = +1e9, zmax = -1e9;
	hitvolume(view, zmin, zmax);
	if (zmin == +1e9 || zmax == -1e9)
		return (Point){0,0,0};

// DDA algorithm，求得射線碰到的體素們。
	Point pmin = eye + (view * zmin);
//	Point pmax = eye + (view * zmax);
//	Point length = pmax - pmin;
	Point length = view * (zmax - zmin);
	int step = ceil(max(
		fabs(length.x),
		fabs(length.y),
		fabs(length.z)
	));
	Point gap = length / step;

for (int i=0; i<step; ++i)
	{
		pmin = pmin + gap;
		int x = round(pmin.x);
		int y = round(pmin.y);
		int z = round(pmin.z);
		if (!involume(x, y, z)) continue;

// 體素小於60，射線就持續穿越體素。
		if (img[z][y][x] < 60) continue;

// 體素大於等於60，即為所求。馬上著色。
		Point normal = normalize(gradient(x,y,z));
		float s = illuminate(normal, normalizeview);
		return (Point){1,1,1} * s;
	}

// 射線穿透整個容積，無色。
	return (Point){0,0,0};
}

設定不同的臨界值，就得到不同器官。皮肉密度低、骨骼密度高，因此臨界值低就得到皮肉，臨界值高就穿越皮肉、得到骨骼。

最困難的地方，就是如何讓電腦自動找到臨界值、讓電腦自動畫出最顯眼的圖片。這屬於影像辨識、人工智慧的領域，讀者可以自行尋找相關資料。

一種臨界值，得到一種器官。採用半透明顏色，就可以同時畫出很多器官。

Ray Marching（Ray Casting）

此處介紹液體、氣體的演算法：

一、體素一共有256種數值（以此模型為例，其他模型不見得皆如此），替每一種數值設定顏色、暨透明程度。

二甲、由遠往近讀取體素，揉合顏色：遠方顏色與當前顏色，依照比例（當前顏色的透明程度）揉合，遞推下去。

二乙、亦得由近往遠讀取體素：透明程度增加一個次方，乘上當前顏色，累加下去。

設定不同的顏色暨透明程度，得到不同的視覺感受。

Ray Marching: Sparse Voxel Octree

若想加速，將所有體素存入「Octree」資料結構，依照區域分類所有體素。適用於體素數量稀少、遠少於容積大小的情況。

Isosurface Rendering

Isosurface（Implicit Surface）

isosurface: f(x,y,z) = t , f(x,y,z) - t = 0 , g(x,y,z) = 0
g(x,y,z) > 0: point (x,y,z) is outside
g(x,y,z) = 0: point (x,y,z) is coincide
g(x,y,z) < 0: point (x,y,z) is inside

「等值表面」和「體素」概念相仿，連續與離散的差別而已。

訂立一個連續函數，讓空間中每個地點皆有數值。再訂立一個臨界值，作為表面。請參考這篇文章。

臨界值移項，使右式為零。左式的正、負、零，恰好代表該地點位於表面的外部、內部、上面。

找到射線與表面的交點，有兩種方式：一、鉤勒直線演算法，一步一步走，隨時判斷內外。二、解聯立方程式，二分法求根。

另外，記得訂立場景範圍，讓射不中物體的射線停止行進。

找到表面的法向量：使用梯度當作法向量，如同體素。

Distance Field（Signed Distance Field）

distance field: f(x,y,z) = signed distance from point (x,y,z) to surface

union        of f and g: min(f, g)
intersection of f and g: max(f, g)
complement   of f      : -f
difference   of f and g: f - g

「距離場」。等值表面究極進化。函數值恰是點到表面的最短距離！請參考這篇文章的參考文獻、這篇文章的模型範例。

當模型是距離場，多個模型的聯集、交集、差集非常好算！聯集：最小值。交集：最大值。補集：負值。差集：減法。

找到射線與表面的交點：當前位置代入函數，得到當前最短距離。走一步，步伐是當前最短距離，頂多碰到表面，而不會穿越表面。就算因為浮點數誤差穿越表面也無妨，下一步會得到負的距離，依舊持續接近表面。反覆行走，直到足夠靠近表面，或者直到超出場景範圍。

找到表面的法向量：使用梯度當作法向量，如同體素。

Ray Marching（Ray Casting）

當模型是距離場，射線與表面的交點，計算過程比較特別。

朝向表面：最短距離越來越小，迅速逼近。偏離表面：起初最短距離越來越小，掠過表面之後，最短距離越來越大，迅速遠離。

當表面平滑柔順，最短距離越來越小，逼近速度大致穩定。

當表面凹凸起伏，最短距離忽大忽小，逼近速度不穩定。

利用逼近速度，可以製做粗糙的ambient occlusion。當逼近速度不穩定，表示該處附近有物體凸出，迫近射線，遮擋交點。

沿著射線小步小步走，累計每一步當中，理想的、實際的最短距離的差值，作為陰影強度。

範例：Blobby Surface

範例：Fractal Surface