円弧の描画方法

　円（または１／８円弧）の描画アルゴリズムとして
ミッチェナーの円のアルゴリズムが知られています。
しかし、任意の円弧（例えば、長方形に内接するような→MFCのライブラリに
あるような指定方法や３点を指定して円弧を描画する）を高速に描画する
方法はありますか？ミッチェナーの方法の変形（制限）でもいいので
教えてください。

No.1 回答者： ametsuchi 回答日時： 2001/05/09 22:23

ミッチェナーの方法の変形をすればいいのでは？

円を描くのも1/8円弧の組み合わせに過ぎない訳です。で、1/8円弧の特徴として、
■単調性
即ち、xが増加した時、yが増加するか、減少するか何れかで、1/8円弧の範囲内で増減することはないことが挙げられます。

ですから、
１）描くべき円弧の中心と、始角・終角を求めておく。
２）描くべき円弧を分割する。1/8円弧に。
３）以下、各1/8円弧毎、
３．１）始角・終角がこの1/8円弧にかからないなら、ミッチェナーの円弧そのまま
３．２）始角・終角がこの1/8円弧にかかるなら、「単調性」を利用して、開始・終了x、yを抑えておき、ミッチェナーの円弧を生成する際、x,y（＝座標値）が、
３．２．１）開始・終了範囲に入るなら、Pixelを描く。
３．２．２）開始・終了範囲に入らないなら、Pixelを描かない。この際、「描く」-->「描かない」になったら、Loopを抜けてよい。

この回答への補足

ミッチェナーの円を具体的にプログラムにすると
//************************************************************
// DrawCircle - ミッチェナーによる円描画
//
//入力値：
// int xo, yo 円の中心
// int r 円の半径
// int color 円周の色
//************************************************************
void DrawCircle(int xo, int yo, int r, int color)
{
int x, y;

x = r;
y = 0;
while (x >= y){
SetPoint(xo + x, yo + y, color);
SetPoint(xo + x, yo - y, color);
SetPoint(xo - x, yo + y, color);
SetPoint(xo - x, yo - y, color);
SetPoint(xo + y, yo + x, color);
SetPoint(xo + y, yo - x, color);
SetPoint(xo - y, yo + x, color);
SetPoint(xo - y, yo - x, color);

r -= (y << 1) - 1;
if (r < 0){
r += (x - 1) << 1;
x--;
}
y++;
}
}
こうなんですが、ここでSetPointで点を打つときに判定して
打つ、打たないをすれば良いわけですが、どうもややこしく
なってしまいそうです。短いプログラムでできますか。

補足日時：2001/05/14 13:48

No.2 回答者： ametsuchi 回答日時： 2001/05/26 19:38

すいません。

質問がでてたんですね。sekinegooさんは、

１）プログラムの簡潔さ
２）処理効率

どちらを優先させようとしているのでしょう？実はsekinegooさんの目的自体イマイチ分かっていません。

私自身は、図形処理自体は職業として20年以上携わっていますが、ベクターデータが殆どです。ラスターデータは、画像データとか、地形標高データ（＝DEM）位のものです。今回の

・ベクター-->ラスター変換

は、昔某テレビ局向けにPolygonデータを放送用の画像データに変換した位で、「ミッチェナーのアルゴリズム」の本の上で知っているに過ぎないのでシロートです。

我々図形処理屋は、ベクターデータが基本ですから、円弧は３点から求めたり、２線に接する半径Rの円弧を求めたり様々です。２D円弧は、
・中心
・半径
・始角（反時計回り）、これを「FROM角」と名付けましょう。
・終角（反時計回り）、これを「TO角」と名付けましょう。
という「標準形」で持つことが多いです。最近流行のNURBSで表現する手もありますが...。

で、効率を無視するなら、

(1)まず円弧を求める。（上の「標準形」で）
(2)円弧を折れ線近似する。（折れ線の精度が問題！！）
(3)折れ線の各線分をＤＤＡ変換

