ブレゼンハムのアルゴリズムをn次元に拡張する

ブレゼンハムのアルゴリズムは2次元ピクセル、3次元ボクセル上で直線を近似する手段としてよく知られていますが、n次元に拡張することができます。ただ、日本語の記事にn次元拡張に関するものが見つからなかったので、備忘録も兼ねて簡単にまとめることにします。

ブレゼンハムのアルゴリズムとは

ブレゼンハムのアルゴリズムは、二点間を結ぶ線分を整数格子上に近似するためのアルゴリズムです。浮動小数点を使わず、かなり低コストでピクセル上に直線を描画できます。
アルゴリズムについてざっと説明すると、（線分の傾きが0以上1未満の場合）

1. (x,y)に線分の始点の座標を代入し、誤差を初期化する
2. xを1増やし、誤差を更新して、誤差が閾値を超えていたらyを1増加させる
3. 座標のピクセルを描画する
4. xが終点の座標を超えるまでこれを2~3を繰り返す

というものです。誤差の初期化や更新で具体的に何を行うのかについては、詳しくはwikipediaをご参照ください。

また、これでどうして直線が描けるのかということについては、過去に記事を書いたのでそちらをご参照ください。
https://qiita.com/JetPenguin1224/items/29e36296b49da626ba24

n次元に拡張する

二次元でのブレゼンハムのアルゴリズムの擬似コードは以下です。

// 傾きが0以上1未満のとき用
void drawLine(int x0, int y0, int x1, int y1){
	int x = x0;
	int y = y0;
	int e = -(x1-x0);　// 誤差の初期化
	while(x<=x1){
		drawPoint(x,y);//(x,y)のピクセルを塗りつぶす関数
		x++;
		e += 2*(y1-y0);
		if(e >= 0){
			y++;
			e -= 2*(x1-x0);
		}
	}
}

これをn次元に拡張することを考えます。この時に注意しなけばいけないのは以下の3つです。

1. 主軸の決定
2. 各次元をインクリメントするかデクリメントするかを決定
3. 各次元の誤差の管理

まず、1についてですが、上記の擬似コードでは、線分の傾きが0以上1未満で固定されています。つまり、始点と終点のy座標の距離よりも、x座標の距離の方が大きいということです。そして、x座標をループのたびに増やすメインの方向とし、y座標は誤差が閾値を超えたときのみ増やすという処理を行っているわけですが、本来、このメインの方向（主軸）には、最も距離が離れている方向を選択する必要があります。
そこで、n次元に拡張する際には、各次元の中で最も距離が離れている軸を主軸として選択するようにします。

次に、2についてですが、上記の擬似コードでは始点の座標より終点の座標が大きいことを前提としていました。そのため、xもyも常に増加させるだけで良かったわけですが、終点の座標の方が小さかった場合は減少させる必要があります。そこで、n次元の場合は各次元の座標差を計算して、増加させるか減少させるかを決定します。

最後に、3についてですが、二次元の時には管理する誤差はy座標の誤差だけでしたが、n次元に拡張した際には主軸以外の誤差を全て管理する必要があります。具体的には誤差の初期化と更新ですね。
上記の擬似コードのx座標の距離は主軸の距離に、y座標の距離は各次元の距離に読み替える必要があります。例えば、上記の擬似コードでは誤差はx座標の距離の-1倍で初期化していますが、主軸の距離の-1倍に変更します。ループのたびに誤差に加える値も、各次元の距離の二倍に変更します。そして、誤差が閾値を超えていたら、その次元の座標をインクリメント（デクリメント）し、誤差から主軸の距離の二倍を引きます。

さて、これに従って書いたn次元版のブレゼンハムのアルゴリズムは以下になります。

擬似コード

void drawLineND(int* p0, int* p1, int n) {
    int d[n], s[n], err[n], max_d;
    
    // 各次元における距離と符号を計算
    for (int i = 0; i < n; i++) {
        d[i] = abs(p1[i] - p0[i]);
        s[i] = (p0[i] < p1[i]) ? 1 : -1;
    }
    
    // 最大の距離を持つ軸を決定
    max_d = 0;
    for (int i = 1; i < n; i++) {
        if (d[i] > d[max_d]) {
            max_d = i;
        }
    }
    
    // 主軸以外の誤差を初期化
    for (int i = 0; i < n; i++) {
        if (i != max_d) {
            err[i] = - d[max_d];
        }
    }
    
    // 線分の描画
    while (p0[max_d] != p1[max_d]) {
        drawPoint(p0); // 点を描画
        
        p0[max_d] += s[max_d];
        
        for (int i = 0; i < n; i++) {
            if (i != max_d) {
                err[i] += 2 * d[i];
                if (err[i] >= 0) {
                    p0[i] += s[i];
                    err[i] -= 2 * d[max_d];
                }
            }
        }
    }
}

余談

使い道はよくわからないです。三次元まではわかりますが、四次元以上で直線を整数座標に近似する機会ってあるんですかね？　ご存知の方がいたら教えていただきたいです。

参考