ド・ブールのアルゴリズム

数値解析の数学分野において、de Boorのアルゴリズム^[1]は、Bスプライン形式のスプライン曲線を評価するための多項式時間で数値的に安定したアルゴリズムです。これは、ベジェ曲線に対するde Casteljauのアルゴリズムの一般化です。このアルゴリズムは、ドイツ系アメリカ人の数学者Carl R. de Boorによって考案されました。de Boorアルゴリズムの単純化された、潜在的に高速な変種が作成されましたが、比較的安定性が低いという問題があります。^{[2] [3]}

B-スプラインの一般的な紹介は、メイン記事で行われています。ここでは、位置 にあるスプライン曲線を評価するための効率的で数値的に安定した方式である de Boor のアルゴリズムについて説明します。曲線は、潜在的にベクトル値の定数 （制御点と呼ばれる） を乗じた B-スプライン関数の合計から構築されます。次数 の B-スプラインは、ノットのグリッド上で定義された次数 の連結された区分多項式関数です（以下では常にゼロベースのインデックスを使用します）。De Boor のアルゴリズムは、スプライン曲線を評価するためにO (p ² ) + O (p) の演算を使用します。注: B-スプラインに関するメイン記事と古典的な出版物^[1]では、異なる表記法が使用されています。B-スプラインはで のようにインデックス付けされます。 ${\displaystyle \mathbf {S} (x)}$${\displaystyle x}$${\displaystyle B_{i,p}(x)}$${\displaystyle \mathbf {c} _{i}}$$${\displaystyle \mathbf {S} (x)=\sum _{i}\mathbf {c} _{i}B_{i,p}(x).}$$${\displaystyle p+1}$${\displaystyle p}$${\displaystyle {t_{0},\dots ,t_{i},\dots ,t_{m}}}$${\displaystyle B_{i,n}(x)}$${\displaystyle n=p+1}$

Bスプラインは局所サポートを持ち、多項式は有限領域でのみ正となり、それ以外の領域ではゼロとなる。Cox-de Boorの再帰式^[4]はこれを示している。 $${\displaystyle B_{i,0}(x):={\begin{cases}1&{\text{if }}\quad t_{i}\leq x<t_{i+1}\\0&{\text{otherwise}}\end{cases}}}$$ $${\displaystyle B_{i,p}(x):={\frac {x-t_{i}}{t_{i+p}-t_{i}}}B_{i,p-1}(x)+{\frac {t_{i+p+1}-x}{t_{i+p+1}-t_{i+1}}}B_{i+1,p-1}(x).}$$

インデックスが位置 を含むノット間隔を定義するとします。再帰式を見ると、このノット間隔では の B スプラインのみが非ゼロであることがわかります。したがって、合計は次のように簡略化されます。 ${\displaystyle k}$${\displaystyle x\in [t_{k},t_{k+1})}$${\displaystyle i=k-p,\dots ,k}$$${\displaystyle \mathbf {S} (x)=\sum _{i=k-p}^{k}\mathbf {c} _{i}B_{i,p}(x).}$$

から、ということがわかります。同様に、再帰では、最も高いクエリされたノット位置はインデックス にあることがわかります。つまり、実際に使用されるノット間隔には、その前後に少なくとも1 個の追加のノットが必要です。コンピュータ プログラムでは、これは通常、最初に使用されたノット位置の回数と最後に使用されたノット位置の回数を繰り返すことによって実現されます。たとえば、および実際のノット位置 の場合、ノット ベクトルを にパディングします。 ${\displaystyle i\geq 0}$${\displaystyle k\geq p}$${\displaystyle k+1+p}$${\displaystyle [t_{k},t_{k+1})}$${\displaystyle p}$${\displaystyle p}$${\displaystyle p=3}$${\displaystyle (0,1,2)}$${\displaystyle (0,0,0,0,1,2,2,2,2)}$

これらの定義により、de Boor のアルゴリズムを説明できるようになりました。このアルゴリズムは、B スプライン関数を直接計算しません。代わりに、同等の再帰式を通じて 評価します。${\displaystyle B_{i,p}(x)}$${\displaystyle \mathbf {S} (x)}$

を の新たな制御点とします。に対して次の再帰が適用されます。 ${\displaystyle \mathbf {d} _{i,r}}$${\displaystyle \mathbf {d} _{i,0}:=\mathbf {c} _{i}}$${\displaystyle i=k-p,\dots ,k}$${\displaystyle r=1,\dots ,p}$$${\displaystyle \mathbf {d} _{i,r}=(1-\alpha _{i,r})\mathbf {d} _{i-1,r-1}+\alpha _{i,r}\mathbf {d} _{i,r-1};\quad i=k-p+r,\dots ,k}$$ $${\displaystyle \alpha _{i,r}={\frac {x-t_{i}}{t_{i+1+p-r}-t_{i}}}.}$$

