ブレゼンハム・アルゴリズムの原典を読み、理解する - WebGLでレトロPCグラフィックスを楽しむ（２３）

公開：2016-05-06 17:34
更新：2020-02-15 04:37
カテゴリ:PCグラフィックス,javascript,webgl,html5,WebGLでレトロPCグラフィックスを楽しむ

さて、ブレゼンハム・アルゴリズムである。

ブレゼンハム・アルゴリズムは整数化・加減算およびシフト演算のみで実装できるということでよく紹介されているけれども、ブレゼンハム・アルゴリズムやDDAの文献を読んでみて、実はそこがポイントではないのではないかと考えるようになった。

点が、ブレゼンハム・アルゴリズムの優れたところなのだろうと今は思っている。思っているだけで、検証できてないけど。。

ブレゼンハムさんが、このアルゴリズムを実装したのは1962年、IBM 1401コンピュータ上だそうだ。ソースコードは下記である。これはIBM 1401のAutoCoderというアセンブラ言語で記述されている。

http://www.mirrorservice.org/sites/www.bitsavers.org/1401/progs/bresenham/bresenham.s

このコード、何とか読めそうで読めない感じ。。

このアルゴリズムは1965年にIBM System Journalで発表された。当時の記事はPDF化されていて読むことが可能だ。このアルゴリズムはプロッタの描画アルゴリズムとして考案された。

Algorithm for computer control of a digital plotter

でこの記事を参考にしながら、ブレゼンハム・アルゴリズムを説明してみることにする。といってもほとんど原文をなぞっているだけだが。。

その前にプロッタというものに触れておこう。 プロッタというのは下の写真のような形をしていた。ペンによってグラフィックを紙に印刷するのである。今のようにディスプレイが高精細でなかったころ、精細なコンピュータ・グラフィックスを表現するにはプロッタに頼るしかなかった。当時は非常に重宝したのだろうし、CADの図面出力用としてはつい最近まで使われていた。現在は大判のインクジェット・プリンターに移行しているものと思うが。

