How to fit an equation (catenary equation in 2D) to a 3D dataset (experimental data of a catenary curve in 3D) ?

21 次查看（过去 30 天）
VIGNESH BALAJI
VIGNESH BALAJI 2023-7-20
评论： Matt J 2023-7-28
I have an experimental data (3D data) of a catenary curve (basically a rope with 2 attached ends and hanging freely, data got from markers on them).
I want to fit this data with a catenary equation y = a*cosh(x/a).
a is parameterized by an equation stating the length of the physical cable, (length/2) - a* sinh(x/a) = 0
I tried using curve fitting tool and it is asking me an equation in terms of x,y and z. If I give such an equation then it is a surface.
Further, I don't know how to give an equation to parameterize parameter a in the curve fitting tool.
My ultimate requirement is, I want to fit a line (2D equation) in 3D space.
Can anyone please let me know the right tool for this ?

请先登录，再回答此问题。

回答（4 个）

Matt J
Matt J 2023-7-20
编辑：Matt J 2023-7-20
In the Examples tab of this FEX page,
https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics
there is a section Fitting a 2D shape to 3D Points. In particular these lines of code coould easily be adapted to your problem,
pFit=planarFit(XYZ0);%Preliminary plane fit
xy0=pFit.project2D(XYZ0); %Map measured 3D samples to 2D
eFit=ellipticalFit(xy0); %Perform ellipse fit in 2D
XYZ=pFit.unproject3D( cell2mat(eFit.sample(1:360)) ); %Post-sample the ellipse fit and map back to 3D
except that instead of the ellipse fit, you would fit with your own model.
  6 个评论
显示 4更早的评论
VIGNESH BALAJI
VIGNESH BALAJI 2023-7-27
Thank you @Matt J it worked. I was trying it with plot3 initially instead of plot and internally the structure of pFit is like this, like the image attached below. How does it plot using this structure, can you explain me how.

请先登录，再进行评论。

