ベジエ曲線に対する最近傍点計算の精度保証

(1)

2008年度修士論文

ベジエ曲線に対する最近傍点計算の精度保証

Precision provement for closest point problem with bezier curve.

提出日:2008年2月4日

指導教授

柏木雅英准教授

早稲田大学大学院理工学研究科情報・ネットワーク専攻

3606U012-4

内村創

Uchimura Hajime

(2)

図目次

iii

(5)

表目次

4.1 実行時間(sec) . . . . 15 4.2 区間半径平均 . . . . 15

iv

(6)

第 1 _{章序論}

1.1 研究の背景と目的

ベジエ曲線(Bezier Curve)はN個の制御点から得られるN−1次曲線である．Pierre E. Bezier によって，1972年に自動車デザインのために考案された．ベジエ曲線はわずか数点のみで制御できるというそのコントロール性の良さから，さまざまなグラフィックツールやコンピュータ支援デ

ザイン(CAD)などで用いられてきた．

さらに現代になってディスプレイ装置が高解像度化し，計算機性能も向上した結果，画面上のさまざまなデータをベジエ曲線集合で表現することによってユーザインタフェースの解像度フリー化が行われるようになるなど，ベジエ曲線の応用は更に広がりを見せ始めている．

しかし，ベジエ曲線は本質的には陰関数であり，パラメータ表記にした際に座標が非線形に移動するなどの数値計算的には扱いづらい性質も持っている．そこで，ベジエ曲線に対する基本的な操作を研究し，より安定な解を得る方法を発見することには一定の価値があると考えた．

本研究においては，N次元のベジエ曲線に対し，ある点Pに対してもっとも近づく点ならびにそのパラメータを得る問題の新しい解法アルゴリズムを提案し、区間演算を用いて精度保証を行う。

1

(7)

1.2本論文の構成 2

1.2 本論文の構成

本研究は全5章で構成される．

第2章では，ベジエ曲線の概要ならびに既存研究の概要と，その問題点を述べる．

第3章では具体的な提案手法に言及する．

第4章では性能評価，ならびに既存研究との性能比較を行い，結果を述べる．

第5章では，本稿のまとめについて述べる．

(8)

第 2 _{章概要}

2.1 はじめに

この章では，ベジエ曲線の概要を説明し，既存研究についても触れ，問題点の調査を行う．

2.2 ベジエ曲線とは

ベジエ曲線(Bezier curve)は，1972年にフランスの自動車会社「ルノー」に所属するエンジニ

ア Pierre E. Bezier によって提案された，パラメトリック曲線であり，いくつかの制御点と呼ば

れる点群とパラメータのみで定義される曲線である．

ベジエ曲線は以下のような，いくつかの非常に有用な特徴を持っている．

• パラメータが(0,1)の範囲ならば微分可能である

• 人間が美しいと感じる曲線を得やすい

• 制御点の凸包から曲線が外に出ない

• 形状を維持したまま再分割をすることができる

ベジエ曲線の持つこれらの特徴は特に設計やイラストレーションなどの分野と親和性が高いため，さまざまなアプリケーションでベジエが利用される結果となった．

以下，ベジエ曲線の具体的な概要について述べる．

2.2.1 定義

ここで，ベジエ曲線Rに対して，N個の制御点をP_k(k: 0..N −1)と置く．このとき，ベジエ曲線RはN + 1次元ベジエ曲線と呼ばれる．多くの場合は平面上の点P_k = (x_k, y_k)を制御点として用いる．

3

(9)

2.2ベジエ曲線とは 4

実際のベジエ曲線は以上の点P_kならびにパラメータtに対して，

R(t) = XN k=0

B_kⁿ(t)P_k

として与えられる．この式中のB(t)は，ベルンシュタイン基底を与える関数である．

2.2.2 ベルンシュタイン基底関数

ベルンシュタイン基底関数の定義を行う．ある実数tと，整数n, rが与えられたとき，

B_rⁿ(t) =_nC_rt^r(1−t)^n−r

という関数を考える．これをベルンシュタイン基底ならびにベルンシュタイン関数と呼ぶ．ここで_nC_rはn個からr個を選ぶ組み合わせの数であって、

B_rⁿ(t) = _nC_rt^r(1−t)^n−r

= n!

r!(n−r)!t^r(1−t)^n−r である．

