スプライン曲線は面白い

　２０１８年１２月２日付のブログで、単純な「平均顔」の画像をデザインし、いくつかの特徴領域のカテゴリを選択し、各特徴領域のデータをランダム化して、多数の変形顔の画像をつくることに言及した。

　変形顔の輪郭は、ランダム・データに基づいているので折れ線となるであろう。この折れ線をなめらかな曲線に変えたいという欲求が生じる。

　そこで、スプラインという技法を学んで、設定されたデータ点をなめらかな曲線で補間する技術を習得することにした。

　スプライン技法によれば、目的とする曲線は、Ｂスプラインと呼ばれる基底関数にその係数を掛けたものの線形結合で表現される。この基底は、少なくとも一つの多項式片を連結したもので表される。２次曲線ならば３個までの２次式片を連結したもの、３次曲線ならば４個までの３次式片を連結したものとなる。

　つまり、スプラインによって表現される曲線は、３層の階層構成をもつアーキテクチャとなっているのである。分かってくると、このアーキテクチャの華麗さとでも言えるものに魅了される。スプライン補間するための計算の手間は半端ではないが、定型計算をライブラリ化したり、数学ツールを使うことによって、その手間を軽減できるであろう。

　基底となるＢスプラインは、Ｎr,i（ｘ）と表される。ｒは階数と呼ばれ、次数＋１となる数である。ｉはその階数に属するＢスプラインの番号を表している。

　例えば、ｒ＝２（１次）のＢスプラインＮ2,i（ｘ）は、ｘ軸上のｘi-2＝＜ｘ＜ｘi，ｘi-2＜ｘi-1＜ｘiの区間に存在し、次の式で表される。
　　　Ｎ2,i（ｘ）＝（ｘ－ｘi-2）／（ｘi-1－ｘi-2）＋（ｘi－ｘ）／（ｘi－ｘi-1）

　Ｎ2,i（ｘ）は、ｘi-2＝＜ｘ＜ｘi-1にある第１項の１次式と区間ｘi-1＝＜ｘ＜ｘiにある第２項の１次式とを連結したものである。ｘ＜ｘi-2，ｘ＝＞ｘiの区間ではＮ2,i（ｘ）＝０になる。ｘi-2，ｘi-1，ｘiを節点という。

　Ｎ2,i-1（ｘ）とは、Ｎ2,i（ｘ）のｉをｉ－１に変えたものであり、Ｎ2,i（ｘ）の存在区間を左方向にシフトしたものであり、次の式で表される。
　　　Ｎ2,i-1（ｘ）＝（ｘ－ｘi-3）／（ｘi-2－ｘi-3）＋（ｘi-1－ｘ）／（ｘi-1－ｘi-2）

　ｘ＜ｘi-3，ｘ＝＞ｘi-1の区間ではＮ2,i-1（ｘ）＝０になる。

　Ｎ3,i（ｘ）は２次のＢスプラインであり、次の漸化式で定義される。
　　　Ｎ3,i（ｘ）＝（ｘ－ｘi-3）Ｎ2,i-1（ｘ）／（ｘi-1－ｘi-3）＋（ｘi－ｘ）Ｎ2,i（ｘ）／（ｘi－ｘi-2）

　Ｎ3,i（ｘ）は、４つの２次式で構成されるが、区間ｘi-2＝＜ｘ＜ｘi-1に存在する２つの２次式は加算によって１つの２次式にまとまる。従って、Ｎ3,i（ｘ）は、３つの２次式を連結したものとなる。ｘ＜ｘi-3，ｘ＝＞ｘiの区間ではＮ3,i（ｘ）＝０になる。

　この漸化式を用いると、関数Ｎ3,i（ｘ）は、ｘ＝ｘi-2およびｘ＝ｘi-1において関数値と一次微分係数が連続な曲線になることが確かめられる。

　Ｎ4,i（ｘ）など高次階のＮr,i（ｘ）も同様の漸化式で定義されるが、詳細は省略する。

　Ｎ3,i（ｘ）のｘi-3をｘi-2の方向に移動させ、ｘi-3＝ｘi-2と重ねると、両節点は多重化される。こうすると、区間ｘi-1＝＜ｘ＜ｘiにある第３の２次式は変わらないが、区間ｘi-3＝ｘi-2＝＜ｘ＜ｘi-1の第２の２次式は更新される。

　更新された第２の２次式をｙ＝ａｘ＾２＋ｂｘ＋ｃとすると、「ｘ＝ｘi-2のときｙ＝０；ｘ＝ｘi-1のときｙは正規化条件から決まる関数値；ｘ＝ｘi-1のときの一次微分係数が第３の２次式の同係数に一致」の３条件からａ，ｂ，ｃが求められる。これによってＮ3,i（ｘ）は２つの２次式を連結したものに更新される。なお、第２の条件によって、ｘ＝ｘi-1のときの関数値は、第３の２次式の関数値に一致することが要請される。

　同様に、Ｎ3,i（ｘ）のｘiをｘi-1の方向に移動させ、節点ｘi-1＝ｘiのように多重化したときも、区間ｘi-3＝＜ｘ＜ｘi-2の第１の２次式は変わらないが、区間ｘi-2＝＜ｘ＜ｘi-1＝ｘiの第２の２次式は更新される。

