Базовая матрица сплайнов

Я реализую алгоритм для генерации базисной функции для bsplines, используя следующую ссылку: http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-coef.html

Я должен получить трехдиагональную систему при решении системы линейных уравнений
однако я получаю нижнюю треугольную матрицу.
Я был бы рад, если бы кто-то мог сказать мне, что я делаю не так? Спасибо.

Я написал computingCoefficients (b_spline_interpolation.cpp), и я попробовал другой базис функции для вычисления коэффициентов. Однако ни один не дает трехдиагональную матрицу. Пожалуйста помоги.

Вот мой код

b_spline_interpolation.cpp

#include "mymain.h"
vector<double> generateParameters(int &nplus1)
{
cout<<"generating parameters.."<<endl;
size_t  i= 0;//unsigned i
double  n = nplus1-1;
vector<double> t(nplus1);
t[0]=0.0;
t[n]=1.0;
for(i=1;i<=n-1;i++)
t[i]=i/n;
cout<<"the parameters are:"<<endl;
for(i=0;i<t.size();i++)//print parameters
cout<<t[i]<<endl;
return   t;
}
void generateKnotsUsingDeBoorFormula(int &nplus1 , int &degree_p , int &mplus1 , vector< double > &knots_u)
{int i,j,jj,p=degree_p,n=nplus1-1,m=mplus1-1;
vector <double> t(nplus1);
t= generateParameters( nplus1 );
for( i=0;i<=p;i++)
knots_u[i]=0;
for( j= 1;j <= n-p;j++)
{
for(i = j; i<=j+p-1 ; i++)
knots_u[ j+p ]+=t[i]/p;
}
for ( jj=0;jj<=p;jj++)
knots_u[m-p+jj]=1;
cout<<"generating knots using DeBoor Formula..."<<endl;
cout<<"the knots are:"<<endl;
for(i=0;i<=m;i++)//print knots
cout<<knots_u[i]<<endl;

}
void generateUniformKnots(int nplus1 , int degree_p , int mplus1 , vector< double > &knots_u)
{   cout<<"generating knots with uniform knot spacing..."<<endl;
int i,j,jj,p=degree_p,n=nplus1-1,m=mplus1-1;
for( i=0;i<=p;i++)
knots_u[i]=0;
for( j= 1;j <= n-p;j++)
knots_u[ j+p ]=((double)j)/(n-p+1);
for ( jj=0;jj<=p;jj++)
knots_u[m-p+jj]=1;

cout<<"the knots are:"<<endl;
for(i=0;i<=m;i++)//print knots
cout<<knots_u[i]<<endl;

}
double calcDist(double x1 , double y1 , double x2 , double y2)
{double d=sqrt( pow( x1 - x2 , 2 )  +pow( y1 - y2 , 2 ));
return d;}
void generateKnotsUsingCentripetalMethod(int nplus1 , int degree_p , int mplus1 , vector< double > &knots_u,MatrixXd temp_matrix)
{   cout<<"generating knots with centripetal method..."<<endl;
int  i , p = degree_p , n=nplus1-1,m=mplus1-1;
vector <double> t(nplus1);
t[ 0 ]= 0;
t[ n ]= 1;
double l = 0 ;
for(i=1;i<=n;i++) //calculating the value of L
{
l=l+calcDist(temp_matrix(i,0),temp_matrix(i,1),temp_matrix(i-1,0),temp_matrix(i-1,1) );
}
for(i=1;i<=n-1;i++)
{
t [i] = ( calcDist(temp_matrix(i,0),temp_matrix(i,1),temp_matrix(i-1,0),temp_matrix(i-1,1) ) )/l;
}for(i=0;i<=m;i++)//print knots
cout<<knots_u[i]<<endl;

}
void computingCoefficients( int nplus1 ,int degree_p , int mplus1 , double u , vector < double >  &knots_u  ,vector <double>  &N )//nplus1 is the number of control points of bspline curve
{

int m = mplus1 -1;
int n = nplus1 -1;

if ( u == knots_u [ 0 ])  // rule out special cases
N[0] = 1.0;
if (u == knots_u [ m ])
N[n] = 1.0 ;

if(u != knots_u[0] && u!=knots_u[m])
{
int  i =0, j =0, ktemp , k =0;

for (  ktemp = 0;ktemp < m;ktemp ++ )//finding k
{
bool b1 = (u >= knots_u [ ktemp ]);
bool b2 = (u < knots_u [ ktemp+1 ]) ;
if (b1&& b2 )
k=ktemp;
}
N[k] = 1.0; // degree 0 coefficient
int degree_d = 0 ;
for (degree_d =1 ;degree_d<=degree_p;degree_d++)  // degree d goes from 1 to p ,outer loop
{
if((k-degree_d>=0)&&(k-degree_d+1>=0))

N[k- degree_d ] = ( ( knots_u[ k + 1 ]-u )/( knots_u[ k + 1 ]-knots_u[ k - degree_d + 1] ) )* N[(k-degree_d)+1]; // right (south-west corner) term onlyfor (  i = k - degree_d + 1 ; i <= k-1 ; i++)  // compute internal terms
if( ( i >= 0 )&&( i + degree_d + 1< N.size())&&( (i+1)<N.size() ) )//negative subcripts of vector give out of range vector
N[ i ] =( ( u - knots_u [ i ] )/( knots_u [ i  + degree_d ] -knots_u [ i ] ) )   * N [ i ] + ( ( knots_u[ i + degree_d + 1  ]- u )/( knots_u[ i + degree_d + 1 ] - knots_u[ i + 1 ] ) ) * N[ i + 1 ];N[ k ] = ( ( u - knots_u [ k  ])/(knots_u[ k + degree_d ]-knots_u[  k  ]) )* N[k]; // let (north-west corner) term only

}
}// end for this if(u != knots_u[0] && u!=knots_u[m])
cout<<"Computing Coefficients.."<<endl;
cout<<"The basis functions are:"<<endl;
for(int i=0; i<nplus1 ; i++)
cout<<N[i]<<endl;

}
void knot( int n,int c, vector<double>  &x)
{cout<<"in Knot function:";
int nplusc,nplus2,i;
nplusc = n + c;
nplus2 = n + 2;
x[1] = 0;
for (i = 2; i <= nplusc; i++)
{
if ( (i > c) && (i < nplus2) )
x[i] = x[i-1] + 1;
else
x[i] = x[i-1];
}
printf("The knot vector is ");
for (i = 1; i <= nplusc; i++){
cout<<x[i]<<"   ";
}
printf("\n");
}

