Commit ee88af43 authored by Nils Rummler's avatar Nils Rummler
Browse files

Splines of Order 3 tested and implemented. Added .m file to compare with...

Splines of Order 3 tested and implemented. Added .m file to compare with matlab (not a knob condition).
parent 44c52249
......@@ -52,6 +52,8 @@ namespace ITABase
//! 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>& vdValues) const;
//! Finds the index of Interval at a given x value
int IntervalIndex(const float xValue) const;
//! Return a piecewise polynomial containing the derivations of the polynomials of this one
......
......@@ -40,6 +40,14 @@ float CPiecewisePolynomial::Value(const float xValue) const
ITA_EXCEPT_INVALID_PARAMETER("currently only splines of order 3");
}
std::vector<float> ITABase::CPiecewisePolynomial::Value(const std::vector<float>& vdValues) const
{
std::vector<float> result(vdValues.size());
for (int i = 0; i < vdValues.size(); i++)
result[i] = Value(vdValues[i]);
return result;
}
int ITABase::CPiecewisePolynomial::IntervalIndex(const float xValue) const
{
if (xValue < vdBreakPoints[0] || xValue > vdBreakPoints[NumIntervals()])
......@@ -160,7 +168,7 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
int iCoefficients = (iPoints - 1) * (iPolynomialOrder + 1);
std::vector<float> vdCoefficients(iCoefficients);
int k = 0;
for (int k = 0; k < iPolynomialOrder + 1;k++) {
for (int k = 0; k < iPoints - 1;k++) {
vdCoefficients[4*k] = vdA[k];
vdCoefficients[4*k + 1] = vdB[k];
vdCoefficients[4*k + 2] = vdC[k];
......
......@@ -24,19 +24,49 @@
#include <iostream>
//using namespace std;
bool MyCompare(std::vector<float> a, std::vector<float> b, float epsilon) {
if (a.size() != b.size()) return false;
bool temp = true;
for (int i = 0; i < a.size(); i++) {
if (std::abs(a[i] - b[i]) > epsilon) {
std::cout << "Value of index " << i << " differs over the tolerance of +-"<<epsilon << std::endl;
temp = false;
}
}
return temp;
}
int main(int iNumInArgs, char* pcInArgs[])
{
/*beispiel von TM interaktiv
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 };
int iPolynomialOrder = 3;
IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly = IW_ITA_BASE_SPLINE::ITABase::Spline(vdSupportingPoints,vdDataPoints,iPolynomialOrder);
std::cout << "Quick and dirty test: "<< myPoly.Value(5)<<" -hier soll grob 10 rauskommen;)";
/*for (int i=0; it < vdDataPoints.size(); i++) {
std::cout << vdSolution[i] - myPoly.get_Coefficient(i);
}
IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly1 = IW_ITA_BASE_SPLINE::ITABase::Spline(vdSupportingPoints,vdDataPoints,iPolynomialOrder);
if MyCompare(std::cout << "Quick and dirty test: "<< myPoly1.Value(5)<<" -hier soll grob 10 rauskommen;)";
*/
float epsilon = 0.1;
//Matlab spline. DataPoints sind aus randi
std::vector<float> vdSupportingPoints{ -10, - 9, - 8, - 7, - 6, - 5, - 4, - 3, - 2, - 1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 };
std::vector<float> vdDataPoints = { 13,14,2,14,10,2,5,9,15,15,3,15,11,4,0,-1,-5,-15,-11,-3,2 };
std::vector<float> vdX = { -10, float(1.5), float(-7.898),3, float(3.3), float(9.9),-2 };
std::vector<float> vdResult = { float(13.0000), float(17.4325), float(2.4282), float(8.0000), float(9.4143), float(11.4819), float(15.0000) };
int iPolynomialOrder = 3;
/*IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly1 = IW_ITA_BASE_SPLINE::ITABase::Spline(vdSupportingPoints, vdDataPoints, iPolynomialOrder);
std::cout << "MATLAB SPLINE TEST 1" << std::endl;
if (MyCompare(vdResult, myPoly1.Value(vdX),epsilon))
std::cout << "Test 1 bestanden" << std::endl;
*/
//Test 2
std::vector<float> vdDataPoints2{ 77, 80, 19, 49, 45, 65, 71, 76, 28, 68, 66, 17,12, 50, 96, 35, 59, 23, 76, 26, 51 };
std::vector<float> vdResult2 = { float(77.0000), float(7.9485), float(18.7512), float(50.0000), float(70.7647), float(36.3051), float(28.0000) };
iPolynomialOrder = 3;
IW_ITA_BASE_SPLINE::ITABase::CPiecewisePolynomial myPoly2 = IW_ITA_BASE_SPLINE::ITABase::Spline(vdSupportingPoints, vdDataPoints2, iPolynomialOrder);
std::cout << "MATLAB SPLINE TEST 2" << std::endl;
if (MyCompare(vdResult2, myPoly2.Value(vdX),epsilon))
std::cout << "Test 2 bestanden" << std::endl;
return 0;
}
%% Create Matlab and ITA::BASE splines
X = [-10:1:10];
Y = [77 80 19 49 45 65 71 76 28 68 66 17 12 50 96 35 59 23 76 26 51]; %randi(100,1,21)
Coeffs = [-24.4868088000000 0 27.4868088000000 77;...%Copy & Paste from the resulting ITA::BASE Spline
58.4340477000000 -73.4604263000000 -45.9736176000000 80;...
-54.2493782000000 101.841713000000 -17.5923386000000 19;...
33.5634651000000 -60.9064217000000 23.3429565000000 49;...
-22.0044861000000 39.7839775000000 2.22050858000000 45;...
16.4544754000000 -26.2294807000000 15.7750053000000 65;...
-30.7892181747901 23.0766833120257 12.6794720000000 71;...
54.7927242498993 -69.2909712123447 -33.5017530375546 76;...
-47.3816788248071 95.0872015373532 -7.70552271254611 28;...
4.73399104932904 -47.0578349370681 40.3238438877390 68;...
23.4457146274910 -32.8558617890810 -39.5898528384100 66;...
-7.51684955929288 37.4812820933919 -34.9644325340990 17;...
5.62168360968059 14.9307334155133 17.4475829748062 12;...
-49.9698848794294 31.7957842445550 64.1741006348744 50;...
79.2578559080372 -118.113870393733 -22.1439855143039 96;...
-75.0615387527194 119.659697330378 -20.5981585776589 35;...
77.5513382000000 -106.515846000000 -7.03549194000000 59;...
-85.7249985000000 126.138168000000 12.5868301000000 23;...
73.3486862000000 -131.036835000000 7.68815231000000 76;...
-29.6697369000000 89.0092087000000 -34.3394737000000 26];
myPoly = mkpp(X,Coeffs);
%% Plot
plot(X,Y,'o',xx,spline(X,Y,xx));
hold on
grid on
title 'Splines: Matlab vs. ITA'
plot(xx,ppval(myPoly2,xx));
plot(xx,spline(X,Y,xx)-ppval(myPoly2,xx));
legend('Data','Matlab','ITA','Diff');
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment