« 2008年06月 | メイン | 2008年08月 »

2008年07月 アーカイブ

2008年07月03日

乗除を使わずに円を描く

三角関数や平方根はもちろん、×や÷をも使わずに円を描くアルゴリズムがあることを知った。通称ブレゼンハムのアルゴリズム(Bresenham's circle algorithm)と呼ばれるらしい。

さぞかし恐ろしく複雑な反復演算をするんだろう、乗除を使った方が速いというオチではないか、と思いながらウェブ上で調べてみると、どのサンプルコードを見てもとてもシンプルだった。ここ「伝説のお茶の間」)とかここ(Wikipedia(英語))とかに図解入りの説明とサンプルコードがある。これは速そうだ。

それらのページの図を見て基本的な原理はすぐに理解できたが、説明とソースコードが理解できなかったので、自分で考えてみた。
Javaのデモページへのリンク
ソースコード

基本的な原理は、例えば格子上のXY平面(X座標とY座標は整数のみ)で(x,y)=(r,0)の点から45°分の円弧を描く場合、同じY座標でX軸方向に複数の点を描く必要は無いので、yを1ずつ動かしながら、適当な時にxを1つ動かして、点を描いていけば良い。45°分の円弧が描ければ、そのX軸対称、Y軸対称や90°回転を使って360°分の円弧が作れる。

yを+1していく時、いつxを-1すれば良いだろうか。
円の式はx2+y2=r2であるが、x,y,rを整数に限ると、もちろんこの等式は成立せず、なるべく誤差の小さいxを選んでいくことになる。ここでは左辺と右辺の差e(x,y)=x2+y2-r2を誤差として考える。
仮にyを固定してxを-1すると、誤差は(x-1)2+y2-r2になるので、誤差の差分は2x-1である。つまり、誤差eが2x-1より大きい時は、xを-1してもまだ(x,y)は円の外側なので、xを-1する必要がある。
yが1増えるとeが2y+1増えるので、eに2y+1を足していって2x-1を超えればxを-1すれば良い、と考えられる。そう考えてとりあえず書いたのがサンプルコードのmethod1()である。

なるほど、こういう手があったか。
要領が解ったので、もう少し真剣に考えてみる。

上の方法では点が外側に行き過ぎないようにしてるだけなので、円の内側に点を取った方が円弧に近い場合も、円弧より遠い外側の点を取ってしまう。
そこで、次に(x,y)と(x-1,y)の中点が円弧の内側にあるか外側にあるかで点を選ぶことが考えられる。そのためには、円弧とその中点との誤差e(x,y)=(x-1/2)2-y2-r2が正か負かを判定すれば良い。
このeはyが+1すれば2y+1増え、xが-1すれば2x-2減る。(r,0)から始めると初期値が-r+1/4なので、eは整数にならない。従って、e=0になることはなく、e>0かe<0である。yを+1していって、e>0の時は中点が円の外側にあるから、内側の点(x-1,y)を取れば(xを-1すれば)良く、逆にe<0の時は外側の点(x,y)を取れば(xをそのままにすれば)良い。
ということで書いたのが、サンプルコードのmethod2()である。
なお、このコードは上記のWikipediaにあるコード(通称Bresenham's circle algorithm)と意味的にはほぼ同じであり、全く同じ円を描くことを確認した。


MSXを使っていた頃、マシン語で書いてもSIN()やCOS()や乗除が遅いのに対して、BASICのCIRCLE文の円てどうやってこんなに速く描いてるんだろう、と不思議に思っていたものだが、内部でこういうのが動いてたんだなと思うと、感慨深い。

続きを読む "乗除を使わずに円を描く" »

2008年07月20日

整数演算のみで楕円を描く

前のエントリーでブレゼンハムのアルゴリズムを使って乗除を使わずに整数演算だけで円を描いてみたが、同じ原理を応用して、乗除を使わずに整数演算で楕円を描けないだろうか、と考えてみた。
Javaのデモページへのリンク
ソースコード

横幅が2a、縦幅が2bの楕円を考え、円の時と同様に、(a,0)や(0,b)から始めて、楕円の外側か内側かを判定してポインターを1ずつ縦横斜めに動かしながら、点を描いていく。第2象限~第4象限(90°~360°)はx,yの符号を変えながら第1象限(0°~90°)と同じ曲線を描けば良いので、第1象限だけを考える。

楕円の式は(x/a)2+(y/b)2=1であるが、整数演算のために割り算を避けて、b2x2+a2y2=a2b2とする。
まず、(x,y)=(a,0)からyを+1していって、xは(x,y)と(x-1,y)の中点が実際の楕円の外側になったら-1する(前のエントリー参照)。それを、xの動きがyの動きより大きくならない所、すなわち接線の傾きが135°になる所まで繰り返す。
次に、(x,y)=(0,b)からxを+1していって、同様に(x,y)と(x,y-1)の中点が楕円の外側になったら-1する、というのを、接線の傾きが135°になる所まで繰り返す。

式で書くと、(x,y)と(x-1,y)の中点は(x-1/2,y)なので、楕円の式に当てはめて、b2(x-1/2)2+a2y2 > a2b2なら中点が外側ということになる。E(x,y)を左辺-右辺とすると、E(a,0)=b2(-a+1/4)であり、E(x,y)はyが+1するとa2(2y+1)増え、xが-1するとb2(2x-2)減る。このE(x,y)の1回の増減分をそれぞれdEdy, dEdxと置くと、dEdyはyが+1すると2a2増え、dEdxはxが-1すると2b2減る。a2やb2は先に計算しておけばいいので、これで乗算を使わずに点を動かしていけることがわかる。
(x,y)と(x,y-1)の中点については、今の話のxとyを置き換えれば良い。

ということで書いたのが、サンプルコードのdraw1()である。コード上では上記のa,bはWx,Wyとしている。whileループの中では乗算は使っていない。
デモでは、赤線で描いているのがそれである。


さて、デモでは、比較のために青線でjava.awt.Graphics#drawArc()で楕円を描いて、その上にdraw1()で同じ大きさの楕円を赤線で描いているが、Java VMにも依るだろうが、青線と赤線が結構ずれている。原因として1つ間違いないのは、X軸、Y軸をy=0, x=0に取って上下対象、左右対称にしているので、楕円の横幅や縦幅が偶数の時も奇数しか書いていないことだ。-1<=y<=1で書くと縦幅は3、-2<=y<=2で書くと縦幅は5である。実は前のエントリーのブレゼンハムの円描画のアルゴリズムでも偶数幅の円が描けないという同じ問題があるのだが、コードの簡潔さを重視して敢えて触れなかった(巷のサンプルコードでも大概無視してるし)。
今回はI/F仕様をjava.awt.Graphics#drawArc()に合わせて楕円のバウンディングボックス(外接矩形)のサイズを入力するようにしているので、それに見合った動作をさせたいものである。

横幅が偶数の場合、-a < x < a-1の範囲で描くとすると、横方向の対称軸はx=-0.5となる。x=0とx=-1とでY座標が同じになるのだが、対称軸がx=0の楕円のx=0.5, x=-0.5の時のyをX方向に-0.5ずらして描くと、yの最大点、最小点が無くなってしまう。実際、a,bをそのままにして(x,y)=(a-1/2,0)や(0,b-1/2)からプロットすると、楕円の上下左右の端が十分に膨らまず、少し四角に近くなってしまった。
そもそも、バウンディングボックスを中心に考えると、x=-a, x=a-1の点はそれぞれ必ずxの最小点と最大点なのであり、中心軸がx=-0.5なので、X軸方向の半径がa-1/2の楕円を描かないといけないのであり、半径がaの円を平行移動する方針ではきれいな楕円にならないようだ。

そこで、曲線描画中は、楕円の横幅や縦幅が偶数の時は径を1/2減らして計算する方向で考えてみた。2aや2bが偶数の場合、楕円の式は(b-1/2)2(x-1/2)2+(a-1/2)2y2 = (a-1/2)2(b-1/2)2となり、(x,y)=(a-1,0)から始めてxを-1、yを+1していく時、(x-1,y)と(x,y)の中点が楕円の外側かどうかを判定する式は(b-1/2)2(x-1)2+(a-1/2)2y2-(a-1/2)2(b-1/2)2 > 0 となる。以下、左辺を4倍して考える。左辺の初期値はy=0なので(2b-1)2(-x+1/4)であり、yが+1すると(2a-1)2(2y+1)増え、xが-1すると(2b-1)2(2x-1)減る。以下略。結局、辛うじて乗算を使わずに点を動かしていける。

ということで書いたのが、サンプルコードのdraw2()である。コード上では上記の2a,2bがWx,Wyに対応する。
デモでは、比較のために青線でjava.awt.Graphics#drawArc()で楕円を描いて、その上にdraw2()で同じ大きさの楕円を緑色の線で描いている。JRE 6.0だと、緑線と青線はほとんど差が無い。

続きを読む "整数演算のみで楕円を描く" »

2008年07月26日

Maximaでドラゴンカーブ

MathematicaとかMapleとかの数式処理ソフトは、何であんなに高いのだろう。2008/7/26現在、個人向けライセンス価格はMathematicaは40万円超え、Mapleも30万円超えだ。教育機関にも政府機関にも属していない、日曜数学ファンの私にはとても手が出せない。
会社でも研究職とかでなく、ソフトウェア開発の付帯業務をやってるので、数学は一切使っておらず、買ってもらえる見込みは無い。バグ収束曲線が云々で必要とか言っても、そんな統計処理ならExcelで十分と言われてしまうであろう。
従って、MathematicaやMapleに触れる機会が無い。

そんな折、Maximaというフリーの強力な数式処理ソフトがあることを知ったので、今回ヤボ用があって、インストールして使ってみた。
とりあえず、テーラー展開させてみる。

 KURO-BOX> maxima

Maxima 5.10.0 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1) taylor(sin(x),x,0,10);
                          3    5      7       9
                         x    x      x       x
