単調3次スプライン補間法 - タンタン的思考

　この2ヶ月間，投稿論文2本の執筆に明け暮れてました．久々の更新．今回は補間について

Wolberg G., Alfy L., "Monotonic Cubic Spline Interpolation," IEEE Int. Conf. Computer Graphics, 1999.

　今更補間もないだろう……とは思うのですが，これはあくまでツール．データが単調に増加していたら，補間曲線も単調増加にしましょう．という要求．3次関数で補完しますが，結果としてはC2滑らかさを放棄します（補間結果の見た目を重視）．

　本論文では，3次スプライン（Cubic Spline）は
$f_k(x)=a_k(x-x_k)^3+b_k(x-x_k)^2+c_k(x-x_k)+d_k$
と定義されるようです．係数の詳しい定義は省略．しかし，媒介変数表示じゃないのか……まぁ，論文の意図としてはこれで正しいけど．

　ここで，簡便のために（あとで利用しやすくするために），k=0の時を考え，さらにx_0=0，x_1=1，x=tとおいて整理する．その結果，
$f(t)=(1-3t^2+2t^3)y_0+(t-2t^2+t^3)y_0'+(3t^2-2t^3)y_1+(-t^2+t^3)y_1'$
と整理できます．実は，始点y_0とその勾配y_0'，終点y_1とその勾配y_1'を用いて曲線を指定していたことになります．さらにこの式はファーガソン曲線（Ferguson Curve）を表しています．行列で表現すると
$\left( \begin{array}{rrrr}1&t&t^2&t^3\end{array} \right) \left( \begin{array}{rrrr}1&0&0&0\\0&0&1&0\\-3&3&-2&-1\\2&-2&1&1\end{array} \right) \left( \begin{array}{r}y_0\\y_1\\y_0'\\y_1'\end{array} \right)$
となります（超重要）．

　ところで，y_0'とy_1'は
$m_0=y_1-y_0 \\ y_0'=\alpha_0 m_0 \\ y_1'=\beta_0 m_0$
と本論文では定義しています．

　さて，与えられた点データを補間するためには，その点を通過して，かつ1次微分が前後で一致する必要があります．また，今回は冒頭のように2次微分は制約条件から除外します．さらに単調な曲線になるように，αとβの値を調整しなければ行けません．これが今回のミソなんですが，結果だけ．

この図において，x軸，y軸がそれぞれα，βを表す．そして，この塗りつぶし領域の中にαとβの組み合わせがあれば，かならず単調な曲線になります．以上が制約条件になります．

　曲線を客観的に評価する方法も提案されています．最小化するべき目的関数は
$\sum_{k=0}^{n-2}f''(x_k^-)-f''(x_k^+)+K$
となります．意図するところは与えられた点における2次微分が近いほど，評価が良くなるということです．ただし，Kは正定数です．

　目的関数と制約条件について，冒頭で求めた媒介変数表示かつα，βを代入してまとめます．
$\min\sum_{k=0}^{n-2}(-6+2\alpha_k+4\beta_k)m_k+(-3+2\alpha_{k+1}+\beta_{k+1})m_{k+1}+K, \\ m_k=y_{k+1}-y_{k}, \\\beta_km_k=\alpha_{k+1}m_{k+1}, \\ (-6+2\alpha_k+4\beta_k)m_k+(-3+2\alpha_{k+1}+\beta_{k+1})m_{k+1}+K\ge0, \\ -\alpha_km_k\le0, \\ -\beta_km_k\le0\\ (\alpha_k-\beta_k-3)m_k\le0, \\ (2\alpha_k+\beta_k-9)m_k\le0, \\ (\alpha_k+2\beta_k-9)m_k\le0, \\ (-\alpha_k+\beta_k-3)m_k\le0.$
だいたいこんな感じです．m_k=0の時は場合分けが必要だったり，添え字の範囲がまずかったりしますが……併せて，スプライン関数も整理しておくと
$f_k(t)=(-2+\alpha_k+\beta_k)m_kt^3+(3-2\alpha_k-\beta_k)m_kt^2+\alpha_km_kt+y_k$
$f_k(t)=y_0(1-3t^2+2t^3)+\alpha_km_k(t-2t^2+t^3)+(y_k+m_k)(3t^2-2t^3)+\beta_km_k(-t^2+t^3)$
となります．

　とりあえず，ここまで公開．ToDoは具体的な計算結果例を挙げることと，glpkを使ってプログラムに落とし込むこと．