【十天自制軟渲染器】DAY 02:畫一條直線(DDA 演算法 & Bresenham’s 演算法)

鹵代烴發表於2021-01-19

寫文不易,懇求各位觀眾老爺 點贊 ?,收藏 ?,評論 ? 三連支援一下!!!謝謝你,這對我真的很重要!


第一天我們搭建了 C++ 的執行環境並畫了一個點,根據 點 → 線 → 面 的順序,今天我們講講如何畫一條直線。

本文主要講解直線繪製演算法的推導和思路(莫擔心,只涉及到一點點的中學數學知識),最後會給出程式碼實現,大家放心的看下去就好。

1.DDA 直線演算法

1.1 簡單實現

我們先來回顧一下中學的幾何知識,如何在二維平面內表示一條直線?最常見的就是斜截式了:

$$y = kx + b$$

其中斜率是 $k$,直線在 $y$ 軸上的截距是 $b$。


斜截式在數學上是沒啥問題的,但是在實際的工程專案中,因為硬體資源是有限的,我們不可能也沒必要表示一條無限長度的直線,現實往往是已知一條線段的起點 $(x_1, y_1)$ 和終點 $(x_2, y_2$),然後把它畫出來。

這時候用兩點式表示一根直線是最方便的,其中 $x_1 \leq x \leq x_2$,$y_1 \leq y \leq y_2$:

$$\frac{x-x_{1}}{x_{2}-x_{1}}=\frac{y-y_{1}}{y_{2}-y_{1}}$$

把上面的式子稍作變形,可以把 $x$ 和 $y$ 用引數 $\lambda$ 表示:

$$\left.\begin{array}{l}x=\lambda\left(x_{2}-x_{1}\right)+x_{1} \\ y=\lambda\left(y_{2}-y_{1}\right)+y_{1}\end{array}\right\} 0 \leq \lambda \leq 1$$

這時候我們只要取不同的 $\lambda$,就可以得出對應的 x 和 y。


按照以上的思路,我們可以用程式碼實現一下。C++ 的實現也很簡單,如下所示(dl 表示 $d \lambda$):

void line(
  int x1, int y1, 
  int x2, int y2, 
  TGAImage &image, TGAColor color) { 
    const float dl = 0.01;
    int dx = x2 - x1;
    int dy = y2 - y1;
    for (float t=0.0; t<1.0; t+=dl) { 
        int x = x1 + dx * t;
        int y = y1 + dy * t;
        image.set(x, y, color);
    } 
}

這個是直線演算法的初步實現,只能說「能用」,地位和排序演算法裡的「氣泡排序」一樣,目的達到了,但是效能不太好:

  • 每畫一個點,都要執行兩次乘法
  • 大量使用浮點運算(眾所周知,$v_{float}$ < $v_{fixed}$ < $v_{int}$)
  • 如果 dl 取的比較小,會導致一個畫素點會被繪製多次,重複計算
  • 如果 dl 取的比較大,會導致直線斷掉

1.2 優化

下面我們就一步一步優化上面的演算法。

首先我們注意到,對於螢幕繪製直線這個場景,理論上是連續的,但實際是離散的

比如說 $x$ 從 $x_1$ 變化到 $x_2$ 時,每次繪製時,$x$ 都是按步長 1 增長的,也就是 $x_{new} = x_{old} + 1$。

這時候 $y_{new} = y_{old} + \frac{y_2 - y_1}{x_{2}-x_{1}} = y_{old} + k$。

我們把上面的公式寫成程式碼,就是下面這個樣子:

void line(
  int x1, int y1, 
  int x2, int y2, 
  TGAImage &image, TGAColor color) { 
  float x = x1;
  float y = y1;
  float step = std::abs(x2 - x1);
  float dlx = (x2 - x1) / step;
  float dly = (y2 - y1) / step;
  
  for (int i=1; i<step; i++) { 
    image.set(x, y, color);
    x = x + dlx;
    y = y + dlx;
  } 
}

這個演算法其實還有一點兒問題,就是繪製斜率大於 1 的直線時,繪製出的直線會斷掉。比如說從 (0, 0) 點繪製到 (2, 4) 點,按照上面的演算法只會繪製兩個點,但是我們期望的是右圖那樣,起碼各個畫素要連線起來:

不連續的線 vs 連續的線

解決方法也很簡單,繪製這種比較「陡峭」的直線時(斜率絕對值大於 1),以 y 的變化為基準,而不是以 x,這樣就可以避免上面直線不連續情況。

最後的直線演算法就是這樣:

void line(
  int x1, int y1, 
  int x2, int y2, 
  TGAImage &image, TGAColor color) { 
  float x = x1;
  float y = y1;
  int dx = x2 - x1;
  int dy = y2 - y1;
  float step;
  float dlx, dly;

  // 根據 dx 和 dy 的長度決定基準
  if (std::abs(dx) >= std::abs(dy)) {
    step = std::abs(dx);
  } else {
    step = std::abs(dy);
  }

  dlx = dx / step;
  dly = dy / step;

  for (int i=1; i<step; i++) {
    image.set(x, y, color);
    x = x + dlx;
    y = y + dly;
  }
}

