May I ask if there are any good methods in MATLAB that can fit broken lines

Is there any good method in MATLAB to fit broken lines

[DGM]

What exactly do you mean by "fit broken lines"? What are the inputs, and what are the expected outputs? If there are relevant example files/data that you can share, you can attach them.

[tabf]

Due to the fact that I have not yet completed the preprocessing of the data, I am unable to provide an accurate data source. However, the desired result is that I have obtained the original point cloud data of the polygonal line type weld seam, and I am currently working hard to extract the weld seam feature points. Finally, I need to perform a polygonal line fitting on the obtained weld seam feature points.

[Image Analyst]

OK, no rush on our part but post your data when you're ready. Will your point set be in 2-D (x,y) or 3-D (x,y,z)?

In a message a few days ago @John D'Errico mentioned the "CRUST" algorithm. It might be applicable to your case.

[tabf]

I need to process 3D data, including xyz. Could you please send me the algorithm link you just mentioned，thanks

[Mathieu NOE]

hello again @tabf

does it has some relationship with your other post ? 3D curve fitting - MATLAB Answers - MATLAB Central (mathworks.com)

[DGM]

This was the thread in question. See the comments following the main answer.

https://www.mathworks.com/matlabcentral/answers/1984884-how-to-compute-intersection-of-union-for-two-polygons

The goal there was to take an unordered set of points and sort them into an ordered polyline/polygon.

[tabf]

Hehe, this is indeed related to my other question. There are a total of three types of line segments that I need to fit in 3-D space: straight lines, S-curves, and polygonal lines. Thank you very much for your help

[tabf]

This is the latest points