という手もあります。それならプログラムは簡単ですが、精度を上げるには折れ線を細かく分割しなくてはならず、三角関数計算を辺数に比例して行なわなくてはならないので処理効率が落ちます。処理効率を上げようとすると、今度は精度が失われます。しかも、ミッチェナーのアルゴリムとは何の関係もなくなっちゃいます。

一方、処理効率受重視なら、プログラムが多少複雑でも、
「SetPointで点を打つときに判定して打つ、打たない...」とし、「ややこしくなって」しまうのはやむを得ないのではないでしょうか？

先ず、「FROM角」、「TO角」から、８つの内そのセルどういうセルか分かります。

イ）全く描かないセル
ロ）開始側の一部だけ描くセル
ハ）終了側の一部だけ描くセル
ニ）全部描くセル

イ）とニ）は問題ないです。問題はロ）とハ）です。これも「FROM角」、「TO角」から分かります。というより、

■「FROM角」、「TO角」から、８つのセルがイ）～ニ）どれに該当し、ロ）とハ）の場合、各セル毎の範囲が分かる。

わけです。できれば三角関数はコストのかかる計算ですから、避けたいですね。whileループになってますが、forループで、整数の範囲ではじければ理想です。No.1でも言ったように、1/8セル内では、円弧は単調増加（または減少）ですから、工夫次第で、角度範囲をforループインデックスの範囲（＝整数値）に置き換えることができるはずです。

ではご健闘を！ウソがあったらごめんなさい！

この回答への補足

わざわざ回答していただき、感謝致します。
三角関数は、使いたくないので、整数計算だけで行いたいと思います。
Ｃマガジンに図形描画の記事があって、バックナンバーをすべて購入して
読みましたが、円はあっても円弧は、ありませんでした。インターネットを
探しても、円はあっても円弧はありません。できれば３点円から描画したい
と思いますが、ミッチェナーのアルゴリズムが高速なので、これの変形（点を打つ／打たないの判定を入れる）として
求めたいのですが、ＭＦＣのように範囲指定になってもかまいません。
虫のいい話ですが、Ｃの関数として、欲しいのです。ですから、簡単に
言えば、ＭＦＣのような引数で、実際のソースがあったらいいなと思って
います。インターネットでもＣマガジンでも円のソースはあっても、円弧の
ソースは、どこにもありません。3点指定の円描画はあきらめるとしても
始点終点指定のプログラムがどこかにないかと探しています。
私が作ると、点描画する／しないの判定が、どうしても冗長になって
しまいます。よろしくお願いします。

補足日時：2001/05/26 20:08

No.3 ベストアンサー 回答者： ametsuchi 回答日時： 2001/05/27 18:20

問題を２つに分けましょう。

（１） ３点から円弧を求める方法：

とりあえず約束として、

P0:円弧始点
P1:円弧通過点
P2:円弧終点
O:求める円弧中心
R:円弧半径

とします。時計回りか、反時計回りかはここでは問いません。
線分P0P1、線分P1P2 にOから下ろした垂線の足は線分の中点ですから、線分長を

L01 = dist(P0, P1)
L12 = dist(P1, P2)
として、

<P1-P0, O> - 0.5*L01^2 = 0
<P2-P1, O> - 0.5*L12^2 = 0

です。ここに、<, >は内積を表わします。

これは、Oの成分、Ox,Oyに関する連立一次方程式になり、

det = (P1x-P0x)*(P2y-P1y) – (P1y-P0y)*(P2x-P1x)

として、

Ox = 0.5* {L01^2*(P2y-P1y) – L12^2*(P1y-P0y)} / det
Oy = 0.5* {-L01^2*(P2x-P1x) + L12^2*(P1x-P0x)} / det……………Equ.1)

が求まります。但し、３点の内２点が重なるとか、３点が一直線上に並ぶ場合はdet=0になるので注意が必要です。「円弧は求まらない」としてエラー処理しなくてはなりません。
半径は以下のように容易に求まります。

R = |P0 - O|…………………………………………………………………Equ.2)

始終角を求めます。先ず、

θ0 = atan2(P0y-Oy, P0x-Ox)
θ1 = atan2(P1y-Oy, P1x-Ox)
θ2 = atan2(P2y-Oy, P2x-Ox)…………………………………………………….Equ.3)

