Debugging error- Vector subscript out of range

[gsaluja]

I am trying to implement B-spline approximation. I am getting an error of vector subscript out of range for Qk.

#include "math.h" 
using namespace std;

class Point2D
{
public:
	double x, y, z;
	double w = 1;
	void addWeight(double weight){
		this->w = weight;
	}
	float distance(Point2D B, float power){
		this->rationalize();
		B.rationalize();
		float dist = 0;
		dist += (this->x - B.x)*(this->x - B.x);
		dist += (this->y - B.y)*(this->y - B.y);
		dist += (this->z - B.z)*(this->z - B.z);
		dist = pow(dist, 0.5);
		dist = pow(dist, power);
		return dist;
	}
	void rationalize(){
		if (this->w != 0){
			this->x = this->x / this->w;
			this->y = this->y / this->w;
			this->z = this->z / this->w;
			this->w = 1;
		}
		else{
			this->x = 0;
			this->y = 0;
			this->w = 1;
		}
	}
    Point2D(){
		this->x = 0;
		this->y = 0;
		this->z = 0;
		this->w = 1;
	}
	Point2D(float x, float y, float z, float w){
		this->x = x;
		this->y = y;
		this->z = z;
		this->w = w;
	}

	Point2D operator*(float temp)
	{
		Point2D point;
		point.x = temp*this->x;
		point.y = temp*this->y;
		point.z = temp*this->z;
		point.w = temp*this->w;
		return point;
	}
	Point2D operator-(Point2D temp)
	{
		Point2D point;
		point.x = this->x-temp.x;
		point.y = this->y-temp.y;
		point.z = this->z-temp.z;
		point.w = this->w-temp.w;
		return point;
	}
	Point2D operator+(float temp)
	{
		Point2D point;
		point.x = temp + this->x;
		point.y = temp + this->y;
		point.z = temp + this->z;
		point.w = temp + this->w;
		return point;
	}
	Point2D operator+(Point2D temp)
	{
		Point2D point;
		point.x = this->x + temp.x;
		point.y = this->y + temp.y;
		point.z = this->z + temp.z;
		point.w = this->w + temp.w;
		return point;
	}
	
};


#endif


//#include <stdio.h>
#include <GL/glut.h>
#include "curve.h"
#include "point2D.h"
#include "color.h"
#include <armadillo>


extern Curve *C;
extern int parametertype;
extern int power;

class Spline
{
private:
    vector<Point2D> B;
	vector<Point2D> D;
	vector<Point2D> Qk;
	vector<Point2D> Q;
    vector<double> m; //knot vector
    int p = 2; //degree of the spline
	vector<double> t; 
public:
    void setDegree(int);
	void Chord_length();
	void Centripetal(float);
	Point2D RationalbSpline(float);
	vector<double> Spline::Uniformly_space_t(int);
	Point2D rationalize(Point2D);
	void addToSpline(Curve *);
	float BasisFun(int, int, float);
	Point2D bSpline(float);
	void getUniformKnots();
	void getClampedKnots();
	void drawRationalBspline(int);
	void drawBspline(int);
	arma::Mat <float> Spline::Nmatrix_approximation();
	void Spline::Qmatrix();
	arma::fmat ConvertoMat(vector<Point2D>);
	vector <Point2D> ConvertoPoint(arma::fmat); 
	void calculateControlPoints_Approximation();

};

#endif /* bspline_h */
#include "bspline.h"
#include "iostream"
#include "curve.h"
#include "color.h"
#include <armadillo>


void Spline::setDegree(int degree){
	this->p = degree;
}

void Spline::addToSpline(Curve *C){
	this->B = C->B;
	this->D = C->B;
	
}
void Spline::getUniformKnots()
{
    	int mLength = this->B.size() + this->p;
	for (int it = 0; it <= mLength ; it++){
		if (this->m.size() == it){
			this->m.push_back(0);
		}
			this->m[it] = it;
	}
}


