computation for coefficients is working. only sorted intervalls accepted

src/ITABase/Spline.cpp

@@ -74,8 +74,8 @@ std::vector<float>::const_iterator ITABase::CPiecewisePolynomial::CoefficientsOf
 CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const int iPolynomialOrder)
 {
 	/*setUp the bandmatrix with the the h values, called M. The vektor with the Koefficents c is called c.
-	The vector with the linear approximations on the right hand side is called g.
-	M * c = g
+	The vector with the linear approximations on the right hand side is called r.
+	M * c = r
 	Solve the bandmatrix M with the thomsas algorithm and reconstruct the other coefficients according to:
 	www.tm-mathe.de/Themen/html/funnatsplines.html
 	out of the c vector.
@@ -84,7 +84,7 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
 	int iPoints = vdSupportingPoints.size();
 	std::vector<float> vdH(iPoints-1);
 	
-	for (int i = 0; i < iPoints-1; i++) {
+	for (int i = 0; i < iPoints-1; i++) {//TODO needs tp be positiv
 		vdH[i] = vdSupportingPoints[i + 1] - vdSupportingPoints[i];
 	}
 
@@ -101,24 +101,24 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
 
 	//excitation vector
 	//TODO only for natural splines with S0=S_iPoints = 0
-	std::vector<float> vdG(iPoints - 2);
+	std::vector<float> vdR(iPoints - 2);
 	for (int i = 0; i < iPoints - 2; i++) {
-		vdG[i] = (vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i];
-		vdG[i] *= 3;
+		vdR[i] = (vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i];
+		vdR[i] *= 3;
 	}
-
+/*
 	// solving the band matrix with Thomas Algorithm
 	//https://www.quantstart.com/articles/Tridiagonal-Matrix-Solver-via-Thomas-Algorithm/
 	std::vector<float> vdC_star(iPoints - 3);
 	std::vector<float> vdG_star(iPoints - 2);
 
 	vdC_star[0] = vdUpper[0] / vdMiddle[0];
-	vdG_star[0] = vdG[0] / vdMiddle[0];
+	vdG_star[0] = vdR[0] / vdMiddle[0];
 	for (int i = 1; i < iPoints - 2; i++) {
 		if(i<iPoints-3){
 			vdC_star[i] = vdUpper[i] / (vdMiddle[i] - vdLower[i - 1] * vdC_star[i - 1]);
 		}
-		vdG_star[i] = (vdG[i] - vdG_star[i - 1] * vdLower[i - 1]) / (vdMiddle[i] - vdC_star[i - 1] * vdLower[i - 1]);
+		vdG_star[i] = (vdR[i] - vdG_star[i - 1] * vdLower[i - 1]) / (vdMiddle[i] - vdC_star[i - 1] * vdLower[i - 1]);
 	}
 	
 	std::vector<float> vdB(iPoints - 1); // coeffzent b and solution of bandmatrix
@@ -127,22 +127,42 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
 		vdB[i] = vdG_star[i] - vdC_star[i-1] * vdB[i + 1];
 		i--;
 	}
+	*/
+	
+	//solving bandmatrix with "numerical recipees in C" page 51
+	//if (vdMiddle[0] == 0)
+	//ITA_EXCEPT_INVALID_PARAMETER("Matrix needs to be semi positiv definit.");
+	std::vector<float> vdGamma(iPoints - 2);
+	float fBeta = vdMiddle[0];
+	std::vector<float> vdB(iPoints - 1); // coeffzent b and solution of bandmatrix
 	
+	vdB[1] = vdR[0]/vdMiddle[0];
+	for (int j = 1; j < iPoints - 2; j++) {
+		vdGamma[j] = vdUpper[j - 1] / fBeta;
+		fBeta = vdMiddle[j] - vdLower[j - 1] * vdGamma[j];
+		vdB[j + 1] = (vdR[j] - vdLower[j - 1] * vdB[j]) / fBeta;
+	}
+	for (int j = iPoints - 3; j > 0; j--) {//backsubstitution
+		vdB[j] -= vdGamma[j] * vdB[j + 1];
+	}
+
+
 	//bandmatrix is solved, now get the other coefficients
 
 	std::vector<float> vdA(iPoints - 1);
 	std::vector<float> vdC(iPoints - 1);
 	std::vector<float> vdD(iPoints - 1);
 	
-	
-	for (int i = 0; i < iPoints - 2; i++) {
-		vdC[i] = (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i] - 1 / 3 * vdH[i] * (2 * vdB[i] + vdB[i + 1]);
+	vdB.push_back(0.0);//for ease of loop
+	for (int i = 0; i < iPoints - 1; i++) {
+		vdC[i] = (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i] - (vdH[i]/3) * (2 * vdB[i] + vdB[i + 1]);
 		vdA[i] = (vdB[i + 1] - vdB[i]) / (3 * vdH[i]);
 	}
-	
+	vdB.pop_back();
 	vdD = vdDataPoints;
 	vdD.pop_back();
 
+	//sort for Constructor CPiecewisePolynomial
 	int iCoefficients = (iPoints - 1) * (iPolynomialOrder + 1);
 	std::vector<float> vdCoefficients(iCoefficients);
 	int k = 0;

tests/SplineTest.cpp

@@ -29,10 +29,11 @@ int main(int iNumInArgs, char* pcInArgs[])
 {
 	std::vector<float> vdSupportingPoints{1,7,12,15,19};
 	std::vector<float> vdDataPoints = {8,10,7,8,7};
-	std::vector<float> vdSolution{-0.01045, 0, 0.7091, 8, 0.0304, -1.88185, -0.4194, 10, -0.05005, 0.268014, -0.02026, 7, 0.01520, -0.182432, 0.2365, 8}
+	//std::vector<float> vdSolution{ -0.01045, 0, 0.7091, 8, 0.0304, -1.88185, -0.4194, 10, -0.05005, 0.268014, -0.02026, 7, 0.01520, -0.182432, 0.2365, 8 };
 	int iPolynomialOrder = 3;
 	IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly = IW_ITA_BASE_SPLINE::ITABase::Spline(vdSupportingPoints,vdDataPoints,iPolynomialOrder);
 	
+
 	/*for (int i=0; it < vdDataPoints.size(); i++) {
 		std::cout << vdSolution[i] - myPoly.get_Coefficient(i);
 	}