ベジエ曲線を整数の加減算だけで描画する方法

　だいぶ前にちらっと書いたことがあるのだが、かれこれ 20 年ほど前に、ベジエ曲線を整数の加減算だけで描画する方法を考えたことがある。こういう会社員時代の仕事を紹介するためには、いろいろと制約があったりするのだが、このアルゴリズムは多分守秘義務には含まれていないし、技術自体がいい加減時代遅れの技術で、多分今はもっといい方法が開発されているはずなので、公開してもどこからも文句は来ないと思う。だから、機会があればいつか紹介したいと思っていた。

　だだ、このアルゴリズムを説明するには、ややこしい数式を並べる必要があって、それが面倒で避けていたところもあった。しかし、今回 Maxima を使って計算したら、わりと簡単に再現できたので、せっかくだから記録にとどめておくことにする。

　ベジエ曲線というのは、空間からいくつかの点を選択したときに、その点との関係で定義される。たとえば、２次のベジエ曲線だったら、P0、P1、P2 の 3 つの点から、以下のように定義される。

この t は 0 <= t <= 1 の範囲にある実数で、この t　が 0 から 1 に変化するにしたがって、Bezier2(t) の表す座標はある曲線に沿って空間中を移動する。その曲線が P0、P1、P2 によって定義される 2 次のベジエ曲線ということになる。ちなみに、このようにパラメータによって定義される曲線を、一般にパラメトリック曲線という。

　普通、ベジエ曲線を描画するには、t を 0　から 1 まで少しずつ変化させながら Bezier2(t) の座標を計算するのだが、t が実数のままでは、明らかに浮動小数点計算が必要だ。だからまず、パラメータを整数に変換する方法を考える。