　同様に、Ｎ3,i（ｘ）のｘi-1，ｘi-2，ｘi-3をｘi-1＝ｘi-2＝ｘi-3のように三重節点とすると、Ｎ3,i（ｘ）の第３の２次式のみが更新されて残る。更新された第３の２次式をｙ＝ａ（ｘ－ｘi）＾２とすると、ｘ＝ｘiのときｙ＝０であるから、ｘ＝ｘi-1のときａ（ｘi-1－ｘi）＾２＝１の正規化条件からａが求められる。Ｎ3,i（ｘ）のｘi，ｘi-1，ｘi-2を三重節点とする場合も同様の類推でｙ＝ａ（ｘ－ｘi-3）＾２が決まる。

　　Ｎ3,i（ｘ）に関し、ｘ軸上の節点をｘ-2＝ｘ-1＝ｘ0＝＜ｘ1＜ｘ2＜．．．＜ｘn＝＜ｘn+1＝ｘn+2＝ｘn+3とするとき、区間［ｘ0，ｘn+1］において、２次のスプラインＳ（ｘ）は２次のＢスプラインの線形結合として表され、
　　　Ｓ（ｘ）＝ｃ1Ｎ3,1（ｘ）＋ｃ2Ｎ3,2（ｘ）＋．．．＋ｃn+3Ｎ3,n+3（ｘ）
となる。ｎは多重節点以外にこの区間の内部に配置する節点の数である。ｃiは、Ｂスプライン係数である。

　データ（Ｘj，Ｙj）（ｊ＝１，２，．．．，Ｎ）が与えられているとき、このデータをスプラインで補間する場合には、補間条件
　　　Ｓ（Ｘj）＝Ｙj
によってＢスプライン係数ｃiを決定する。ここで、データの数（Ｎ）はパラメータの数ｎ＋３に一致させる必要がある。

　Ｂスプラインの局所性により、区間［ｘ0，ｘn+1］内の［ｘ0，ｘ1］など各小区間には３個のＮ3,i（ｘ）のみが実際の計算に関与する。

　例として、ｎ＝３とし、ｘ0＝０，ｘ1＝１，ｘ2＝２，ｘ3＝３，ｘ4＝４とする。Ｘ1＝０，Ｘ2＝０．５，Ｘ3＝１．５，Ｘ4＝２．５，Ｘ5＝３．５，Ｘ6＝４と設定すると、ｃ1＝Ｙ1，ｃ6＝Ｙ6が決定する。

　残りのｃ2～ｃ5を決定するためのＳ（Ｘj）に関する連立方程式は次の通りである。
　　　 ｃ1Ｎ3,1（Ｘ2）＋ｃ3Ｎ3,3（Ｘ2）＋ｃ5Ｎ3,5（Ｘ2）＝Ｙ2
　　　 ｃ2Ｎ3,2（Ｘ3）＋ｃ3Ｎ3,3（Ｘ3）＋ｃ5Ｎ3,5（Ｘ3）＝Ｙ3
　　　 ｃ2Ｎ3,2（Ｘ4）＋ｃ4Ｎ3,4（Ｘ4）＋ｃ5Ｎ3,5（Ｘ4）＝Ｙ4
　　　ｃ2Ｎ3,2（Ｘ5）＋ｃ4Ｎ3,4（Ｘ5）＋ｃ6Ｎ3,6（Ｘ5）＝Ｙ5

　ここで、Ｎ3,1（ｘ）は（０，１）；（１，０）を境界とする三重節点をもつＢスプライン、Ｎ3,2（ｘ）は（１，０）；（４，０）を境界とする３連のＢスプライン、Ｎ3,3（ｘ）は（０，０）；（２，０）を境界とする２連、Ｎ3,4（ｘ）は（２，０）；（４，０）を境界とする２連、Ｎ3,5（ｘ）は（０，０）；（３，０）を境界とする３連のＢスプライン、Ｎ3,6（ｘ）は（３，０）；（４，０）を境界とする三重節点をもつＢスプラインである。

　例えば、小区間［１，２］において、Ｎ3,2（ｘ）＋Ｎ3,3（ｘ）＋Ｎ3,5（ｘ）＝１であることが確かめられる。すなわち、区間［０，４］においてＮ3,1（ｘ）＋．．．＋Ｎ3,6（ｘ）＝１に正規化されていることを知る。

　この方程式を解いてｃ2～ｃ5が求められたら、各小区間について上式のＸjの代わりにｘの刻みを代入してｙの値を計算する。

　当然のことながら、Ｓ（ｘ＝Ｘj）＝Ｙjとならなければならないから、この計算は、ｃ2～ｃ5の計算過程およびＮ3,1（ｘ）～Ｎ3,6（ｘ）の計算式のチェックにもなる。

　下図は、Ｙ1＝２．１，Ｙ2＝２．３，Ｙ3＝１．９，Ｙ4＝１．７，Ｙ5＝２．４，Ｙ6＝２．６と設定し、Ｅｘｃｅｌを用いてｙ＝Ｓ（ｘ）を計算し、グラフ表示したものである。



　（Ｘ1，Ｙ1）～（Ｘ6，Ｙ6）をマーカー表示とし、スプライン補間した部分と区別した方が分かりやすくなるが、一部のデータ点（Ｘj，Ｙj）のみをマーカー表示する方法が見当たらないので、すべてのデータについてマーカー表示なしとした。

　参考文献
　吉本富士市著「インターネット時代の数学―スプライン」（共立出版）