 |
|
|
|
|
|
| |
| |
|
|
|
|
| |
| |
|
|
"jr" <cre### [at] gmail com> wrote:
> I'll reply tomorrow (Sunday), via email, on topic.
Well, before I go changing the documentation for POV-Ray, I'd like to point out
that currently I have only worked out some form of the Catmull-Rom spline that
seems to work in a predictable way.
I'm currently at the state where I'm looking at the source code in sor.cpp to
try and figure out if the 2 things match in some way.
https://github.com/POV-Ray/povray/blob/master/source/core/shape/sor.cpp
1. The source uses the x and y values of the control points separately in the
terms of the equations used to solve for the coefficients. I do not.
2. The source uses powers of the control point y values in the matrix, whereas I
do not.
3. The source performs a matrix inversion before calculating the coefficients.
No idea why, or what that does (yet).
4. Then there are other adjustments after that,
c[0] = 3.0 * A;
c[1] = 2.0 * B;
c[2] = C;
5. and then they somehow need to "solve a polynomial" - presumably to find a
root.
n = Solve_Polynomial(2, c, r, false, 0.0, stats);
6. then this n value is used as the iterator in the interpolation
while (n--)
{
if ((r[n] >= y[0]) && (r[n] <= y[1]))
{
x[n] = sqrt(r[n] * (r[n] * (r[n] * A + B) + C) + D);
}
}
Which all seems bewilderingly complex to me, but seems to jive with current
documentation.
HOW that all works is currently beyond me.
It's just a giant big black box.
Somewhere around 2000, it seemed to be beyond Warp as well.
I have some thing to do today and tomorrow, but Maybe I can spin my version
around the y-axis and compare it the actual sor object and see if they're even
close.
- BW
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
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
|
|
| |
| |
|
|
|
|
| |
| |
|
|
"Bald Eagle" <cre### [at] netscapenet> wrote:
> 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.
Sent values of A-D for each segment to the debug stream, and there are
definitely minus signs.
> 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.
This apparently has to do with the value of "alpha" in the interpolation, with 0
being Uniform, 1/2 being Centripetal, and 1 being Chordal.
I'm not really sure where "r" comes from in the source, or how that while loop
interpolates the spline, but I just ignored all of that, and things worked out
anyway. :D
After fixing some boneheaded visualization bugs, the math seemed to be working
out with the spline, so I took my SDL version of the SOR spline and used the X
value (minus the minor radius) as the major radius of a torus, and "plotted" a
torus at each interpolation point. (RED)
Seems to overlay the sor {} object (BLUE) very closely. Errors are probably due
to my janky simulation method, and coincident surfaces.
I'll try to work out the details as I write things up, and I'll probably make a
parametric {} and isosurface {} object to further test how well things match up.
- BE
Post a reply to this message
Attachments:
Download 'surface of revolution documentation.png' (255 KB)
Preview of image 'surface of revolution documentation.png'
|
|
| |
| |
|
|
|
|
| |
| |
|
|
"Bald Eagle" <cre### [at] netscapenet> wrote:
> ... I'll probably make a
> parametric {} and isosurface {} object to further test how well things match up.
So this is the isosurface version, overlaid on the POV-Ray sor {} object
object {SOR scale 0.999999}
(It's actually one isosurface for every individual segment of the CR spline)
If I scale to 5 nines, I get a LOT more coincident surface artefacts, so I'd say
that everything is working very nicely.
I would still like to see that original 1989 paper, and perhaps I can work out
the square root thing with some level of rigor.
I'm pretty happy with what I've got, and will get working on some sort of
write-up.
If anyone has any questions, or would like to see a specific type of render,
animation, or mathematical derivation, just let me know, so I can include that.
- BE
Post a reply to this message
Attachments:
Download 'surface of revolution documentation.png' (80 KB)
Preview of image 'surface of revolution documentation.png'
|
|
| |
| |
|
|
|
|
| |
|
|
