Bresenhamの線分描画アルゴリズム
目次
ベースの考え方となるアルゴリズム（傾き$ 0 \leq \theta \leq \pi/4のとき） 最適化
線分の傾きに関して一般化
具体的な実装例
付録. $ -1 \leq \Delta y / \Delta x \leq 1の時のアルゴリズムをしっかり導出する https://scrapbox.io/files/685469fcca71c3a0aa813de6.png
目的：$ (x_0, y_0), (x_1, y_2)の2点を結ぶ線分$ lについて次の制約を満たしながら描く。 $ x_0 < x_1, 0 \leq \Delta y / \Delta x\leq 1（$ \Delta y = y_1 - y_0, \Delta x = x_1 - x_0）
線分を描く区間$ x_0 \leq x \leq x_1について、任意の$ xに対して塗られる画素$ (x, y)がただ一つ存在する。
$ xにおいて塗った画素の中心点$ (x, y_p)と、実際に直線が通る点$ (x, y)の誤差$ eを考える。
この誤差が全体において最小になるように、一つずつ点を打っていく。
ベースの考え方となるアルゴリズム（傾き$ 0 \leq \theta \leq \pi/4のとき） $ y_p \leftarrow y_0
for $ x \in \{x_0, x_0 + 1, \cdots, x_1\}
$ xに対して、実際の線分$ l上の点$ (x, y)を求める。
現状$ y_pに対する誤差$ e_p := y - y_pを考える。
誤差$ e_pが$ 0.5を超えたとき、$ y_pに$ 1足す。
（そうすれば、誤差が最小の点$ (x, y_p)を$ xに対して選べる。）
$ (x, y_p)に点を描く。
最適化
最適化1
（$ e_pを一から計算せず、前の計算結果を使ってもっと手早く求める）
$ y_p \leftarrow y_0
$ e_p \leftarrow 0
for $ x \in \{x_0, x_0 + 1, \cdots, x_1\}
$ e_p+= $ \Delta y / \Delta x
if $ e_p > 0.5:
$ y_p += $ 1
$ e_p -= 1
$ (x, y_p)に点を描く。
最適化2.の準備
（分数を扱う処理をなくすべく、$ e_pの計算を全てスケールする）
そのために、誤差をまず0.5減らした上で扱う
その上で、$ 2\Delta x倍にスケールする。
$ f(e) = 2\Delta x( e - 0.5)とする。
関数$ fに通した上で、不等号の条件$ e_p > 0.5を考えても、それは全く同じ意味合いを持つ。
すなわち、$ e_p > 0.5 \Leftrightarrow f(e_p) > f(0.5) = 0
このため、上記における演算はそれぞれ
$ e_p \leftarrow 0を
$ f(e_p) \leftarrow f(0) = -\Delta xに
$ e_p += $ \Delta y/\Delta xを
これは$ e_p \leftarrow e_p + \Delta y / \Delta xと解釈できて
$ f(e_p) \leftarrow f(e_p + \Delta y / \Delta x)
$ = 2\Delta x( e_p + \Delta y /\Delta x - 0.5 )
$ = 2\Delta x( e_p - 0.5 ) + 2 \Delta y
$ = f(e_p) + 2 \Delta y
従って、$ f(e_p) += $ 2\Delta yに
$ e_p -= $ 1を
これは$ e_p \leftarrow e_p - 1と解釈できて
$ f(e_p) \leftarrow f(e_p - 1)
$ = 2\Delta x( e_p - 1- 0.5 )
$ = 2\Delta x( e_p - 0.5 ) - 2 \Delta x
$ = f(e_p) - 2 \Delta x
従って、$ f(e_p) -= $ 2\Delta x
に置き換えられる
最適化2.
$ y_p \leftarrow y_0
$ f(e_p) \leftarrow -\Delta x（ $ f(e_p)を変数として扱っている）
for $ x \in \{x_0, x_0 + 1, \cdots, x_1\}
$ f(e_p)+= $ 2\Delta y
if $ e'_p > f(0.5) = 0:
$ y_p += $ 1
$ f(e_p) -= $ 2 \Delta x
$ (x, y_p)に点を描く。
線分の傾きに関して一般化
ここまでは$ 0 \leq \Delta y/\Delta x \leq 1 \land \Delta x \geq 0のとき、では一般のケースに使えるようにするためには？
まず$ -1 \leq \Delta y / \Delta x \leq 1に拡張してみる。
ざっくりと
誤差$ e_pの性質として、（$ e_p := y - y_p）
$ e_pは、$ y_pの値が変化しない限り、$ \Delta yの符号の方向に進んでいく。
$ |e_p| > 0.5であった時に、上/下のピクセルに移動させることで誤差を$ 1減少させられる。
ということで、$ \Delta yの符号の方向に合うように$ e_pか$ f(e_p)の取り方を反転させれば同じアルゴリズム（誤差の蓄積の計算方法）を使い回すことができる。 つまり、$ \Delta y > 0のときは$ e_p := y - y_p、$ \Delta y < 0のときは$ e_p := y_p - y
（この誤差はいずれの場合も負の値を取りうることに留意する）
$ f(e_p)の初期値として同じく$ -dxを選ぶ。
$ f(e_p)は必ず増える方向に足し合わせる（+=$ |2\Delta y|）。
$ y_p値を更新するときは、傾きが正なら += 1, 負なら -1
さすれば誤差$ f(e_p)はいかなる時でも$ 1減る。
$ \Delta x < 0
1ステップあたり$ |\Delta y|を誤差として足す。累積誤差が0を上回ったら$ |\Delta x|を引く。
どの方向に誤差が増えるかを意識すれば実装できる。
2点をスワップして対応する。
これをすると、$ (x_0, y_0), (x_1, y_1)の順に辿りたい時に困るappbird.icon
スワップするのではなく、$ xを$ \Delta xの符号の向きに$ 1ずつ進めたい
その方向に$ \Delta xが進むと誤差は絶対に増える
ここからは残ったケースである$ \Delta x < 0 \land \Delta y / \Delta x \notin [-1 , 1\rbrackを求めていく。
$ \Delta y / \Delta x < -1 \lor 1 < \Delta y/\Delta xのとき
$ x座標と$ y座標をスワップする（$ y = xに関して描く線分を対称移動してアルゴリズムを走らせる）
ただし、点を打つときは$ (x, y_p)ではなく$ (y_p, x)に点を打つ。
$ B((y_0, x_0), (y_1, x_1))
これにより全てのケースに対応することができた。
なんかもっといい方法はないかな？appbird.icon
コードが一番簡単になるのはこれだと思うappbird.icon
具体的な実装例
code:rs
// color色の点を(x, y)に打つ関数
fn draw_pixel(buffer: &mut Vec<u32>, x: usize, y: usize, color: &Color);
use crate::{snapshot, util::Point2};
/*
fn ceil_div(a:i32, b:i32) -> i32 {
return (a + b - 1) / b;
}
*/
/**
* 与えられたある2点の間を太さ1の線分で描いたときに通る画素の点を列挙するイテレータ。
* 列挙方法はBresenhamの線分描画アルゴリズムに基づく。
* 線分の傾きをkとしたとき、-1 < k < 1の範囲で動作する。
*/
pub struct BresenhamLineIterAlongX {
p1:Point2,
/** このイテレータがどの点を指しているかを示す。 */
current:Point2,
/** p2 - p1 */
offset: Point2,
/** {i32::signum(self.offset.x), i32::signum(self.offset.y)} */
sign_offset:Point2,
/** 現在描こうとしている点currentと本来描こうとしている線分Y方向における累積誤差を表す。
* より正確には、この値に対してスケーリングと並行移動を施したものであり、こうすることで整数範囲内での演算に収まるように工夫されている。
* この値が0を超えるとy値にsが加算される。 */
err:i32,
}
impl BresenhamLineIterAlongX {
pub fn new(p1:&Point2, p2:&Point2) -> Self {
let offset = p2 - p1;
BresenhamLineIterAlongX {
p1: p1.clone(),
current: p1.clone(),
offset,
sign_offset: Point2::new(i32::signum(offset.x), i32::signum(offset.y)),
err: -offset.x,
}
}
fn err_per_x_step(&self) -> i32 {
return 2 * self.offset.y * self.sign_offset.y;
}
fn err_per_y_step(&self) -> i32 {
return -2*self.offset.x *self.sign_offset.x;
}
/*
pub fn y(&self, x:i32) -> i32 {
let n = x - self.current.x;
let err_per_x = self.err_per_x_step();
let err_per_y = self.err_per_y_step();
let required_step = (n*err_per_x - self.offset.x)/(-err_per_y);
return self.current.y + required_step;
}
pub fn x(&self, y:i32) -> i32 {
let m = y - self.current.y;
let err_per_x = self.err_per_x_step();
let err_per_y = self.err_per_y_step();
let required_step = ceil_div(m*-err_per_y - self.offset.x, err_per_x);
return self.current.x + required_step;
}
*/
}
impl Iterator for BresenhamLineIterAlongX {
type Item = Point2;
fn next(&mut self) -> Option<Self::Item> {
if self.current.x == self.p1.x + self.offset.x { return None; }
self.current.x += self.sign_offset.x;
self.err += self.err_per_x_step();
if self.err > 0 {
self.current.y += self.sign_offset.y;
self.err += self.err_per_y_step();
}
return Some(self.current.clone());
}
}
pub struct BresenhamLine{
base_iter:BresenhamLineIterAlongX,
flipped: bool
}
impl BresenhamLine {
pub fn trace(p1:&Point2, p2:&Point2) -> BresenhamLine {
let offset = p2 - p1;
let (base_iter, flipped) =
if offset.x.abs() >= offset.y.abs() {
(BresenhamLineIterAlongX::new(p1, p2), false)
} else {
(BresenhamLineIterAlongX::new(&p1.flipped(), &p2.flipped()), true)
};
BresenhamLine { base_iter, flipped }
}
/*
pub fn x(&self, y:i32) -> i32 {
if self.flipped { self.base_iter.y(y) } else { self.base_iter.x(y) }
}
pub fn y(&self, x:i32) -> i32 {
if self.flipped { self.base_iter.x(x) } else { self.base_iter.y(x) }
}
*/
}
impl Iterator for BresenhamLine {
type Item = Point2;
fn next(&mut self) -> Option<Self::Item> {
let next= self.base_iter.next();
if self.flipped {
match next {
Some(x) => Some(x.flipped()),
None => None
}
} else {
next
}
}
}
参考文献
付録. $ -1 \leq \Delta y / \Delta x \leq 1の時のアルゴリズムをしっかり導出する 最適化1.のアルゴリズムは$ \Delta y < 0について次のように書き換えられる。（$ e_p := y - y_pの定義をとって考える。）
$ y_p \leftarrow y_0
$ e_p \leftarrow 0
for $ x \in \{x_0, x_0 + 1, \cdots, x_1\}
$ e_p+= $ \Delta y / \Delta x
if $ e_p < -0.5:
$ y_p += $ -1
$ e_p -= $ -1
（つまりこう）
$ g(e_p) = -2\Delta x(e_p + 0.5)と定めれば、最適化2.のアルゴリズムも同じように書き換えられる。
$ y_p \leftarrow y_0
$ g(e_p) \leftarrow -\Delta x
for $ x \in \{x_0, x_0 + 1, \cdots, x_1\}
$ g(e_p)+= $ -2\Delta y
if $ g(e_p) > 0: （符号が入れ替わることに注意する）
$ y_p += $ -1
$ g(e_p) -= $ 2 \Delta x
$ (x, y_p)に点を描く。
総括すると、次のように書ける。$ s = \operatorname{sign}(\Delta y) = (\Delta y > 0)\; ?\; 1:-1としよう。
$ y_p \leftarrow y_0
$ e_p \leftarrow -\Delta x
for $ x \in \{x_0, x_0 + 1, \cdots, x_1\}
$ e_p+= $ 2|\Delta y| (= 2s\Delta y)
if $ e_p > 0:
$ y_p += $ s
$ e_p -= $ 2 \Delta x
$ (x, y_p)に点を描く。