[计算机图形学经典算法] 直线段和圆弧在屏幕上的绘制 （附matlab代码）

刚学习了计算机图形学这门课程，为奠定根基的算法所倾倒，特此记录一二。

直线—中点 Bresenham 算法

DDA算法在效率上较低的原因是需要计算 k，并以之作为累加项。一个直观的改进方式，是在整个运算过程中将涉及到的数值乘以 dx (或dy)，转化为整型进行运算。

中点 Bresenham 算法采用一种不同的观点来解决这一问题－判别式。我们先考虑一般的直线方程，在图形学中我们一般给出的已知条件是两个端点：

此判别方程的意义，在于直线上的所有点均使 F(x,y)=0。

现在，我们考虑离散的情形，在|k| < 1时，为判断如下图所示的与P相邻的像素点Pd和Pu中的哪一个归属于直线，我们取两者的中点M作为判决依据。将其坐标代入判别式
(1) F(M) < 0, 取上方点 Pu
(2) F(M) > 0, 取下方点 Pd

因此，若(Xi, Yi)为直线上的点，在判定下一个点(Xi＋1,Yi+1)并在(Xi＋1,Yi+1)和(Xi＋1,Yi)中进行选择时，可以根据中点坐标(Xi＋1,Yi+0.5)代入判别式进行判定。
F(Xi,Yi) = a Yi + b Xi + c = 0
F(Xi＋1,Yi+0.5) = a Yi + 0.5a + b Xi + b +c
= 0.5a + b
(1) F(M) = 0.5a + b < 0, 取上方点 Pu (Xi＋1,Yi+1)
(2) F(M) = 0.5a + b > 0, 取下方点 Pd (Xi＋1,Yi)

两个问题：

(1) F(Xi,Yi)=0 不一定总能满足，如何解决？
(2) 尽管 a, b, c 均为整数(?)，仍有一个0.5为小数。
问题(2)容易解决，仅需在计算过程中乘2即可；
问题(1)则较为麻烦。我们知道在起始点时，F(Xi,Yi)=0 总是成立的，在后续计算中我们将偏差 F(Xi,Yi)=di 累计表达为 F(Xi,Yi)-di=0，再累积起来通过一个变量来考虑的。然而，因为前面的推导基于 +1(x方向) 和 +0.5(y方向)，所以对该变量要进行调整。具体通过递推完成。

递推式

因为不能保证各像素点处的F(Xi,Yi)=0，因此只能从起始点处考虑 F(X0,Y0)=0，然后依次考虑F(X0＋1,Y0＋0.5), F(X0+2,Y0+?),… 此处我们可以看出，必须要分两种情况。
当 F(X0＋1,Y0＋0.5) < 0，选择(Xi＋1,Yi+1)时，再判断的下一个中点应为(Xi＋2,Yi+1.5)；
当 F(X0＋1,Y0＋0.5) > 0，选择(Xi＋1,Yi)时，再判断的下一个中点应为(Xi＋2,Yi+0.5)；
注意每一步只有两种可能，所以这个模式是重复出现的，换句话说可用迭代的方式实现两个中点间的递推关系。

累积误差项的递推（di+1 < 0）

累积误差项的递推关系，反映了相邻中点间的递推关系。

补充细节:

算法整体描述:

• 计算初始值△x、△y、d=△x-2△y、x=x0、y=y0。
• 绘制点(x,y)。判断d的符号。若d<0，则(x,y)更新为(x+1,y+1)，d更新为d+2△x-2△y；否则(x,y)更新为(x+1,y), d更新为d-2△y。
  (这种递推关系类似数学归纳法，相邻中点间构建关系)
• 当直线没有画完时，重复上一步骤，否则结束。

例题：

matlab代码

% P,Q为直线段端点，pixs_line为计算出的像素点
function pixs_line = bresenham_line(P,Q)

% 判断斜率 k 是否绝对值小于1
% 如果斜率大于1, X,Y坐标互换
if abs(Q(2)-P(2)) > abs(Q(1)-P(1))  
    P = [P(2) P(1)];  Q = [Q(2),Q(1)];
end

% 初始化变量
n = abs(Q(1)-P(1)) + 1; 
pixs_line = zeros(n,2);
pixs_line(1,1:2) = [P(1) P(2)];

%中点 Bresenham 算法扫描转换
dx = Q(1) - P(1);  dy = Q(2) - P(2);
d = dx - 2*dy;  % 计算中点
for i = 2:n
    pixs_line(i,1) = P(1) + i - 1;
    if d < 0
        pixs_line(i,2) = pixs_line(i-1,2) + 1;
        d = d + 2*dx - 2*dy;
    else 
        pixs_line(i,2) = pixs_line(i-1,2);
        d = d - 2*dy;
    end
