Bresenham直线插补算法

Bresenham直线插补算法原理与matlab实现

问题说明：近期研究双振镜（四轴）激光加工系统。由于振镜有最小的脉冲当量（即最小的移动单位），需要使用到插补算法对期望直线进行处理。以下是对Bresenham直线插补算法的学习笔记。

参考资料

clever101博客

算法原理

情景：

加工直线起始点$(x_1,y_1)$，终止点$(x_2,y_2)$

理论部分

Bresenham算法是DDA算法画线算法的一种改进算法。本质上它也是采取了步进的思想。不过它比DDA算法作了优化，避免了步进时浮点数运算，同时为选取符合直线方程的点提供了一个好思路。首先通过直线的斜率确定了在x方向进行单位步进还是y方向进行单位步进：当斜率k的绝对值|k|<1时，在x方向进行单位步进；当斜率k的绝对值|k|>1时，在y方向进行单位步进。

令$dy=x_2-x_1$，$dx=y_2-y_1$
该算法是将直线方程转换为一个迭代的过程。迭代过程是：在计长方向上（设x方向），每次坐标必变化一个单位步长（通常单位步长取1），即$x_i=x_{i-1}+1$ (加、减1视dx的正负而定)。
在非计长方向上（设y方向），每次坐标值是否变化一个单位步长视误差项d而定。
$d_i>=0$时，即$y_i=y_{i-1}+1$(加、减1视dy的正负而定)，产生新的误差$d_{i+1}=d_i+2(dy-dx)$
$d_i<0$时，即$y_i=y_{i-1}$，产生新的误差$d_{i+1}=d_i+2dy$。
误差项d是在同计长方向垂直的方向上，是直线的精确路径同实际产生点之间的距离，如下图所示。

误差项d的推导过程如下：
设$dx>dy$（即斜率k<1,x方向上的总行驶路径要大于y方向）且dx,dy>0，单位步长为t（设t=1），当前点$(x_{i-1},y_{i-1})$是离直线最近的点。
计算下一点的坐标值，由下图可知

当$S_i>=T_i$，下一点$(x_{i-1}+1,y_{i-1}+1)$更接近直线。
当$S_i<T_i$，下一点$(x_{i-1}+1,y_{i-1})$更接近直线。

所以
当$S_i>=T_i$，下一点为$(x_{i-1}+1,y_{i-1}+1)$
当$S_i<T_i$，下一点为$(x_{i-1}+1,y_{i-1}+1)$
由相似三角形得：
$S_i/t=dy/dx$
$S_i=(dy/dx)t$
$T_i=t-S_i$
$T_i=t-(dy/dx)t$
因此$S_i-T_i=(2dy/dx)t-t$
可以看出$(S_i-T_i)$的正负决定了坐标值的变化，而与$S_i-T_i$的大小无关。因此$dx(S_i-T_i)/t$并不影响$S_i-T_i$的正负值。
令 $d_i=dx(S_i-T_i)/t=2dy-dx$
则当$d_i>=0$时，走步到点
$x_i=x_{i-1}+1$
$y_i=y_{i-1}+1$
此时产生新的误差项d_{i+1}，如下图所示：

由于 $d_{i+1}=dx(S_{i+1}-T_{i+1})/t$
从相似三角形的关系得：
$(S_{i+1}+t)/2t=dy/dx$
$S_{i+1}=dy/dx\ast 2t-t$
$T_{i+1}=t-S_i+1$
$T_{i+1}=2t-dy/dx\ast 2t$
$S_{i+1}-T_{i+1}=dy/dx\ast 2t-t-2t+dy/dx\ast 2t$
因此$d_{i+1}=2dy-dx-2dx+2dy=d_{i+1}\ast (dy-dx)$
当$d_i<0$时，走步到点x_i=x_{i-1}+1,y_i=y_{i-1}
此时产生新的误差项$d_{i+1}$
$S_{i+1}/2t=dy/dx$
$S_{i+1}=dy/dx\ast 2t$
$T_{i+1}=t-S_{i+1}$
$T_{i+1}=t-dy/dx\ast 2t$
故$d_{i+1}=2dy-dx+2dy=d_i+2dy$

结论：当i=1时，$d_i$的初值为
$d_1=2dy-dx$
当$d_i>=0$时，x,y方向坐标值都将会发生改变
$x_i=x_{i-1}+-1$
$y_i=y_{i-1}+-1$
加减1根据dx,dy的正负而确定。
产生新的误差项$d_{i+1}=d_i+2(dy-dx)$
当$d_i<0$时，y方向坐标值不发生改变
$x_i=x_{i-1}+-1$
$y_i=y_{i-1}$
产生新误差项
$d_{i+1}=d_i+2dy$

可以将上述结论推广为适应于$dy>dx$不同走向的直线。此时以y方向为计长方向，每次坐标值变化一个单位步长（加减视dy的正负而定）
在x方向，坐标值是否变化一个单位步长，视误差项d而定。
当$d_i>=0$时，
$x_i=x_{i-1}+-1$
$y_i=y_{i-1}+-1$
产生新的误差
$d_{i+1}=d_i+2(dx-dy)$
当$d_i<0$时
$x_i=x_{i+1}$
$y_i=y_{i-1}+-1$
产生新误差
$d_{i+1}=d_i+2dx$
当i=1
当i=1时, d的初值为$d_1=2dx-dy$。从该式中可看出，当dx＜dy时，只要将x与y的坐标值进行交换，就可直接使用dx＞dy时的结论来画直线。同时只要用|dy|与|dx|就可以画出不同走向的直线。
下图为用Bresenham算法画直线的流程图。

