How to automatically identify and fit segments of a force-distance curve in MATLAB?

Dear All,

Background:

• The curve has multiple distinct segments.
• The goal is to automatically identify the breakpoints and fit a function to each segment.
• Manual identification is time-consuming and inaccurate.
• The data points contain force (F) and distance (Z) values.

Desired Solution:

• A program that automatically identifies the breakpoints.
• Function fitting for each segment.

Additional Information:

• A TXT file is attached with the curve data.

Question:

Does anyone have any ideas for a solution?

[Matthew Sisson]

I'd suggest using a highpass filter on the data first, then differentiating it using the "diff" function to get a plot of slope. Once you have that slope plot, you can use that to determine when the slope changes from the shallow curve to the steep curve. Once you know that index, you can split the original data into the two curves and then do a linear fit on them.

There are many ways to do this. I'm just providing one example, assuming your data is fairly consistant.

Answer 1

[Mathieu NOE]

hello

maybe this ?

I simply wanted stright lines (polyfit order 1)

I suggest x0,y0 and x1,y1 as possible candidates for the inflexion point - maybe you could use simply the median point of both

code :

fileDir = pwd; % current directory (or specify which one is the working directory)
S = dir(fullfile(fileDir,'*++.txt')); % get list of data files in directory

% % optionnal 
% S = natsortfiles(S); % sort file names into natural order , see : 
%            %(https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort)

%% main loop
for k = 1:numel(S)
    filename = S(k).name % to actually show filenames are sorted (see command window)
    out = readmatrix( fullfile(fileDir, filename)); % 
    force = out(:,1);
    dist = out(:,2);
    
    [N,EDGES] = histcounts(force,100);
    center_values = (EDGES(1:end-1) + EDGES(2:end))/2; % center value between two edges

    
    %% first segment selection
    % find force upper limit in histogram corresponding to 1 percent of max histogram value
    % after the histogram peak
    [m,id] = max(N);
    ind = find(N(id:end)<0.01*m);
    ind = ind(1) + id - 1;
    force_up = center_values(ind);
    ind = (force<=force_up);
    dist_seg1 = dist(ind);
    force_seg1 = force(ind);
    x0 = dist_seg1(end);

     %% second segment selection
     ind = (dist>x0);
    dist_seg2 = dist(ind);
    force_seg2 = force(ind);
    
    
    %% make polyfit fits 
    % segemnt 1 is easy 
    p = polyfit(dist_seg1,force_seg1,1);
    yfit1 = polyval(p,dist_seg1);
    y0 = yfit1(end);
    
    % segment 2 requires a bit more attention 
    % we remove the first 30% of the data to not be influenced by the
    % lower kinck zone
    stop = numel(dist_seg2);
    start = 1+ round(0.3*stop);
    
    p = polyfit(dist_seg2(start:stop),force_seg2(start:stop),1);
    yfit2 = polyval(p,dist_seg2);
    ind = find(yfit2>=y0);
    xfit2 = dist_seg2(ind);
    yfit2 = yfit2(ind);
    x1 = xfit2(1);
    y1 = yfit2(1);  
    
       % Plot results
        figure;
        plot(dist_seg1,force_seg1,'b',dist_seg2,force_seg2,'r');
        hold on
        plot(dist_seg1,yfit1,'g',xfit2,yfit2,'m')
        plot(x0,y0,'c+','LineWidth',1.5,'MarkerSize',17)
        plot(x1,y1,'k+','LineWidth',1.5,'MarkerSize',17)
        grid on; xlabel('X'); ylabel('Y'); 
        tit = strrep(filename,'_',' ');
        
        title(tit)
        legend('Data segment 1','Data segment 2','Fit segment 1','Fit segment 2','x0,y0','x1,y1','Location','North')

end

filename = 'D3EC_240208_0003Force_Ext++.txt'

filename = 'D3EC_240208_0004Force_Ext++.txt'

filename = 'D3EC_240208_0005Force_Ext++.txt'

[Mathieu NOE]

as the amount of data is quite large , here a slight modification of the code that allows you to introduce a downsampling factor (here 10, no impact on the results)

fileDir = pwd; % current directory (or specify which one is the working directory)
S = dir(fullfile(fileDir,'*++.txt')); % get list of data files in directory

% % optionnal 
% S = natsortfiles(S); % sort file names into natural order , see : 
%            %(https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort)

downsampling_factor = 10;
%% main loop
for k = 1:numel(S)
    filename = S(k).name % to actually show filenames are sorted (see command window)
    out = readmatrix( fullfile(fileDir, filename)); % 
    [m,n] = size(out);
    % extract dist / force data with (optionnal) downsampling
    
    force = out(1:downsampling_factor:m,1);
    dist = out(1:downsampling_factor:m,2);
    
    [N,EDGES] = histcounts(force,100);
    center_values = (EDGES(1:end-1) + EDGES(2:end))/2; % center value between two edges
    
    %% first segment selection
    % find force upper limit in histogram corresponding to 1 percent of max histogram value
    % after the histogram peak
    [m,id] = max(N);
    ind = find(N(id:end)<0.01*m);
    ind = ind(1) + id - 1;
    force_up = center_values(ind);
    ind = (force<=force_up);
    dist_seg1 = dist(ind);
    force_seg1 = force(ind);
    x0 = dist_seg1(end);

     %% second segment selection
     ind = (dist>x0);
    dist_seg2 = dist(ind);
    force_seg2 = force(ind);
    
    
    %% make polyfit fits 
    % segemnt 1 is easy 
    p = polyfit(dist_seg1,force_seg1,1);
    yfit1 = polyval(p,dist_seg1);
    y0 = yfit1(end);
    
    % segment 2 requires a bit more attention 
    % we remove the first 30% of the data to not be influenced by the
    % lower kinck zone
    stop = numel(dist_seg2);
    start = 1+ round(0.3*stop);
    
    p = polyfit(dist_seg2(start:stop),force_seg2(start:stop),1);
    yfit2 = polyval(p,dist_seg2);
    ind = find(yfit2>=y0);
    xfit2 = dist_seg2(ind);
    yfit2 = yfit2(ind);
    x1 = xfit2(1);
    y1 = yfit2(1);  
    
       % Plot results
        figure;
        plot(dist_seg1,force_seg1,'b',dist_seg2,force_seg2,'r');
        hold on
        plot(dist_seg1,yfit1,'g',xfit2,yfit2,'m')
        plot(x0,y0,'c+','LineWidth',1.5,'MarkerSize',17)
        plot(x1,y1,'k+','LineWidth',1.5,'MarkerSize',17)
        grid on; xlabel('X'); ylabel('Y'); 
        tit = strrep(filename,'_',' ');
        
        title(tit)
        legend('Data segment 1','Data segment 2','Fit segment 1','Fit segment 2','x0,y0','x1,y1','Location','North')

end

[Mathieu NOE]

hello

problem solved ?

do you mind accepting my answer ? Tx !!