反復が完了すると、 が得られ、これが望ましい結果であることを意味します。 ${\displaystyle \mathbf {S} (x)=\mathbf {d} _{k,p}}$${\displaystyle \mathbf {d} _{k,p}}$

De Boor のアルゴリズムは、必ずゼロが乗算される項を計算しないため、Cox-de Boor 再帰式を使用した B スプラインの明示的な計算よりも効率的です。${\displaystyle B_{i,p}(x)}$

上記のアルゴリズムは、コンピュータでの実装に最適化されていません。一時的な制御点用のメモリが必要です。各一時的な制御点は、正確に 1 回書き込まれ、2 回読み取られます。 の反復を逆にすると(カウントアップではなくカウントダウン)、 のメモリを再利用することで、一時的な制御点用のメモリのみを使用してアルゴリズムを実行できます。同様に、各ステップで使用される の値は 1 つだけなので、メモリも再利用できます。 ${\displaystyle (p+1)+p+\dots +1=(p+1)(p+2)/2}$${\displaystyle \mathbf {d} _{i,r}}$${\displaystyle i}$${\displaystyle p+1}$${\displaystyle \mathbf {d} _{i,r}}$${\displaystyle \mathbf {d} _{i,r-1}}$${\displaystyle \alpha }$

さらに、一時的な制御点にはゼロベースのインデックスを使用する方が便利です。前のインデックスとの関係は です。したがって、改善されたアルゴリズムが得られます。 ${\displaystyle j=0,\dots ,p}$${\displaystyle i=j+k-p}$

を とします。を について繰り返します。j はカウントダウンする必要がある ことに注意してください。反復が完了すると、結果は になります。 ${\displaystyle \mathbf {d} _{j}:=\mathbf {c} _{j+k-p}}$${\displaystyle j=0,\dots ,p}$${\displaystyle r=1,\dots ,p}$$${\displaystyle \mathbf {d} _{j}:=(1-\alpha _{j})\mathbf {d} _{j-1}+\alpha _{j}\mathbf {d} _{j};\quad j=p,\dots ,r\quad }$$ $${\displaystyle \alpha _{j}:={\frac {x-t_{j+k-p}}{t_{j+1+k-r}-t_{j+k-p}}}.}$$${\displaystyle \mathbf {S} (x)=\mathbf {d} _{p}}$

Python プログラミング言語による次のコードは、最適化されたアルゴリズムの単純な実装です。

 deBoor の定義( k :  int 、 x :  int 、 t 、 c 、 p :  int ):
    """S(x)を評価します。

    引数
    ---------
    k: x を含むノット間隔のインデックス。
    x: 位置。
    t: ノット位置の配列。上記のようにパディングする必要があります。
    c: 制御点の配列。
    p: Bスプラインの次数。
    「」
    d  =  [ c [ j  +  k  -  p ] 、  jは 範囲( 0 、p + 1 )内]      

     rが 範囲( 1 , p + 1 )の 場合:   
         jが 範囲( p , r - 1 , - 1 )の 場合:    
            アルファ =  ( x  -  t [ j  +  k  -  p ])  /  ( t [ j  +  1  +  k  -  r ]  -  t [ j  +  k  -  p ]) 
            d [ j ]  =  ( 1.0  - アルファ)  *  d [ j  -  1 ]  + アルファ *  d [ j ]

     d [ p ]を返す

• ドゥ・カステルジョのアルゴリズム
• ベジェ曲線
• 非一様有理Bスプライン

• ド・ブールのアルゴリズム
• デブール・コックス計算

• PPPACK: Fortranのスプラインアルゴリズムを多数収録
• GNU 科学ライブラリ: C ライブラリ、PPPACKから移植されたスプラインのサブライブラリを含む
• SciPy: Python ライブラリ。FITPACKに基づくスプライン関数を備えたサブライブラリscipy.interpolateが含まれています。
• TinySpline: C++ ラッパーと C#、Java、Lua、PHP、Python、Ruby のバインディングを備えたスプライン用の C ライブラリ
• Einspline: Fortran ラッパーを使用した 1、2、3 次元のスプライン用の C ライブラリ

引用文献

• Carl de Boor (2003)。『スプライン実践ガイド』改訂版。Springer- Verlag。ISBN 0-387-95366-3。