Mathieu NOE
Mathieu NOE 2023-7-20
This is a lazy (tired) guy answer
i figure out I could easily transform the 3D array into 2D then do a parabola fit (yes I know it's only an approximation for the catenary cosh model)
load('export_dataset1.mat')
x = cat1_2_x';
y = cat1_2_y';
z = cat1_2_z';
m = numel(x);
figure(1),
plot3(x,y,z,'*-');
xlabel('X')
ylabel('Y')
zlabel('Z')
P = [x y z]; % create array of x,y,z coordinates
V = diff(P,1); % array of vectors
d = sqrt(sum(V.^2,2)); % array of vectors length
% check distance between first and last point (debug)
dmax = P(m,:)-P(1,:);
dmax3D = sqrt(sum(dmax.^2))
% compute raltive angle between successive vectors
ang(1) = 0;
for k = 1:m-2
ang(k+1) = acos(dot(V(k,:)/norm(V(k,:)),V(k+1,:)/norm(V(k+1,:))));
end
angc = cumsum(ang);
angc = angc - angc(end)/2; % so the 2plot shows a balanced (symetrical) start and end angles
% create the 2D points by adding vector coordinates
x(1)= 0;
y(1) = 0;
for k = 1:m-1
x(k+1) = x(k)+d(k)*cos(angc(k));
y(k+1) = y(k)+d(k)*sin(angc(k));
end
% check distance between first and last point (debug)
P = [x y];
dmax = P(m,:)-P(1,:);
dmax2D = sqrt(sum(dmax.^2))
% this is the same as the original data
% now you can fit your 2D data
% as I am lazy and it's getting late here I am doing a parabola fit (and
% not a true cosh model fit
%Set up the appropriate matrix A to find the best-fit parabola of the form y=C+Dx+Ex^2. The
%first column of A will contain all 1's, using the ones() command. The second column of A
%contains x values that are stored in X. The third column of A contains the squared x values
%that are stored in X. Elementwise multiplication of X by itself, using .* operator, will
%produce the desired values for the third column.
A = [ones(m,1) x x.*x];
A_transposeA = A.' * A;
A_transposeY = A.' * y;
%backslash operation to solve the overdetermined system.
Soln2 = A_transposeA\A_transposeY;
%x values to use for plotting the best-fit parabola.
xfit = linspace(min(x),max(x),100);
yfit = Soln2(1) + Soln2(2)*xfit + Soln2(3)*xfit.^2; %
figure(2),
plot(xfit, yfit, x, y, 'r*','Markersize',15);
legend('fit','data');
xlabel('X')
ylabel('Y')
  3 个评论
显示 1更早的评论
Mathieu NOE
Mathieu NOE 2023-7-20
well
that was much faster than I thougt
second catenary results
load('export_dataset2.mat')
x0 = cat6_1_x';
y0 = cat6_1_y';
z0 = cat6_1_z';
figure(1),
plot3(x0,y0,z0,'*-');
xlabel('X')
ylabel('Y')
zlabel('Z')
hold on
% isolate one catenary data (the one which is parrallel to the Y
% direction)
n= find(y0>-0.5 & y0<-0.4);
x = x0(n);
y = y0(n);
z = z0(n);
% sort / unique in x ascending order
[x,ind,~] = unique(x);
y = y(ind);
z = z(ind);
m = numel(x);
plot3(x,y,z,'*r','Markersize',15);
xlabel('X')
ylabel('Y')
zlabel('Z')
hold off
P = [x y z]; % create array of x,y,z coordinates
V = diff(P,1); % array of vectors
d = sqrt(sum(V.^2,2)); % array of vectors length
% check distance between first and last point (debug)
dmax = P(m,:)-P(1,:);
dmax3D = sqrt(sum(dmax.^2))
% compute raltive angle between successive vectors
ang(1) = 0;
for k = 1:m-2
ang(k+1) = acos(dot(V(k,:)/norm(V(k,:)),V(k+1,:)/norm(V(k+1,:))));
end
angc = cumsum(ang);
angc = angc - angc(end)/2; % so the 2plot shows a balanced (symetrical) start and end angles
% create the 2D points by adding vector coordinates
x(1)= 0;
y(1) = 0;
for k = 1:m-1
x(k+1) = x(k)+d(k)*cos(angc(k));
y(k+1) = y(k)+d(k)*sin(angc(k));
end
% check distance between first and last point (debug)
P = [x y];
dmax = P(m,:)-P(1,:);
dmax2D = sqrt(sum(dmax.^2))
% this is the same as the original data
% now you can fit your 2D data
% as I am lazy and it's getting late here I am doing a parabola fit (and
% not a true cosh model fit
%Set up the appropriate matrix A to find the best-fit parabola of the form y=C+Dx+Ex^2. The
%first column of A will contain all 1's, using the ones() command. The second column of A
%contains x values that are stored in X. The third column of A contains the squared x values
%that are stored in X. Elementwise multiplication of X by itself, using .* operator, will
%produce the desired values for the third column.
A = [ones(m,1) x x.*x];
A_transposeA = A.' * A;
A_transposeY = A.' * y;
%backslash operation to solve the overdetermined system.
Soln2 = A_transposeA\A_transposeY;
%x values to use for plotting the best-fit parabola.
xfit = linspace(min(x),max(x),100);
yfit = Soln2(1) + Soln2(2)*xfit + Soln2(3)*xfit.^2; %
figure(2),
plot(xfit, yfit, x, y, 'r*','Markersize',15);
legend('fit','data');
xlabel('X')
ylabel('Y')
VIGNESH BALAJI
VIGNESH BALAJI 2023-7-20
@Mathieu NOE If I am correct you first found a plane where the 2D data lie in the 3D plane and then fitted the 2D data.
If I am correct these lines of code below is same as plane fit -
P = [x y z]; % create array of x,y,z coordinates
V = diff(P,1); % array of vectors
d = sqrt(sum(V.^2,2)); % array of vectors length
% check distance between first and last point (debug)
dmax = P(m,:)-P(1,:);
dmax3D = sqrt(sum(dmax.^2))
% compute raltive angle between successive vectors
ang(1) = 0;
for k = 1:m-2
ang(k+1) = acos(dot(V(k,:)/norm(V(k,:)),V(k+1,:)/norm(V(k+1,:))));
end
angc = cumsum(ang);
angc = angc - angc(end)/2; % so the 2plot shows a balanced (symetrical) start and end angles
% create the 2D points by adding vector coordinates
x(1)= 0;
y(1) = 0;
for k = 1:m-1
x(k+1) = x(k)+d(k)*cos(angc(k));
y(k+1) = y(k)+d(k)*sin(angc(k));
end

请先登录，再进行评论。

John D'Errico
John D'Errico 2023-7-20
编辑：John D'Errico 2023-7-20
  1. How is this a 3-d probem? The curve will lie in a plane. So it is just a 2-d problem. That you have measurements in x,y,z is irrelevant. So reduce it to TWO dimensions first, along the line in the (x,y) plane between the two endpoints.
  2. Then fit the curve, as a function of the new variable along that line between the endpoints. WTP?
  3. Finally, Since you have the equation of the line between the endpoints, convert it back into 3-d, if that is what you really wanted.
  3 个评论
显示 1更早的评论
VIGNESH BALAJI
VIGNESH BALAJI 2023-7-20
Sure, I am sharing the datasets. I am also adding an image of the datasets, you can load the mat files and get their x,y and z co-ordinates.
Dataset 1 conatins one catenary
Dataset 2 contains 2 catenaries
I haven't tried fitting a plane. I was just thinking about it. What could be a good tool to fit a plane, You mean PlanarFit or some other tool to fit a plane ?

请先登录，再进行评论。

VIGNESH BALAJI
VIGNESH BALAJI 2023-7-20
编辑：VIGNESH BALAJI 2023-7-20
@Mathieu NOE thanks a lot for trying this problem even when you are tired.
I also wrote a script in the meantime where I seperated each catenary curve and fitted and fitted a 2d curve,
then i manually rotated and level shiffted the plane to fit the 3d data and to show the results in a 3d fit format.
I have put the code below. This process was very manually intensive as I have to manuallty check the fit.
I think your work showed me an inspiration to find a plane for 3d data and to curve fit in the 2d plane, then I can rotate this plane back to 3d original space to show the fit in 3d space.
I have put down my code for catenary fit (attached images) -
%% curve fitting with catenary equation
% data
xco = cat1_2_y; % data rearranged and not necessary
yco = cat1_2_z;
zco = cat1_2_x;
%length of catenary used in actual setup
l = 1.76; % measured from actual cable used in experiment
span = (abs(cat1_2_y(:,1)) + abs(cat1_2_y(:,6))); % span tells how wide the curve(catenary) is opened up
% finding paramegter "a" to find the sag of the curve
fun = @(a) (- l/2 + a * sinh(span/(2*a))); x0 = 0.15;
a = abs(fzero(fun,x0));
xm = linspace(xco(:,1), xco(:,6),100); % divide x data into 100 divisions to get a smooth curve
size_xm = size(xm)
ym = zeros(1,size_xm(2));
for k=1:size_xm(2)
ym(1,k) = a*cosh((xm(1,k)/a)) + (0.4965-a); % catenary equation in 2d, the plus term is a scaling factor as the curve starts from zero height and actual data is at a particular height
end
zm = linspace(zco(:,1), zco(:,6),100); % divide z data into 100 divisions to get a smooth curve
zm = zm + 0.04; % 4 cm z error correction, mismatch between experiment and analytical plotting
%% rotate this generated 2d catenary curve in 3d space to fit the 3d experimental data
new = zeros(3,size_xm(2));
for i=1:size_xm(2)
new(:,i) = rotx(-1.9)*roty(0)*rotz(0)*[xm(:,i); ym(:,i); zm(:,i)]; % x rotation, done manually by seeing plots
end
x_new = new(1,:);
y_new = new(2,:);
z_new = new(3,:);
%% plotting results for comparisoon pf data and fit
plot3(xco(:), yco(:), zco(:), 'o'); % experimental data
hold on
plot3(xm(:), ym(:), zm(:)) % generated catenary before rotation
xlabel('x axis (m)')
ylabel('y axis (m)')
zlabel('z axis (m)')
plot3(x_new(:), y_new(:), z_new(:), 'g') % generated catenary after rotation
legend('Experimental data', 'Analytical expression before x rotation and static error', 'Analytical expression after x rotation and static error')
length = 2* a*sinh(span/(2*a)); % check length of curve is same as length in experiment
%% fit evaluation, do meansquare error to see goodness of the fit
% data rearranging
exp_data = [xco;yco;zco];
ana_x = [x_new(:,1), x_new(:,9), x_new(:,23), x_new(:,49), x_new(:,92), x_new(:,100)];
ana_y = [y_new(:,1), y_new(:,9), y_new(:,23), y_new(:,49), y_new(:,92), y_new(:,100)];
ana_z = [z_new(:,1), z_new(:,9), z_new(:,23), z_new(:,49), z_new(:,92), z_new(:,100)];
ana_data = [ana_x; ana_y; ana_z];
error_x = immse(xco, ana_x)
error_y = immse(yco, ana_y)
error_z = immse(zco, ana_z)
% E = rmse( exp_data, new, "all");
%disp(err)
  4 个评论
显示 2更早的评论
Mathieu NOE
Mathieu NOE 2023-7-21
well , I don't know howyu could make all this "automatic"
myself I had to do a few manual operations especially for the two catenary data set , to extract each one
then today I have no time for the 3D "backwards" projection topic (sorry)
VIGNESH BALAJI
VIGNESH BALAJI 2023-7-21
@Mathieu NOE I think I didn't explain it clearly. I will explain it now.
I have a script to extract catenaries (that is fine).
I want the fitting part of catenaries in 3D space as automatic instead of manual.
The automatic part of fitting could be a tool or a standard one time written script.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Delaunay Triangulation 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by