例えば具体的に，n= 3のベルンシュタイン関数は B₀³(t) = (1−t)³

B₁³(t) = 3t(1−t)² B₂³(t) = 3t²(1−t) B₃³(t) = t³

となる．

このベルンシュタイン関数は

F(t) = ((1−t) +t)ⁿ

を二項展開した際に得られる各項を示している．

(10)

2.2.3 ド・カステリョのアルゴリズム

前節において，ベルンシュタイン既定関数を定義したことによって，ベジエ曲線を多項式として表記することができるようになった．具体的に，n= 3のベジエ曲線つまり4次のベジエ曲線を多項式表記してみると，以下のようになる．

ベジエ曲線の定義が R(t) =

XN k=0

B_kⁿ(t)P_k

= XN k=0

N!

k!(N−k)!t^k(1−t)^N−k であるから，

R(t) = B₀³(t)P₀+B₁³(t)P₁+B₂³(t)P₂+B₃³(t)P₃

= (1−t)³P₀+ 3t(1−t)²P₁+ 3t²(1−t)P₂+t³P₃ となる．

この式は

R(t) = (1−t)³P₀+ 3t(1−t)²P₁+ 3t²(1−t)P₂+t³P₃

= (1−t)(1−t)(1−t)P₀+ 3t(1−t)(1−t)P₁+ 3tt(1−t)P₂+tttP₃

= (1−t)³P₀+t(1−t)²P₁+t(1−t)²P₁+t(1−t)²P₁ + t²(1−t)P₂+t²(1−t)P₂+t²(1−t)P₂+t³P₃

= (1−t)²((1−t)P₀+tP₁) +t(1−t)((1−t)P₁+tP₂) + t(1−t)²P₁+t²(1−t)P₂+t²(1−t)P₂+t³P₃

= (1−t)²((1−t)P₀+tP₁) +t(1−t)((1−t)P₁+tP₂) + t²((1−t)P₂+tP₃) +t(1−t)((1−t)P₁+tP₂)

= (1−t)((1−t)((1−t)P₀+tP₁) +t((1−t)((1−t)P₁+tP₂))) + t((1−t)((1−t)P₁+tP₂) +t((1−t)P₂+tP₃))

(11)

として展開し，変形することができる．

ここでこの式を良く見てみると，tA+ (1−t)Bの形になっている部分が多い．そこで，実数t に対して

l_t(A, B) =tA+ (1−t)B

という関数¹を定義すると，ベジエ曲線R(t)は以下のように書き換えることができる．

R(t) = (1−t)((1−t)((1−t)P₀+tP₁) +t((1−t)((1−t)P₁+tP₂))) + t((1−t)((1−t)P₁+tP₂) +t((1−t)P₂+tP₃))

= l_t(l_t(l_t(P₀, P₁), l_t(P₁, P₂)), l_t(l_t(P₁, P₂), l_t(P₂, P₃))) この式から，あるベジエ曲線R(t) 上の点C は

step.1 制御点P_k(k: 0..N)を線形補間し，次の制御点P_k²(k: 0..N−1)を求める step.2 制御点P_k²(k: 0..N−1)を線形補間し，次の制御点P_k³(k: 0..N−2)を求める step. ...

step.N-1 C= (P₀^N−1t+P₁^N−1(1−t))を計算するとして得られることがわかる．

この計算方法をド・カステリョ(de Castljau)のアルゴリズムと呼ぶ．

2.2.4 次数上げ

さて、前節においてド・カステリョのアルゴリズムによって、ベジエ曲線R(t)のパラメータt における点P =R(t)を計算することができた。さらに、このベジエ曲線は点P で分割することができる。

ベジエ曲線Rならびにその制御点P_kを、パラメータα,1−αで分割する場合、ド・カステリョのアルゴリズムの、各ステップにおける制御点C_k^N を

P_kⁿ = αP_kⁿ⁻¹+ (1−α)P_k+1ⁿ⁻¹

として計算することができ、新たなベジエ曲線R₁, R₂ならびにその制御点Q¹_k, Q²_kをそれぞれ Q¹_k = P₁^k

Q²_k = P_n−k^k

1線形補間関数:linear interpolation method

(12)

と与えれば、ベジエ曲線を新たなベジエ曲線二つで表すことができる。なお、このベジエ曲線は分割元のベジエ曲線と形状が完全に一致する。

2.2.5 既存研究

