グラフィックス系のスプライン

参考図書

3Dグラフィックスのための数学入門:クォータニオン・スプライン曲線の基礎

この文章ではノットやインデックス、区間の判定（特に大なり小なり演算子にイコールが付く付かない）などが間違っている可能性があります。参考文献をよく読んで理解し、わたくしに噛んで含めるように教えてください。

前回

標準Ｂスプライン（新版）

indexをちゃんとしました。

List<PVector> PxSpline(List<PVector>spline_cp, int degree, double resolution)
{
 //int n = spline_cp.size();
 int d = degree;
 //int m = n+d+1;
 
 var knots = make_knots_b(spline_cp, degree);
 return PxSpline(spline_cp, knots, degree, resolution);
}

List<PVector> PxSpline(
  List<PVector>spline_cp, double[] knots, int degree, double resolution)
{
 int n = spline_cp.size();
 int d = degree;
 int k = d+1;
 int m = n+d+1;//ノット個数
 
 var ret = new ArrayList<PVector>();  
 for(double t = knots[d]; t<knots[m-d-1]; t+=resolution)
 {
   PVector v = Bx(spline_cp, knots, k, t);    
   ret.add(v);
 }
 
 //ノットの終端値t= knots[m-d-1]は0を吐く
 //PVector v = Bx(spline_cp, knots, k, knots[m-d-1]);     
 //ret.add(v);
 
 //制御点をそのままコピーする場合
 //ただしこの処理は次数が高いと無意味（端点にたどり着かないから）
 //var v = spline_cp.get(spline_cp.size()-1).copy();
 //ret.add(v);

 //このあたりの処理を再帰関数内部で済ませたい場合は後述
 //開スプラインの項のNikO関数を使う
 
 return ret;
}

//knots生成
double[] make_knots_b(List<PVector> spline_cp, int degree) 
{
 int n = spline_cp.size();
 int d = degree;
 int m = n+d+1;
 
 double[] knots = new double[m];
 //knots作成
 for(int i = 0; i<knots.length; i++)
 {
   knots[i]=i;
 }
 return knots; 
}

PVector Bx(List<PVector>spline_cp, double[] knots, int k, double t)
{
   //int degree = k-1;
   PVector sum = new PVector(0,0);
   for(int i = 0; i<spline_cp.size(); i++)
   {
     PVector v = spline_cp.get(i).copy();
     float val = (float)Nik(knots, i, k, t);
     v.mult(val);
     sum.add(v);
   }
   return sum;
}

double Nik(double[] knots, int i, int k, double t)
{
 if(k==1)
 {
   //特別に最終区間の処理をする場合はここでする
   //開スプラインのNikO関数を参照

   if(knots[i]<=t&&t<knots[i+1]){return 1;}
   else{return 0;}
 }
 
 double al = 0;
 if(knots[i+k-1]-knots[i]!=0)
 {
   al = (t-knots[i])/(knots[i+k-1]-knots[i])*Nik(knots, i,k-1,t);
 }
 double ar = 0;
 if(knots[i+k]-knots[i+1]!=0)
 {
   ar = (knots[i+k]-t)/(knots[i+k]-knots[i+1])*Nik(knots, i+1,k-1,t);
 }
 return al+ar;
}

標準Ｂスプライン（旧版）

(両端点を通らない、届かない）

//森北
//標準Ｂスプライン:knots自動生成
//補間が制御点末端まで届かない
//p157,p158,p163
ArrayList<PVector> PxSpline(ArrayList<PVector>spline_cp, int degree, double resolution)
{
 int n = spline_cp.size()-1;
 int d = degree;
 int m = n+d+1;
 
 double[] knots = make_knots_b(spline_cp, degree);
 
 //double[] knots = new double[m+1];
 ////knots作成
 //for(int i = 0; i<knots.length; i++)
 //{
 //  knots[i]=i;
 //}
 
 return PxSpline(spline_cp, knots, degree, resolution);
}

//森北
//標準Ｂスプライン:knots引数入力
//p157,p158
ArrayList<PVector> PxSpline(ArrayList<PVector>spline_cp, double[] knots, int degree, double resolution)
{
 //accuracy:v0,v1の区間の分割精度

 int n = spline_cp.size()-1;
 int d = degree;
 int m = n+d+1;
 
 ArrayList<PVector> ret = new ArrayList<PVector>();  
 for(double t = knots[d]; t<knots[m-d]; t+=resolution)
 {
   PVector v = Bx(spline_cp, knots, d, t);    
   ret.add(v);
 }
 PVector v = Bx(spline_cp, knots, d, knots[m-d]);    
 ret.add(v);
 return ret;
}

//knots生成
double[] make_knots_b(ArrayList<PVector> spline_cp, int degree) 
{
 int n = spline_cp.size()-1;
 int d = degree;
 int m = n+d+1;
 
 double[] knots = new double[m+1];
 //knots作成
 for(int i = 0; i<knots.length; i++)
 {
   knots[i]=i;
 }
 return knots; 
}

//Bx(t)一回分:knots自動生成
PVector Bx(ArrayList<PVector>spline_cp, int degree, double t)
{
 double[] knots = make_knots_b(spline_cp, degree);
 return Bx(spline_cp, knots, degree, t);
}

//Bx(t)一回分:knots引数入力
//森北p157//入力(t_d <=t<= t_{m-d})
PVector Bx(ArrayList<PVector>spline_cp, double[] knots, int degree, double t)
{
 PVector sum = new PVector(0,0);
 for(int i = 0; i<spline_cp.size(); i++)
 {
   PVector v = spline_cp.get(i).copy();
   v.mult((float)Nik(knots, i, degree, t));
   sum.add(v);
 }
 return sum;
}

//森北p157
double Nik(double[] knots, int i, int k, double t)
{
 if(k==0)
 {
   if(knots[i]<=t&&t<knots[i+1]){return 1;}
   else{return 0;}
 }
 
 double al = 0;
 if(knots[i+k]-knots[i]!=0)
 {
   al = (t-knots[i])/(knots[i+k]-knots[i])*Nik(knots, i,k-1,t);
 }
 double ar = 0;
 if(knots[i+k+1]-knots[i+1]!=0)
 {
   ar = (knots[i+k+1]-t)/(knots[i+k+1]-knots[i+1])*Nik(knots, i+1,k-1,t);
 }
 return al+ar;
}

De Boor-Cox（旧版）

両端点を通らない（届かない）De Boor-Coxのアルゴリズム
De Boor-Coxの漸化式とは多分別。制御点＊基底関数の所が制御点からそのまま求められる形になっていると思われる。

//ドボアコックス:Knots自動生成
ArrayList<PVector> PxtSpline(ArrayList<PVector>spline_cp, int degree, double accuracy)
{
 int n = spline_cp.size()-1;
 int d = degree;
 int m = n+d+1;
 
 double[] knots = new double[m+1];
 //knots作成
 for(int i = 0; i<knots.length; i++)
 {
   knots[i]=i;
 }

 return PxtSpline(spline_cp, knots, degree, accuracy);
}

//ドボアコックス:Knots引数入力
ArrayList<PVector> PxtSpline(ArrayList<PVector>spline_cp, double[] knots, int degree, double accuracy)
{
 //accuracy:v0,v1の区間の分割精度
 
 int n = spline_cp.size()-1;
 int d = degree;
 int m = n+d+1;
 
 ArrayList<PVector> ret = new ArrayList<PVector>();
 
 //入力tに応じて、そのtに応じた区間のノットインデックスを取得する必要がある
 int r = d;
 for(double t = knots[d]; t<knots[m-d]; t+=accuracy)
 {
   if(knots[r+1]<t){r++;}    //ただしノットがすっ飛ばされる可能性を考慮に入れていない
   //再帰の場合これで良いと思われる
   int i = r;
   PVector v = Pxt(spline_cp, knots, d, d, i, t);    
   ret.add(v);
 }
 double t = knots[m-d];
 PVector v = Pxt(spline_cp, knots, d, d, r, t);    
 ret.add(v);
 
 return ret;
}

//ドボアコックスPx(t):１回分
PVector Pxt(ArrayList<PVector>spline_cp, double[] knots, int degree, int j, int i, double t)
{
 //p163
 //degreeは設定された曲線の次数
 //j:現在のステップの次数
 //r:入力tを囲むノット区間のインデックス
 //i:制御点のインデックス
 
 int n = spline_cp.size()-1;
 int d = degree;
 int m = n+d+1;

 if(i==-1)
 {
   return new PVector(0,0);
 }
 if(i==m-d+1)
 {
   return new PVector(0,0);
 }
 //次数０
 if(j==0)
 {
   return spline_cp.get(i);    
 }
 
 PVector L = PVector.mult(Pxt(spline_cp,knots,d,j-1,i-1,t),(float)(1-weight(knots,d,j,i,t)));
 PVector R = PVector.mult(Pxt(spline_cp,knots,d,j-1,i,t),(float)weight(knots,d,j,i,t));
 return PVector.add(L,R);
}

double weight(double[] knots, int degree, int j, int i, double t)
{
 return (t-knots[i])/(knots[i+degree-j+1]-knots[i]);
}

開スプライン

両端点を通る開スプライン(端点を通過させるために両端のノットを重ねる）

//開スプライン
//補間が制御点末端まで届く
ArrayList<PVector> CalcuOpenSpline(ArrayList<PVector>spline_cp, int degree, double accuracy)
{  
 double[] knots = make_knots_o(spline_cp,degree); 
 return OSpline(spline_cp, knots, degree, accuracy);
}

//開スプライン用ノット
double[] make_knots_o(ArrayList<PVector> spline_cp, int degree)
{
 return make_knots_o(spline_cp.size(), degree);
}

double[] make_knots_o(int cp_count, int degree)
{
 int n = cp_count-1;//p175:中段
 int d = degree;
 int m = n+d+1;//p175:中段
 
 double[] q = new double[m+1];    
 
 //森北
 //p175
 for(int i = 0; i<=d; i++)
 {
   q[i] = 0;
 }
 for(int i = 0; i<=n-d; i++)
 {
   q[d+1+i] = i+1;
 }
 for(int i = m-d; i<=m; i++)
 {
   //ここでノットが圧縮されるのがキーか
   q[i] = n-d+1;
 }
 return q;
}//make_knots_o

//開スプラインp174
ArrayList<PVector> OSpline(ArrayList<PVector>spline_cp, double[] knots, int degree, double accuracy)
{
 //森北
 //p157,p158
 //accuracy:v0,v1の区間の分割精度
 
 ArrayList<PVector> ret = new ArrayList<PVector>();
 
 //int n = vectors.size()-1;
 int m = knots.length-1;
 int d = degree;

 //p175
 for(double i = knots[d]; i<knots[m-d]; i+=accuracy)
 {
   PVector v = Ox(spline_cp, knots, d, i);    
   ret.add(v);
 }
 PVector v = Ox(spline_cp, knots, d, knots[m-d]);    
 ret.add(v);
 return ret;
}


//開スプラインO(x)(knot省略)
PVector Ox(ArrayList<PVector>spline_cp, int degree, double t)
{
 //tは0から（制御点の個数-次数)まで
 //例えば制御点３、次数２ならばtは0から1
 //例えば制御点４、次数２ならばtは0から2
 double[] knots = make_knots_o(spline_cp,degree); 
 return Ox(spline_cp, knots, degree, t);
}

//森北p157,p175
PVector Ox(ArrayList<PVector>spline_cp, double[] knots, int degree, double t)
{
 PVector sum = new PVector(0,0);
 for(int i = 0; i<spline_cp.size(); i++)
 {
   PVector v = spline_cp.get(i).copy();
   v.mult((float)NikO(knots, degree, i, degree, t));
   sum.add(v);
 }
 return sum;
}

//森北p175
double NikO(double[] knots,int degree, int i, int k, double t)
{
 int m = knots.length-1;
 //int d = k;//再帰用次数ステップ
 int d = degree;
 
 if(k==0)
 {
   //最終区間p175
   if(i==m-d-1)
   {
     if(knots[m-d-1]<=t&&t<=knots[m-d])
     {
       return 1;
     }
     else
     {
       return 0;
     }      
   }
   //通常区間p157
   else if(knots[i]<=t&&t<knots[i+1]){return 1;}
   else{return 0;}    
 }
 
 double al = 0;
 if(knots[i+k]-knots[i]!=0)
 {
   al = (t-knots[i])/(knots[i+k]-knots[i])*NikO(knots, degree, i, k-1, t);
 }
 double ar = 0;
 if(knots[i+k+1]-knots[i+1]!=0)
 {
   ar = (knots[i+k+1]-t)/(knots[i+k+1]-knots[i+1])*NikO(knots, degree, i+1, k-1, t);
 }
 return al+ar;
}

周期的スプライン

//周期的スプライン
ArrayList<PVector> CalcuCyclicSpline(ArrayList<PVector>spline_cp, int degree, double accuracy)
{  
 int d = degree;
 int n = spline_cp.size()-1;
 
 ArrayList<PVector>new_vectors = SetUPCyclicVectors(spline_cp, degree);
 double[] knots = make_knots_cyclic(spline_cp, degree);  
 return PxSpline(new_vectors, knots, degree, accuracy);
}

//一時的に座標を追加。あくまで計算用であり、元のベクトル群に反映されない。
ArrayList<PVector> SetUPCyclicVectors(ArrayList<PVector>spline_cp, int degree)
{
 int n = spline_cp.size()-1;
 int d = degree;
 
 ArrayList<PVector> ret = new ArrayList<PVector>();
 //copy
 for(int i = 0; i<spline_cp.size(); i++)
 {
   ret.add(spline_cp.get(i).copy());
 }
 //追加分
 for(int i = 0; i<d; i++)
 {
   ret.add(ret.get(i).copy());
 }
 
 return ret;
}

double[] make_knots_cyclic(ArrayList<PVector>old_vectos, int degree)  
{
 int n = old_vectos.size()-1;
 int d = degree;
     
 double[] knots = new double[n+2*d+2];
   
 for(int i = 0; i<=n+1; i++)
 {
   knots[i] = i;
 }
 
 for(int i = 1; i<=2*d; i++)
 {
   knots[n+1+i]=knots[n+i]+(knots[i]-knots[i-1]);    
 }
 
 return knots;
}

//周期スプラインC(x)(knot省略)
PVector Cx(ArrayList<PVector>spline_cp, int degree, double t)
{
 double[] knots = make_knots_cyclic(spline_cp,degree); 
 return Ox(spline_cp, knots, degree, t);
}

制御点逆算

補間された密なスプライン座標点からコントロールポイントを逆算する

ArrayList<PVector> InvOpenSpline(ArrayList<PVector>spline, int degree, int cp_count)
{
 double[] knots = make_knots_o(cp_count,degree);
   
 PVector P0 = spline.get(0);
 PVector P3 = spline.get(spline.size()-1);
 
 int n = spline.size();
 double[] x = new double[n];
 double[] y = new double[n];
 double[][] B = new double[n][cp_count];
 
 for(int i = 0; i<spline.size(); i++)
 {
   PVector v = spline.get(i);
   //index位置にある座標のスプライン曲線上におけるt値
   //ただしこの関数は近似に過ぎないため、精度を落とす原因となっている
   double t = GetSpline_t(spline, i, cp_count, degree);
   
   x[i] = v.x;
   y[i] = v.y;    
   for(int j = 0; j<cp_count; j++)
   {
     //制御点が4つならば、
     //ある座標に対して４つの基底関数出力値が関わる
     double b = NikO(knots, degree, j, degree, t);    
     B[i][j]=b;
   }   
 }
 
 double[] Px = svd_solve(B,x);
 double[] Py = svd_solve(B,y);
 
 ArrayList<PVector> ret = new ArrayList<PVector>();
 for(int i = 0; i<Px.length; i++)
 {
   if(cp_count-1<i){break;}
   ret.add(new PVector((float)Px[i],(float)Py[i]));
 }
 return ret;
 
}

double GetSpline_t(ArrayList<PVector>spline, int v_index, int cp_count, int degree)
{
 int t_max = cp_count-degree;
 return t_max*GetSplineLengthRatio(spline, v_index);
}

//0-1の比で返す
double GetSplineLengthRatio(ArrayList<PVector>spline, int index)
{
 double to_v = 0;//目的のindexの座標に至るまでの長さ
 double to_ev = 0;//splineの全長
 
 for(int i = 0; i<spline.size()-1; i++)
 {
   PVector v0 = spline.get(i);
   PVector v1 = spline.get(i+1);
   PVector lv = v1.copy().sub(v0);    
   
   float len = lv.mag();
   if(i<index)
   {
     to_v+=len;
   }      
   to_ev+=len;    
 }  
 return to_v/to_ev;
}

//ソルバー使用、CommonMath
import org.apache.commons.math3.linear.*; 
import org.apache.commons.math3.exception.*;//例外
import org.apache.commons.math3.linear.RealMatrix;

double[] svd_solve(double[][]A, double[]b)
{
 double[]x = new double[b.length]; 
 RealMatrix solver_U = exe_svd(A, b);
 if(solver_U == null)    {      return null;    }
 for(int i = 0; i<x.length; i++)
 {
   if(A[0].length-1<i){break;}
   float u = (float)solver_U.getEntry(i, 0); 
   x[i] = u;
 }
 return x;    
}

boolean is_singlar = false;  
RealMatrix exe_svd(double[][]A, double[]b)
{
 RealMatrix solver_A = new Array2DRowRealMatrix(A);
 RealMatrix solver_b = new Array2DRowRealMatrix(b);      
 try
 {
   //SVD
   SingularValueDecomposition svd = new SingularValueDecomposition(solver_A);    
   RealMatrix solver_u = svd.getSolver().solve(solver_b);
   is_singlar = false;
   
   return solver_u;      
 }
 catch(SingularMatrixException sme)
 {
   is_singlar = true;
 }
 catch(NullPointerException npe)
 {
   is_singlar = true;
 }    
 catch(MathArithmeticException mae)
 {
   is_singlar = true;      
 }
 catch(MaxCountExceededException mcee)
 {
   is_singlar = true;          
 }
 return null;  
}

テストコード

  ArrayList<PVector> control_points = new ArrayList<PVector>();
 PVector cp0 = new PVector(100,100);
 PVector cp1 = new PVector(200,200);
 PVector cp2 = new PVector(300,100);
 PVector cp3 = new PVector(400,200);
 control_points.add(cp0);control_points.add(cp1);control_points.add(cp2);control_points.add(cp3);
 
 DrawLiner(control_points);
 DrawPoints(control_points);
 
 ArrayList<PVector> spline = new ArrayList<PVector>();
 spline = CalcuOpenSpline(control_points,3,0.01);  
 DrawLiner(spline);
 
 ArrayList<PVector> Inv_CP = InvOpenSpline(spline,3,4);  
 stroke(255,0,0);
 DrawPoints(Inv_CP);
 DrawSpline(Inv_CP,3,0.01);
 stroke(0);