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

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 ?

Answer 1

[Matt J]

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.

[VIGNESH BALAJI]

@Matt J I looked into the FEX examples and the example code you mentioned.

It looks like in this tool, it first fits the shape to get preliminary xyz co-ordinates and generates a plane out of these co-ordinates. Then used to fit a 3D shape. How do I fit a plane directly for given 3D data irrespective of the shape of the data ?

I have attached the extension of the code example, you showed me which shows the need to find the shape before the plane.

Can you let me know, how do I fit a plane directly for given 3D data irrespective of the shape of the data ?

XYZ0=ellipticalFit.xysim([x0,y0],[a,b], 0, thetaSamps); %Generate noisy measured samples
 XYZ0(3,:)=0;
 XYZ0=R(1:3,1:3)*XYZ0+t + 0.3*randn(size(XYZ0));

pFit=planarFit(XYZ0);%Preliminary plane fit

  xy0=pFit.project2D(XYZ0); %Map measured 3D samples to 2D

eFit=ellipticalFit(xy0); %Perform ellipse fit in 2D

[Matt J]

I have attached the extension of the code example, you showed me which shows the need to find the shape before the plane.

No, the extra code you've added merely generates simulated input data. It is not part of the fitting process. The plane fit is done in the suitably commented line,

pFit=planarFit(XYZ0);%Preliminary plane fit

[VIGNESH BALAJI]

@Matt J Thank you. This helped me to progress in my analysis. Can you please let me know, how do I plot the plane that fitted my data in this framework ?

I looked at the code of how the plane fit is done. I am trying to understand the plane fit code, any help in that direction is also welcome

        function obj=planarFit(XYZ)
        %planarFit constructor
        %
        %    obj=planarFit(xyz)                    
        %                            
        %IN:                                                    
        %    xyz: a 3xN  matrix whose columns are sample coordinates on a plane to be fitted.      
           
            if ~nargin, XYZ=[]; end

            if isempty(XYZ), return; end
            
            
            XYZ0=XYZ;
            
            [XYZ,T]=conicFit.homogNorm(XYZ);
            
            normal=conicFit.mostnull(XYZ.');
            normal=T(1:3,1:3)\normal;
             s=sign(normal(find(normal,1,'last')));
            normal = s*normal/norm(normal);
            
            d=dot([normal;0], T\[0;0;0;1]);
            
            R=[normal,null(normal.')];

            R(:,3)=cross(R(:,1),R(:,2));

            obj.XYZ=XYZ0;
            obj.R=R;
            obj.d=d;
 
        end

[Matt J]

plot(pFit)

[VIGNESH BALAJI]

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.

[Matt J]

scatter3 is used to plot the data points and this is then overlaid with a plot of the fittedplane using fimplicit3.

Answer 2

[Mathieu NOE]

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')

[Mathieu NOE]

this is the same idea for the second data set

I had to do some more work to "isolate" one catenary and apply the same logic

tomorrow I will try to finish my work for the second catenary (shoudn't be too complicated !!)

 
  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(x0>-0.4 & x0<-0.3);
  x = x0(n);
  y = y0(n);
  z = z0(n);
  
  % sort / unique in y ascending order
  [y,ind,~] = unique(y);
  x = x(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')
  

[Mathieu NOE]

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]

@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

Answer 3

[John D'Errico]

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.

[VIGNESH BALAJI]

@John D'Errico It is a 3D problem because although the data is of a curve, it is slightly tilted along a plane (Basically a 2D curve rotated in 3D plane, Hence it is a 3D problem). I have attached an image for clarity.

Can you please let me know any tools to solve this to fit to a known equation (In my case 2D catenary equation rotated in 3D plane)

[Mathieu NOE]

it would be more efficient to share the data

has you tried to find (fit) the plane to which the data belong ?

[VIGNESH BALAJI]

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 ?

Answer 4

[VIGNESH BALAJI]

@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)

[Mathieu NOE]

ok , sounds good but in your original request , I don't see the need to plot the 2D fitted curve back in the 3D plot

so I skipped this and focus solely on 2d curve fitting

[VIGNESH BALAJI]

@Mathieu NOE thanks, I have this already but my main problem is for each dataset, I have to do a lot of manual checking. I want to fit the 2d curves automatically without manual tunning and rotating or manual 2d fitting. Can you please let me know any tools or scripts used to do it ?

[Mathieu NOE]

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]

@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.