ベジエ曲線P =R(t)と点Cの最小距離を計算する方法をもっともシンプルに考えた場合、ユークリッド距離を

d(t) = |R(t)−C|

と表記して、これを最適化する問題だと考えることができ、n次元ベクトル空間でのN次元ベジエ曲線と点Cからの距離は

d(t) = |R(t)−C|

= vu utXⁿ

i=1

((R_i(t)−C_i)²)

= vu utXⁿ

i=1

((

XN k=0

B_k^N(t)P_k_i−C_i)²)

として表せるので、この多項式多変数に対しての最適化問題を例えばニュートン法などを用いて直接解くことで解を得ることはできる。しかし、この手法はベジエ曲線が制御点の凸包(Convex

hull)からはみ出ないなどの性質を全く利用しないため、非効率であった。

そこで、西田らによって提案されたのがBezier clippingである。Bezier clippingはベジエ曲線が制御点凸包からはみ出ないという性質に着目した手法である。

ベジエ曲線Rの持つ制御点P_k(k: 0..N−1)から d_k = (|P_k−C|, k

N)

として新たなベジエ曲線S(t)ならびにその制御点d_k(k: 0..N)を計算する。このベジエ曲線は第二軸が単調増加なので極値を得られ、解を得ることができる。しかしこのBezier clippingはそもそも高速処理を求めて作られた手法であり、区間演算などを用いて精度保証を行いにくいという問題点を抱えている。

(13)

2.3まとめ 8

2.3 まとめ

本節においては、ベジエ曲線の概要ならびにその数学的性質と、最近傍点問題の既存の解法について述べた。

2.4 区間演算

本節においては、精度保証ツールとして、区間演算について説明する。

2.4.1 定義

区間演算では、ある実数Xを誤差を含めたその上限x_H と下限x_Lで包み込み、

X= [x_L, x_H] ={x_L≤x≤x_H}

として表す。

このとき、四則演算を以下のように定義する。

X+Y = [x_L+y_L, x_H +y_H] X−Y = [x_L−y_L, x_H −y_H]

X·Y = [min(x_L·y_L, x_H ·y_L, x_L·y_H, x_H·y_H), max(x_L·y_L, x_H ·y_L, x_L·y_H, x_H·y_H)]

X/Y = [min(x_L/y_L, x_H/y_L, x_L/y_H, x_H/y_H), max(x_L/y_L, x_H/y_L, x_L/y_H, x_H/y_H)]

本研究においては、すべての数値演算が四則演算のみで十分なので、それ以外の演算子、例えばtanやsin等は不要である。

(14)

第 3 _{章提案手法}

3.1 はじめに

本節では、ベジエ曲線に対して最近傍点問題を解く新たな手法を提案する。

3.1.1 距離区間推定

ベジエ曲線に対する最近傍点問題における最大の問題は重解ならびに局所解の存在である。そこで、まずベジエ曲線と最近傍点との距離を関数とした場合にこの関数がパラメータ区間t∈[0,1]

に対してとり得る最大・最小値を見積もることが重要となる。この最大・最小値の組み合わせをここでは距離区間と呼ぶ。

最大値

ベジエ曲線が制御点の凸包内に必ず収まるという性質から、N次元のベジエ曲線R(t)ならびに制御点P_kと任意の点Cに対して、

|R(t)−C| ≤max(|P₀−C|,|P₁−C|, . . . ,|P_N₋₁−C|)

が必ず成立するので、距離区間の上界は

max(|R(t)−C|)≤max(|P₀−C|,|P₁−C|, . . . ,|P_N−1−C|)

と表せる。

最小値

同様に、N次元のベジエ曲線R(t)ならびに制御点P_kに対して任意の点Cからの距離がとり得る最小値を考える。

9

(15)

3.2最近傍点探索アルゴリズム 10

しかし厳密な最小値を得る事はすなわち最小距離問題を解くことそのものなので、ここでは最小値の下界を見積もることを考える。

ベジエ曲線は制御点の凸包内に必ず収まるので、点Cが制御点凸包内に含まれる場合は、最小値の下界は0となる。

もしベジエ曲線が制御点の凸包内に含まれない場合は、点と制御点凸包の辺との距離が下界となる。実際には制御点から任意の二点の組み合わせを全て取り出し、二点P_i, P_j を通る直線と任意の点Cとの距離を計算すれば良い。

