Think Stitch
PRINCIPIA  最近の更新


円の描画

線分描画のアルゴリズムを調べるために奥村晴彦先生のアルゴリズム辞典を見ていたら、暗号のようなプログラムが出ていました。 円を描画するプログラムです。

void gr_circle(int xc, int yc, int r, int color)
{
    int x, y;

    x = r; y = 0;
    while (x >= y) {
        gr_dot(xc + x, yc + y, color);
        gr_dot(xc + x, yc - y, color);
        gr_dot(xc - x, yc + y, color);
        gr_dot(xc - x, yc - y, color);
        gr_dot(xc + y, yc + x, color);
        gr_dot(xc + y, yc - x, color);
        gr_dot(xc - y, yc + x, color);
        gr_dot(xc - y, yc - x, color);
        if ((r -= (y++ << 1) - 1) < 0)
            r += (x-- - 1) << 1;
    }
}

解析

どうしてこれで円が描けるのかとても不思議なので,コードを解析してみます。 奥村先生のプログラムに合わせるために、座標平面上の 1/8 象限、x >=0, y >=0, y <= x の範囲で考えます。 x を y の関数 x(y) として考えることにします。 y 座標を 0 から 1 ずつ変えていくと、x は半径 r から単調に減少します。 このとき、曲線の傾きの絶対値は1よりも小さいので、y が1ピクセル進むときに x は変わらないか、あるいは1ピクセルだけ下がります。

ピクセル線分の方程式と同じように考えると、ピクセル (x, y) を描画する条件は次のようになります。

したがって 1/8 円弧のピクセル方程式は次のようになります。 これはついでに書いただけで、今回は使いません。

x(y) は単調に減少するので、式 (1) の左側の不等式だけが問題になります。 y を 1 増やしたときに、この不等式が成り立つならば、x は変わらないということです。 もし成り立たない場合は x は 1 だけ減ります。

平方根は扱いにくいので、代わりに次の関数 h(y) を考えます。 2乗して移項しただけです。

すると x の変換を次のように整理することができます。

まとめると、式 (3) と (4) を使えば、円を書くことができます。 平方根も切り捨て関数も消えてしまいました。 ちょっと不思議です。

h(y) の差分化

式 (3) には 1/2 とか2乗が入っているので、まだ計算が面倒です。 そこでちょっと差分を計算してみます。

x(y+1) と x(y) は式 (4) で結ばれていますから、これを使うと次のように整理できます。

これを見ると、y をインクリメントする毎回のステップで、h(y) を更新するには2乗も 1/2 も要らないことがわかります。

残る問題は h(y) の初期値です。 式 (3) で y=0, x=r とおくと、次のようになります。

残念ながら、ここに 1/4 という半端な数が出てしまいました。惜しい。 そこでえいやっと英断します。この 1/4 を捨ててしまうのです。 すると h(y) の初期値は r です。 さらに式 (4) と (6) には r が入っていません。 r を使うのは x と h の初期化だけです。 そこで、x の初期化が済んだら変数 r を h として流用してしまいましょう。

式 (3) を見ると、一番最初のステップ、つまり y=0 では、次の y の2乗、つまり 1 を h から引かなければなりません。 ところが h の更新ステップでは共通部分として必ず 2y+1 を引きます。 最初は y=0 ですから、この共通部分の更新を先にやってしまえば、初期値の計算を省略できます。 その結果、h つまり変数 r が正または 0 であれば x は変わりません。 r が負の場合は x をデクリメントして、さらに r に 2x-2 を加えます。

ではもう1度奥村先生のプログラムを見てみます。 y は毎回インクリメントするので、y++ となっています。 次に 2y+1 を先に引きます。 2倍はシフトにして (y++ << 1) + 1 ・・・あれ?符号が合いません。 ま、いっか。

        if ((r -= (y++ << 1) - 1) < 0)
            r += (x-- - 1) << 1;

もし r が負だったら x は1つ減ります。これが x-- の部分です。 さらに r を更新します。2x-2 を足す代わりに (x-1)<<1 を足しているだけです。

これでプログラムの謎が解けました。 めでたし、めでたし(?)。

1/4 はどうした?

1/4 を捨てたということは、h(y) に 1/4 を加えているということです。 すると式 (3) から、このプログラムは半径 r ではなく、次の値を半径とする円を書いているということです。 小さい値なので気がつかないんですね。

符号の食い違いは?

上の文章を書いている途中で符号が違うことに気がつきました。 もしやと思って検索したらありました。どうやら間違いだったようです。

こんどこそほんとにめでたし、めでたしでした。

2013/09/02
© 2013,2014,2015 PRINCIPIA Limited