として、円弧中心から見た３点の角度が－π～＋πの範囲で求まります。
これを

θ0’ = 0
θ1’ =θ1 -θ0
θ2’ =θ2 -θ0………………………………………………………………………Equ.4)

と始点基準にします。

これを、０～２πの値に変換して、

・θ2’ < θ1’ なら反時計回り
・θ1’ < θ2’ なら時計回り

なのです。あと処理の都合のため、時計回り回りなら始終点を入れ替わると、必ず反時計回りになります。以下、「円弧」と言うと、

・中心（上のOですが、あなたのプログラムでは、「x0, y0」です）
・半径（上のRですが、あなたのプログラムでは「ｒ」です）
・始角（反時計回り）
・終角（反時計回り）

で表現されるものとします。ただ、Equ.3)で分かるとおり、既に、コストのかかる逆三角関数を３回も用いています。

（２） 円弧ベクター－>ラスター変換：

先ず、「ミッチェナーのアルゴリズム」ですが、ご存知のように、

１） 円を８分割し、対称性を利用して基本１パタン（東->北、0°～45°）を利用。
２） 基本パタンでは、y：必ず１Pixelずつ増加、x:０または１Pixel減少。
３） x: ０または－１を判断するのに、「真円からのx方向差分」を巧妙な計算で求める。
ただし、この差分は、x: ０または－１で異なる。

という非常に巧妙な仕掛けになっている訳です。三角関数は勿論、sqrt()も一度も使われていません。しかも近似計算ではなく、誤差は1/2Pixel以内というか、可能な限りでの正確な計算になっています。

本題に入ると、whileループに入る前に、（１）で求めた始終角から、

0) この1/8セルの処理区分， 開始y, 終了y
1) この1/8セルの処理区分， 開始y, 終了y 時計回り
2) この1/8セルの処理区分， 開始y, 終了y 時計回り
3) この1/8セルの処理区分， 開始y, 終了y
4) この1/8セルの処理区分， 開始x, 終了x 時計回り
5) この1/8セルの処理区分， 開始x, 終了x
6) この1/8セルの処理区分， 開始x, 終了x
7) この1/8セルの処理区分， 開始x, 終了x 時計回り

なる８つの行からなるテーブルを作ったら如何でしょうか？8個固定ですから配列でなくてもいいです。

・ 1/8セルの処理区分：
(1) 全域無効、描かない、
(2) 全域有効、描く。
(3) 一部有効、開始～終了範囲だけ描く。

・開始・終了y(x)：必ず開始<終了になります。
０）～７）何れも、最大で＝ｒ、最小で=0.707＊ｒです。ここは気を付けて下さい。

No.2で「forループにして、forループインデックスの範囲」てなこと言いましたが、whieループでいいと思います。「開始終了ｙ（ｘ）」を角度から求めるには始終点位置の角度に関しては残念ながら、三角関数を一度ずつ使わなくてはいけないでしょう。それ以外の場合は、ｒ0.707＊ｒです。三角関数計算は不要です。

No.4 回答者： ametsuchi 回答日時： 2001/05/27 18:31

訂正します。

先程、始終角位置では、三角関数を一度づつ使わなくてはいけないと言いましたが、元々、３点が分かっているので、三角関数は使わなくて済みます。失礼しました。

逆三角関数（＝atan2）を使わなくて出来るかどうか分かりません。上手く場合分けすると、逆三角関数を使わなくて済むかも知れません。

申し訳ないが、後はご自分で考えて下さい。

No.5 回答者： ametsuchi 回答日時： 2001/05/27 18:48

くどくて申し訳ない。

三角関数も使わずに済みそうです。

１）円の中心と半径は「No.3」のとおり求める。

２）「No.3」では、始終角を求めていたが、これの替わりに、３点が８セクタの内どこに落ちるか判定する。これは

・P0-O
・P1-O
・P2-O

のx成分、ｙ成分の大小関係などを調べれば容易に分かる。「O」は円弧中心座標。

３）一番難しいのは、円弧の有効側を調べ、８セクタ各々、「No.3」での３分類をすること。