void basis(int c , double t , int npts , vector < double > &x , vector < double >  &n )

{
int nplusc;
int i,k;
double  d,e;
vector <double>  temp(36);

nplusc = npts + c;

printf("knot vector is \n");
for (i = 1; i <= nplusc; i++){
cout<<" "<<x[i]<<endl;
}
printf("t is %f \n", t);/* calculate the first order basis functions n[i][1]    */

for (i = 1; i<= nplusc-1; i++){
if (( t >= x[i]) && (t < x[i+1]))
temp[i] = 1;
else
temp[i] = 0;
}

/* calculate the higher order basis functions */

for (k = 2; k <= c; k++){
for (i = 1; i <= nplusc-k; i++){
if (temp[i] != 0)    /* if the lower order basis function is zero skip the calculation */
d = ((t-x[i])*temp[i])/(x[i+k-1]-x[i]);
else
d = 0;

if (temp[i+1] != 0)     /* if the lower order basis function is zero skip the calculation */
e = ((x[i+k]-t)*temp[i+1])/(x[i+k]-x[i+1]);
else
e = 0;

temp[i] = d + e;
}
}

if (t == (double )x[nplusc]){       /*    pick up last point    */
temp[npts] = 1;
}

/* put in n array   */

for (i = 1; i < npts; i++) {
n[i] = temp[i];
}
}

MatrixXd getDataPoints( int &nplus1,MatrixXd &temp_matrix)//the second argument is required for centripetal spacing
{
int n=nplus1-1;
cout<<"Getting Data Points...."<<endl;
MatrixXd D=MatrixXd :: Zero( nplus1 , 2);//second argument is 2 for x,y annd 3 for x,y,z .currently in 2D
int i=0 ,j=0,k=0;

cout<< "Enter the data points:"<<endl;
for (i=0;i<=n;i++)
{ double temp[2] = {0};
for (j=0;j<1;j++)
{
cout<<"x:";
cin>>temp[j];
cout<<"y:";
cin>>temp[j+1];
cout<<endl;
D(i,j)=temp[j];
D(i,j+1)=temp[j+1];
}
}
temp_matrix =D ;//copy the matrix D into a temporary global matrix.Use this matrix in calculating knots for centripetal method
return D;
}
MatrixXd solveLinearEquation(MatrixXd &N,MatrixXd &D)
{
cout<<"solving System of Linear Equations using LU decomposition..."<<endl;

MatrixXd A =N;
MatrixXd B =D ;
MatrixXd X;
cout << "Here is the matrix A:\n" << A << endl;
cout << "Here is the right hand side b:\n" << B << endl;
X = A.ldlt().solve(B);
cout << "The solution is:\n" << X << endl;
cout << "Here is the matrix A:" << endl << A << endl;
cout << "Here is the matrix B:" << endl << B << endl;
X = A.fullPivLu().solve(B);
if((A*X).isApprox(B))
{
cout << "Here is a solution X to the equation AX = B:" << endl << X << endl;
}
else
cout << "The equation AX = B does not have any solution." << endl;

return X;}

