| 9 | the other segments are in C[1], C[2], ... C[n-1] */ |
| 10 | |
| 11 | CubicVal* CubicSpline2D::SolveCubic(std::vector<float> vals) |
| 12 | { |
| 13 | int n = (int) vals.size() - 1; |
| 14 | const float* x = &vals[0]; |
| 15 | |
| 16 | float* gamma = new float[n+1]; |
| 17 | float* delta = new float[n+1]; |
| 18 | float* D = new float[n+1]; |
| 19 | int i; |
| 20 | /* We solve the equation |
| 21 | [2 1 ] [D[0]] [3(x[1] - x[0]) ] |
| 22 | |1 4 1 | |D[1]| |3(x[2] - x[0]) | |
| 23 | | 1 4 1 | | . | = | . | |
| 24 | | ..... | | . | | . | |
| 25 | | 1 4 1| | . | |3(x[n] - x[n-2])| |
| 26 | [ 1 2] [D[n]] [3(x[n] - x[n-1])] |
| 27 | |
| 28 | by using row operations to convert the matrix to upper triangular |
| 29 | and then back sustitution. The D[i] are the derivatives at the knots. |
| 30 | */ |
| 31 | |
| 32 | gamma[0] = 1.0f/2.0f; |
| 33 | for ( i = 1; i < n; i++) |
| 34 | gamma[i] = 1/(4-gamma[i-1]); |
| 35 | gamma[n] = 1/(2-gamma[n-1]); |
| 36 | |
| 37 | delta[0] = 3*(x[1]-x[0])*gamma[0]; |
| 38 | for ( i = 1; i < n; i++) |
| 39 | delta[i] = (3*(x[i+1]-x[i-1])-delta[i-1])*gamma[i]; |
| 40 | delta[n] = (3*(x[n]-x[n-1])-delta[n-1])*gamma[n]; |
| 41 | |
| 42 | D[n] = delta[n]; |
| 43 | for ( i = n-1; i >= 0; i--) |
| 44 | D[i] = delta[i] - gamma[i]*D[i+1]; |
| 45 | |
| 46 | /* now compute the coefficients of the cubics */ |
| 47 | CubicVal* C = new CubicVal[n]; |
| 48 | for ( i = 0; i < n; i++) |
| 49 | { |
| 50 | C[i].Set((float)x[i], D[i], 3*(x[i+1] - x[i]) - 2*D[i] - D[i+1], |
| 51 | 2*(x[i] - x[i+1]) + D[i] + D[i+1]); |
| 52 | } |
| 53 | return C; |
| 54 | } |
| 55 | |
| 56 | CubicSpline2D::CubicSpline2D() |
| 57 | { |