documentation

include/ITABase/Spline.h

@@ -81,7 +81,7 @@ namespace ITABase
 	//! Uses splines of given polynomial order to create a piecewise polynomial for the given data pairs
 	ITA_BASE_API CPiecewisePolynomial Spline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const int iPolynomialOrder);
 
-	//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs
+	//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs. vdSupportingPoints need to be sorted in ascending order
 	ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints) { return Spline(vdSupportingPoints, vdDataPoints, 3); };
 }
 

src/ITABase/Spline.cpp

@@ -27,6 +27,7 @@ CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vdBreakPoin
 float CPiecewisePolynomial::Value(const float xValue) const
 {
 	//TODO: Check if x-value is inside boundaries
+
 	//TODO: Find index of interval, get corresponding coefficients and calculate values
 
 	//Something like this might be necessary
@@ -65,7 +66,6 @@ std::vector<float>::const_iterator ITABase::CPiecewisePolynomial::CoefficientsOf
 	if (idxInterval < 0 || idxInterval >= NumIntervals())
 		ITA_EXCEPT_INVALID_PARAMETER("Index for interval is out of bounds.");
 
-	//TODO: Check if this is correct
 	return vdCoefficients.begin() + idxInterval * NumCoefficients();
 }
 
@@ -73,19 +73,23 @@ 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.
+	/*setUp the bandmatrix M with the the h values. The vector with the coefficent b is called b.
 	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.
+	M * b = r
+	Solve the bandmatrix M with the thomas algorithm and reconstruct the other coefficients according to:
+	https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwjcofWO_IrqAhUux4sKHT52AigQFjABegQIBRAB&url=http%3A%2F
+	out of the b vector.
 	*/
+
 	//Calculate h_i
 	int iPoints = vdSupportingPoints.size();
 	std::vector<float> vdH(iPoints-1);
 	
 	for (int i = 0; i < iPoints-1; i++) {//TODO needs tp be positiv
 		vdH[i] = vdSupportingPoints[i + 1] - vdSupportingPoints[i];
+		if (vdH[i] < 0) {
+			ITA_EXCEPT_INVALID_PARAMETER("vdSupportingPoints needs to be sorted in ascending order");
+		}
 	}
 
 	//SetUp the band matrix with lower, middle and upper diagonals. momentary only natural splines
@@ -106,28 +110,6 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
 		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] = 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] = (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
-	vdB[iPoints - 2] = vdG_star[iPoints - 3];
-	for (int i = iPoints - 3; i > 0;) {
-		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)