QuickDrawはどのように素早く円を描いていたのか?
かつてのMac OS9までの描画エンジンの主役はQuickDrawが担っていた。GUIなOSでは、文字も含めてすべてをグラフィックとして扱うので、画面に見えているすべてのもの*1はQuickDrawによって描かれていたことになる。描画エンジンは、GUIなOS開発の要となる技術である。その出来が、GUIなOS開発の成否を分けるとも言える。
そして、最初期のQuickDrawは、ビル・アトキンソンがたった一人で開発したそうである。
当時(25年以上前)のCPUは、動作クロックが8MHzという性能だった。(現在は2GHz=2000MHzかつ、複数コアが当たり前)
そのような性能であっても、違和感なくマウスで操作できるOS環境にするために、斬新な発想や試行錯誤を重ね、相当な努力の末に開発されたのがLisaやMacintoshであった。
Amazon.co.jp: レボリューション・イン・ザ・バレーの中で、QuickDrawを開発する逸話が少し紹介されている。
当時、研究所の中で開発されていたGUIなOS「Alto」のデモを見たスティーブ・ジョブズとビル・アトキンソンは、その斬新さに感銘を受け、自身でも開発を始めた。
Altoは、実際には専用に開発されたハードウェアの性能に助けられて実現していた環境だった。
にもかかわらず、ビル・アトキンソンは、Altoのすべてはソフトウェアで実現されていると勝手に思い込んで、そのまま開発を続けて、最終的にはQuickDrawを完成させてしまった。
とにかく、素早く描画するために、斬新な発想で様々な工夫がされた。
例えば、円(本当は楕円ルーチンだけどシンプルに円で考える)を描くには普通に考えれば、sin・cosあるいは平方根の計算が必要になる。
x^2 + y^2 = r^2 y^2 = r^2 - x^2 y = √(r^2 - x^2)
しかし、小数を扱うと、浮動小数点計算ユニットを持たない当時のCPUではスピードが問題になった。
そこで、奇数の数列の和が、二乗の数列になる(1 + 3 = 2^2、1 + 3 + 5 = 3^2、1 + 3 + 5 + 7 = 4^2、... )ことを利用して、座標までの距離を整数値を概算しながら描画することを思い付いたそうである。
√(r^2 - x^2)の小数点以下の数値が分からなくても、ピクセルは整数で表現されるので、概算で十分なのである。
賢い!本当にそのアルゴリズムで円を描けるのか、自分でも試してみたくなった。早速やってみる。
基本
- 円を描くには、グラフィック環境が必要になるが、実装を手抜きするため、ターミナル環境にアスキーアートで描くことにする。
- ターミナル環境の縦横比を合わせるために、環境設定... >> 設定 >> テキスト >> フォントの変更...ボタンで、以下の設定にしておいた。
- CourierNew 12pt・文字間隔 1.0・行間隔 0.5
- 表示倍率を適当に縮小して、縦横比が等倍になりそうな大きさにした。
- 言語も手軽さを狙ってRubyを使用する。よって、アルゴリズムの確認は出来るが、実用上のメリットは全くない。
- また、たとえC言語で書いたとしても、高速な浮動小数点計算ユニットを持つ最新のCPUでは普通に小数を含めて計算した方が速いかもしれない。
まずは、最も基本的な円の公式をそのまま計算して円を描いてみた。
module Graph def self.round(r) half = [] quarter = [] (0..r).each do |y| x = Math::sqrt(r*r - y*y).round quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end end puts Graph.round(20)
奇数の数列で√計算
次に、ビルの技で描いてみる。
module Calc def self.root(sqr) n = 0 sum = 0 while sum < sqr n += 1 sum += 2*n - 1 end n end end module Graph def self.round_bill(r) half = [] quarter = [] (0..r).each do |y| x = Calc.root(r*r - y*y) quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end end puts Graph.round_bill(20)
計測
しかし、ここでベンチマークを計測してみると...
user system total real 0.570000 0.010000 0.580000 ( 0.584603) round 7.550000 0.040000 7.590000 ( 7.662687) round_bill
ビルの技はめちゃくちゃ遅い...。
require 'benchmark' ...(中略)... n = 1000 Benchmark.bm do |x| x.report { n.times do; Graph.round(100); end } x.report { n.times do; Graph.round_bill(100); end } end
最適化
いくら最新の浮動小数点計算ユニットが高速であっても、ここまで差がつくのは、たぶん自分の実装が悪いはず。もう少し、工夫する必要がある。
- 平方根を計算する時に、境界値として前回の値も一緒に渡すようにした。
- 境界値から処理することで、1から順に調べる方式より、かなりの無駄なループを削減できる。
- また、ループ中の半径rの2乗は、ループ前に計算しておくことで、無駄な重複計算を減らせる。
require 'benchmark' module Calc def self.root(sqr) n = 0 sum = 0 while sum < sqr n += 1 sum += 2*n - 1 end n end def self.root2(sqr, limit) n = limit sum = limit * limit while sum >= sqr && n > -1 n -= 1 sum -= 2*n + 1 end n + 1 end end module Graph def self.round_bill(r) half = [] quarter = [] (0..r).each do |y| x = Calc.root(r*r - y*y) quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end def self.round_bill2(r) half = [] quarter = [] x = r rr = r*r (0..r).each do |y| x = Calc.root2(rr - y*y, x - 1) quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end end n = 1000 Benchmark.bm do |x| x.report { n.times do; Graph.round(100); end } x.report { n.times do; Graph.round_bill(100); end } x.report { n.times do; Graph.round_bill2(100); end } end
- ベンチマークの結果を見ると、かなりの改善が見られた!(10倍速)
user system total real 0.570000 0.020000 0.590000 ( 0.593901) round 7.510000 0.030000 7.540000 ( 7.556229) round_bill 0.730000 0.020000 0.750000 ( 0.753958) round_bill2
- しかし、それでもなお小数点以下を含めて計算した方が早い...。
円を正確に描く
- ここで、半径rが1から4まで、描いてみる。
- 上がroundで描画した結果
- 下がround_bill2で描画した結果
- 結果が微妙に異なるのだ...。
- 自分の感覚で言えば、以下のように描画される方が良いと思った。(それぞれの一番上が理想の描画)
- OSX標準でインストールされているアプリケーション >> ユーティリティ >> Grapherで描いて確認してみた。
- 例:半径r = 4の場合
- 自分自身の人間的感覚でマーキングすれば、赤枠のピクセルが円周を辿るラインになる、と考えた。
- 今まで、無意識にy=0の時、x=4で、x軸直下の4pxを塗りつぶすと思い込んでいたが、
- 厳密に考えれば、y=0.5の時、x=?と考える必要があったのではないか?(以降、y=1.5、2.5、3.5に対応するx座標を求めるべき)
- yが0から4まで変化する時、境界値は5つあるのに、その間の画素は4つしかない。
- そのジレンマを解消するためには、y=0.5、1.5、2.5、3.5を考えて処理する必要があったのだ。
- たとえば、画素の中心点(2.5, 2.5)は円の内側に存在するが、画素の中心点(3.5, 2.5)は外側にある。
- このように、すべての画素を中心点の座標で判断することで、より正確な円を描画できると考えた。
- しかし、このままでは少数を含む計算をする必要がある。
- せっかくの高速化アルゴリズムも、無駄になってしまう...。
- そこで、半径r=4なら、2倍してr=8と考えて処理すれば良いことに気付いた。
- y=0.5、1.5、2.5、3.5は、y=1、3、5、7に対応する。
- y=5では、x=√(8^2 - 5^2) = √39 = 6.24...、切り捨てて6。
- 半径r=4の円に戻すには、上記値に1を加算して、2で割り、切り捨てて3。
- つまり、y=2.5のピクセルを水平に3つ分塗り潰すのがベストということになるのだ。
以上の処理を実装してみると...
- Calc.root3を追加した
- Graph.round_exactを追加した
#!/usr/bin/ruby require 'benchmark' module Calc # 切り上げて整数にする(遅い) def self.root(sqr) n = 0 sum = 0 while sum < sqr n += 1 sum += 2*n - 1 end n end # 切り上げて整数にする(早い) def self.root2(sqr, limit) n = limit sum = limit * limit while sum >= sqr && n > -1 n -= 1 sum -= 2*n + 1 end n + 1 end # 切り捨てて整数にする(早い) def self.root3(sqr, limit) n = limit sum = limit * limit while sum > sqr n -= 1 sum -= 2*n + 1 end n end end module Graph # 円の公式に忠実バージョン def self.round(r) half = [] quarter = [] (0..r).each do |y| x = Math::sqrt(r*r - y*y).round quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end # ビルの技 低速バージョン def self.round_bill(r) half = [] quarter = [] (0..r).each do |y| x = Calc.root(r*r - y*y) quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end # ビルの技 高速バージョン def self.round_bill2(r) half = [] quarter = [] x = r rr = r*r (0..r).each do |y| x = Calc.root2(rr - y*y, x - 1) quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end # 画素の中心点で正確に計算する & ビルの技 高速バージョン def self.round_exact(r) half = [] quarter = [] r2 = r*2 x2 = r2 rr2 = r2*r2 1.step(r2 - 1, 2) do |y2| x2 = Calc.root3(rr2 - y2*y2, x2) x = ((x2 + 1)/2) quarter << "o" * x + " " * (r - x) end quarter.each {|x| half << x.reverse + x} half.reverse + half end end 0.upto(10) do |i| p i puts Graph.round(i) #puts Graph.round_bill(i) puts Graph.round_bill2(i) puts Graph.round_exact(i) puts "\n"*4 end n = 1000 Benchmark.bm do |x| x.report { n.times do; Graph.round(100); end } x.report { n.times do; Graph.round_bill(100); end } x.report { n.times do; Graph.round_bill2(100); end } x.report { n.times do; Graph.round_exact(100); end } end
試してみる!
- それぞれ3番目が、画素の中心点で正確に計算したバージョン
- かなり良い感じになった気がする。真円に近付いただろうか?
ベンチマークも計測してみた。
- 中心点を考えないときと比べて少し遅くなったけど、それほど悪い結果ではない。(4番目のround_exact)
user system total real 0.570000 0.020000 0.590000 ( 0.605281) round 7.480000 0.030000 7.510000 ( 7.554570) round_bill 0.740000 0.030000 0.770000 ( 0.770146) round_bill2 0.820000 0.030000 0.850000 ( 0.870190) round_exact
ビルを超える技
円を描く方法を調べていると、感動的にシンプルな方法に出会った。
- 円を描く (1)円弧描画のアルゴリズム
- 伝説のお茶の間 No007-09(1) 円の描画(1) MichenerとBresenham
- Midpoint circle algorithm - Wikipedia, the free encyclopedia
ブレゼンハムのアルゴリズムと言うらしい。
- 3時の位置から、1/8の円弧を描くことを考えた場合...
- 始点Pを(x, y)と表現してみる。
- 1/8の円弧の傾きは1以上である。(最終地点で傾きがちょうど1になる)
- よって、上図の矢印のように、次に進むべき座標は真下Pd(x, y + 1)か、左下(x + 1, y + 1)のどちらかになるはず。
- どちらの点に進むべきかは、円の中心からの距離を計算して、より半径rに近い点を選択すべきである。
- より半径Rに近い点を繰り返し選択していくと、最終的にそれは円弧になる。
- 以上の仕組みでコードを実装してみた。
基本
module Graph ...(中略)... def self.round_bresenham(r) half = [] g1 = [] g2 = [] x = r y = 1 while (x >= y) do g1 << "o" * x + " " * (r - x) d0 = r**2 - x**2 - y**2 d1 = r**2 - (x - 1)**2 - y**2 if d0**2 > d1**2 then x -= 1 g2 << "o" * y + " " * (r - y) end y += 1 end g2.pop if g1.size + g2.size > r quarter = g1 + g2.reverse quarter.each {|x| half << x.reverse + x} half.reverse + half end end 1.step(10,1) do |i| p i puts Graph.round_exact(i) puts puts Graph.round_bresenham(i) puts "\n"*4 end
- 早速、実行してみる。
- 上がround_exactで描画したもの
- 下がround_bresenhamで描画したもの
- round_exactは、より正確な円を描くと思いきや...
- 半径7の円については、round_bresenhamの方が、より丸く見える。
最速
上記の基本のままでは、平方根の計算は不要なものの、掛け算使いまくりであまり感動はない。
ところが、可能な限り数式をシンプルにしていくと、以下のコードでも円を描けてしまうのだ!
- 何と!コードの中から掛け算が消えてしまった!
- 加減算とシフト操作だけで円を描いているのだ。
- 2進数で動いているCPUにとって、2倍、4倍という計算は単なるシフト操作になるので、加減算と同等のスピードで処理できる。
- アスキーアートを生成する文字列の掛け算記号は、文字を繰り返し生成する記号で、数値の掛け算とは意味が違う。
module Graph ...(中略)... def self.round_bresenham2(r) half = [] g1 = [] g2 = [] x = r y = 0 f = -2*r + 3 while (x > y) do g1 << "o" * x + " " * (r - x) if f > 0 then x -= 1 f -= 4*x g2 << "o" * (y + 1) + " " * (r - (y + 1)) end y += 1 f += 4*y + 2 end g2.pop if g1.size + g2.size > r quarter = g1 + g2.reverse quarter.each {|x| half << x.reverse + x} half.reverse + half end end 1.step(10,1) do |i| p i puts Graph.round_exact(i) puts puts Graph.round_bresenham(i) puts puts Graph.round_bresenham2(i) puts "\n"*4 end n = 1000 Benchmark.bm do |x| x.report { n.times do; Graph.round(100); end } #x.report { n.times do; Graph.round_bill(100); end } x.report { n.times do; Graph.round_bill2(100); end } x.report { n.times do; Graph.round_exact(100); end } x.report { n.times do; Graph.round_bresenham(100); end } x.report { n.times do; Graph.round_bresenham2(100); end }
- ベンチマークも計測してみると、これまでで最速である。
user system total real 0.570000 0.020000 0.590000 ( 0.587615) round 0.730000 0.020000 0.750000 ( 0.757364) round_bill2 0.810000 0.030000 0.840000 ( 0.832937) round_exact 0.660000 0.020000 0.680000 ( 0.687514) round_bresenham 0.550000 0.020000 0.570000 ( 0.578169) round_bresenham2
掛け算が、加減算とシフト操作になる理由
- 原点O(0, 0)を中心とする半径Rの円を想像する。
- その円周上の3時の位置を始点P(r, 0)と考える。
- yを0、1、2、3...と変化させた時、最も円周に近いxを計算してみる。
- xとyは必ず整数値を取る。
始点P(r, 0)に対してyを1増やした時、
- xはそのままの点P0(x, y)か、または-1減らした点P1(x - 1, y)のどちらかを選択する。
- どちらの点がが円周に近いかは、円の中心からの距離が半径とどれだけ差分があるかを求めて、比較する。
- 半径Rとの差分が小さいほど、円周に近い点ということになる。
- 点P0の差分
d0 = x^2 +y^2 -R^2
- 点P1の差分
d1 = (x - 1)^2 +y^2 -R^2 = (x^2 - 2x + 1) +y^2 -R^2 = x^2 +y^2 -R^2 -2x +1
E = x^2 +y^2 -R^2 とすると、以下のように定義できる。......定義A
d0 = E d1 = E -2x +1
- 差分は正負の整数として計算されるが、絶対値で比較するため、さらにd0、d1を2乗してから比較する。
- d0^2 - d1^2 <= 0 なら点P0を選択し、d1^2 - d1^2 > 0 なら点P1を選択するべきである。
- 因数分解すると、(d0 + d1)(d0 - d1) <= 0 と表現できる。
- 上記数式に定義Aを代入してみる。
(E + (E -2x +1))(E - (E -2x +1)) <= 0 (E +E -2x +1)(E -E +2x -1) <= 0 (2E -2x +1)(2x -1) <= 0
- 両辺を(2x -1)で割る。
2E -2x +1 <= 0
- Eを展開しておく。
2x^2 +2y^2 -2R^2 -2x +1 <= 0
- ところで、上記左辺を関数F(x, y)と定義すると...
F(x, y) = 2x^2 +2y^2 -2R^2 -2x +1 F(x, y + 1) = 2x^2 +2(y + 1)^2 -2R^2 -2x +1 = 2x^2 +2y^2 +4y +2 -2R^2 -2x +1 = 2x^2 +2y^2 -2R^2 -2x +1 +4y +2 F(x - 1, y) = 2(x - 1)^2 +2y^2 -2R^2 + 2(x - 1) +1 = 2x^2 -4x +2 +2y^2 -2R^2 +2x -2 +1 = 2x^2 +2y^2 -2R^2 -2x +1 -4x
- 上記の関数より、yが1増える・あるいはxが1減少すると、関数Fの結果がどう変化するのか計算してみる。
F(x, y + 1) - F(x, y) = 4y + 2 F(x - 1, y) - F(x, y) = -4x
- そして、3時の始点からの最初の分岐点は、関数F(R, 1)を計算すると求められる。
F(R, 1) = 2R^2 +2*1^2 -2R^2 -2R +1 = -2R + 3
- 上記のように、関数F(x, y)を毎回計算すれば、結果を求められる。
- しかし、先に求めたように、以下のように連動することがわかっている。
- yが1増加すると、関数Fの値も 4y + 2 増加する。
- xが1減少すると、関数Fの値も -4x 減少する。
- つまり、初期値 -2R +3 を計算して、そこに 4y +2 または -4x を加減算しても、関数Fの値は求められるのだ。
- 関数Fの結果が0より小さければ点P0、0より大きければ点P1と判定していくと、ほら、立派な円が描けている!
円を描くための、複雑な掛け算や平方根の計算が綺麗になくなってしまった。数学ってすごい!何だか、しみじみ思った。
逸話の続き
QuickDrawの開発逸話には続きがあって、
- 奇数の数列を利用して高速に円を描く発想を得たビル・アトキンソンは、高速に丸や四角を描くデモを見せびらかして喜んでいた。
- しかし、そのデモを見たスティーブ・ジョブズは、全く誉めもせずに「角が丸い四角は描けないのか?」と聞いてきた。
- それを聞いたビル・アトキンソンは、不機嫌そうに「そんなことするのは面倒だし、必要だとも思わない」と答えた。
- その答えにスティーブジョブズは、部屋中のものを指差し角が丸いことを指摘した。
- さらには、ビルを外に連れ出し、街中に角が丸い四角が溢れていることを指摘した。
- スティーブ・ジョブズの説得にギブアップしたビル・アトキンソンは、角が丸い四角を描画するルーチンの開発を約束した。
- 次の日、またしてもビルはデモを見せびらかしていた。
- そのデモは、ほとんど普通の四角と変わらない、猛烈な速さで角が丸い四角を描いていた。
- ビル・アトキンソンはそれをRoudRectと名付けて、QuickDrawに追加した。
- その後、RoundRectはGUIの様々な部分で利用され、必要不可欠な存在になったということだ。
たぶん、ビルはブレゼンハムのアルゴリズムにも気付いたのかもしれない。スティーブ・ジョブズのこだわりが、QuickDrawをより高いレベルに引き上げてしまった。今でこそ角の丸い四角はGUIとして一般的になったが、25年以上も前からRoundRect一発で描けるQuickDrawって素晴らしい!
ブックマークコメント御中
ビルアトキンソンの話だと、MacPaintのソースはKnuth先生のおかげで公開されている。http://slashdot.jp/apple/05/09/15/1331242.shtml
すごいことだと思って、探してみました。しかし、自分では以下の経過までしか見つけることが出来ませんでした...。できることなら、見てみたいです!
ビル・アトキンソンは今は写真家やってるんだっけ。
少し前の記事ですが、写真家やってるみたいですね。
一方FM-7はCIRCLE命令で正64角形を描いた
何と!円に見えていたのは64角形だったとは...。遥か昔の少年時代、自分もnew 7で試行錯誤した記憶が甦りました。
プログラマー現役続行 (技評SE新書) http://bit.ly/bQ1uwl
面白そうです!
ビルのが遅い原因はRubyのメソッド呼び出し? Cで書くとちがいそう。
Math::sqrtもCalcモジュールの中のメソッドに組み入れて、実験してみました。(2番目の結果)
user system total real 0.570000 0.020000 0.590000 ( 0.595689) Math::sqrt 直接実行するコード 0.610000 0.010000 0.620000 ( 0.624908) Math::sqrt 別モジュールのメソッドに含めたコード 0.730000 0.020000 0.750000 ( 0.748067) round_bill2 無駄な繰り返し処理を無くしたコード 0.800000 0.010000 0.810000 ( 0.827443) round_exact 画素の中心座標で比較して、より正確な円を描くコード 0.750000 0.030000 0.780000 ( 0.775632) round_bresenham 二乗演算使いまくりのブレゼンハムのアルゴリズム 0.550000 0.020000 0.570000 ( 0.575092) round_bresenham2 整数の加減算のみのブレゼンハムのアルゴリズム
少しだけ遅くなりましたが、それでもビルの技よりは速いですね。Cではまだ試してません。
もう一つの実装
ブレゼンハムのアルゴリズムを少し視点を変えて実装しているページに出会った。(たいへん参考になります。感謝です。)
- レンダリングサンプル7 - ブレゼンハムのアルゴリズム(円) - Windows Live
- レンダリングサンプル6 - 誤差の最小化を基準とする基本的なアルゴリズムで円を描く - Windows Live
整数値の座標 = 画素の中心 と捉えて実装しているようだ。(自分の場合は、整数値の座標 = 画素と画素の境界)
module Graph ...(中略)... def self.round_bresenham3(r) half = [] g1 = [] g2 = [] x = r y = 0 f = -2*r + 2 while (x > y) do if f > -x then f += -2*x + 1 x -= 1 g2 << "o" * y + " " * (r - y) end #if f < y then f += 2*y + 1 y += 1 #end g1 << "o" * x + " " * (r - x) end g2.pop if g1.size + g2.size > r quarter = g1 + g2.reverse quarter.each {|x| half << x.reverse + x} half.reverse + half end end
- 半径によっては、こちらの方が良い結果(自分の感覚では)になった。(下の方が上記コードで描画した円)
- 数学的には途切れなく連続した状態のアナログな円を、画素で構成するデジタルな円に変換しているのだから、何を優先するかで描画結果は変わってくるのだ。
- ちなみに、半径を指定して円を描く方式(これまでの実験方式)では、その直径は決して奇数にはならない。
- しかし、QuickDraw本来の実装では、四角形の座標を渡して、その四角形に内接する円・楕円を描く仕様だ。
- 四角形の大きさが奇数の場合、円の中心が画素の中心と考えて処理する必要があるのだ。
デジタルな円を描く奥深さを実感している。