|
 |
At the moment, everything is looking to more or less ok, there probably only
needs to be a minor amount of clarification and explanation to make the
documentation sufficiently user-friendly.
I added some comments to the source code for sor.cpp
and hopefully I can make a little more progress on comparing my own version
written from scratch to what's going on in source.
There are normal, chordal, and centripetal forms.
I suspect mine is the normal form, and the one in source is centripetal.
(I will need to check source against the following, which I stumbled upon)
https://splines.readthedocs.io/en/latest/euclidean/catmull-rom-properties.html
So the first thing we need to do is set up our system of equations that relates
a cubic interpolation between the START (P1) and END (P2) points of the spline
to the constraints for the data, namely that the spline has to pass _through_
the endpoints, and the tangents at the endpoints for the spline segments are
defined using the additional spline points P0 and P3.
Following along with the comments that summarize the matrix values, it will be
fairly easy to match this up with what is going on in this video:
https://www.youtube.com/watch?v=DLsqkWV6Cag
/* Use cubic interpolation. */
k[0] = P[i+1][X] * P[i+1][X]; // P1.x^2, actual spline start
k[1] = P[i+2][X] * P[i+2][X]; // P2.x^2, actual spline end
k[2] = (P[i+2][X] - P[i][X]) / (P[i+2][Y] - P[i][Y]); // df/dt at
spline start
k[3] = (P[i+3][X] - P[i+1][X]) / (P[i+3][Y] - P[i+1][Y]); // df/dt at
spline end (tangent)
k[2] *= 2.0 * P[i+1][X]; // = 2 P1.x * (P2.x - P0.x)
----------------------
(P2.y - P0.y)
k[3] *= 2.0 * P[i+2][X]; // = 2 P2.x * (P3.x - P1.x)
----------------------
(P3.y - P1.y)
w = P[i+1][Y];
Mat[0][0] = w * w * w; // P1.y^3
Mat[0][1] = w * w; // P1.y^2
Mat[0][2] = w; // P1.y
Mat[0][3] = 1.0; // 1
Mat[2][0] = 3.0 * w * w; // 3*P1.y^2
Mat[2][1] = 2.0 * w; // 2*P1.y
Mat[2][2] = 1.0; // 1
Mat[2][3] = 0.0; // 0
w = P[i+2][Y];
Mat[1][0] = w * w * w; // P2.y^3
Mat[1][1] = w * w; // P2.y^2
Mat[1][2] = w; // P2.y
Mat[1][3] = 1.0; // 1
Mat[3][0] = 3.0 * w * w; // 3*P2.y^2
Mat[3][1] = 2.0 * w; // 2*P2.y
Mat[3][2] = 1.0; // 1
Mat[3][3] = 0.0; // 0
-------------------------------------------------------------------------
One can solve a system of linear equations in a matrix using a number of
different methods.
The source uses the method of finding the matrix inverse, and then multiplying
both sides of the matrix equation by the matrix inverse.
(This is where the over-brief equation in the docs is confusing, without any
other context)
-------------------------------------------------------------------------
// Now we find the matrix inverse so that we can solve for the coefficients
/*
ax = b
(1/a)*ax = (1/a)*b
(a^-1)ax = (a^-1)b
(a^-1)*x = 1, so
x = (a^-1) b
*/
MInvers(Mat, Mat); // in math/matrix.cpp
// Once we have the matrix inverse, we take the data constraints and the matrix
inverse and
// multiply them to calculate the coefficients.
/* Calculate coefficients of cubic patch. */
// k[] is the column vector of the constraints
// and Mat[] is the *matrix_inverse* M
// so multiplying together equals the column vector of the coefficients
A = k[0] * Mat[0][0] + k[1] * Mat[0][1] + k[2] * Mat[0][2] + k[3] *
Mat[0][3];
B = k[0] * Mat[1][0] + k[1] * Mat[1][1] + k[2] * Mat[1][2] + k[3] *
Mat[1][3];
C = k[0] * Mat[2][0] + k[1] * Mat[2][1] + k[2] * Mat[2][2] + k[3] *
Mat[2][3];
D = k[0] * Mat[3][0] + k[1] * Mat[3][1] + k[2] * Mat[3][2] + k[3] *
Mat[3][3];
I still feel like there should be some minus signs in there to make the
interpolation work - but maybe those pop out the other side of the MInvers () .
.. . I'll have to write that out in SDL and run it to find out.
I also don't quite understand why we need to take the square root of the result
in line 1082
x[n] = sqrt(r[n] * (r[n] * (r[n] * A + B) + C) + D);
to get the interpolated spline values, but perhaps that may become clear as I
find more bits of time here and the to properly focus on this.
Just wanted to issue a progress report. :)
- BE
Post a reply to this message
|
|