程序流程图

为了简化Bresenham算法，可把对误差项d的计算与判别式转变为不断比较直线所对的水平距离上的测试值(x，y)，并重复对该测试值进行增加并测试。每当x，y的值增加到一个可能的新点时，就在屏幕上画出这个点。测试值(x，y)的增量与dx，dy有关。下面是简化后的Bresenham算法流程图

设从点$(x_1,y_1)$到点 $(x_2,y_2)$则
$dx=x_2-x_1$,$ix=|dx|$
$dy=y_2-y_1$,$iy=|dy|$
设x,y为测试值，初值为0,0；plot为逻辑量。当plot为真时，画新点；plot为假时，不画新点。(plotx，ploty)为标识一个新点的坐标，初值为$plotx＝x_1$，$ploty＝y_1$。$inc＝max(|dx|,|dy|)$为终点判断条件。依据两种不同的Bresenham算法画直线(2，7）到(7，3)时所生成的数值，以及由直线算法选择的像素所绘制的直线，在所走路径上有些小小的区别，如下图所示。

总结

1. 输入线段的起点和终点。

2. 判断线段的斜率是否存在（即起点和终点的x坐标是否相同），若相同，即斜率不存在，
   只需计算y方向的单位步进（△Y+1次），x方向的坐标保持不变即可绘制直线。

3. 计算线段的斜率k，分为下面几种情况处理

   a. k等于0，即线段平行于x轴，即程序只需计算x方向的单位步进，y方向的值不变

   b. |k|等于1，即线段的x方向的单位步进和y方向的单位步进一样，皆为1。直接循环△X次计算x和y坐标。

4. 根据输入的起点和终点的x、y坐标值的大小决定x方向和y方向的单位步进是1还是-1

5. 画出第一个点。

6. 若|k| <1，设m =0,计算P0，如果Pm>0，下一个要绘制的点为（Xm+单位步进，Ym），
   Pm+1 = Pm -2△Y;
   否则要绘制的点为（Xm+单位步进，Ym+单位步进）
   Pm+1 = Pm+2△X-2*△Y；

7. 重复执行第七步△X-1次；

8. 若|k| >1，设m =0,计算Q0，如果Qm>0，下一个要绘制的点为（Xm，Ym+单位步进），
   Pm+1 = Pm -2△X;
   否则要绘制的点为（Xm+单位步进，Ym+单位步进）
   Pm+1 = Pm+2△Y-2*△X；

matlab实现

% AUTHOR    :CHANDAN KUMAR 
% TITLE     :BRESENHAM LINE ALGORITHM
% SUBJECT   :COMPUTER GRAPHICS AND SOLID MODELLING
% DISCLAIMER:CODE DRAWS A LINE USING BRESENHAM LINE ALGORITHM.

function bresenham_line()
clc
clear all
point = input('Give coord[ x0 y0 x1 y1]: ');

if (abs(point(4)-point(2)) > abs(point(3)-point(1)))       % If the line is steep                                
    x0 = point(2);y0 = point(1); x1 = point(4);y1=point(3);% then it would be converted to 
    token =1;                                              % non steep by changing coordinate
else
    x0 = point(1);y0 = point(2); x1 = point(3);y1=point(4);
    token = 0; 
end
if(x0 >x1)
    temp1 = x0; x0 = x1; x1 = temp1;
    temp2 = y0; y0 = y1; y1 = temp2;
end
dx = abs(x1 - x0) ;                              % Distance to travel in x-direction
dy = abs(y1 - y0);                               % Distance to travel in y-direction
sx = sign(x1 - x0);                              % sx indicates direction of travel in X-dir
sy = sign(y1 - y0);                              % Ensures positive slope line

clf, subplot(2,1,1) ,hold on , grid on ,axis([x0-1 x1+1 y0-1 y1+1]);
title('Bresenham Line Generation Algorithm: Point form')
x = x0; y = y0;                                  % Initialization of line
param = 2*dy - dx ;                              % Initialization of error parameter
for i = 0:dx-1                                   % FOR loop to travel along X
    x_coord(i+1) = x;                            % Saving in matrix form for plot
    y_coord(i+1) = y;
    if (token ==0)                               % Plotting in point form 
        plot(x,y,'r*');                          % For steep line coordinate is again
    else                                         % converted for actual line drawing.
        plot(y,x,'r*');
    end
    param = param + 2*dy;                        % parameter value is modified
    if (param >0)                                % if parameter value is exceeded
        y = y +1*sy;                             % then y coordinate is increased
        param = param - 2*(dx );                 % and parameter value is decreased
        
    end
    x = x + 1*sx;                                % X-coordinate is increased for next point
end
subplot(2,1,2)                                   % Plotting in line fragment form
if (token ==0)
    plot(x_coord,y_coord,'r-','LineWidth',2);
else
    plot(y_coord,x_coord,'r-','LineWidth',2);
end
grid on
axis([x0-1 x1+1 y0-1 y1+1]);
title('Bresenham Line Generation Algorithm: Line fragment form')