void Spline::getClampedKnots()
{ 
	{
		int mLength = this->B.size() + this->p;
		for (int i = 0; i <= mLength; i++)
		{
			if (this->m.size() == i){
				this->m.push_back(0);
			}
				if (i <= this->p){
					this->m[i] = 0;
				}
				else if (mLength - i <= this->p){
					this->m[i] = 1;
				}
				else{
					double den = 1 + mLength + 1 - 2 * (this->p + 1);
					this->m[i] = this->m[i - 1] + 1/den ;
				}	
		}
	}
}

float Spline::BasisFun(int i, int pw, float t)
{
	if (pw == 0)
	{
		if (t >= this->m[i] && t <= this->m[i + 1]) 
		{ return 1; }
		else
		{ return 0; }
	}
	else{
		float first = (t - this->m[i]) / (this->m[i + pw] - this->m[i]); 
		float second = (this->m[i + pw + 1] - t) / (this->m[i + pw + 1] - this->m[i + 1]);
		if (this->m[i + pw] == this->m[i]){
			first = 0;
		}
		if (this->m[i + pw + 1] == this->m[i + 1]){
			second = 0;
		}
		float BasisFu = first*BasisFun(i, pw - 1, t) + second*BasisFun(i + 1, pw - 1, t);
		return BasisFu;
	}
}



arma::Mat <float> Spline::Nmatrix_approximation()
{
	arma::fmat N(this->D.size()-2, this->B.size()-2);
	for (int k = 0; k < this->D.size()-2; k++){
		for (int i=0; i < this->D.size()-2; i++){
			float temp = this->BasisFun(i, this->p, this->t[k]);
			N(k, i) = temp;
		}
	}
	return N;
}

void Spline::Qmatrix()
{
	for (int i = 0; i < this->B.size()-2; i++)
	{
		for (int k = 0; k < this->D.size()-2; k++)
		{
			float temp = this->BasisFun(0, this->p, this->t[k]);
			float temp2 = this->BasisFun(this->B.size()-1, this->p, this->t[k]);
			Qk.push_back( (this->D[k] -(this->D[0] * temp) - (this->D[D.size()-1] * temp2)));
			Point2D temp3 =  (this->Qk[k])* this->BasisFun(i, this->p, this->t[k]);
			Point2D temp4 = temp4 + temp3;
			Q.push_back(temp4);
		}
	}
}


void Spline:: calculateControlPoints_Approximation()
{

	using namespace arma;
	Spline::Qmatrix();
	arma::fmat N_a = this->Nmatrix_approximation();
	arma::fmat::iterator it = N_a.begin();
	arma::fmat::iterator it_end = N_a.end();
	arma::fmat Nt = arma::trans(N_a);
	arma::fmat M = Nt*N_a;
	arma::fmat Mv = arma::inv(M);
	arma::fmat P = Mv*this->ConvertoMat(this->Q);
	for (int i = 0; i < this->D.size()-2; i++)
	{	
		this->B[i] = Point2D(P(i, 0), P(i, 1), P(i, 2), P(i, 3));
	}
}


arma::fmat Spline:: ConvertoMat(vector <Point2D> x)
{
	arma::fmat result(x.size(), 4);
	for (int i = 0; i < x.size(); i++){
		result(i, 0) = x[i].x;
		result(i, 1) = x[i].y;
		result(i, 2) = x[i].z;
		result(i, 3) = x[i].w;	
	}
	return result;
}


Point2D Spline::RationalbSpline(float u)
{
	Point2D Sp{};
	Sp.w = 0;
		int n = this->B.size() - 1;
		//loop
		for (int i = 0; i <= n; i++)
		{
			if (u >= this->m[i] && u < this->m[i + this->p + 1]){
				float basicFunctionVal = BasisFun(i, this->p, u);
				Sp.x = Sp.x + basicFunctionVal*this->B[i].x*this->B[i].w;
				Sp.y = Sp.y + basicFunctionVal*this->B[i].y*this->B[i].w;
				Sp.z = Sp.z + basicFunctionVal*this->B[i].z*this->B[i].w;
				Sp.w = Sp.w + basicFunctionVal*this->B[i].w;
			}
		}
		return Sp;
}

