Included linear extraploation for values outside the defined piecewise polonomial

include/ITABase/Math/PiecewisePolynomial.h

@@ -54,12 +54,12 @@ namespace ITABase
 			//! Returns the Coefficients of that Polynom. Sorted according to ascending intervals [a0, b0, c0, d0, a1, b1, c1, d1]
 			inline const std::vector<float>& Coefficients() const { return vdCoefficients; };
 
-			//! Evaluate the piecewise polynomial at a given x value.
-			float Value(const float& xValue) const;
-			//! Evaluates the piecewise polynomial at all given x values.
-			std::vector<float> Value(const std::vector<float>& vdXValues) const;
-			//! Finds the index of Interval at a given x value
-			int IntervalIndex(const float xValue) const;
+			//! Evaluate the piecewise polynomial at a given x value. If desired, xValues from outsite the breakpoints can be extrapolated linearly.
+			float Value(const float& xValue, bool bExtrapolate = false) const;
+			//! Evaluates the piecewise polynomial at all given x values. If desired, xValues from outsite the breakpoints can be extrapolated linearly.
+			std::vector<float> Value(const std::vector<float>& vdXValues, bool bExtrapolate = false) const;
+			//! Finds the index of Interval at a given x value. If extrapolation is desired, it will return -1 for xValue outsite the breakpoints
+			int IntervalIndex(const float xValue, bool bExtrapolate = false) const;
 			//! Return a piecewise polynomial containing the derivations of the polynomials of this one
 			CPiecewisePolynomial Derivation() const;
 

src/ITABase/Math/PiecewisePolynomial.cpp

@@ -27,41 +27,52 @@ CPiecewisePolynomial::CPiecewisePolynomial(const std::vector<float>& vdBreakPoin
 
 }
 
-float CPiecewisePolynomial::Value(const float & xValue) const
+float CPiecewisePolynomial::Value(const float & xValue, bool bExtrapolate) const
 {
 	//a(x-x0)^(n) + b(x-x0)^(n-1) + c(x-x0)^(n-2) +...
-	int iInterval = IntervalIndex(xValue);
+	int iInterval = IntervalIndex(xValue, bExtrapolate);
+	
+
+	if (iInterval == -1) {//out of bounds
+		CPiecewisePolynomial ppDerivative = this->Derivation();
+		if (xValue < vdBreakPoints[0]) iInterval = 0;
+		else iInterval = NumIntervals();
+		
+		float fX0 = vdBreakPoints[iInterval];
+		return ppDerivative.Value(fX0) * ( xValue - fX0) + this->Value(fX0); //y=m(x-x0)+d
+	}
+	//inside bounds
 	std::vector<float>::const_iterator itCoef = CoefficientsOfInterval(iInterval);
 	const float x0 = vdBreakPoints[iInterval];
 	float fResult = 0;
-
-	for (int idx = 0; idx <= iOrder;  idx++) {
+	for (int idx = 0; idx <= iOrder;  idx++) { //adding up the parts of the polynomial
 		fResult += itCoef[idx] * pow(xValue - x0, iOrder - idx);
 	}
 	return fResult;
 }
 
-std::vector<float> CPiecewisePolynomial::Value(const std::vector<float>& vdXValues) const
+std::vector<float> CPiecewisePolynomial::Value(const std::vector<float>& vdXValues, bool bExtrapolate) const
 {
 	std::vector<float> result(vdXValues.size());
 	for (int i = 0; i < vdXValues.size(); i++)
-		result[i] = Value(vdXValues[i]);
+		result[i] = Value(vdXValues[i], bExtrapolate);
 	return result;
 }
 
-int CPiecewisePolynomial::IntervalIndex(const float xValue) const
+int CPiecewisePolynomial::IntervalIndex(const float xValue, bool bExtrapolate) const
 {
 	if (xValue < vdBreakPoints[0] || xValue > vdBreakPoints[NumIntervals()])
-		ITA_EXCEPT_INVALID_PARAMETER("xValue is out of bounds");
+		if (bExtrapolate) return -1;
+		else ITA_EXCEPT_INVALID_PARAMETER("xValue is out of bounds");
 	for (int i = 1; i < vdBreakPoints.size(); i++) {
 		if (vdBreakPoints[i] > xValue)
 			return i - 1;
 	}
+	return NumIntervals()-1; //last possibility is xValue == last BreakPoint --> xValue in last interval
 }
 
 CPiecewisePolynomial CPiecewisePolynomial::Derivation() const
 {
-	//TODO: Check if this works correctly
 	const int iNewOrder = std::max(0, iOrder - 1);						// Minimum order is zero
 	const int iNewNumCoefficients = std::max(1, NumCoefficients() - 1);	// At least one coefficient
 

tests/SplineTest.cpp

@@ -97,22 +97,26 @@ bool GeneralTest(){
 	//should get specific Value
 	std::cout << "At random Point x=0.5 the spline should interpolate ~40.9: " << myPoly.Value(0.5) << std::endl;
 	
-	//deviation test
+	
 	std::vector<float> vdBreakPoints = { 0, 1 };
 	std::vector<float> vdCoeff = { 1,1,1 };
 	int iOrder = 2;
 	CPiecewisePolynomial myPoly2 = CPiecewisePolynomial(vdBreakPoints, vdCoeff, iOrder);
-	myPoly2 = myPoly2.Derivation();
-	std::cout << "The deviation of x^2+x+1 should be: " << myPoly2.Coefficients()[0] << " * x + " << myPoly2.Coefficients()[1] << std::endl;
-	std::cout << "so the new polynoial only got " << myPoly2.NumCoefficients() << " coefficients." << std::endl;
+	
+	//deviation test
+	CPiecewisePolynomial myPoly3 = myPoly2.Derivation();
+	std::cout << "The deviation of x^2+x+1 should be: " << myPoly3.Coefficients()[0] << " * x + " << myPoly3.Coefficients()[1] << std::endl;
+	std::cout << "so the new polynoial only got " << myPoly3.NumCoefficients() << " coefficients." << std::endl;
 		
 	//exceptions for Piecewise Polynoms
-	std::cout << "Testing the Value function with Argument out of it's range:" << std::endl;
-	try { myPoly2.Value(2); }
+	std::cout << "Testing the Value function with Argument out of it's range, without extrapolation:" << std::endl;
+	try { myPoly2.Value(-2); }
 	catch (ITAException& err) {
 		std::cout << err << std::endl;
 	}
-	
+	std::cout << "Testing the Value function with Arguments out of it's range, with extrapolation:" << std::endl;
+	std::cout << "At -2, the calculated value should be -1: " << myPoly2.Value(-2,true) << std::endl;
+	std::cout << "At 2, the calculated value should be 6: " << myPoly2.Value(2, true) << std::endl;
 	//exceptions for spline
 	std::cout<<"Testing spline with to few Points:" << std::endl;
 	try {