図形処理やN次元データの解析などにいわゆる幾何的曲線当てはめが有用なので実装してみましょう。
・近似する対象は点列とする。
点群から点列への構成が必要な場合は独自に取得してください。
・近似曲線はNUBSとする。(without Rational)
もちろんBezier曲線だと曲線セグメントが複数必要になり扱いづらい、
また、NUBSのようにノットベクトルで曲率が変化した方が便利です。
・端点条件は近似曲線の出発点、到着点を、点列の開始点と最終点に重ねる条件と、
端点条件無しの全点列を距離近似の対象にする条件の二つを用意したい。
近似対象の点列 ： $P_i(i=0,\cdots,np-1)$ npは点の数
階数 ： m
コントールポイント数 ： n
ノットベクトル ： $T=\left[t_0,\cdots,t_{m+n-1}\right]$
近似の為のパラメータ ： $\bar{t_i}\quad i=0,\cdots,np-1$
$\bar{t_i}$はコード長の比やコード長の平方根の比を利用します。
$\bar{t_i}=[0,\bar{t}_{i-1}+delta,\cdots,1]$
$delta=\frac{|P_i-P_{i-1}|}{d},\quad d=\sum_{i=1}^{np-1}|P_i-P_{i-1}|$
平方根の比の場合は、 $delta=\frac{\sqrt{|P_i-P_{i-1}|}}{d},\quad d=\sum_{i=1}^{np-1}\sqrt{|P_i-P_{i-1}|}$
ノットベクトルは$\bar{t_i}$を反映して設定します。
$t_0～t_{m-1}=0$
$t_m～t_{n+m-1}=1$
$t_{j+m-1}=(1-a)\bar{t}_{i-1}+a\bar{t}_i,\quad j=1,\cdots,n-m$
それぞれ、
$i=int(jd)$
$a=jd-i$
$d=\frac{np}{n-m+1}$
近似対象の点列$P_i$とパラメータ$\bar{t_i}$における曲線上の点$P(\bar{t_i})$との距離$d_i$の二乗和$f$が最少になるようにコントロールポイント$q_i(i=1\cdots,n-2)$を求めます。
$f=\displaystyle\sum_{i=1}^{np-2}d_i^2=\displaystyle\sum_{i=1}^{np-2}\left|P_i-P(\bar{t_i})\right|^2=\displaystyle\sum_{i=1}^{np-2}\left|P_i-\sum_{j=0}^{n-1}N_{j,m}(\bar{t_i})q_j\right|^2$
$f$を最小化する
$\dfrac {\partial f}{\partial q_k}=\displaystyle\sum_{i=1}^{np-2}\{-2N_{k,m}(\bar{t_i})R_i+2N_{k,m}(\bar{t_i})\sum_{j=1}^{n-2}N_{j,m}(\bar{t_i})q_j\}=0$
それでは具体的にはどんな演算をしているのかを例として、
近似対象の点列中の点毎にどのコントロールポイントセグメントが曲線描画の影響範囲にあるか抽出しておいて、
コントロールポイント毎に微分係数を抽出していきます。画像では$q_1$のみについて
↑ 端点条件無しの場合は、全点列を距離近似の対象にして微分係数をピックアップしていけばいいです。
※参考書物 実践NURBS 著者：三浦曜、望月一正 発行所：工業調査会