[John D'Errico]

This is your data?

So a broken path in 3-d, but also very noisy, where you don't know the break points?

In this case, the points all seem to fall exactly in a plane, but an inclined plane.

No, CRUST would not apply there. You want to find the 4 line segments essentially. Will the points ALWAYS lie in a perfect plane? If they do, then the first thing I would do is reduce it to a planar problem. Even then though, there is no simple tool that will extract the 4 line segments, and decide where the break points lie. It could be done of course, but the code would need to be written.

Is the variability always in a VERY specific direction, as it would appear form that data. (That alone seems a bit strange.)

Answer 1

[Mathieu NOE]

hello again @tabf

this would be my suggestion - of course it has some ressemblance to my answer in your other post.

when looking at your data , we can easily see that z is a linear function of both x and y , so a simple order 1 polynomial fit is enough here.

it's still using the excellent polyfitn for the heavy lifting tasks and also the smoothing of the noisy y data is now performed with the super smoother ! you can have it for free here

Supersmoother - File Exchange - MATLAB Central (mathworks.com)

(but you can get reasonnably good results also with the regular smoothdata)

hope it helps

N = readmatrix('zx.txt');
x = N(:, 1);
y = N(:, 2);
z = N(:, 3);

order = 1;
p = polyfitn([x,y],z,order);
% % FEX : https://fr.mathworks.com/matlabcentral/fileexchange/34765-polyfitn?s_tid=ta_fx_results
% % The result can be converted into a symbolic form to view the model more simply. 
% % Here I'll use the sympoly toolbox, but there is also a polyn2sym function provided.
% % FEX : https://fr.mathworks.com/matlabcentral/fileexchange/9577-symbolic-polynomial-manipulation?s_tid=srchtitle

% polyn2sympoly(p)
% ans =
%     0.40796*X1 - 0.95562*X2 + 429.1794
    
pC = p.Coefficients; % get the polynomial coefficients
pTerms = p.ModelTerms;
% create clean smooth x,y data 
xf = linspace(min(x),max(x),200);
yf = interp1(x,y,xf);

% smooth a bit yf
% pick the method / solution you prefer
% yf = smoothdata(yf,'loess',25);
% yf = smoothdata(yf,'gaussian',15);
yf = supsmu((1:numel(yf)),yf); % % Fex : https://fr.mathworks.com/matlabcentral/fileexchange/17986-supersmoother?s_tid=ta_fx_results

% create the polynomial model (z = f(x,y))
zf = 0;
for k = 1:numel(pC)
    zf = zf + pC(k)*(xf.^pTerms(k,1)).*(yf.^pTerms(k,2));
end

  plot3(x,y,z,'r*',xf,yf,zf,'k','linewidth',2)
  xlabel('X');
  ylabel('Y');
  zlabel('Z');
  legend('raw data','fitted line');
 
  axis tight square

[Mathieu NOE]

hello again

this is a second iteration, if you prefer having y as straight segments and no smoothed corners as in my first answer

here I am suing lsq_lut_piecewise.m function available here :

Piecewise linear least square fit - File Exchange - MATLAB Central (mathworks.com)

and extrema.m is attached (there are many alternatives to choose from)

N = readmatrix('zx.txt');
x = N(:, 1);
y = N(:, 2);
z = N(:, 3);

order = 1;
p = polyfitn([x,y],z,order);
% % FEX : https://fr.mathworks.com/matlabcentral/fileexchange/34765-polyfitn?s_tid=ta_fx_results
% % The result can be converted into a symbolic form to view the model more simply. 
% % Here I'll use the sympoly toolbox, but there is also a polyn2sym function provided.
% % FEX : https://fr.mathworks.com/matlabcentral/fileexchange/9577-symbolic-polynomial-manipulation?s_tid=srchtitle

% polyn2sympoly(p)
% ans =
%     0.40796*X1 - 0.95562*X2 + 429.1794

pC = p.Coefficients; % get the polynomial coefficients
pTerms = p.ModelTerms;
% create clean smooth x,y data 
xf = linspace(min(x),max(x),numel(x));
yf = interp1(x,y,xf);

% smooth a bit yf
% pick the method / solution you prefer
% yf = smoothdata(yf,'sgolay',25);
% yf = smoothdata(yf,'gaussian',15);
ys = supsmu((1:numel(yf)),yf); % % Fex : https://fr.mathworks.com/matlabcentral/fileexchange/17986-supersmoother?s_tid=ta_fx_results

% create some linear fit segments if you don't want y to be a smooth curve
[xmax,imax,xmin,imin] = extrema(ys);
s = sort([imax imin]);
XI = xf(s);
% % obtain vector of 1-D look-up table "y" points
YI = lsq_lut_piecewise( xf(:), yf(:), XI ); % you can choose to make the fit on the raw data (yf) or the smooted data (ys) 
yy = interp1(XI,YI,xf);

figure(1),
plot(x,y,'r*',xf,ys,'b',xf,yy,'m');


% create the polynomial model (z = f(x,y))
zf = 0;
for k = 1:numel(pC)
    %     zf = zf + pC(k)*(xf.^pTerms(k,1)).*(ys.^pTerms(k,2)); % for smoothed y display
    zf = zf + pC(k)*(xf.^pTerms(k,1)).*(yy.^pTerms(k,2)); % for straight lines y display
end

figure(2),
%   plot3(x,y,z,'r*',xf,ys,zf,'k','linewidth',2); % for smoothed y display
plot3(x,y,z,'r*',xf,yy,zf,'k','linewidth',2);  % for straight lines y display
xlabel('X');
ylabel('Y');
zlabel('Z');
legend('raw data','fitted line');

axis tight square

Answer 2

[Bruno Luong]

Use my FEX BSFK

xyz=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1417239/zx.txt')

xyz0 = mean(xyz,1);
xyzc = xyz-xyz0;
[~,~,V] = svd(xyzc);
xyzr = xyzc*V;
pp = BSFK(xyzr(:,1),xyzr(:,2),2);

% Generate broken lines
s = pp.breaks.';
t = ppval(pp,s);
u = 0*t; % One can do better here if the data is not on a plane
xyzl = [s,t,u]*V'+xyz0;

% Graphical check
plot3(xyz(:,1),xyz(:,2),xyz(:,3),'.r', ...
      xyzl(:,1),xyzl(:,2),xyzl(:,3),'-b','LineWidth', 2)

[tabf]

Thank you very much for your reply, but after obtaining the fitted polyline, my idea is to perform interpolation processing as I need to send this path to the robot, so I need to know the expressions of these line segments clearly. Can you provide further feedback?

[Bruno Luong]

xyz=readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1417239/zx.txt')

xyz0 = mean(xyz,1);
xyzc = xyz-xyz0;
[~,~,V] = svd(xyzc);
xyzr = xyzc*V;
pp = BSFK(xyzr(:,1),xyzr(:,2),2);

% Generate broken lines
s = pp.breaks.';
t = ppval(pp,s);
u = 0*t; % One can do better here if the data is not on a plane
xyzl = [s,t,u]*V'+xyz0;

% Graphical check
plot3(xyz(:,1),xyz(:,2),xyz(:,3),'.r', ...
      xyzl(:,1),xyzl(:,2),xyzl(:,3),'-b','LineWidth', 2); hold on

The 3D coordinates of the 5 break points of the broken lines are

xyzl

xyzl =

  -22.4698   15.0194  405.6614
    1.0716   17.2526  413.1280
   23.8660   15.0848  424.5020
   47.6580   17.2998  432.0883
   69.4099   15.1267  443.0419

And the interpolation of coordinates (here 100 points) are:

t = linspace(0,1,100);
tb = linspace(0,1,size(xyzl,1));
xyzi = interp1(tb, xyzl, t),
plot3(xyzi(:,1),xyzi(:,2),xyzi(:,3),'b+')
legend('data','broken line', '100 waitpoints');

ans =

  -22.4698   15.0194  405.6614
  -21.5186   15.1096  405.9631
  -20.5674   15.1998  406.2648
  -19.6163   15.2901  406.5664
  -18.6651   15.3803  406.8681
  -17.7139   15.4705  407.1698
  -16.7628   15.5608  407.4715
  -15.8116   15.6510  407.7732
  -14.8604   15.7412  408.0749
  -13.9093   15.8315  408.3765
  -12.9581   15.9217  408.6782
  -12.0069   16.0119  408.9799
  -11.0558   16.1021  409.2816
  -10.1046   16.1924  409.5833
   -9.1534   16.2826  409.8849
   -8.2023   16.3728  410.1866
   -7.2511   16.4631  410.4883
   -6.3000   16.5533  410.7900
   -5.3488   16.6435  411.0917
   -4.3976   16.7338  411.3934
   -3.4465   16.8240  411.6950
   -2.4953   16.9142  411.9967
   -1.5441   17.0045  412.2984
   -0.5930   17.0947  412.6001
    0.3582   17.1849  412.9018
    1.3018   17.2307  413.2429
    2.2228   17.1431  413.7025
    3.1438   17.0555  414.1620
    4.0648   16.9679  414.6216
    4.9858   16.8803  415.0811
    5.9068   16.7928  415.5407
    6.8277   16.7052  416.0002
    7.7487   16.6176  416.4598
    8.6697   16.5300  416.9193
    9.5907   16.4424  417.3789
   10.5117   16.3548  417.8384
   11.4327   16.2672  418.2980
   12.3537   16.1796  418.7575
   13.2747   16.0921  419.2171
   14.1957   16.0045  419.6766
   15.1166   15.9169  420.1362
   16.0376   15.8293  420.5958
   16.9586   15.7417  421.0553
   17.8796   15.6541  421.5149
   18.8006   15.5665  421.9744
   19.7216   15.4789  422.4340
   20.6426   15.3914  422.8935
   21.5636   15.3038  423.3531
   22.4846   15.2162  423.8126
   23.4055   15.1286  424.2722
   24.3467   15.1296  424.6552
   25.3080   15.2190  424.9617
   26.2693   15.3085  425.2682
   27.2306   15.3980  425.5748
   28.1919   15.4875  425.8813
   29.1531   15.5770  426.1878
   30.1144   15.6665  426.4943
   31.0757   15.7560  426.8008
   32.0370   15.8455  427.1073
   32.9983   15.9350  427.4139
   33.9596   16.0245  427.7204
   34.9209   16.1140  428.0269
   35.8822   16.2035  428.3334
   36.8435   16.2930  428.6399
   37.8048   16.3825  428.9465
   38.7661   16.4720  429.2530
   39.7274   16.5615  429.5595
   40.6886   16.6510  429.8660
   41.6499   16.7405  430.1725
   42.6112   16.8300  430.4790
   43.5725   16.9195  430.7856
   44.5338   17.0090  431.0921
   45.4951   17.0985  431.3986
   46.4564   17.1880  431.7051
   47.4177   17.2775  432.0116
   48.3172   17.2340  432.4202
   49.1960   17.1462  432.8628
   50.0749   17.0584  433.3053
   50.9538   16.9706  433.7479
   51.8326   16.8828  434.1905
   52.7115   16.7950  434.6330
   53.5903   16.7072  435.0756
   54.4692   16.6194  435.5182
   55.3481   16.5316  435.9608
   56.2269   16.4438  436.4033
   57.1058   16.3560  436.8459
   57.9847   16.2682  437.2885
   58.8635   16.1804  437.7310
   59.7424   16.0926  438.1736
   60.6213   16.0048  438.6162
   61.5001   15.9170  439.0588
   62.3790   15.8292  439.5013
   63.2579   15.7414  439.9439
   64.1367   15.6536  440.3865
   65.0156   15.5658  440.8291
   65.8944   15.4780  441.2716
   66.7733   15.3901  441.7142
   67.6522   15.3023  442.1568
   68.5310   15.2145  442.5993
   69.4099   15.1267  443.0419

[tabf]

Thank you very much for your reply. The method you proposed can indeed be applied to my current situation.