Getting equation for a interpolated spline using csaps
Mostrar comentarios más antiguos
Hello good people
I have set of points (x,y,z) which form closed curves. I interpolated the points using MATLAB csaps to get the closed splines. I need the interpolated points to be at a certain equal angular distance suppose (360 degrees/N , N= being the the number of interpolated points on the closed curve) from the centroid of each spline.
- My first question is how do I do that?
- Can I get the spline equation from Matlab?
- For my reference I am attaching an image
Respuesta aceptada
John D'Errico
el 1 de Mzo. de 2023
Editada: John D'Errico
el 1 de Mzo. de 2023
I'll post this as a separate answer, since I did not use a spline at all to build the model, and so it does not get lost in the comments to my other answer.
As I said, the result will be hundreds of coefficients if you built this as a spline model. So nothing at all you can write down. I've attached your data as a .mat file to this comment.
load data.mat
plot3(xyz(:,1),xyz(:,2),xyz(:,3),'-o')
box on
grid on
So we see a simple circular looking curve. But in fact, the points seem to lie exactly in a plane, although they are not in fact in a perfect circle in that plane. The curve is a bit egg shaped. That is, if I transform that 3-d set of points into a 2-d plane, we would see they lie almost perfectly in a palne, but not exactly so. The failure to lie exactly in a plane probably arises because you have written down the points with only 5 significant digits, and that is roughly the amount of error I see as a deviation from the plane.
xyzbar = mean(xyz,1);
svd(xyz - xyzbar)
So the smallest singular value there is roughly 5 powers of 10 smaller than the largest singular value. Again, that is good evidence the deviation from a perfect plane lies in having rounded off the 5th decimal place in your numbers.
Anyway, we can now try to build a model for your data that will be roughly as accurate as your data. I could transform the problem into a 2-dimensional one. Or, simpler yet to follow, we can just work in cylindrical coordinates. Since you wanted points that were at a set of uniform angles, that may be best.
I notice that your first and last points were the same. That was done to fit the spline as a periodic function, I presume.
x = xyz(2:end,1); xbar = mean(x);
y = xyz(2:end,2); ybar = mean(y);
z = xyz(2:end,3);
[theta,R] = cart2pol(x-xbar,y-ybar);
plot(theta,[x,y,z])
grid on
legend('x','y','z')
xlabel 'Polar angle theta (in radians)'
Simplest now might be to just build a simple fourier series model for each of x, y, and z. Alternatively, I could have continurd to work in cylindrical coordinates, but the simple Fourier model for x(theta), y(theta), z(theta) will work nicely, and the added layer of complexity of working in cylindrical coordinates is not worth the effort.
mdl = fittype('b0+a1*sin(t)+b1*cos(t)+a2*sin(2*t)+b2*cos(2*t)+a3*sin(3*t)+b3*cos(3*t)+a4*sin(4*t)+b4*cos(4*t)+a5*sin(5*t)+b5*cos(5*t)+a6*sin(6*t)+b6*cos(6*t)','indep','t')
[mdlx,goodnessx] = fit(theta,x,mdl)
The model has an RMSE of 0.0797, which is not quite as small as the noise in your data due to rounding. An extra trm or two in the model would decrease that error considerably. I'm not sure it is warranted. That would be your choice. The noise in your data due to rounding was on the oder of 10% or less of that number. Still it is probably adequate.
plot(mdlx,theta,x)
As you can see, the model fits quite nicely. As such, it gives you a nice function you can write down. As well, you can evaluate that model at a list of points in polar angle.
Next, build similar models for y and z.
mdly = fit(theta,y,mdl);
mdlz = fit(theta,z,mdl);
plot(mdly,theta,y)
plot(mdlz,theta,z)
Again, the y and z models seem entirely adequate for this data. Or you could add an extra pair of terms if you wanted more.
Each of these models are now easily written down, since they only required 13 terms to fit as Fourier series. And, we can now evaluate them easily enough.
T = (0:30:330)'; % degrees, to make it easier to follow
Trad = deg2rad(T);
table(T,mdlx(Trad),mdly(Trad),mdlz(Trad))
Finally, we can now return and plot the three curves in 3-d.
fplot3(@(t) mdlx(t),@(t) mdly(t),@(t) mdlz(t))
hold on
plot3(x,y,z,'ro')
1 comentario
mehlil ahmed
el 1 de Mzo. de 2023
Más respuestas (1)
John D'Errico
el 1 de Mzo. de 2023
Editada: John D'Errico
el 1 de Mzo. de 2023
0 votos
Can you get the spline equation? No. There is no simple equation you can write down. Or, if you could write it down, it would involves at least many dozens of coefficients. This is perhaps the most common misunderstanding made by users about splines. They seem to think that all functions have some nice, simple form they can write down in a paper or report. No, you can't, at least not in general with a spline. (In a very simple case, where the spline was built from only 3 points or so, well, yes. Then you could. But that never actually happens. You would not be bothering to use a spline then anyway.)
Can you get a set of points at specific degree increments? Not directly, no, at least, not easily. (You COULD do it, but that would essentially involve calls to a tool like fzero to generate each point on the curve.) If you had generated the curves as a function of angle, then yes, you can do that trivially, because then the curves would be represented as a function of polar angle. And there should be no reason you could not have done that, but I gather you just started with (x,y,z) coordinates from your comments. But you could easily convert (x,y) into polar coordinates, then build the splines r(theta), z(theta). Now it would be trivial to compute points on the curve at a list of angles. The point being, you generated your splines the wrong way to achieve what you seem to want.
Anyway, we don't have your data, so I cannot give an example of what you might have done differently.
1 comentario
mehlil ahmed
el 1 de Mzo. de 2023
Categorías
Más información sobre Spline Postprocessing en Centro de ayuda y File Exchange.
el 1 de Mzo. de 2023
el 1 de Mzo. de 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
