[計算機圖形學]光柵化演算法：DDA和Bresenham演算法
一、DDA
DDA演算法是最簡單的直線繪製演算法。主要思想是利用直線的斜截式：\(y=kx+b\)
對於一條直線的繪製，往往會給定兩個端點：\(P_A = (0,0)\)和\(P_B = (60,60)\)
然後呼叫函式：OLED_DrawLine(0, 0, 60, 60);
首先，我們來看一下繪製直線都可以用哪些方法。
確定了兩個端點，那麼兩個端點之間的點如何確定？
第一個可以想到：知道了兩個點，可以用兩點式來確定直線方程，然後帶入點不就行了！
calc line expression f(x) from two point for x in range(x_start, x_end): y = f(x) draw(x, y)
這是一個方法，不過這個方法可能要反覆進行復雜數學運算（畢竟不可能所有的直線都是：\(y=x、y=x+1、...\)這類吧，但凡帶點乘除、浮點的，統統都是複雜運算）
這個方法貌似不可行，不僅要確定直線的表示式，還要進行多次複雜的浮點數乘法運算。還記得在之前的文章（將三角形重心座標時）簡單介紹過線性組合...
巧了，這裡正好有兩個點，妥妥的線性組合哇。
線性組合公式：\(P' = tP_A + (1-t)P_B, 其中 t \in [0,1]\)
for t in range(0, accu, x num):
(x, y) = tPa + (1-t)Pb
draw(x, y)
雖然不用先計算直線的表示式，但是需要2N次的浮點數乘法。有點得不償失。
這個方法也不得行。
那麼可不可以每次讓\(P_A\)加上一個值，\(P_B\)減去一個值，然後利用迭代的思想，逐個求出每一個點呢？
答案是可以的。
不妨用下面的形式進行稍微改進：
知道兩個點，那麼必然會知道某一個方向的步進次數（如果你樂意，甚至可以隨便選取一個方向，我們這裡選\(x\)方向），那麼另一個方向（\(y\)方向）的步長\(\delta\)就知道了，雖然說確定步長會涉及到浮點除法運算，但是畢竟只用計算一次，還是可以接受的。然後根據步進次數，\(y=y'+\delta\)，\(y'\)是上一個點的值（注意，雖然繪製到螢幕前要取整，但這裡儲存上一次值時保留取整前的。
calc y step value dy y = y start loop x num: y += dy draw(x, y)
可以很明顯的察覺這個方法相比之前的暴力計算，計算量少了許多：整體來看，此方法用到一次浮點數除法，N次浮點數加法，相比N次浮點數乘法（暫且認為乘除是一樣的，實際這種認為有失偏頗，畢竟大多時候能不用除法就不用）運算量降低了許多。
第三種方法便是DDA演算法的前身。
但是DDA演算法給出了更為明確的流程。
設當前點：\((x_i, y_i)\)
根據兩個端點便可計算出\(\delta x\)和\(\delta y\)
則下一點：
\(x_{i+1} = x_i + \delta x\)
\(y_{i+1} = y_i + \delta y\)
根據前面提到的，需要一個方向的\(\Delta\)值為1，即每次步進一個畫素
那麼如何確定這個步進方向？
DDA演算法給出明確的方法。
分別計算\(x\)和\(y\)方向上的差值\(\Delta x\)和\(\Delta y\)
- 如果\(\Delta x > \Delta y\)，說明\(x\)軸變化的大，所以把\(x\)方向作為主步進方向，即：\(\delta x = 1, \delta y = \frac{\Delta y}{\Delta x} = k\)
- 如果\(\Delta y > \Delta x\)，說明\(y\)軸變化的大，所以把\(y\)方向作為主步進方向，即：\(\delta y = 1, \delta x = \frac{\Delta x}{\Delta y} = \frac{1}{k}\)
仍然通過迭代的方式，即可求出每一個點。
可以看到，DDA演算法去掉了浮點數乘法運算，仍需要多次浮點數加法運算和浮點數取整。因此還是有優化空間的。
二、Bresenham
Bresenham演算法是一種基於誤差判別式來生成直線的方法。
同樣採用步進的思想，令每次最大變化方向的座標步進一個畫素單位（與DDA演算法相同），另一個方向根據誤差判別式的符號決定是否也要步進一個畫素單位（與DDA演算法不同）。
從Bresenham演算法的思想描述中可以看出，本質上和DDA沒有太大區別，只不過是另一個方向的步進值的確定產生了變化。
為什麼在另一個方向上每次最大隻步進一個畫素？
這一點很好解釋：
因為DDA演算法和Bresenham演算法都選取最大變化方向為主步進方向，這也就意味著，另一個方向的步進值無論是\(\delta y = \frac{\Delta y}{\Delta x}\)還是\(\delta x = \frac{\Delta x}{\Delta y}\)都必然小於等於1 。
另外，Bresenhan演算法誤差判別過程如下圖所示。
那麼，Bresenham演算法和DDA演算法區別在哪？就一個步進值麼？
不是的，Bresenham演算法和DDA演算法的區別在於最後的光柵化過程（就是望螢幕上繪製的時候）。至於這個步進值的差異，不是很關鍵。
- DDA演算法光柵化時，使用了一個浮點數轉化巨集：
#define FloatToInteger(fn) ((fn>0)?(fn+0.5):(fn-0.5)) - 而Bresenham演算法光柵化時，使用的誤差判別式
可以看到DDA始終於偏向選擇無腦“步進”
如果Bresenham演算法就到這，那也太low了：不僅沒有去掉DDA演算法中的浮點數加法運算，僅僅是為了讓步進更絲滑而引入誤差判別式，結果又引入了新的浮點數運算，圖啥？
對此，改進版的Brehensam演算法應運而生。
實際上，我沒有去考證Bresenham在1965年發表Brehensam演算法的論文，所以也不清楚那篇論文中就是改進後的。完了，不嚴謹了。
圖片來自https://www.cnblogs.com/LiveForGame/p/11706904.html
\(d_1 = y - y_i = k(x_i + 1) + b - y_i\)
\(d_2 = y_{i+1} - y = (y_i + 1) - [k(x_i + 1) + b]\)
兩式相減，得：\(d_1 - d_2 = 2k(x_i + 1) - 2y_i + 2b - 1\)
因為：\(k = \frac{\Delta y}{\Delta x}\)
所以：\(\Delta x (d_1 - d_2) = 2 \Delta y x_i + 2 \Delta y - 2y_i \Delta x + 2b \Delta x - \Delta x\)
又因為：\(\Delta y、\Delta x、b\)對於一條指向來說，是常量
所以：\(\Delta x (d_1 - d_2) = 2 \Delta y x_i - 2 \Delta x y_i + c\)
令：\(\Delta x (d_1 - d_2) = e_i\)，\(e_i\)稱為誤差測量引數
若\(e_i > 0\)，即：\(d_1 - d_2 > 0\)，則實際點更靠近右上方的點（應選用右上方的點）
若\(e_i < 0\)，即：\(d_1 - d_2 < 0\)，則實際點更靠近右側的點（應選用右側的點）
若\(e_i = 0\)，即：\(d_1 - d_2 = 0\)，則隨緣。實際不容易遇到，畢竟\(d_1、d_2\)都是浮點數，相等太難了（這一點參考浮點數的編碼方式就知道了）
現在通過判斷\(e_i\)的符號就可以判斷下一個點是否需要步進了。
那麼，如何去掉判別時的浮點運算呢？即如何確定\(d_1\)和\(d_2\)的值？
不忙，繼續推導。
當前點的誤差測量引數：\(e_i = 2 \Delta y x_i - 2 \Delta x y_i + c\)
下一點的誤差測量引數：\(e_{i+1} = 2 \Delta y x_{i+1} - 2 \Delta x y_{i+1} + c\)
兩式相減，得：\(e_{i+1} - e_i = 2 \Delta y x_{i+1} - 2 \Delta x y_{i+1} - [2 \Delta y x_i - 2 \Delta x y_i]\)
整理，得：\(e_{i+1} - e_i = 2 \Delta y (x_{i+1} - x_i) - 2 \Delta x (y_{i+1} - y_i)\)
又因為：\(x_{i+1} - x_i = 1\)
所以：\(e_{i+1} - e_i = 2 \Delta y - 2 \Delta x (y_{i+1} - y_i)\)
所以：
當選擇右側的點時：\(e_{i+1} = e_i + 2 \Delta y\)
選擇右上角的點時：\(e_{i+1} = e_i + 2 \Delta y - 2 \Delta x\)
可以發現，並不需要確定\(d_1\)和\(d_2\)的值。
根據\(e_i\)的符號可以遞推出下一點的誤差判別引數\(e_{i+1}\)，反過來根據這個新得到的誤差判別引數，可以繼續確定下下一點的誤差判別引數...
遞迴，完美。
但是，初始的\(e_0\)怎麼確定？
對於初始點：
因為：\(\Delta x (d_1 - d_2) = e_i\)，所以\(\frac{e_0}{\Delta x} = d_1 - d_2\)
又因為：\(d_1 - d_2 = 2k(x_0 + 1) - 2y_0 + 2b - 1 = 2kx_0 + 2k - 2y_0 + 2b - 1\)
又因為：\(y_0 = kx_0 + b\)
所以：\(d_1 - d_2 = 2(kx_0 + b) + 2k - 2y_0 - 1 = 2\frac{\Delta y}{\Delta x} - 1 = \frac{e_0}{\Delta x}\)
所以：\(e_0 = 2 \Delta y - \Delta x\)
好了，初始點有了，遞推公式也有了，剩下的就是寫程式了。
至此，改進版的Brehenham演算法全部推導完成。
後面會附上Brehensam演算法繪製直線的C語言程式，可能和這裡的推導過程由出入，但演算法的核心是一樣的。
三、繪製圖形
1. 繪製直線
對於水平直線和垂直直線，大可不必通過演算法去求解，畢竟這兩類直線只在一個方向有步進，而另一個方向步進值始終為0。因此，對於這兩種情況，可以單獨討論。
/**
* @brief :畫線(畫素座標，左上為基點，右下增)
* @note :--
* @param :xStart, 行起始座標(0~127)
yStart, 列起始座標(0~63)
xEnd , 行終止座標(0~127)
yEnd , 列終止座標(0~63)
* @return :void
*
* @date :2016/09/09
* @design :
**/
void OLED_DrawLine(uint32_t xStart, uint32_t yStart, uint32_t xEnd, uint32_t yEnd)
{
int8_t x_width; //x軸寬度
int8_t y_height;//y軸高度
int8_t x_inc; //x方向自增標記
int8_t y_inc; //y方向自增標記
int8_t rem; //current remainder
uint8_t start, end;
uint8_t i;
if(yStart == yEnd)//繪製水平線，horizon line
{
if(xStart > xEnd)
{
start = xEnd;
end = xStart;
}else{
start = xStart;
end = xEnd;
}
for(i=start; i<=end; i++){
OLED_DrawPixelPoint(i, yStart, 1);
}
}else if(xStart == xEnd){//繪製垂直線，vertical line
if(yStart > yEnd)
{
start = yEnd;
end = yStart;
}else{
start = yStart;
end = yEnd;
}
for(i=start; i<=end; i++){
OLED_DrawPixelPoint(xStart, i, 1);
}
}else{//繪製任意直線
x_width = xEnd - xStart;
y_height = yEnd - yStart;
if(x_width < 0) x_width = 0 - x_width;
if(y_height < 0) y_height = 0 - y_height;
x_inc = (xEnd > xStart) ? 1 : -1;
y_inc = (yEnd > yStart) ? 1 : -1;
if(x_width >= y_height)
{
rem = x_width/2;
for(; xStart!=xEnd; xStart+=x_inc)
{
OLED_DrawPixelPoint(xStart, yStart, 1);
rem += y_height;
if(rem >= x_width)
{
rem -= x_width;
yStart += y_inc;
}
}
}else{
rem = y_height/2;
for(; yStart!=yEnd; yStart+=y_inc)
{
OLED_DrawPixelPoint(xStart, yStart, 1);
rem += x_width;
if(rem >= y_height)
{
rem -= y_height;
xStart += x_inc;
}
}
}
}
}
2. 繪製圓
沒有什麼特別的，主要注意利用圓的八分對稱性，可以減少數學運算的次數。
同時使用改進版本，避免了浮點運算。
/**
* @brief :八分對稱法(畫素座標)
* @note :--畫出給定點的八分對稱點(畫圓基礎演算法)
* @param :xc, 圓心行座標
yc, 圓心列座標
x , 給定點
y , 給定點
* @return :void
*
* @date :2017/01/02
* @design :
**/
static void Circle8Point(uint32_t xc, uint32_t yc, uint32_t x, uint32_t y)
{
//直角座標系第一象限x軸開始，逆時針旋轉！
OLED_DrawPixelPoint((xc+x), (yc+y), 1);//1
OLED_DrawPixelPoint((xc+y), (yc+x), 1);//2
OLED_DrawPixelPoint((xc-y), (yc+x), 1);//3
OLED_DrawPixelPoint((xc-x), (yc+y), 1);//4
OLED_DrawPixelPoint((xc-x), (yc-y), 1);//5
OLED_DrawPixelPoint((xc-y), (yc-x), 1);//6
OLED_DrawPixelPoint((xc+y), (yc-x), 1);//7
OLED_DrawPixelPoint((xc+x), (yc-y), 1);//8
}
/**
* @brief :改進畫圓(畫素座標)
* @note :--避免浮點運算(軸上點不突進！)！
* @param :xc, 圓心行座標
yc, 圓心列座標
r , 半徑
* @return :void
*
* @date :2017/01/02
* @design :
**/
void OLED_DrawCircle(uint32_t xc, uint32_t yc, uint32_t r)
{
uint32_t x, y;
int32_t d;//改進，避免浮點運算！
x = 0;
y = r;
d = 3-2*r;
Circle8Point(xc ,yc, x, y);
while(x < y)
{
if(d < 0)
{
d += 4*x+6;
}else{
d += 4*(x-y)+10;
--y;
}
++x;
Circle8Point(xc, yc, x, y);
}
}
3. 繪製橢圓
和圓繪製過程類似，同樣利用了橢圓的對稱性。
/**
* @brief :四分對稱法(畫素座標)
* @note :--畫出給定點的四分對稱點(畫橢圓基礎演算法)
* @param :xc, 橢圓中心行座標
yc, 橢圓中心列座標
x , 給定點
y , 給定點
* @return :void
*
* @date :2017/01/04
* @design :
**/
static void Ellipse4Point(uint32_t xc, uint32_t yc, uint32_t x, uint32_t y)
{
//直角座標系第一象限開始，逆時針旋轉！
OLED_DrawPixelPoint((xc+x), (yc+y), 1);//1
OLED_DrawPixelPoint((xc-x), (yc+y), 1);//2
OLED_DrawPixelPoint((xc-x), (yc-y), 1);//3
OLED_DrawPixelPoint((xc+x), (yc-y), 1);//4
}
/**
* @brief :畫橢圓(畫素座標)
* @note :--
* @param :xc, 橢圓中心行座標
yc, 橢圓中心列座標
a , 半長軸長度
b , 半短軸長度
* @return :void
*
* @date :2017/01/04
* @design :
**/
void OLED_DrawEllipse(uint32_t xc, uint32_t yc, uint32_t a, uint32_t b)
{
int32_t x=0;
int32_t y=b;
int32_t b2=(int32_t)b;
float sqa=a*a;
float sqb=b*b;
float d=sqb+sqa*(-b2+0.25f);
Ellipse4Point(xc, yc, x, y);
while((sqb*(x+1)) < (sqa*(y-0.5f)))
{
if(d < 0)
{
d += sqb*(2*x+3);
}else{
d += sqb*(2*x+3)+sqa*(-2*y+2);
--y;
}
++x;
Ellipse4Point(xc, yc, x, y);
}
d = (b*(x+0.5))*2 + (a*(y-1))*2 - (a*b)*2;
while(y > 0)
{
if(d < 0)
{
d += sqb*(2*x+2)+sqa*(-2*y+3);
++x;
}else{
d += sqa*(-2*y+3);
}
--y;
Ellipse4Point(xc, yc, x, y);
}
}