Point2D Spline::bSpline(float u)
{
	Point2D Sp{};
	int n = this->B.size() - 1;
	//loop
	for (int i = 0; i <= n; i++)
	{
		if (u >= this->m[i] && u < this->m[i + this->p + 1]){
			float basicFunctionVal = BasisFun(i, this->p, u);
			Sp.x = Sp.x + basicFunctionVal*this->B[i].x;
			Sp.y = Sp.y + basicFunctionVal*this->B[i].y;
			Sp.z = Sp.z + basicFunctionVal*this->B[i].z;
		}
	}
	return Sp;
}

void Spline::drawRationalBspline(int showBspline)
{
	if (1 < this->D.size() && showBspline && this->p < this->D.size()) {
	Spline::getClampedKnots();
	
	Spline::calculateControlPoints_Approximation();

		float u = this->m[this->p];
		while (u <= this->m[this->B.size()]-0.01)
		{
			glLineWidth(1);
			glMaterialfv(GL_FRONT, GL_DIFFUSE, BLUE);
			glColor3f(1.0, 0.0, 0.0);
			glBegin(GL_LINES);
			Point2D bzDecast = this->RationalbSpline(u);
			bzDecast.rationalize();
			glVertex3f(bzDecast.x, bzDecast.y, bzDecast.z);
			u += 0.01;
			bzDecast = this->RationalbSpline(u);
			bzDecast.rationalize();
			glVertex3f(bzDecast.x, bzDecast.y, bzDecast.z);
			glEnd();
		};
		
	}
}

[lastchance]

Point2D temp4 = temp4 + temp3;
Are you sure that is legal?

[gsaluja]

Yes i did an operator overloading for that.

[tpb]

Even if it's "legal", it's unlikely to work. The only sense I can make out of it is:

Point2D temp4 = Point2D() + temp3;   // add temp3 to default Point2D 

