Bresenham演算法

original author : http://liar.pangwa.com/2009/04/12/bresenham/

在看一本書《Windows遊戲編程大師技巧》 (Tricks of Windows Game Programming Gurus). 這次繼續書裏的內容: 直線光柵化的Bresenham演算法. 書上講的比較含糊, 沒有講演算法的推導過程, 更沒講演算法是怎麼想出來的. 所以我們只好自己動手, …

直線光柵化

直線光柵化是指用圖元點來類比直線. 比如下圖中用藍色的圖元點來類比紅色的直線. 圖中坐標系是顯示器上的坐標系: x軸向右, y軸向下.

 

 

設deltaX = endX – startX, deltaY = endY – startY. 那麼斜率為k = deltaY / deltaX. 我們先考慮簡單的情況: 當 0 < k < 1即直線更貼近x軸. 在這種情況下deltaY < deltaX, 所以在光柵化的過程中, 在y軸上描的點比在x軸上描點少. 那麼就有一個很直觀的光柵化演算法:

line_bresenham(startX, startY, endX, endY)

{

    deltaX = endX - startX;

    deltaY = endY - startY;

    k = deltaY / deltaX;

    for (x = startX, y = startY; x <= endX; ++x)

    {

        if (滿足一定條件)

        {

            ++y;

        }

        drawPixel(x, y);

    }

}

基於斜率 / 距離的兩個簡單直線光柵化演算法

好了，貌似很簡單, 就剩一個問題: “滿足一定條件”是什麼? 可以用斜率判斷, 也可以用上圖中直線與光柵線交點 (紅點) 與光柵點 (藍點) 的距離來判斷. 繼續用偽代碼說話:

// 演算法1: 用斜率判斷

void line_bresenham_k(startX, startY, endX, endY)

{

    deltaX = endX - startX;

    deltaY = endY - startY;

    k = deltaY / deltaX;

    for (x = startX, y = startY; x <= endX; ++x)

    {

        if (x - startX != 0)

        {

            // 計算當前斜率

            currentK = (y - startY) / (x - startX);

            // 如果當前斜率 < k, 則增加y座標

            if (currentK < k)

            {

                ++y

            }

        }

        drawPixel(x, y);

    }

}

// 演算法2: 用距離判斷. 計算直線與光柵線交點y座標我們需要用到

// 直線方程 y = k (x - startX) + startY

line_bresenham_dist(startX, startY, endX, endY)

{

    deltaX = endX - startX;

    deltaY = endY - startY;

    k = deltaY / deltaX;

    for (x = startX, y = startY; x <= endX; ++x)

    {

        // 計算直線與光柵線交點的y座標, 以及與光柵點的距離

        ptY = k * (x - startX) + startY;

        dist = ptY - y;

        // 如果距離 > 0.5或者 < -0.5, 說明我們需要增加y以

        // 將距離的絕對值控制在0.5之類

        if (dist > 0.5 || dist < -0.5)

        {

            ++y;

        }

        drawPixel(x, y);

    }

}

 

消滅浮點數!

以上都是很直觀的演算法, 下面不直觀的來了 – 上面的演算法都需要在循環體內執行乘法, 準確的說, 是進行浮點數的乘法. 我們怎麼能減少這些浮點數的乘法開銷呢? 以基於距離的演算法2為例: 首先, k是一個浮點數, 0.5也是浮點數. 我們可以通過將這些運算式都乘以2 * deltaX (整數) 來解決浮點數的問題. 偽代碼:

// 演算法3: 在演算法2的基礎上消滅浮點數!

line_bresenham_dist(startX, startY, endX, endY)

{

    deltaX = endX - startX;

    deltaY = endY - startY;

    for (x = startX, y = startY; x <= endX; ++x)

    {

        // 計算直線與光柵線交點的y座標, 以及與光柵點的距離

        ptY1 = deltaY * (x - startX) + startY * deltaX;

        dist1 = ptY1 - y * deltaX;

        dist1 = dist1 << 1; // dist1 = dist1 * 2

        // 如果距離 > 0.5或者 < -0.5, 說明我們需要增加y以

        // 將距離的絕對值控制在0.5之類

        if (dist1 > deltaX || dist < -deltaX)

        {

            ++y;

        }

        drawPixel(x, y);

    }

}

 

消滅乘法!

圓滿解決浮點數運算問題! 不過…乘法運算還在. 消滅乘法問題的辦法比較不直觀, 讓我們想一想: 還有什麼辦法能簡化運算. 直線方程已經不能再簡化, 所以唯一的突破口就是能不能利用遞推 / 用上一次迴圈的計算結果推導下一次迴圈的計算結果.