　具体的には、区間 [0, 1] を N 等分し、t がそのうちの何番目まで進んだかをパラメータにする。言い換えれば t = i/N (i、N は整数、i=0..N）という変換をほどこすことにする。すると、上の式は以下のようになる。

　これで、P0、P1、P2 が整数座標であれば、この式に出てくる項はすべて整数になるので、iBezier2(i) は一応整数の加減乗除だけで計算できる。ただ、これではまだまだ計算量が多すぎるので、ここからさらに計算量を減らすことを考えよう。

　ここで発想のヒントになったのは、直線描画のためのブレゼンハムのアルゴリズムである。ブレゼンハムのアルゴリズムでは、直線上の各点を最初から計算する代わりに、隣の点の座標との差を計算し、その差を加算していくことにより、計算量を減らしている。

　同じように、iBezier2(i) も、隣の点同士の差を計算すれば簡単になるのかもしれない。この iBezier2(i) のように、整数に対して定義される関数（数列とも言う）の隣り合った点同士の差を計算することを、差分をとると言うが、ここでまず、iBezier2(i) の差分を計算してみる。

　さっきよりは大分すっきりしたが、まだまだ計算量が多い。そこで、もう一回差分をとってみる。

　この式の右辺を見ると、なんと、変数の i がなくなって定数項だけになっていることがわかる。これは実は偶然ではまったくない。元の式が i の 2 次関数だったことを思い出してほしい。2 次関数を 2 回微分すると定数項だけになるのと同じ原理で、2 回差分をとってもやはり定数項だけになるのである。ただし、微分のときほど単純にはならないが。

　さて、このことを利用すると、かなり計算量を減らすことができる。実際に描画する際には、i を 1 ずつ増やしながら iBezier(i) を計算していくことになるわけだが、まず i=0 のときの初期値として iBezier2(0)、dBezier2(0)、d2Bezier2(0) を計算しておく。すると、その次の i=1 の場合、iBezier2(1) は iBezier2(0) に dBezier2(0) を足せば計算できるし、dBezier2(1) は dBezier2(0) に d2Bezier2(0) を足せば計算できる。そして、d2Bezier2(i) は上の通り i に依存しない定数なのだから、d2Bezier2(1) は d2Bezier2(0) と同じである。これを i=2,3,…と続けていけば、定数をどんどん足していくだけで iBezier(i) のすべての点が計算できるはずである。

（実はこの方法は、オイラー法やルンゲクッタ法のような微分方程式の数値解法にも似たところがある。ただ違うのは、あっちは未知関数を微分で近似的に計算しているのにたいし、こっちは既知関数を差分で厳密に計算しているという点である。）

　ただ、ここで気になるのは、d2Bezier2(i) の分母に N の 2 乗があって値が有理数になっていることである。もし、この割り算の結果を切り捨てて整数にして、それを加算していったら、打切り誤差が集積するのでとても実用にならない。

　ところが、上の式をよく見ると、iBezier2(i)、dBezier2(i)、d2Bezier2(i) の分母はすべてが同じ N の 2 乗である。ということは、全体に N の 2 乗を掛けておけば、整数の加減算だけで計算できそうである。そこで、あらためて以下のように定義し直す。

　これで、NBezier2(i) の各点は、整数の加減算だけで計算できる。そしてそれを N の 2 乗で割れば、iBezier2(i)、つまりベジエ曲線上の各点の座標になる。この方法なら、打切り誤差が出るのは最後の割り算のときだけで、誤差が集積することはない。

　また、この N は勝手に決めた定数であるから、誤差が出ないだけの十分な細かさがあれば、どんな整数でもよい。したがって、N が 2 の整数乗になるようにすれば、最後の割り算も、より計算量の少ないビットシフトで計算できる。通常、画面の解像度は 1000～2000 ドット程度であるから、N = 2^10 = 1024 とか N = 2^11 = 2048 ぐらいにしておけば十分であろう。

　改めてアルゴリズムとしてまとめると、NBezier2(i)、dNBezier2(i)、d2NBezier2(i) の初期値が、

で、NBezier2(i)、dNBezier2(i)、d2NBezier2(i) の間の関係が、

であることを利用すれば、NBezier2(i) が逐次的に計算できる。それを ビットシフトすれば各点の座標が求まるというわけだ。

　一応、C++ 風の擬似コードで書くと、こんな感じになる。

（Point クラスの演算子関数は適宜インラインで定義されているものとする。ちなみに、このコードは Point の次元数に依存していないので、Point が何次元でもかまわないし、テンプレート化して任意の次元数に対応することすら可能である。）

void drawBezier2( const Point& p0, const Point& p1, const Point& p2)

{

const int nbit = 10;

const int n2bit = nbit + nbit;

const int N = 1 << nbit;

const int N2 = 1 << n2bit;

Point NBezier2 = p0 * N2;

Point dNBezier2 = p2 + ( 2 * N - 2 ) * p1 + ( 1 - 2 * N ) * p0;

Point d2NBezier2 = 2 * p2 - 4 * p1 + 2 * p0;

for( int i = 0; i <= N; i++ )

{

draw( NBezier2 >> n2bit );

NBezier2 += dNBezier2;

dNBezier2 += d2NBezier2;

}

}

　ループの中が足し算とシフトだけになっているのがおわかりいただけると思う。

　注意してほしいのは、このアルゴリズムは、決して大雑把な近似値を計算しているわけではなく、少なくとも、浮動小数点を使ったアルゴリズムと同等の精度を持っているということである。数値計算にうとい人はピンと来ないかもしれないが、浮動小数点というのは、決して「実数」ではなくあくまで「有理数」に過ぎない。つまり、整数／整数という形式の一つの表現に過ぎないので、整数演算だけでも同等の精度を出すことは可能なのである。

　まあ、今みたいに CPU の速度が速くなってコプロセッサや GPU が進化した時代には無用の技術だろうが、ぼくがこれを開発したのは、8086 とか 80286 とかの時代だったので、このアルゴリズムで描画速度は劇的に速くなった。少なくとも当時ぼくが読んだどの文献にも出ていなかった、オリジナルのアルゴリズムでもある(^^)。また、このアルゴリズムは、2 次のベジエだけではなく、3 次のベジエや B スプラインにも応用できるので、物好きな人は考えてみてほしい。

　ぼくが当時このアイデアを思いついてから、実際に描画処理をインプリメントするまでには相当時間がかかった。特に、ぼくは計算が下手なので、上のゴチャゴチャした差分計算をするのに、紙を何枚も使って何時間も掛けて計算したものだ。

　それが今 Maxima を使えば一瞬で計算できるのだから、隔世の感がある。もし当時こんなソフトがあったら、ぼくももっと斬新なアルゴリズムをバンバン開発できていたに違いないと思うとちょっとくやしい(^^)。そういう意味でも、今のエンジニアは恵まれていると思う。

追記：　今、自分で読み直してみて気づいたけど、このビットシフトで割り算を回避しているところは、もっとスマートな方法があるし、そうすればさらに高速化できますね(^^)。今度ヒマがあったら試してみようと思います。ぼくも 20 年たって少しは頭がよくなってるようで、安心しましたね。みなさんは気づきましたか？