前回は、真円を描くアルゴリズムを紹介しましたが、これを少し一般化して楕円を描けるようにします。文献に書かれていた内容では、座標の算出を真円のつもりで行い、点を打つ直前に偏平率を掛けることで対応していますが、ここではまじめにアルゴリズムから検討していきたいと思います。
まず、楕円を表す方程式は次式で表すことができます。
(但し a>0, b>0)
または極座標を使って
x | = | a'Rcosθ + x0 |
y | = | b'Rsinθ + y0 |
(但し a'>0, b'>0)
整数演算のみを使って楕円を描くアルゴリズムは真円の場合とほとんど変わりませんが、楕円は直線 y = x, y = -xに対しては対称とならないので、描画範囲を円の 4分の 1(90度分)として考える必要があります。下図のように、原点(0,0)を中心とする楕円 ax2 + by2 = R2を、点(R/√a, 0)から点(0, R/√b)まで時計回りに描いていく場合、描画途中のある点 P(x,y)の次に描くべき点は
真下の点 PV | ( x, y+1 ) |
左下の点 PD | ( x-1, y+1 ) |
真左の点 PL | ( x-1, y ) |
のいずれかになり、この三点のうち真の楕円に近い点を選択すればよいことになります。
点 PV, PD, PLのどれが一番真の円弧に近いかを知るため、真円の場合と同様に、それぞれの点について、円周からの距離の平方と半径 Rの平方の差 Eを求めます。
E = ax2 + by2 - R2 より
PV( x, y+1 )の誤差 | EV = E + 2by + b |
PD( x-1, y+1 )の誤差 | ED = E - 2ax + 2by + ( a + b ) |
PL( x-1, y )の誤差 | EL = E - 2ax + a |
以上3つをを比較して
|EV| - |ED| < 0 | かつ | |EV| - |EL| < 0 | なら | 点 PV | ・・・ (1) |
|ED| - |EL| < 0 | かつ | |EV| - |ED| ≧ 0 | なら | 点 PD | ・・・ (2) |
|EV| - |EL| ≧ 0 | かつ | |ED| - |EL| ≧ 0 | なら | 点 PL | ・・・ (3) |
を選択すればよいことになります。
(1)の場合、不等式は
EV2 - ED2 < 0 | ・・・ (a) |
EV2 - EL2 < 0 | ・・・ (b) |
と等価であるので、それぞれの式を展開すると
a( 2x - 1 )( 2E - 2ax + 4by + a + 2b ) < 0 | ・・・ (a') |
{ a( 2x - 1 ) + b( 2y + 1 ) }( 2E - 2ax + 2by + a + b ) < 0 | ・・・ (b') |
よって
F = 2E - 2ax + 4by + a + 2b | |
G = 2E - 2ax + 2by + a + b |
としたとき、x ≧ 1, y ≧ 0 なら F < 0, G < 0が成り立てば (1)の式を満たします。
同様に、(2)の場合は
ED2 - EL2 < 0 | ・・・ (c) |
EV2 - ED2 ≧ 0 | ・・・ (d) |
と等価であるので、それぞれの式を展開すると
b( 2y + 1 )( 2E - 4ax + 2by + 2a + b ) < 0 | ・・・ (c') |
a( 2x - 1 )F ≧ 0 | ・・・ (d') |
よって
としたとき、x ≧ 1, y ≧ 0 なら F ≧ 0, H < 0が成り立てば (2)の式を満たします。
(3)の場合は
EV2 - EL2 ≧ 0 | ・・・ (e) |
ED2 - EL2 ≧ 0 | ・・・ (f) |
と等価であるので、それぞれの式を展開すると
{ a( 2x - 1 ) + b( 2y + 1 ) }G ≧ 0 | ・・・ (e') |
b( 2y + 1 )H ≧ 0 | ・・・ (f') |
よって x ≧ 1, y ≧ 0 なら G ≧ 0, H ≧ 0が成り立てば(3)の式を満たします。
以上まとめると、
F | = | 2E - 2ax + 4by + a + 2b |
G | = | 2E - 2ax + 2by + a + b |
H | = | 2E - 4ax + 2by + 2a + b |
のとき
F < 0, | G < 0 | ならば | (1) |
F ≧ 0, | H < 0 | ならば | (2) |
G ≧ 0, | H ≧ 0 | ならば | (3) |
となります。
さらに、x ≧ 1, y ≧ 0の条件下では
であることから H < G < Fが必ず成り立つので
F < 0 | ならば | (1) |
F ≧ 0, H < 0 | ならば | (2) |
H ≧ 0 | ならば | (3) |
と条件分岐することで真円の場合と同様に描画が可能となります。(1),(2)の場合、PV, PDを選択したことになるため yが 1増加し、(2),(3)の場合、PD, PLを選択したことになって、xが 1減少します。よって、上記場合分けは次のように変更することができます(F<0ならば H<0が、H≧0ならば F>0が常に成り立つことに注意)。
F ≧ 0 | ならば | xが 1減少 |
H < 0 | ならば | yが 1増加 |
Fは、xが 1減少したとき
F(x-1,y) - F(x,y) | = | ( 2E(x-1,y) - 2a( x - 1 ) + 4by + a + 2b ) - ( 2E(x,y) - 2ax + 4by + a + 2b ) |
= | 2( a( x - 1 )2 + by2 - R2 ) - 2( ax2 + by2 - R2 ) + 2a | |
= | -4ax + 4a |
増加し、yが 1増加すると
F(x,y+1) - F(x,y) | = | ( 2E(x,y+1) - 2ax + 4b( y + 1 ) + a + 2b ) - ( 2E(x,y) - 2ax + 4by + a + 2b ) |
= | 2( ax2 + b( y + 1 )2 - R2 ) - 2( ax2 + by2 - R2 ) + 4b | |
= | 4by + 6b |
増加します。
また Hは、xが 1減少したとき
H(x-1,y) - H(x,y) | = | ( 2E(x-1,y) - 4a( x - 1 ) + 2by + 2a + b ) - ( 2E(x,y) - 4ax + 2by + 2a + b ) |
= | 2( a( x - 1 )2 + by2 - R2 ) - 2( ax2 + by2 - R2 ) + 4a | |
= | -4ax + 6a |
増加し、yが 1増加すると
H(x,y+1) - H(x,y) | = | ( 2E(x,y+1) - 4ax + 2b( y + 1 ) + 2a + b ) - ( 2E(x,y) - 4ax + 2by + 2a + b ) |
= | 2( ax2 + b( y + 1 )2 - R2 ) - 2( ax2 + by2 - R2 ) + 2b | |
= | 4by + 4b |
増加するので
F = F - 4a( x - 1 )
H = H - { 4a( x - 1 ) - 2a }
F = F + { 4b( y + 1 ) + 2b }
H = H + 4b( y + 1 )
のように変化量を加えることで、F, Hを順次求めることができます。
F, Hの初期値 F0, H0は座標値( R/√a, 0 )のときで
となります。
以上をまとめると、楕円 a( x - x0 )2 + b( y - y0 )2 = R2の 4分円を描くアルゴリズムは
1) | x = R/√a, y = 0, F = -2R√a + a + 2b, H = -4R√a + 2a + bで初期化する |
2) | (x,y)に点を打つ |
3) | F ≧ 0であれば x = x - 1, F = F - 4ax, H = H - ( 4ax - 2a ) |
4) | H < 0であれば y = y + 1, F = F + ( 4by + 2b ), H = H + 4by |
5) | x > 0なら 2)に戻る |
楕円は X軸、Y軸に対して対象な図形ですから、あとは点(x,y)と同時に(-x,y),(x,-y),(-x,-y)にも点を打てば楕円全体が完成します。また、中心座標を加えて平行移動させてやれば、任意の点を中心とする楕円が描けるようになります。
以下に、楕円描画のためのサンプル・プログラムを示します。
/* 楕円を描く */ void ellipse( int x0, int y0, int r, int a, int b, int col ) { int x = (int)( (double)r / sqrt( (double)a ) ); int y = 0; double d = sqrt( (double)a ) * (double)r; int F = (int)( -2.0 * d ) + a + 2 * b; int H = (int)( -4.0 * d ) + 2 * a + b; while ( x >= 0 ) { pset( x0 + x, y0 + y, col ); pset( x0 - x, y0 + y, col ); pset( x0 + x, y0 - y, col ); pset( x0 - x, y0 - y, col ); int hIsNeg = ( H < 0 ); if ( F >= 0 ) { --x; F -= 4 * a * x; H -= 4 * a * x - 2 * a; } if ( hIsNeg ) { ++y; F += 4 * b * y + 2 * b; H += 4 * b * y; } } }
F が 0 以上だったとき、H の値が変化するため、先に H の判定結果を変数 hIsNeg へ格納しています ( 以前はこの箇所が誤っており、ご指摘がありました。ありがとうございました )。
◆◇◆更新履歴◆◇◆
◎ 以前のサンプルコードの中に、一部間違いがあったため訂正しました。御指摘どうもありがとうございました(アルゴリズムの説明書いてからソースコード作って間違いに気づき、慌てて説明文を直したんですがまた間違っていたというアホなことをやってしまいました)。
◎ 内容を大幅に見直しました。サンプル・プログラムも大幅に変更しています。
◎ サンプル・プログラムの F, H の条件判定に誤りがあり、修正をしました。ご指摘、ありがとうございました。
前に戻る | タイトルに戻る |