(But that's not what the code-as-written does.)

[lastchance]

But how do you work out the RHS (which involves temp4) when temp4 hasn't yet been declared (and given a value)?

You haven't given remotely like enough code to sort through this.

Please sort out your indentation (LOOK AT YOUR POST!) and supply the definition of the spline class (presumably in bspline.h).

Why do you need to keep writing this-> ?

This line doesn't look very likely either: this->D = C->B;

It's a while since I coded up a cubic spline, but why do you need both Qk and Q?

[tpb]

Yes i did an operator overloading for that.

You're not understanding the problem lastchance has discovered (well-spotted, BTW, lastchance!).

Look at this:

int y = 1;
int x = x + y;

What is the value of the x on the right-hand-side? IT HAS NO DEFINITE VALUE! It's basically a random number. So adding y to it just gives you garbage.

[gsaluja]

@tpb.. thank you, i understood now

[gsaluja]

http://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/CURVE-APP-global.html
I am trying to implement this algorithm. I was successful in implementing the interpolation part.. having trouble with approximation

[lastchance]

@gsaluja,

You don't need to write this->, since you are within the spline definition and as far as I can see there is no chance of a name clash. It would make your code more readable if you removed it.

Your code would also be more readable if you used consistent indentation.

The second of these lines is almost certainly wrong:

this->B = C->B;
this->D = C->B;

Your routine void Spline::Qmatrix() is a bit of a mess (including the glaring error we have already pointed out). Do you have a link to the equations you are actually coding up here?

[gsaluja]

http://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/INT-APP/CURVE-APP-global.html

[lastchance]

OK, it look like your routine void Spline::Qmatrix() is intended to do the following bit

Input: n+1 data points D0, D1, ..., Dn, a degree p, and the desired number of control points h+1;

....

Algorithm:

    Obtain a set of parameters t0, ..., tn and a knot vector U;
    Let P0 = D0 and Ph = Dn;
    for k := 1 to n-1 do  // i.e. all except end points
        Compute Qk with the following formula:
                Qk = Dk - N0,p(tk)D0-Nh,p(tk)Dn

    for i := 1 to h-1 do  // i.e. all except end points
        Compute the following and save it to the ith row of matrix Q:
                Sum(from k=1 to n-1) Ni,p(tk)Qk

(Rest of algorithm probably OK ... famous last words!)
This is your current code of this bit (AND IS WRONG, but I need them next to each other to see what is going on). Also removed this-> pointer

void Spline::Qmatrix()
{
	for (int i = 0; i < B.size()-2; i++)
	{
		for (int k = 0; k < D.size()-2; k++)
		{
			float temp = BasisFun(0, p, t[k]);
			float temp2 = BasisFun(B.size()-1, p, t[k]);
			Qk.push_back( (D[k] -(D[0] * temp) - (D[D.size()-1] * temp2)));
			Point2D temp3 =  (Qk[k])*BasisFun(i, p, t[k]);
			Point2D temp4 = temp4 + temp3;  // YUK!!!
			Q.push_back(temp4);
		}
	}
}

Your code seems to have a nested loop. Aren't you supposed to calculate all the Qk first, then put weighted sums into the rows of Q?

I think you would be better declaring the (2-d) sizes of Q and the 1-d size of Qk at the start of this routine. Then you can refer to elements by index, rather than using push_back().

[gsaluja]

The problem i had with that is multiplication of a point to a float value

[lastchance]

The following code looks closer to your algorithm. Note, I've left some indices and sizes as ??? because I don't quite know how your code corresponds to that in the web link.

Also, the web link suggests that Q is a 2-d matrix rather than vector? Not sure - this is going to take me a long time to thrash through the maths.

void Spline::Qmatrix()
{
   Qk = vector<Point2D>( ??? );     // you sort out the size
   Q = vector<Point2D>( ??? );

   for ( int k = ???; k < ???; k++)
   {
      Qk[k] = D[k] - D[0] * BasisFun(0, p, t[k]) - D[D.size()-1] * BasisFun(B.size()-1, p, t[k]);
   }
 
   for ( int i = ???; i < ???; i++ )        // can't guarantee that number of rows      
   {
      Q[i] = 0.0;
      for (int k = ???; k < ???; k++)
      {
         Q[i] += Qk[k] * BasisFun(B.size()-1, p, t[k]);      // VERY unsure about this line
      }
   }
}

The writer of that web link has made a hell of a dog's breakfast of least-squares fitting.

[gsaluja]

I have a code written to convert a vector to a matrix. I could use that to convert vector Q to matrix Q.
I do not understand what you are doing here:
Q = vector( ???, vector<Point2D>( ??? );

[lastchance]

Sorry @gsaluja, I was reading that web link as a "matrix" when it now looks more like a vector of Point2Ds - I've edited my code.

I'll see if I can thrash through that maths tomorrow - I've still got work to finish tonight.

[lastchance]

Hello @gsaluja.

I've got absolutely no way of testing this, as I don't have your other libraries.

I've introduced variables h and n to make it look more like the description in the given algorithm.

Note, if you are using Q to manufacture a matrix, the "working" rows of Q are 1,...,h-1, whereas arrays in C++ start from 0: be careful what you do about this. I am VERY sure that it will affect what you do in routine ConvertoMat(Q);

I haven't checked any of the other routines in your code.

void Spline::Qmatrix()
{
   int h = B.size() - 1;              // Control points are P0,...,Ph
   int n = D.size() - 1;              // Data points are D0,...,Dn

   vector<Point2D> Qk( n );           // Note: Qk[0] won't be used; last element is Qk[n-1], as per algorithm
                                      // Declare Qk here as a local variable
   Q = vector<Point2D>( h );          // Note: Q [0] won't be used; last element is Q [h-1], as per algorithm

   for ( int k = 1; k <= n - 1; k++)
   {
      Qk[k] = D[k] - D[0] * BasisFun( 0, p, t[k] ) - D[n] * BasisFun( h, p, t[k] );
   }
 
   for ( int i = 1; i <= h - 1; i++ )  
   {
      Q[i] = Point2D(0,0,0,1);        // probably unnecessary        
      for (int k = 1; k <= n - 1; k++)
      {
         Q[i] = Q[i] + Qk[k] * BasisFun( i, p, t[k] );
      }
   }
}