然後我們用這個演算法測試一下不同起點不同斜率的直線,看效果執行良好:


這個演算法就是經典的 DDA (Digital differential analyzer) 演算法),他比我們一開始的程式碼要高效的多:

  • 消除了迴圈內的乘法運算
  • 避免了重複的繪製運算
  • 保證線段連續不會斷掉

但是它還有個很耗效能的問題:計算過程中涉及大量的浮點運算

作為渲染器最底層的演算法,我們肯定希望是越快越好。下面我們就來學習一下,消除浮點運算的 Bresenham’s 直線演算法。

2.Bresenham’s 直線演算法

2.1 初步實現

本節內容不會從一開始就講完善版的 Bresenham’s 演算法,我們先從一個小節開始推導,最後推匯出完善的演算法。

最一開始,我們先考慮所有直線裡的一個子集,即斜率範圍在 $[0, 1]$ 之間的直線:$0 \leq k \leq 1$。

上一小節裡我們說過,對於螢幕繪製直線這個場景,理論上是連續的,但實際是離散的。我們先假設已經繪製了一個點 $(x,\ y)$,那麼在畫素螢幕上,下一個新點的位置,只可能有兩種情況:

  • $(x + 1,\ y)$
  • $(x + 1,\ y + 1)$

那麼問題就轉化為,下一個新點的位置該如何選擇?

我想大家應該都想到方案了,大體思路如下

  • 先把 $x_{new} = x + 1$ 這個值帶入直線方程裡,算出來 $y_{new}$ 的值
  • 然後比較 $y_{new}$ 和 $y + 0.5$ 的大小

    • $y_{new} \leq y + 0.5$,選點 $(x + 1,\ y)$
    • $y_{new} > y + 0.5$,選點 $(x + 1,\ y + 1)$

我們再把思路完善一下,把每次取捨時的誤差考慮進去:

day2_Bresenham_line

如上圖所示,實際上繪製的點的位置是 $(x, y)$,理論上點位置是 $(x,\ y + \epsilon)$。

當點從 $x$ 移動到 $x + 1$ 時,理論上新點的位置應該是 $(x + 1,\ y + \epsilon + k)$,其中 k 是直線的斜率。

實際繪製時,要比較 $y + \epsilon + k$ 和 $y + 0.5$ 的大小:

  • $y + \epsilon + k \leq y + 0.5$,選點 $(x + 1,\ y)$
  • $y + \epsilon + k > y + 0.5$,選點 $(x + 1,\ y + 1)$

對於下一個新點 $x + 2$,我們可以按照下式更新誤差 $\epsilon$:

  • 若前一個點選擇的是 $(x + 1,\ y)$,則 $\epsilon_{new} = (y + \epsilon + k) - y = \epsilon + k$
  • 若前一個點選擇的是 $(x + 1,\ y + 1)$,則 $\epsilon_{new} = (y + \epsilon + k) - (y + 1) = \epsilon + k - 1$

把上面的思考過程用虛擬碼表示一下:

$
begin{aligned}
& epsilon leftarrow 0, quad y leftarrow y_{1} \
& pmb{for} x leftarrow x_{1} pmb{to} x_{2} pmb{do} \
& quad pmb{draw} pmb{point} pmb{at} (x, y) \
& quad pmb{if} ( (epsilon + k) < 0.5) \
& quad quad epsilon leftarrow epsilon + k \
& quad pmb{else} \
& quad quad y leftarrow y + 1 \
& quad quad epsilon leftarrow epsilon + k - 1 \
& quad pmb{end} pmb{if} \
& pmb{end} pmb{for}
end{aligned}
$

2.2 消除浮點運算

觀察上面的虛擬碼,我們可以發現這裡面出現了 0.5,也就是說存在浮點運算。下面我們就通過一些等價的數學變換消除浮點數。

首先對於不等式 $\epsilon + k < 0.5$,我們給它不等號左右兩邊同時乘以 2 倍的 $\Delta x$,這樣就可以同時消除斜率除法和常量 0.5 帶來的浮點運算:

$$\epsilon + \Delta y / \Delta x < 0.5$$

$$2 \epsilon \Delta x + 2 \Delta y <\Delta x$$

然後用 $\epsilon^{\prime}$ 表示 $\epsilon\Delta x$,上式可以轉換為 $$2(\epsilon^{\prime} + \Delta y)< \Delta x$$

同樣的,我們在更新 $\epsilon$ 時,把它也替換為 $\epsilon^{\prime}$ ,也就是對於下面兩式:

  • $\epsilon = \epsilon + m$
  • $\epsilon = \epsilon + m - 1$

等號兩邊同時乘以 $\Delta x$,有:

  • $\epsilon \Delta x = \epsilon \Delta x+\Delta y$
  • $\epsilon \Delta x = \epsilon \Delta x+\Delta y-\Delta x$