end

% k 的绝对值大于1时的特殊处理
if abs(Q(2)-P(2)) > abs(Q(1)-P(1))
    pixs_line = fliplr(pixs_line);
end

圆形—中点Bresenham算法

中点Bresenham画圆算法的主要想法，仍然是构造判别式，然后从起始点开始，依次代入邻近点的的像素中点，并根据判别式的值判定邻近点的归属。

下面列出中点 Bresenham 画圆算法的几个要点：

• 构造判别式；
• 记录一个误差项 di；
• 建立递推式，分别针对取上、下两个像素的不同情形；
• (运算过程的整数化，其它一些算法细节)

构造判别式

误差项 di

构造递推式

算法描述

1. 输入圆的半径 R。
2. 计算初始值d=1-R、x=0、y=R。
3. 绘制点(x,y)及其在八分圆中的另外七个对称点。
4. 判断d的符号。若d < 0，则先将d更新为d+2x+3，再将(x,y)更新为(x+1,y)；否则先将d更新为d+2(x-y)+5，再将(x,y)更新为(x+1,y-1)。
5. 当x < y时，重复步骤3和4。否则结束。

例题

function [X,Y]=Bresenhamcircle(x0,y0,r)
X=ones(1,1000);
Y=X;                                      %坐标向量，用于存储绘制的点的坐标
x1=X;x2=X;x3=X;x4=X;x5=X;x6=X;x7=X;x8=X;
y1=Y;y2=Y;y3=Y;y4=Y;y5=Y;y6=Y;y7=Y;y8=Y;  %将圆对称划分为8个部分，分别用xi，yi记录坐标
D=ones(1,1002);                           %判别向量
i=1;
x1(1)=x0;x2(1)=x0+r;x3(1)=x2(1);x4(1)=x0;x5(1)=x4(1);x6(1)=x0-r;x7(1)=x6(1);x8(1)=x1(1);
y1(1)=y0+r;y2(1)=y0;y3(1)=y2(1);y4(1)=y0-r;y5(1)=y4(1);y6(1)=y0;y7(1)=y6(1);y8(1)=y1(1);     %初始条件
xd=2^(1/2)/2*r+x0;
while x1(i)<xd
    i=i+1;
    x1(i)=x1(i-1)+1;y2(i)=y2(i-1)+1;y3(i)=y3(i-1)-1; x4(i)=x4(i-1)+1; x5(i)=x5(i-1)-1;y6(i)=y6(i-1)-1;y7(i)=y7(i-1)+1; x8(i)=x8(i-1)-1;         %计长方向总加1
    D(i)=2*(y1(i-1)-y0-(r^2-(x1(i)-x0)^2)^(1/2))-1;
    if D(i)<0
        y1(i)=y1(i-1);x2(i)=x2(i-1);x3(i)=x3(i-1);y4(i)=y4(i-1);y5(i)=y5(i-1);x6(i)=x6(i-1);x7(i)=x7(i-1); y8(i)=y8(i-1);
    end
    if D(i)>=0
        y1(i)=y1(i-1)-1;x2(i)=x2(i-1)-1;x3(i)=x3(i-1)-1;y4(i)=y4(i-1)+1;y5(i)=y5(i-1)+1;x6(i)=x6(i-1)+1;x7(i)=x7(i-1)+1; y8(i)=y8(i-1)-1;
    end
end
X(1:i)=x1(1:i);X(i+1:2*i)=x2(1:i);X(2*i+1:3*i)=x3(1:i);X(3*i+1:4*i)=x4(1:i);X(4*i+1:5*i)=x5(1:i);X(5*i+1:6*i)=x6(1:i);X(6*i+1:7*i)=x7(1:i);X(7*i+1:8*i)=x8(1:i);
Y(1:i)=y1(1:i);Y(i+1:2*i)=y2(1:i);Y(2*i+1:3*i)=y3(1:i);Y(3*i+1:4*i)=y4(1:i);Y(4*i+1:5*i)=y5(1:i);Y(5*i+1:6*i)=y6(1:i);Y(6*i+1:7*i)=y7(1:i);Y(7*i+1:8*i)=y8(1:i);
X=X(1:8*i);Y=Y(1:8*i);
a=0:pi/2000:2*pi;
x=x0+r*cos(a);
y=y0+r*sin(a);
plot(x,y,'k');
legend('实际图像');hold on;
plot(X,Y,'r.');
hold off;
end