distance=|(P_j−P_i)×(P_j−C)

|P_j−P_i| |

こうして得られた最小値の下界を、距離区間の下端として用いる。

3.2 最近傍点探索アルゴリズム

以上から、N次元のベジエ曲線B(t)ならびにその制御点P_kに対して任意の点Cがとり得る最小距離を得るアルゴリズムは以下のようになる。

step.1

もしベジエ曲線Bが十分に小さいならば、P_N

k=0|P_k−C|を距離として返す。

step.2

ベジエ曲線B(t)をt= 0.5で再分割する。分割でできたベジエ曲線をそれぞれL, Rとする。

step.3

ベジエ曲線L, Rそれぞれと点Cのとり得る距離区間D_l, D_rを計算する。

step.3.1 もし距離区間D_l, D_rが交差部分を持つならば重解をとり得るので、ベジエ曲線L, R 両方に対して再帰的にアルゴリズムを適用し、より距離の近かったベジエ曲線の解を返す。

step.3.2 もし距離区間D_l, D_rが交差部分を持たないならば重解をとり得ないのでより距離区間の最大値が小さかったベジエ曲線を新たなベジエ曲線Bとしてstep.1から繰り返しアルゴリズムを適用し、解を返す。

なお、このアルゴリズムはstep.2でベジエ曲線を分割する際に分割元の制御点凸包よりも小さい面積の凸包に分割されるので、必ず収束する。また、重解を持つ場合には構造上必ず最小比較を行うので、解が最小であることが保証される。

(16)

3.2.1 実装

コードは以下のようになった。

// 実際は T にはtemplateでdoubleかintervalが入る。

T close_point( const Vector<T> &p, Vector<T> &cp, T &u ){

return close_point_sub( p, cp, 0.0, 1.0, u ); // 初期区間[0,1]

}

T close_point_sub( const Vector<T> &p, Vector<T> &cp, T lo, T hi, T &u ){

Bezier<T> *bez, l, r;

Vector<T> v, c, cp_l, cp_r;

T d, dll, dlh, drl, drh, ul, ur;

bez = this;

for( int iteration = 0; iteration < LOOP ;iteration++ ){

bez->subdivide( 1 / 2.0 , l, r, c );

l.approx_distance( p , dll, dlh ); // 左側(0,0.5) r.approx_distance( p , drl, drh ); // 右側(0.5,1)

if( (dlh.getInf() - dll.getSup()) < EPSILON ||

(drh.getInf() - drl.getSup()) < EPSILON ){

// 区間が十分に小さくなったら解. cp = c;

u = (lo+hi) / 2.0;

return (dlh+dll+drh+drl)/4.0;

}else if((dll.getSup() < drl.getInf() && drl.getSup()<dlh.getInf()) ||

(dll.getSup() < drh.getInf() && drh.getSup()<dlh.getInf()) ||

(17)

(drl.getSup() < dll.getInf() && dll.getSup()<drh.getInf()) ||

(drl.getSup() < dlh.getInf() && dlh.getSup()<drh.getInf()) ){

// 区間が重なりを持つならば分解. T mid = (lo+hi)/2.0;

dll = l.close_point_sub( p, cp_l , lo, mid, ul ); // 左側再帰 drl = r.close_point_sub( p, cp_r , mid, hi, ur ); // 右側再帰

// 小さいほうの解を採用

if( dll.getSup() < drl.getSup() ){

cp = cp_l; u = ul; return dll;

}else{

cp = cp_r; u = ur; return drl;

} } else {

// 解区間が明快に分離->小さいほうだけで再ループ if( dlh.getSup() <= drh.getSup() ){

d = dll; bez = &l; hi = (lo+hi)/2.0;

}else{

d = drl; bez = &r; lo = (lo+hi)/2.0;

} } }

cp = c;

u = (lo+hi)/2.0;

return d;

}

(18)

3.3まとめ 13

3.3 まとめ

本節では、距離関数の区間推定を用いたベジエ曲線と点の最近傍点問題の解法アルゴリズムを提案した。

(19)

第 4 _{章性能評価}

4.1 はじめに

本節では、前節までに述べたアルゴリズムの性能評価を行う。前節のアルゴリズムを、C++言語を用いて実装した。評価に用いたマシンはCPU:Core Duo 1.7GHz、メモリ2GBの構成である。

4.2 実験 1