(%o1)/T/             x - -- + --- - ---- + ------ + . . .
                         6    120   5040   362880
(%i2) taylor(sqrt(r^2-y^2),y,0,10);
                    2      4      6         8       10
                   y      y      y       5 y     7 y
(%o2)/T/       r - --- - ---- - ----- - ------ - ------ + . . .
                   2 r      3       5        7        9
                         8 r    16 r    128 r    256 r
(%i3)
まあ言いたいことは十分わかるが、表示がなんか残念だ。
そこで、満足のいく表示を得るためにTeXmacsを使ってみる。以下の画面はMaxima 5.11.0とTeXmacs 1.0.6.11を使って実際に表示されたものである。%i??という青い部分が入力した部分だ。(表示環境はFreeBSD 6.3+Xorg 1.4.0)

積分もしてみる。

e^(-x^2)を-∞から+∞まで積分すると√πになるのは、手計算でやると面積分にして極座標変換してと技を尽くして偶然解ける感じの有名な例題だが、Maximaは区分求積による近似値計算ではなく、あっさり解析的に計算したようだ。

せっかくTeXmacsを使ってるので、見た目が派手な連分数を表示してみる。

cf(continued fraction)の出力は連分数のリスト表記になっていて、それをcfdisrepすると連分数表記になる。それをratで実数にすると(無限の連分数でないので元々実数だが)約分され、floatで小数にすると元の√2に限りなく近いことがわかる。
蛇足だが、cfdisrepしないと実数や小数に変換できない。

√3についてもやってみる。

方程式を解いてみる。solveで方程式を解析的に解ける。

解けない場合は、xの解にxが混ざったりする。(例は省略)

allrootsを使うと、多項式の方程式を数値的に計算できる。

指数の方程式を渡すと文句を言われる。

find_rootを使うと、(おそらくニュートン法による)近似計算で、解が1つだけ含まれる範囲の解を計算できる。

x^3=3^xのe~100の間の解は当然3なのだが、誤差が出た。


さて、せっかくフリーソフトなので、参考資料も無料のものを、と思って、図書館でMaximaの本を探したら、無かった。その代わり、あるMathematicaの本で、わずか十数行のコードでドラゴンカーブ(ドラゴン曲線)を描くMathematicaのコードを見つけた。マンデルブロー集合とかペアノ曲線とかのフラクタル図形の簡単なやつだ。Maximaの練習とMathematicaとの比較を兼ねて、その本を参考にMaximaで似たようなコードを書いてみた。

dragon(segments_):= block([p, q, r, s],
 p: map(
  lambda([x], [
   [x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]],
   [x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]]
  ]), segments_),
 q: map(
  lambda([x], [
   {[x[1][1], x[1][2]], [x[1][2], x[1][3]]},
   {[x[2][1], x[2][2]], [x[2][2], x[2][3]]}
  ]), p),
 r: map(lambda([x], listify(x)), flatten(q)),
 s: map(lambda([x], [realpart(x), imagpart(x)]), flatten(r)),
 plot2d([discrete, float(s)],
  [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
 r
)$

p: [[[0, 1/2-%i/2], [1/2-%i/2, 1]]]$
for i:1 thru 10 do p:dragon(p)$
gnuplotが使える環境でこれを実行すると、以下のようなドラゴンカーブが描かれる。TeXmacsは必要ない。順に、L字型に配置した2本の線から始めて、全体を90°回転させて先端に繋げる、というのを1回、2回、…10回繰り返したものである。
アルゴリズムとしては、全ての線をL字に折り曲げる、というのを繰り返せばドラゴンカーブと同じ形になるというのを使っている。必要最小限にコードを簡略化して難解になったし、同様のアルゴリズムの説明は多くあると思うので、ここではコードの説明は省く。

続きを読む "Maximaでドラゴンカーブ" »

About 2008年07月

2008年07月にブログ「Weblog on mebius.tokaichiba.jp」に投稿されたすべてのエントリーです。過去のものから新しいものへ順番に並んでいます。

前のアーカイブは2008年06月です。

次のアーカイブは2008年08月です。

他にも多くのエントリーがあります。メインページアーカイブページも見てください。

Powered by
Movable Type 3.35