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的軟體已經發展成熟。然而本篇文章不談如何製作軟體，也不談如何使用軟體。

像素

const int X = 1024, Y = 768; struct Color {unsigned char r, g, b;}; Color image[Y][X]; void drawPixel(int x, int y, Color c = (Color){0,0,0}) { image[y][x] = c; } bool onImage(int x, int y) { return x >= 0 && x < X && y >= 0 && y < Y; } void initImage() { memset(image, 0xff, sizeof(image)); // 白色 }

直線

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

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

void drawSegment(Point p1, Point p2) { int xmin = max(0, (int)ceil (min(p1.x, p2.x))); int xmax = min(X-1, (int)floor(max(p1.x, p2.x))); for (int x = xmin; x <= xmax; ++x) { float y = (p2.x == p1.x ? p1.y : (x - p1.x) / (p2.x - p1.x) * (p2.y - p1.y) + p1.y); drawPixel(x, round(y)); } }

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

void drawSegment(Point p1, Point p2) { Point length = p2 - p1; int step = ceil(max( fabs(length.x), fabs(length.y), fabs(length.z) )); Point gap = length / step; for (int i=0; i<step+1; ++i) { int x = round(p1.x); // 四捨五入 int y = round(p1.y); if (!onImage(x, y)) continue; drawPixel(x, y); p1 = p1 + gap; } }

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

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

曲線

大家習慣使用Bézier curve。

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

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

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

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

void drawCurve(Point p1, Point p2, Point p3, Point p4) { if (fabs(p1.x - p4.x) < 1.0 && fabs(p1.y - p4.y) < 1.0) return; Point a1 = (p1 + p2) / 2.0; Point a2 = (p2 + p3) / 2.0; Point a3 = (p3 + p4) / 2.0; Point b1 = (a1 + a2) / 2.0; Point b2 = (a2 + a3) / 2.0; Point c = (b1 + b2) / 2.0; drawPixel(round(c.x), round(c.y)); drawCurve(p1, a1, b1, c ); drawCurve(c , b2, a3, p4); }

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

手繪線

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

可視作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合在一起時，書寫感受是筆劃太多，需要更簡潔一點，於是連合成新字形。解法同前。

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

未涉及的主題

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

一、輪廓字的編輯介面。

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

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

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

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

六、壓縮。減少字型檔案大小。例如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

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

model

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

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

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

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

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

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

struct Triangle {Point foreColor, backColor, vertex[3], vnormal[3], normal;}; vector<Triangle> triangleList; bool LoadFile(const char* fileName) { ifstream fin(fileName); if (!fin) return false; string s; while (fin >> s) { /* 讀檔 */ Triangle t; fin >> t.foreColor >> t.backColor; for (int i = 0; i < 3; ++i) fin >> t.vertex[i] >> t.vnormal[i]; /* 額外處理 */ // OpenGL的色彩值範圍是0...1，而不是0...255。 t.foreColor = t.foreColor / 256; t.backColor = t.backColor / 256; // 向量預先正規化，往後點積運算就不用除以向量長度。 for (int i = 0; i < 3; ++i) t.vnormal[i] = normalize(t.vnormal[i]); // 三角形法向量（有許多種求法，使用其中一種即可。） t.normal = cross(t.vertex[1] - t.vertex[0], t.vertex[2] - t.vertex[0]); // t.normal = cross(t.vertex[1] - t.vertex[0], t.vertex[2] - t.vertex[1]); t.normal = normalize(t.normal); if (fin) triangleList.push_back(t); } return true; } 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中心點，倒退一萬步，作為焦點的位置。

// 圖片大小 const int X = 400, Y = 300; struct Camera { float radius, depth; Point dirx, diry, dirz, center, eye; } camera; void Camera::init() { // 倒退的距離不要太近也不要太遠，不然看不見整個物體。 radius = 1000; // 圖片與焦點的距離（必須比倒退距離少） depth = 300; // 長的方向、寬的方向、面對方向（單位向量） dirx = (Point){1,0,0}; diry = (Point){0,0,1}; dirz = (Point){0,1,0}; // mesh中心點 center = mesh_center(); // 焦點的位置，從物體中心倒退一萬步。 eye = center - (dirz * radius); } // 求得mesh中心點 Point Camera::mesh_center() { Point a = {+1e9, +1e9, +1e9}; Point b = {-1e9, -1e9, -1e9}; for (int i=0; i<(int)triangleList.size(); ++i) for (int j=0; j<3; ++j) { a = min(a, triangleList[i].vertex[j]); b = max(b, triangleList[i].vertex[j]); } return (a + b) / 2; }

ray casting

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

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

Hit Camera::cast(int x, int y) { // 焦點往像素的射線。（隱性轉型成float） Point view = (dirx * (x - X/2)) + (diry * (y - Y/2)) + (dirz * depth); /* // 加上0.5，對準像素中心。此處不採用，便宜行事。 Point view = (dirx * (x - X/2 + 0.5)) + (diry * (y - Y/2 + 0.5)) + (dirz * depth); */ }

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

for (int i=0; i<(int)triangleList.size(); ++i) { Triangle& t = triangleList[i]; // 平面與直線交點 float dv = dot(view, t.normal); float dt = dot(t.vertex[0] - eye, t.normal); // 三角形與視線平行，距離無限大。 // 其中三角形與視線共面時，當作是太薄了看不見。 if (dv == 0) return 1e9; // 以點積求得dirz方向的距離。dirz已經是單位向量。 float z = dot(view * dt / dv, dirz);

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

float z = dt / dv * depth;

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

float z = dt / dv;

以倍率輕鬆求得交點。

float z = dt / dv; Point hit = eye + (view * z);

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

// 擊中的三角形、交點、深度、擊中正面或反面。 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; }

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

Point color[X][Y]; void make_image() { // 調整迴圈 for (int y = Y-1, i = 0; y >= 0; --y, ++i) for (int x = 0, j = 0; x < X; ++x, ++j) { Hit h = camera.cast(x, y); image[i][j] = h.side ? h.t.foreColor : h.t.backColor; } }

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的數學模型

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

Point reflect(Point v, Point hit, Point normal) { return normal * (dot(v, normal) * -2.0) + v; } Point refract(Point v, Point hit, Point normal, float index1, float index2) { Point nv = normalize(v); float r = index1 / index2; // ratio of refractive index float c = -dot(nv, normal); // cos(θ1) // hit backface (can intergrate to the last line) if (c < 0) {c = -c; normal = -normal;} float d = 1.0 - r * r * (1.0 - c*c); // cos²(θ2) // total internal reflection if (d < 0) return (Point){0,0,0}; return (nv * r) + (normal * (r * c - sqrt(c2))); }

二、考慮各種反射方向、折射方向。此模型稱作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步的結果。

// 模型全部都是球，不是mesh。 http://kevinbeason.com/smallpt/

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

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

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

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」，三塊小三角形的面積，就是內插比重。

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

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

Point Camera::interpolate(Triangle t, Point hit) { // 點到三角形頂點的向量 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); float c = c0 + c1 + c2; // 一次內插。以三角形三個頂點的法向量，得到交點的法向量。 // 由於檔案裡面已經提供頂點法向量，就直接拿來用了。 Point n = (t.vnormal[0] * (c0 / c)) + (t.vnormal[1] * (c1 / c)) + (t.vnormal[2] * (c2 / c)); return n; }

normal encoding

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

perspective projection

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

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

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

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

Point Camera::project(Point p) { // 運用點積求得真實距離 Point view = p - eye; Point d = { dot(view, dirx), dot(view, diry), dot(view, dirz) }; // 相似三角形邊長成比例： // 縮小為 depth/z 倍，求得投影距離。 if (d.z == 0) return d; // 避免除以0！ d = d / fabs(d.z) * depth; // 求得圖片座標 return (Point){X/2 + d.x, Y/2 + d.y, d.z}; }

繪製空心三角形！

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

窮舉每個三角形，把三角形三個頂點投射到圖片上，求得圖片座標。運用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方向的距離、也就是深度。

float Camera::zvalue(int x, int y, Triangle t) { Point view = (dirx * (x - X/2)) + (diry * (y - Y/2)) + (dirz * depth); float dv = dot(view, t.normal); float dt = dot(t.vertex[0] - eye, t.normal); if (dv == 0) return 1e9; // 三角形與視線平行 return dot(view * dt / dv, dirz); }

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

return dt / dv * depth;

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

return dt / dv;

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

float zvalue[X][Y]; // z-buffer，初始化為無限大 Triangle* triangle[X][Y]; // 每個像素對應的三角形 void FillTriangle(Triangle& t) { ...... 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) { // 判斷深度，取深度最小者 float z = camera.zvalue(x, y, t); if (z < zvalue[x][y]) { zvalue[x][y] = z; triangle[x][y] = &t; } } } }

古人曾經使用另一種演算法：運用「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: 求三角形每一點的顏色，法向量內插。

// Phong shading Point Camera::color(int x, int y, Triangle t, float zvalue) { // 從焦點往圖片像素的射線。 Point v = (dirx * (x - X/2)) + (diry * (y - Y/2)) + (dirz * depth); // 取得三角形顏色。判斷焦點面對三角形的正面還是背面。 Point c = (dot(v, t.normal) < 0) ? t.foreColor : t.backColor; // 顏色×亮暗程度=呈現顏色 return c * illuminate(t, eye + (v * zvalue), v); }

繪製實心三角形！

void display() { glClear(GL_COLOR_BUFFER_BIT); // 初始化z-buffer for (int x = 0; x < X; ++x) for (int y = 0; y < Y; ++y) zvalue[x][y] = +1e9; // 填充每個三角形 for (int i=0; i<(int)triangleList.size(); ++i) FillTriangle(triangleList[i]); // 畫出所有像素 for (int x = 0; x < X; ++x) for (int y = 0; y < Y; ++y) { if (zvalue[x][y] == +1e9) continue; Point c = camera.color(x, y, *triangle[x][y], zvalue[x][y]); glColor3f(c.x, c.y, c.z); glBegin(GL_POINTS); glVertex2f(x, y); glEnd(); } glutSwapBuffers(); }

culling

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

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

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

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

bool Camera::toward(Triangle t) { return dot(t.vertex[0] - eye, t.normal) < 0; } void display() { ...... for (int i=0; i<(int)triangleList.size(); ++i) if (camera.toward(triangleList[i])) FillTriangle(triangleList[i]); ...... }

UVa 12628

clipping

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

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

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

void FillEdge(Point p1, Point p2) { // 左邊界和右邊界。太左、太右就不畫。 int xmini = max(20, (int)ceil (min(p1.x, p2.x))); int xmaxi = min(80, (int)floor(max(p1.x, p2.x))); ...... } void FillTriangle(Triangle& t) { ...... // 左邊界和右邊界。太左、太右就不畫。 int xmini = max(20, (int)ceil (xmin)); int xmaxi = min(80, (int)floor(xmax)); for (int x = xmini; x <= xmaxi; ++x) { ymin[x] = +1e9; ymax[x] = -1e9; } ...... for (int x = xmini; x <= xmaxi; ++x) { // 下邊界和上邊界。太低、太高就不畫。 int ymini = max(30, (int)ceil (ymin[x])); int ymaxi = min(70, (int)floor(ymax[x])); for (int y = ymini; y <= ymaxi; ++y) { float z = camera.zvalue(x, y, t); // 近邊界和遠邊界。太近、太遠就不畫。 // if (z < camera.depth) continue; if (z < 1.0) continue; if (z > 5.0) continue; if (z < zvalue[x][y]) { zvalue[x][y] = z; triangle[x][y] = &t; } } } }

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

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

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

省略

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

3D transformation

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

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

struct Matrix {float a[3][3];}; // 矩陣乘矩陣 Matrix operator*(Matrix m1, Matrix m2) { Matrix m = {{{0,0,0},{0,0,0},{0,0,0}}}; for (int i=0; i<3; ++i) for (int j=0; j<3; ++j) for (int k=0; k<3; ++k) m.a[i][j] += m1.a[i][k] * m2.a[k][j]; return m; } // 矩陣乘向量 Point operator*(Matrix m, Point p) { return (Point){ m.a[0][0] * p.x + m.a[0][1] * p.y + m.a[0][2] * p.z, m.a[1][0] * p.x + m.a[1][1] * p.y + m.a[1][2] * p.z, m.a[2][0] * p.x + m.a[2][1] * p.y + m.a[2][2] * p.z }; } Point angle; void Camera::init() { radius = 1000; depth = 300; center = mesh_center(); angle = (Point){0,0,0}; // 一開始沒有旋轉 rotate(); } void Camera::rotate() { // 製作三維旋轉矩陣 float c, s; c = cos(angle.x); s = sin(angle.x); Matrix rx = {{{1,0,0},{0,c,s},{0,-s,c}}}; c = cos(angle.y); s = sin(angle.y); Matrix ry = {{{c,0,-s},{0,1,0},{s,0,c}}}; c = cos(angle.z); s = sin(angle.z); Matrix rz = {{{c,s,0},{-s,c,0},{0,0,1}}}; // 實施旋轉 Matrix r = rx * ry * rz; dirx = r * (Point){1,0,0}; diry = r * (Point){0,0,1}; dirz = r * (Point){0,1,0}; light = r * (Point){1,0,0}; // 光源跟著轉 // 倒退一萬步 eye = center - (dirz * radius); }

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

物體旋轉通常用於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; }

voxel rendering（volume rendering）

voxel

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

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

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

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

const int N = 99; const int VX = 256; const int VY = 256; const int VZ = N * 2 - 1; short img[VZ][VY][VX]; void LoadFile() { // 運用OpenCV依序讀取99張圖片。 for (int z=0; z<N; ++z) { char name[50]; sprintf(name, "cthead-8bit\\cthead-8bit%03d.tif", z+1); cv::Mat image = cv::imread(name); for (int x=0; x<VX; ++x) for (int y=0; y<VY; ++y) { cv::Vec3b intensity = image.at<cv::Vec3b>(y, x); img[z*2][y][x] = (intensity.val[0] + intensity.val[1] + intensity.val[2]) / 3; } } // 所有相鄰圖片之間，一次內插一張圖片， // 讓頭殼看起來不會太扁。 for (int z=0; z<N-1; ++z) for (int x=0; x<VX; ++x) for (int y=0; y<VY; ++y) img[z*2+1][y][x] = (img[z*2][y][x] + img[z*2+2][y][x]) / 2; }

normal vector

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

bool involume(int x, int y, int z) { return x >= 0 && x < VX && y >= 0 && y < VY && z >= 0 && z < VZ; } Point gradient(int x, int y, int z) { if (involume(x+1, y+1, z+1) && involume(x-1, y-1, z-1)) // 稍後會計算單位向量，此處不必除以2。 return (Point){ img[z][y][x+1] - img[z][y][x-1], img[z][y+1][x] - img[z][y-1][x], img[z+1][y][x] - img[z-1][y][x] }; else return (Point){0, 0, 0}; }

ray marching（ray casting）

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

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

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

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

void Camera::hitvolume(Point view, float& zmin, float& zmax) { // 容積的六個面（各面只取一個頂點） static Point vertex[6] = { {0,0,0},{0,0,0},{0,0,0}, {VX,0,0},{0,VY,0},{0,0,VZ} }; // 容積的六個面的法向量 static Point normal[6] = { {1,0,0},{0,1,0},{0,0,1}, {1,0,0},{0,1,0},{0,0,1} }; // 求得射線打中容積的哪一個表面、求得z-value。 zmin = +1e9, zmax = -1e9; for (int i=0; i<6; ++i) { // ray-plane intersection float dv = dot(view, normal[i]); float dt = dot(vertex[i] - eye, normal[i]); if (dv == 0) continue; float z = dt / dv; if (!(z > 1.0 && z < +1e9)) continue; // 由於浮點數誤差， // 我們只能用複雜的方式，判斷射線是否打中表面。 Point hit = eye + (view * z); // if (involume(hit.x, hit.y, hit.z)) if (( (hit.x >= 0 && hit.x < VX) || (normal[i].x == 1.0) ) && ( (hit.y >= 0 && hit.y < VY) || (normal[i].y == 1.0) ) && ( (hit.z >= 0 && hit.z < VZ) || (normal[i].z == 1.0) )) { // 如果射線打中表面，就更新z-value。 if (z < zmin) zmin = z; if (z > zmax) zmax = z; } } } 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種數值（以此模型為例，其他模型不見得皆如此），替每一種數值設定顏色、暨透明程度。

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

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

void InitColor() { // 初始化顏色 for (int i=0; i<255; ++i) color_func[i] = (Point){0,0,0}; // 皮肉顏色 for (int i=60; i<85; ++i) color_func[i] = (Point){1,0.7,0}; // 骨骼顏色 for (int i=100; i<130; ++i) color_func[i] = (Point){1,1,1}; // 初始化為完全透明 for (int i=0; i<255; ++i) alpha_func[i] = 0; // 皮肉不透明程度（皮肉顏色密度） for (int i=60; i<85; ++i) alpha_func[i] = 0.2; // 骨骼不透明程度（骨骼顏色密度） for (int i=100; i<130; ++i) alpha_func[i] = 0.9; } Point Camera::color(int x, int y) { ...... Point color = {0,0,0}; for (int i=0; i<step; ++i) { // 由遠往近讀取體素 pmax = pmax - gap; int x = round(pmax.x); int y = round(pmax.y); int z = round(pmax.z); if (!involume(x, y, z)) continue; // 著色 Point normal = normalize(gradient(x,y,z)); float s = illuminate(normal, normalizeview); Point c = color_func[img[z][y][x]] * s; // 揉合顏色 float a = alpha_func[img[z][y][x]]; color = color * (1.0 - a) + c * a; } return color; }

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

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