ベジエ曲線の制御点としてそれぞれ

• P₀(10,10)

• P₁(90,90)

• P₂(90,10)

• P₃(10,90) の四点を設定した。

この上に1000回、距離の許容最大誤差と最大ループ数をそれぞれ変化させながら(0,100)の範囲の乱数で点を発生させ、合計実行時間ならびに区間半径の累積を測定した。

4.3 結果

4.3.1 許容誤差、最大ループ数の関係

距離の最大許容誤差epsilonと、最大ループ繰り返し数loopを変化させて実行時間を計測した結果、以上のようになった。なおこのときそれぞれの値は区間半径が1e-12程度で精度保証された。

14

(20)

4.4実験2 15

表4.1: 実行時間(sec)

epsilon 1e-1 1e-2 1e-3 1e-4 1e-5 1e-6 1e-7

loop 2 0.2190 0.2500 0.2500 0.2810 0.2960 0.2960 0.3120 loop 4 0.5310 0.6560 0.7970 0.9060 1.0310 1.1400 1.2190 loop 8 0.6720 0.9529 1.2030 1.4539 1.7030 1.9540 2.1870 loop 16 0.6879 1.030 1.4530 1.9210 2.3750 2.8430 3.2500 loop 32 0.6870 1.046 1.4530 1.9690 2.5470 3.2030 3.9070 loop 64 0.6879 1.046 1.4530 1.9690 2.5470 3.2040 3.9209

4.4 実験 2

ランダムに作成した高次元のベジエ曲線に対して、画像上の格子点それぞれから最近傍点を求め、垂線を描画する。解が安定でなかったり、局所解などに落ちている場合は垂線がうまく描画されないことになる。

4.4.1 結果1

ベジエ曲線の制御点としてそれぞれ

• P₀(10,10)

• P₁(90,90)

• P₂(90,50)

• P₃(90,10)

• P₄(10,90) の四点を設定した。

(21)

4.5まとめ 16

4.4.2 結果1

ベジエ曲線の制御点としてランダムに生成した8点を用意した。

どちらも綺麗に垂線が描画できていることがわかる。

4.5 まとめ

本節においては、前節で述べたアルゴリズムの性能評価として、条件を変化させた場合の実行時間ならびに複雑なベジエに対する垂線描画問題を解いた。

(22)

第 5 _{章結論}

5.1 結論

以上のように、高次で複雑なベジエ曲線に対しても安定に解を精度保証しながら解くことができるアルゴリズムを示した。

現在のディスプレイ装置などは紙などのデバイスに比べると解像度が低く、アンチエリアシングを行わなくては曲線を綺麗にレンダリングすることができない。であるならば、曲線と点の距離を直接解くことができれば、より高品質なアンチエリアシングのためのヒントとしても応用することができる。

また、非有利ベジエ曲線やB-スプライン曲線などのベジエ曲線以外のパラメトリック曲線についても適用範囲を広げていきたい。

17

(23)

参考文献

Bezier clipping関連

[1] T.Nishita, T.Sederberg, M.Kakimoto, ”Ray Tracing Trimmed Rational Surface Patches,”

Computer Graphics, Vol.24, No.4, 1990-8, pp.337-345.T.Nishita, T.Sederberg, M.Kakimoto,

”Ray Tracing Trimmed Rational Surface Patches,” Computer Graphics, Vol.24, No.4, 1990-8, pp.337-345.

[2] T.Sederberg, T.Nishita, ”Curve Intersection using Bezier Clipping,”, CAD, Vol.22, No.9, 1990-11, pp.337-345.

参考にしたサイトなど

[3] http://markun.cs.shinshu-u.ac.jp/learn/cg/cg4/index3.html

[4] http://www.sra.co.jp/people/aoki/Kjo/Manuals/Kjo/Curve/Bezier/Bezier.htm [5] http://www.wakayama-u.ac.jp/ wuhy/GSS/10.html

[6] http://d.hatena.ne.jp/nitoyon/20070919/bezier 2

18

(24)

謝辞

遅々として進まない論文執筆を生暖かく見守ってくださった、柏木研究室の皆様、特に柏木雅英助教授に感謝いたします。

また、何だかんだ言いながらも支えてくれた柏木研究室M2、柏木敬一郎氏にも感謝いたします。

2008年2月内村創

19

ベジエ曲線に対する最近傍点計算の精度保証