首先我們來看看在演算法2的基礎上 (因為演算法2計算紅點藍點之間的距離, 比較直觀), 怎麼通過第n – 1次迴圈計算出的dist值 (設為d1) 來推導出第n次迴圈的dist值 (設為d2). 先回顧一下: dist = 直線與光柵線交點的y座標 – 相應光柵點的y座標. 我們從幾何上直觀地考慮: 在第n次迴圈中, 我們先根據上一次迴圈所計算出來的d1, 暫時令d2 = d1 + k, 因為我們要保證-0.5 < d2 < 0.5, 而d1 + k滿足d1 + k > –0.5, 所以我們只需要考慮當d1 + k > 0.5時, 我們需要將光柵點y座標增加1, 並且將d2減去1. 顯然, 設y1是第n – 1次迴圈中光柵點的y座標, y2是第n次迴圈中光柵點的y座標. 我們有
1) d2 = d1 + k –  (y2 – y1)
2) 當d1 + k > 0.5時y2 = y1 + 1, 否則y2 = y1
   我們已經能根據上面的兩個關係式寫出演算法, 不過為了消除乘法和浮點數, 我們將這兩個關係式兩端同時乘以2 * deltaX,   並且設e = 2 * deltaX * d, 則我們有
3) e2 =  e1 + 2 * deltaY – 2 * deltaX * (y2 – y1)
4) 當e1 + 2 * deltaY > deltaX時y2 = y1 + 1, 否則y2 = y1
  終於, 沒有了乘法 (2 * deltaY在循環體外計算且被簡化為左移一位的運算), 沒有了浮點數, 根據關係式3) 和 4), 寫出演算法:

// 演算法4: 在演算法2, 3的基礎上利用遞推消滅乘法和浮點數!

line_bresenham(startX, startY, endX, endY)

{

    deltaX = endX - startX;

    deltaY = endY - startY;

    e = 0;

    deltaX2 = deltaX << 1;

    deltaY2 = deltaY << 1;

    drawPixel(startX, startY);

    for (x = startX + 1, y = startY; x <= endX; ++x)

    {

        // 關係式3) e2 =  e1 + 2 * deltaY – 2 * deltaX * (y2 – y1)

        // 關係式4) 當e1 + 2 * deltaY > deltaX時y2 = y1 + 1, 否則y2 = y1

        e += deltaY2;

        if (e > deltaX)

        {

            e -= deltaX2;

            ++y;

        }

        drawPixel(x, y);

    }

}

 

消滅浮點數! 代數推導

上面遞推關係的推導過程是從圖形上”直觀”地分析得來的, 但是不嚴密. 我們能不能形式化地證明關係式1), 2), 3), 4)呢? 因為關係式3), 4)和1), 2)能互相推導, 我們只證明3), 4)如下:

在演算法3的基礎上設第n – 1次迴圈計算出的dist1值為e1, 對應的y值為y1, 第n次迴圈計算出的dist1值為e2, 對應的y值為y2. 根據演算法3,
dist1 = 2 * deltaY * (x – startX) + 2 * startY * deltaX – 2 * y * deltaX, 則
e2 – e1
= 2 * deltaY * (x – startX) + 2 * startY * deltaX – 2 * y2 * deltaX – [2 * deltaY * (x – 1 – startX) + 2 * startY * deltaX – 2 * y1 * deltaX ]
=  – 2 * y2 * deltaX + 2 * deltaY + 2 * y1 * deltaX
= 2 * deltaY – 2 * deltaX * (y2 – y1)
所以e2 = e1 + deltaY – deltaX * (y2 – y1). 所以我們有關係式
1) e2 = e1 + 2 * deltaY – 2 * deltaX * (y2 – y1)
2) –deltaX <  e1 < deltaX
3) –deltaX < e2 < deltaX
4)  y2 – y1 = 0 或者 1 
我們根據e1 + 2 * deltaY的取值範圍進行討論. 首先, 因為不等式6), 我們有
2 * deltaY – deltaX < e1 + 2 * deltaY < 2 * deltaY + deltaX

情況1: 如果2 * deltaY – deltaX < e1 + 2 * deltaY < deltaX, 則
2 * deltaY – deltaX – 2 * deltaX * (y2 – y1) < e2 < deltaX– 2 * deltaX * (y2 – y1) 
反證: 若y2 – y1 = 1, 則 2 * deltaY – deltaX – 2 * deltaX < e2 < deltaX – 2 * deltaX = -deltaX, 所以y2 – y1 = 1不成立. 即情況1中y2 = y1.

情況2:  如果 deltaX < e1 + 2 * deltaY < 2 * deltaY + deltaX, 則
deltaX – 2 * deltaX * (y2 – y1) < e2 < 2 * deltaY + deltaX – 2 * deltaX * (y2 – y1)
反證: 若y2 – y1 = 0, 則 deltaX < e2 < 2 * deltaY + deltaX 所以y2 – y1 = 0不成立. 即情況2中y2 = y1 + 1.

打了這麼多字, 累…以上就是當0 < k < 1的情況, 剩餘幾種情況 (k > 1, –1 < k < 0, k < –1. 不要挑剔我不用”>=”這種符號…) 都可以通過簡單的x, y交換以及正負交換來搞定.

 

The other line:http://www.cnblogs.com/soroman/archive/2006/07/27/509602.html