java用中点画圆法_图形学笔记: 中点画圆法

正圆

首先我们先研究一下圆本身的性质.圆是高度对称的, 在我们实际画图的时候, 我们只需要计算其中八分之一个圆, 然后把另外八分之七由对称性推出来. 我们选择斜率在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

= 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 | |_|

| |

+--------> +-------->

为了把这一点求出來, 我们需要做一点偏导:

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;

}

其实, 文章所描述的算法仍然没有彻底优化, 在很多细节方面还有很多可以修补的地方, 读者可以慢慢思考一下.