(Wikipediaより）

当時のプロッタは8方向に線形に移動可能で、精度は1/100inchだった。つまりは100DPIか。

プロッタは真ん中の点を現在位置とするとM1～M8の8方向に移動可能である。プロッタの位置指定は1/100inch単位でx,yの直交座標で指定する。当時のプロッターは曲線は線分で近似していた。曲線上に細かく点を配置して、その間を線分で結ぶのである。

続いて図の曲線中の線分D1D2の描画について見ていく。

図２の点D1(x1,y1),D2(x2,y2)を拡大表示したのが下の図３となる。升目の一ますは、1/100 inchである。

プロッターは8方向にしか動けないので、D1-D2間を結ぶ線分通りには描画できず、太線のように8方向の動きで近似することになる。例えば点Qと点Rどちらを選ぶかといえば、点Qのほうがより理想的な線分に近いので、点Qまでの線分を描画する。 で、気が付くのはD1D2の線分が、第一Octantに属する場合、D1からD2への描画はM1,M2の動きのみとなっている点だ。Octantというのは以下の図で表される区画のことである。適切な日本語がわからないのでOctantという言葉のまま進めることにする。 で話を戻すと、続いて線分D1D2にできる限り近い形にM1,M2を選んでいく方法を考えていく。

方法を考える前に、わかりやすくするため下の図４のように線分D1D2をD1を原点としたときのa,b座標系で表すことにする。この座標系ではD2座標は (Δa,Δb) = (x2 - x1,y2 - y1) となる。

で、図4中の P_{i-1} 点が処理済みであるとして、次の動きを考えると、

  1. もし r_i < q_i ならば M1( R_i を選択)
  2. もし r_i \geq q_i ならば M2( Q_i を選択)

の2つのうちのいずれかである。これを少し変形すると

  1. もし r_i - q_i \geq 0 ならば、M1( R_i を選択)
  2. もし r_i - q_i < 0 ならば、M2( Q_i を選択)

となる。

図を見てもらうとわかるけど、線分D1D2と q_i,q_i' で囲まれた三角形と線分D1D2と r_i,r_i' で囲まれた三角形は相似であるので、 r_i' - q_i' の符号は r_i - q_i の符号と同じになる。つまり r_i' - q_i' が0以上か0未満かでM1かM2を選んでもよいということになる。 さらに線分D1D2が第一Octant内の場合、 {\Delta}a \ > \ 0 なので、 \bigtriangledown_i = (r_i' - q_i'){\Delta}a は r_i - q_i と同じ符号である。つまり \bigtriangledown_i が0以上か0未満かでM1かM2を選んでもよいということになる。 なんか突然この式が出てきたけど（原文でもいきなり登場する）、この式はM1かM2のいずれかの動きを選択するための式となる。 しかしこの \bigtriangledown_i = (r_i' - q_i') はまあわかるとして、なぜこの式に \Delta{a} をかけるのだろうかね？それは後程わかることになる。

そして、 \bigtriangledown_i は次の漸化式(1)を満たす。なぜ満たすのかは後ほど検証する。

\left .\begin{align} \bigtriangledown_1 \\ &= 2{\Delta}b - {\Delta}a\\ \bigtriangledown_{i+1} = \left\{ \begin{array}{I} \bigtriangledown_i + 2{\Delta}b - 2{\Delta}a \\ & {\rm if} \\ & \bigtriangledown_i \geq 0 \\ \bigtriangledown_i + 2{\Delta}b \\ & {\rm if} \\ & \bigtriangledown_i < 0 \end{array} \right. \end{align} \right\} \quad \quad (1) \\ {\rm ただし} \\ {\Delta}a = x_2 - x_1 , {\Delta}b = y_2 - y_1 \quad \quad (2)

(1),(2)式によって計算される \bigtriangledown_i の値は、プロッターの移動先を決定するために使える。

{\rm if} \quad \left\{ \begin{align} \bigtriangledown_i \ < \ 0 \ , \quad {\rm execute} \ m_1 \\ \bigtriangledown_i \ \geq 0 \ , \quad {\rm execute } \ m_2 \end{align} \right\} , \quad i = 1 , \cdots , {\Delta}a, \quad (3) \\ {\rm ただし} \\ m_1 = M_1 \quad {\rm and } \quad m_2 = M_2.\quad(4)

結論として、(1),(2),(3),(4)式は第一Octantのためのアルゴリズムを構成する。が、他のOctantについては(2)と(4)式の右辺を修正しなければならない。 このことについては後述するとして、その前に先ほどの漸化式の検証を行うことにする。

まず図4の表記について説明する。

続いて漸化式(1)の検証をする。 先ほどの図4の表記の内容に基づいて \bigtriangledown_i の式を書き換えることにする。

\begin{align} \bigtriangledown_i &= (r_i' - q_i')\Delta{a} \\ &= \{(b_i - {\rm floor}(b_i)) - ({\rm ceil}(b_i) - b_i)\}\Delta{a} \\ &= 2b_i\Delta{a} - ({\rm floor}(b_i) + ({\rm ceil}(b_i))\Delta{a}. \end{align}

ここで b_i = (\frac{\Delta{b}}{\Delta{a}})a_i なので、

\begin{align} \bigtriangledown_i &= 2(\frac{\Delta{b}}{\Delta{a}})a_i\Delta{a} - ({\rm floor}(b_i) + ({\rm ceil}(b_i))\Delta{a} \\ &= 2a_i\Delta{b} - ({\rm floor}(b_i)+{\rm ceil}(b_i))\Delta{a}. \end{align}

上の式でなぜ \bigtriangledown_i = (r_i' - q_i'){\Delta}a としたのかがわかるだろう。 {\Delta}a を掛けておくことによって、分母が払われ、整数化できるからである。 そうしても結果には影響がない（符号に変化がない）ことは前に説明したとおりである。

D1D2は第一Octantに属するため、 a_i = a_{i-1} + 1 である。 さらに {\rm floor}(b_i) = {\dot{b}}_{i-1} + 1 、{\rm ceil}(b_i) = {\dot{b}}_{i-1} である。

これを \bigtriangledown_i の式に適用すると、 a_i,b_i を式から取り除くことができる。

\begin{align} \bigtriangledown_i &= 2(a_{i-1} + 1)\Delta{b} - ({\dot{b}}_{i-1} + 1 + {\dot{b}}_{i-1})\Delta{a} \\ &= 2a_{i-1}\Delta{b} - 2{\dot{b}}_{i-1}\Delta{a} + 2\Delta{b} - \Delta{a}. \end{align}

P_0 の座標の初期状態として、 a_0 = 0,{\dot{b}}_0 = 0 を適用する。

\begin{align} \bigtriangledown_1 &= 2a_0\Delta{b} - 2{\dot{b}}_0\Delta{a} + 2\Delta{b} - \Delta{a} \\ &= 2\Delta{b} - \Delta{a} \\ \end{align} {\rm If}\quad \bigtriangledown_i \geq 0, \\ {\dot{b}}_i = {\dot{b}}_{i-1} + 1, \\ {\rm よって } \\ \begin{align} \bigtriangledown_{i+1} &= 2(a_{i-1} + 1)\Delta{b} - 2({\dot{b}}_{i-1} + 1)\Delta{a} + 2\Delta{b} - \Delta{a} \\ &= (2a_{i-1}\Delta{b} - 2{\dot{b}}_{i-1}\Delta{a} + 2\Delta{b} - \Delta{a}) + 2\Delta{b} - 2\Delta{a} \\ &= \bigtriangledown_i + 2\Delta{b} - 2\Delta{a}. \end{align}

{\rm If}\quad \bigtriangledown_i < 0, \\ {\dot{b}}_i = {\dot{b}}_{i-1}, \\ {\rm よって } \\ \begin{align} \bigtriangledown_{i+1} &= 2(a_{i-1} + 1)\Delta{b} - 2{{\dot{b}}_i}\Delta{a} + 2\Delta{b} - \Delta{a} \\ &= (2a_{i-1}\Delta{b} - 2{{\dot{b}}_{i-1}}\Delta{a} + 2\Delta{b} - \Delta{a}) + 2\Delta{b} \\ &= \bigtriangledown_i + 2\Delta{b}. \end{align}

ということで、(1)の漸化式が成り立つことが検証できた。

続いて、先ほど他のOctantについては(2)と(4)式の右辺を修正しなければならないと書いたが、そのことについて検討する。

今までは最初の点から見て2つ目の点が0第一Octantに属するときのみ考えただけだった。もし最初の点から見た2つ目の点が他のOctantに属する場合、図５で示すように、最初のデータ点からの原点と、各Octantで個別に向いた軸に基づいてab直交座標が再選択される。プロッターの動作 m_1 と m_2 に関連した方向が示される。この情報はm1とm2の動作の割り当てと一緒に表１の列の左側にまとめている。

この表によって、式(1)と(4)の変形は8つのOctantそれぞれについて明示され、前述した(2)と(3)も併せて、汎用的な正しい公式を構成できていることを確認することができる。 （うーむ。この部分意味不明かもしれない。原文を読んで言わんとすることはなんとなくわかるんだけど、うまく翻訳できないんだよねぇ。。）

あと(2)と(4)の適切な式を決定するために任意のデータ点のペアについて適切なoctantを決定するための、コンピュータによる計算可能な手続きが、このアルゴリズムを完成させるために必要である。

(2)の式は |\Delta{x}|-|\Delta{y}| の符号に依存する。 (4)の式を見つけるためには、 \Delta{x},\Delta{y},|\Delta{x}| - |\Delta{y}| に対応したブール値 X,Y,Z を導入する。表１で示すように、これらの変数は負の数であるかそうでないかによって0,1の値をとる。

m_1 の割り当てを決定するためには表１を検査することによって見つけることができる

F(X,Y,Z) = (XZ,Y\bar{Z},\bar{X}Z,\bar{Y}\bar{Z}) を導入する。

同様に、関数

G(X,Y) = (XY,\bar{X}Y,\bar{X}\bar{Y},X\bar{Y})

が表１の G と m_2 列とを併せて、 m_2 に適切に割り当てるために使用される。

まとめると、 \Delta{a},\Delta{b} は |\Delta{x}|-|\Delta{y}| の符号によって \Delta{x},\Delta{y} のいずれかの値をとり、 F,G によって m_1,m_2 の取るべき動作が決まる。これによって式(2),(4)が決定され、(1),(3)の式によって線分を描画することができるようになるというわけである。

ということで、ここまでが原文に書かれてあった内容なのだが、最後のほうはほとんどへたくそな翻訳となってしまった。。 私はこのアルゴリズムを正確に理解しようと頑張っているが、まだちょっと理解度が低いような気がする。。原文もちょっと端折り気味のような気もする。掲載できるページ数の関係だろうか。。

F と G って必要かな？と思っている。 X,Y,Z を3ビットの2進数値とすれば、0-7の値でどのoctantに属するのかがわかる（変換が必要だが）し、Octantがわかれば \Delta{a},\Delta{b},m_1,m_2 のとるべき値は決まるからね。

IBM1401が十進コンピュータであることと何か関係があるのだろうか。それともプロッターの制御APIと何か関係があるのだろうか。。

しかしまあ、これがブレゼンハム・アルゴリズムの原点だとすると、巷のブレゼンハム・アルゴリズムは亜種というか、改良版がはびこっているような気もするな。原文を読んでから、それに近い実装を探してみたけど、ちょっと違うんだよね。やはり先ほどのアセンブラコードを読まなくてはならないのだろうかね。。でもよくよく考えると、これは当然かも。巷のものはビットマップ・ディスプレイ用だからね。

P.S.

この記事の数式をTexで書くのにMathJaxというのを使ってみた。

https://www.mathjax.org/

これ結構すごいんだけど、Texの文法は知らないし、記事中でMathJaxで数式を書くにはTexそのままで書いてもエラーになってしまったりしてそれで四苦八苦してしまった。この記事を書くための時間の大半はTexとMathJaxとの格闘に費やされてしまった。。

ちなみにはてな記法でも数式は書けるけど、これも裏でMathJaxを使っているそうだ。