然後用 $\epsilon^{\prime}$ 表示 $\epsilon\Delta x$,可以得到:

  • $\epsilon^{\prime} = \epsilon^{\prime}+\Delta y$
  • $\epsilon^{\prime} = \epsilon^{\prime}+\Delta y-\Delta x$

這時候我們就可以得到一個去掉浮點數運算的虛擬碼:

$
begin{aligned}
& epsilon^{prime} leftarrow 0, quad y leftarrow y_{1} \
& pmb{for} x leftarrow x_{1} pmb{to} x_{2} pmb{do} \
& quad pmb{draw} pmb{point} pmb{at} (x, y) \
& quad pmb{if} ( 2 (epsilon^{prime} + Delta y) < Delta x) \
& quad quad epsilon^{prime} leftarrow epsilon^{prime} + Delta y \
& quad pmb{else} \
& quad quad y leftarrow y + 1 \
& quad quad epsilon^{prime} leftarrow epsilon^{prime} + Delta y - Delta x \
& quad pmb{end} pmb{if} \
& pmb{end} pmb{for}
end{aligned}
$


C++ 實現如下:

void line(Screen &s,
  int x1, int y1,
  int x2, int y2,
  TGAImage &image, TGAColor color) {
  int y = y1;
  int eps = 0;
  int dx = x2 - x1;
  int dy = y2 - y1;

  for (int x = x1; x <= x2; x++) {
    image.set(x, y, color);
    eps += dy;
    // 這裡用位運算 <<1 代替 *2
    if((eps << 1) >= dx)  {
      y++;
      eps -= dx;
    }
  }
}

這樣我們就實現了斜率在 $[0, 1]$ 區間的高效演算法。也就是說,現在我們可以繪製 1/8 個象限的直線了。剩下範圍的直線,可以通過交換 xy 等方式實現繪製。具體的實現都是些髒活累活,就不擺出來了,感興趣的可以去 GitHub 上看程式碼的完整實現


3.繪製模型

這一部分可以結合原英文教程學習,我只做一些細節上的補充。

前面兩個小節都是演算法基礎學習,本小節開始載入一個非洲人的 .obj 模型,然後把模型上每個三角形面的點連線起來。

OBJ 檔案是一種被廣泛使用的 3D 模型檔案格式(obj 為字尾名),用來描述一個三維模型。模型關鍵字較為繁瑣,限於篇幅本文暫不展開,大家可以自行搜尋學習。

這一節的流程也很清楚:從磁碟上載入 .obj 檔案 → 按行分析 .obj 檔案 → 構建 model → 迴圈 model 中的每個三角形 → 連線三角形的三條邊 → 渲染出圖


上訴流程的前三步已經被原作者封裝好了,我們直接把原始碼裡的 model.hmodel.cpp 拖到主工程裡就可以了,感興趣的人可以看一下原始碼實現,非常簡單,在一個 while 迴圈裡一直 readline 就可以了,因為和圖形學關係不大,我這裡就略過了。

最後的畫三角形的程式碼如下,關鍵步驟我已經用註釋標註了:

// 例項化模型
model = new Model("obj/african_head.obj");

// 迴圈模型裡的所有三角形
for (int i = 0; i < model->nfaces(); i++) {
  std::vector<int> face = model->face(i);

  // 迴圈三角形三個頂點,每兩個頂點連一條線
  for (int j = 0; j < 3; j++) {
    Vec3f v0 = model->vert(face[j]);
    Vec3f v1 = model->vert(face[(j + 1) % 3]);
    
    // 因為模型空間取值範圍是 [-1, 1]^3,我們要把模型座標平移到螢幕座標中
    // 下面 (point + 1) * width(height) / 2 的操作學名為視口變換(Viewport Transformation)
    int x0 = (v0.x + 1.) * width / 2.;
    int y0 = (v0.y + 1.) * height / 2.;
    int x1 = (v1.x + 1.) * width / 2.;
    int y1 = (v1.y + 1.) * height / 2.;
    
    // 畫線
    line(x0, y0, x1, y1, image, white);
  }
}

最後渲染出的影像如下:

toyrenderer_day02_obj



今天學習瞭如何畫一條線,明天我們學習如何畫一個三角形


部落格原文,閱讀體驗更佳(公式沒有錯亂)

大家對圖形學感興趣的話,可以關注 ?️ 號「滷蛋實驗室」,後臺回覆「圖形學」獲取經典教材《虎書4》和《Real Time Rendering 4》

寫文不易,懇求各位觀眾老爺 點贊 ?,收藏 ?,評論 ? 三連支援一下!!!謝謝你,這對我真的很重要!


參考連線:

Line Drawing on Raster Displays

The Bresenham Line-Drawing Algorithm

DDA Line Drawing Algorithm - Computer Graphics

Bresenham's Line Drawing Algorithm



相關文章