正圆
首先我们先研究一下圆本身的性质.圆是高度对称的, 在我们实际画图的时候, 我们只需要计算其中八分之一个圆, 然后把另外八分之七由对称性推出来. 我们选择斜率在0到-1之间的一段. 有了之前画线段的思路, 我们便可以很显然地想到方法: 利用判别式决定画哪个点. 我们希望判别式的作用是这样的:
f(x, y) > 0: (x, y) is out of the circle
f(x, y) < 0: (x, y) is inside the circle
f(x, y) = 0: (x, y) is on the circle
,--- (x, y)
A V
| _ _
| |_|_| ,___ What if we test where (x+1, y-1/2) is?
| |_| `
|
+-------->
测试两个候选点(x+1, y)和(x+1, y-1)的中点, 如果中点在圆内, 则说明圆弧经过上面的点更多, 反则经过下面的点更多. 现在我们希望得到这样的判别式. 大多数人都应该已经想到了, 这个判别式就是 f(x, y) = x^2 + y^2 - r^2. 这样, 我们就得到了第一条程序:
void circle_1(int x0, int y0, int r, const color& c, image& img)
{
auto f = [&](float x, float y) {
x -= x0; y -= y0;
return x * x + y * y - r * r;
};
for(int x = x0, y = y0 + r; x - x0 <= y - y0; x++) {
img(x, y) = c;
img(y - y0 + x0, x - x0 + y0) = c;
if(f(x + 1, y - 0.5) > 0) y--;
}
}
正圆优化时间
上面的程序是正确的. 但是我们希望能用加减法來减少运算. 思路跟优化画线判别式是非常相似的, 就是找出f(x+1, y-1/2)与f(x, y)之间的递推关系. 现在, 我们静下心來做些数学推导:
f(x + 1, y - 1/2) = x^2 + 2x + 1 + y^2 - y + 1/4 - r^2
= f(x, y) + 2x + 1 - y + 1/4
这样就可以由 上一个点 f(x, y)得出当前的判别式了. 但是, 当我希望继续递推下去的时候, 我发现我们并没有得到判别式所需的 当前的点 的判别式f(x+1, y)或f(x+1, y-1). 所以, 我们要把它求出来:
f(x + 1, y) = f(x + 1, y - 1/2) + y - 1/4 <------- Not integral
= f(x, y) + 2x + 1
f(x + 1, y - 1) = f(x + 1, y - 1/2) - y - 1 - 1/4
= f(x, y) + 2x - 2y
有了这些式子, 于是我们的算法就出来了:
-
f(0, r) = 0是第一个值. - 求
f(1, r - 1/2), 是大于0还是小于0. - 如果大于0, 那么画下面的点, 否则画上面的点.
- 如果画的是下面的点, 我们就应该将
f由f(1, r - 1/2)修改成f(x+1, y-1), 否则修改成f(x+1, y) - 现在可以返回第二步, 來处理下一个点的事情了
等等! 在我们开始写代码之前, 还有一个问题要处理: 1/4 不是整数. 所幸, 这几乎影响不到什么. 1/4是一个不足1的数, 也就是说, 它丝毫影响不到f的正负性(如果f是-1, 加上1/4仍然小于0, 如果f是0, 加上1/4就大于0了). 到了后面, 1/4 恰好都是要被减去的. 所以, 直接当它不存在就好了:
void circle_2(int x0, int y0, int r, const color& c, image& img)
{
for(int x = x0, y = y0 + r, f = 0; x - x0 <= y - y0; x++) {
int dx = x - x0, dy = y - y0;
img(x, y) = c;
img(dy + x0, dx + y0) = c;
f += dx + dx - dy + 1;
if(f > 0) {
f += -dy - 1;
y--;
} else f += dy;
}
}
正圆渲染结果
int main()
{
image img(800, 600);
circle_1(300, 300, 250, color(uint32_t(0x00000000)), img);
circle_2(300, 300, 220, color(uint32_t(0x00000000)), img);
ofstream f("3-5-circle.ppm");
f << img;
f.close();
}
椭圆
顺着这个思路, 我们可以想想椭圆该怎么画了. 圆可以八分对称, 但是椭圆只可以四分对称. 所以我们必须将4个象限的内容都算出来. 但是同时, 我们发现原先的中点测试不管用了, 因为它只能应对斜率在[0, -1]部分. 因此, 我们需要另一个迭代, 來应对斜率在[-1, -inf]的部分. 在此之前, 我们需要先知道椭圆的公式和判别式f(x, y)是怎样的.
2 2
x y
F = --- + --- - 1 = 0
2 2
a b
f(x, y) = b^2 x^2 + a^2 y^2 - a^2 b^2,
跟圆几乎没有差别, f(x, y)的性质也是一样的. 现在我们來讨论椭圆在一个象限里面的上半支和下半支. 比如说上半支, 我希望它从(0, y0)出发, 到某一点停止, 而下半支, 则从(x0. 0)出发, 到同一点停止. 这时, 求出这一点就变得相当关键. 考虑这一点是什么 -- 没错, 这一点就是斜率为-1的那点.
A
|(0, y0) ,----- dy/dx = -1
+--------____ |
| ``--.V
| #
| `\
| \
| |
+---------------------+------->
(x0, 0)
,--- (x, y)
A V A ,---- test this while
| _ _ | _V_ dy/dx < -1
| |_|_| ,_ test this while | |_|_|
| |_| ` dy/dx > -1 | |_| <--- (x, y)
| |
+--------> +-------->
为了把这一点求出來, 我们需要做一点偏导:
dy ∂F/∂y b^2 x
---- = ------- = - ------- <= -1
dx ∂F/∂x a^2 y
x <= a^2/b^2 y
y <= b^2/a^2 x
如此以来, 我们便有了起始点的凭据, 可以书写一段代码了:
void eclipse_1(int x0, int y0, int a, int b, const color& c, image& img)
{
const double bsq = b*b, asq = a*a, absq = bsq * asq;
const double sepx = asq/bsq, sepy = bsq/asq;
auto f = [&](float x, float y) {
x -= x0; y -= y0;
return bsq * x * x + asq * y * y - absq;
};
for(int x = x0, y = y0 + b; x - x0 <= sepx * (y - y0); x++) {
img(x, y) = c;
if(f(x + 1, y - 0.5) > 0) y--;
}
for(int y = y0, x = x0 + a; y - y0 <= sepy * (x - x0); y++) {
img(x, y) = c;
if(f(x - 0.5, y + 1) > 0) x--;
}
}
椭圆优化时间
相似地, 我们优化时间需要有下列这些式子: 上一点到中点的关系式, 中点到当前点的关系式(这里有两条), 所以总共是6条:
f(x + 1, y - 1/2) = f(x, y) + 2 b^2 x + b^2 - a^2 y + a^2 /4 (i
f(x - 1/2, y + 1) = f(x, y) - b^2 x - b^2 /4 + 2 a^2 y + a^2 (ii
f(x + 1, y - 1) = f(x, y) + 2 b^2 x + b^2 - 2 a^2 y + a^2
= f(x + 1, y - 1/2) - a^2 y - a^2 / 4 + a^2; (iii
f(x - 1, y + 1) = f(x, y) - 2 b^2 x + b^2 + 2 a^2 y + a^2
= f(x - 1/2, y + 1) - b^2 x - b^2 / 4 + b^2; (iv
f(x + 1, y) = f(x, y) + 2 b^2 x + b^2
= f(x + 1, y - 1/2) + a^2 y - a^2 / 4; (v
f(x, y + 1) = f(x, y) + 2 a^2 y + a^2
= f(x - 1/2, y + 1) + b^2 x + b^2 / 4; (vi
同样地, a^2/4之类的部分, 因为小数部分对f(x+1, y-1/2)的正负性没有影响, 并且在后继的当前点的判别值的时候会被减去抵消, 所以并不必特意去校正它. 这样, 我们可以得到最终的程序是:
void eclipse_2(int x0, int y0, int a, int b, const color& c, image& img)
{
long long int asq = a * a, bsq = b * b,
asq_4 = asq / 4, bsq_4 = bsq / 4,
asq_y, bsq_x, f;
f = 0, asq_y = asq * b, bsq_x = 0;
for(int x = 0, y = b; bsq_x <= asq_y; x++, bsq_x += bsq) {
img(x + x0, y + y0) = c;
f += bsq_x + bsq_x + bsq - asq_y + asq_4; // (i
if(f < 0) f += asq_y - asq_4; // (v
else {
f += - asq_y - asq_4 + asq; // (iii
y--;
asq_y -= asq;
}
}
f = 0, asq_y = 0, bsq_x = bsq * a;
for(int x = a, y = 0; bsq_x >= asq_y; y++, asq_y += asq) {
img(x + x0, y + y0) = c;
f += asq_y + asq_y + asq - bsq_x + bsq_4; // (ii
if(f < 0) f += bsq_x - bsq_4; // (vi
else {
f += - bsq_x - bsq_4 + bsq; // (iv
x--;
bsq_x -= bsq;
}
}
}
椭圆渲染结果
int main()
{
image img(800, 600);
eclipse_1(300, 300, 250, 200, color(uint32_t(0)), img);
eclipse_2(300, 300, 270, 220, color(uint32_t(0)), img);
ofstream file("3-6-eclipse.ppm");
file << img;
file.close();
return 0;
}
其实, 文章所描述的算法仍然没有彻底优化, 在很多细节方面还有很多可以修补的地方, 读者可以慢慢思考一下.
