Commit 418b41be authored by Nils Rummler's avatar Nils Rummler
Browse files

Overloaded function CubicSplines

parent 6454c5ce
......@@ -82,12 +82,10 @@ namespace ITABase
std::vector<float> vdCoefficients;
};
//! 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. 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); };
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs vdSupportingPoints and vdDataPoints. vdSupportingPoints need to be sorted in ascending order
ITA_BASE_API CPiecewisePolynomial CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints);
//! Uses cubic splines to create a piecewise polynomial (order 3) for the given data pairs vdSupportingPoints and vdDataPoints. Evaluates that Polynom at vdCostumPoints. vdSupportingPoints need to be sorted in ascending order.
ITA_BASE_API std::vector<float> CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const std::vector<float>& vdCostumPoints);
}
#endif // IW_ITA_BASE_SPLINE
......@@ -88,26 +88,27 @@ std::vector<float>::const_iterator ITABase::CPiecewisePolynomial::CoefficientsOf
return vdCoefficients.begin() + idxInterval * NumCoefficients();
}
CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const int iPolynomialOrder)
ITA_BASE_API CPiecewisePolynomial ITABase::CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints)
{
/*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 * 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.
*/
if (vdSupportingPoints.size()!= vdDataPoints.size())
The vector with the linear approximations on the right hand side is called r.
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.
*/
if (vdSupportingPoints.size() != vdDataPoints.size())
ITA_EXCEPT_INVALID_PARAMETER("vdSupportingPoints and vdDataPoints need to have the same size");
//TODO andere als Ordnung 3
if (vdSupportingPoints.size() < 2)
ITA_EXCEPT_INVALID_PARAMETER("vdSupportingPoints and vdDataPoints need to contain at least 3 elements");
//Calculate h_i, formular 1.27b
int iPolynomialOrder = 3;
int iPoints = vdSupportingPoints.size();
std::vector<float> vdH(iPoints-1);
for (int i = 0; i < iPoints-1; i++) {
std::vector<float> vdH(iPoints - 1);
for (int i = 0; i < iPoints - 1; i++) {
vdH[i] = vdSupportingPoints[i + 1] - vdSupportingPoints[i];
if (vdH[i] < 0) {
ITA_EXCEPT_INVALID_PARAMETER("vdSupportingPoints needs to be sorted in ascending order");
......@@ -116,32 +117,32 @@ CPiecewisePolynomial ITABase::Spline(const std::vector<float>& vdSupportingPoint
//SetUp the band matrix with lower, middle and upper diagonals. momentary only natural splines
std::vector<float> vdLower(iPoints - 3), vdMiddle(iPoints - 2), vdUpper(iPoints - 3);
for (int i = 0; i < iPoints - 3;i++) {
for (int i = 0; i < iPoints - 3; i++) {
vdLower[i] = vdH[i + 1];
}
vdUpper = vdLower;//TODO: non natural splines behave different -> put it in the loop
for (int i = 0; i < iPoints - 2;i++) {
for (int i = 0; i < iPoints - 2; i++) {
vdMiddle[i] = 2 * (vdH[i] + vdH[i + 1]);
}
//excitation vector
//TODO only for natural splines with S0=S_iPoints = 0
std::vector<float> vdR(iPoints - 2);
for (int i = 0; i < iPoints - 2; i++)
vdR[i] = 3 *( (vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i]);
for (int i = 0; i < iPoints - 2; i++)
vdR[i] = 3 * ((vdDataPoints[i + 2] - vdDataPoints[i + 1]) / vdH[i + 1] - (vdDataPoints[i + 1] - vdDataPoints[i]) / vdH[i]);
//solve bandmatrix and get the other coefficients
std::vector<float> vdB = Math::BandmatrixSolver(vdLower,vdMiddle,vdUpper,vdR);
std::vector<float> vdB = Math::BandmatrixSolver(vdLower, vdMiddle, vdUpper, vdR);
vdB.insert(vdB.begin(), 0.0);
std::vector<float> vdA(iPoints - 1);
std::vector<float> vdC(iPoints - 1);
std::vector<float> vdD(iPoints - 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]);
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();
......@@ -152,14 +153,19 @@ 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 < iPoints - 1;k++) {
vdCoefficients[4*k] = vdA[k];
vdCoefficients[4*k + 1] = vdB[k];
vdCoefficients[4*k + 2] = vdC[k];
vdCoefficients[4*k + 3] = vdD[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];
vdCoefficients[4 * k + 3] = vdD[k];
}
return CPiecewisePolynomial(vdSupportingPoints, vdCoefficients, iPolynomialOrder);
}
ITA_BASE_API std::vector<float> ITABase::CubicSpline(const std::vector<float>& vdSupportingPoints, const std::vector<float>& vdDataPoints, const std::vector<float>& vdCostumPoints)
{
CPiecewisePolynomial myPoly = CubicSpline(vdSupportingPoints, vdDataPoints);
return myPoly.Value(vdCostumPoints);
}
\ No newline at end of file
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