void computeControlPoints( int nplus1 ,int degree_p , int mplus1  , vector< double >  &knots_u , int counter ,vector < double > &N ,MatrixXd &N_matrix , MatrixXd &D )
{

int n= nplus1 - 1;

//Input: n+1 n data points D0, ..... Dn and a degree p
//Output: A B-spline curve of degree p that contains all data points in the given order
int i=0,j=0,k=0;

{
for( j = 0 ; j< N_matrix.cols(); j++ )
{if (counter<N_matrix.rows())

N_matrix(counter, j ) = N[j];
}
}
if(counter==n)
{
cout << "Basis Matrix is :" << endl <<N_matrix << endl;
cout<<"Computing Control Points..."<<endl;
MatrixXd P = solveLinearEquation( N_matrix, D);
}
//Evaluate Nj,p(ti) into row i and column j of matrix N;
/* matrix N is available */
//for (i = 0 ;  i< = n ; i++ )
//Place data point Di on row i of matrix D;
/* matrix D is constructed */
//Use a linear system solver to solve for P from D = N.P
//Row i of P is control point Pi;
//Control points P0, ..., Pn, knot vector U, and degree p determine an

interpolating B-spline curve;

}

mymain.h

#include <iostream>
#include<GL/glut.h>
#include<Eigen/Dense>
#include <vector>
#include<cmath>

using namespace std;  // for iostream library
using namespace Eigen; // for Eigen Libary
//Function declarations
double calcDist(double x1 , double y1 , double x2 , double y2);void knot( int n,int c, vector<double>  &x);
void computeControlPoints( int nplus1 ,int degree_p , int mplus1  , vector< double >  &knots_u , int counter ,vector < double > &N ,MatrixXd &N_matrix , MatrixXd &D );

void computingCoefficients( int nplus1 ,int degree_p , int mplus1 ,double u, vector < double >  &knots_u  ,vector <double>  &N );

void drawBsplineCurve(int nplus1,int mplus1 ,MatrixXd P , vector <double> knots_u );

void generateKnotsUsingCentripetalMethod(int nplus1 , int degree_p , int mplus1 , vector< double > &knots_u,MatrixXd temp_matrix);

MatrixXd getDataPoints( int &nplus1,MatrixXd &temp_matrix);

MatrixXd solveLinearEquation(MatrixXd &N,MatrixXd &D) ;

void generateKnotsUsingDeBoorFormula(int &nplus1 , int &degree_p , int &mplus1 , vector< double > &knots_u);

vector<double> generateParameters(int &nplus1);

void generateUniformKnots(int nplus1 , int degree_p , int mplus1 , vector< double > &knots_u);

void generateKnotsUsingDeBoorFormula(int &nplus1 , int &degree_p , int &mplus1 , vector< double > &knots_u);
void basis(int c,double t,int npts,vector<double> &x,vector<double>  &n);//cubic_spline_interpolation
int spline (int n, int end1, int end2,   double slope1, double slope2,
double x[], double y[],  double b[], double c[], double d[],
int *iflag );

double seval (int n, double u,
double x[], double y[],
double b[], double c[], double d[],
int *last);
double sinteg (int n, double u,
double x[], double y[],
double b[], double c[], double d[],
int *last);
double deriv (int n, double u,
double x[],
double b[], double c[], double d[],
int *last);
//Display
void DisplayText(double positionx,double positiony ,int choice );

void Draw2DOrigin( double xorigin ,double yorigin );void Draw3DOrigin(double xorigin ,double yorigin ,double zorigin);

void SelectObject(int choice);
//linear_interpolation file
void linearInterpolation(double theta1,double tx1,double ty1,double theta2,double tx2,double ty2 ,double timeparameter);
//optimization calcuation
double  optimization(int choice, double a,double b , VectorXd weightsi ,VectorXd weightsj, MatrixXd A ,VectorXd q);

Я прошу прощения за длинный пост. С нетерпением жду ответа!

PS: я использую Visual Studio 2008 Express Edition.

Я использую Eigen Library для матриц: http://eigen.tuxfamily.org/dox/