B spline Basis Matrix Program

Question

I am implementing an algorithm to generate basis function for bsplines using the following link: http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-coef.html

I should get a tridiagonal system on solving the system of linear equations however I get a lower triangular matrix. I would be glad if someone could tell me what am I am doing wrong? Thank you.

I have written computingCoefficients (b_spline_interpolation.cpp) function and I tried another function basis to compute the coeffiicents. However neither gives a tridiagonal matrix. Please help.

Here is my code

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 only 


            for (  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);

I am sorry for the long post. Eagerly awaiting reply!

PS: I am using Visual Studio 2008 Express Edition.

I am using Eigen Library for matrices: http://eigen.tuxfamily.org/dox/