matlab - 在 MATLAB 中查找非线性图的最长线性部分

Question

很抱歉发了这么长的帖子，但这需要一点解释。我正在尝试制作一个脚本来找到情节的最长线性部分。示例数据位于 csv 文件中 here ，它是用于计算 3D 打印 sample 的剪切模量的应力和应变数据。到目前为止，我的代码如下:

x_data = [];
y_data = [];
x_data = Data(:,1);
y_data = Data(:,2);
plot(x_data,y_data);
grid on;

answer1 = questdlg('Would you like to load last attempt''s numbers?');
switch answer1
    case 'Yes'
        [sim_slopes,reg_data] = regr_and_longest_part(new_x_data,new_y_data,str2num(answer2{3}),str2num(answer2{2}),K);
    case 'No'
        disp('Take a look at the plot, find a range estimate, and press any button to continue');
        pause;

        prompt = {'Eliminate values ABOVE this x-value:','Eliminate values BELOW this x-value:','Size of divisions on x-axis:','Factor for similarity of slopes:'};
        dlg_title = 'Point elimination';
        num_lines = 1;
        defaultans = {'0','0','0','0.1'};
        if isempty(answer2) < 1
            defaultans = {answer2{1},answer2{2},answer2{3},answer2{4}};
        end
        answer2 = inputdlg(prompt,dlg_title,num_lines,defaultans);
        uv_of_x_range = str2num(answer2{1});
        lv_of_x_range = str2num(answer2{2});
        x_div_size = str2num(answer2{3});
        K = str2num(answer2{4});
        close all;

        iB = find(x_data > str2num(answer2{1}),1,'first');
        iS = find(x_data > str2num(answer2{2}),1,'first');

        new_x_data = x_data(iS:iB);
        new_y_data = y_data(iS:iB);

        [sim_slopes, reg_data] = regr_and_longest_part(new_x_data,new_y_data,str2num(answer2{3}),str2num(answer2{2}),K);
end

[longest_section0, Midx]= max(sim_slopes(:,4)-sim_slopes(:,3));
longest_section=1+longest_section0;

long_sec_x_data_start = x_div_size*(sim_slopes(Midx,3)-1)+lv_of_x_range;
long_sec_x_data_end = x_div_size*(sim_slopes(Midx,4)-1)+lv_of_x_range;

long_sec_x_data_start_idx=find(new_x_data >= long_sec_x_data_start,1,'first');
long_sec_x_data_end_idx=find(new_x_data >= long_sec_x_data_end,1,'first');

long_sec_x_data = new_x_data(long_sec_x_data_start_idx:long_sec_x_data_end_idx);
long_sec_y_data = new_y_data(long_sec_x_data_start_idx:long_sec_x_data_end_idx);

[b_long_sec, longes_section_reg_data] = robustfit(long_sec_x_data,long_sec_y_data);
plot(long_sec_x_data,b_long_sec(1)+b_long_sec(2)*long_sec_x_data,'LineWidth',3,'LineStyle',':','Color','k');

function [sim_slopes,reg_data] = regr_and_longest_part(x_points,y_points,x_div,lv,K)
reg_data = cell(1,3);

scatter(x_points,y_points,'.');
grid on;
hold on;

uv = lv+x_div;

ii=0;
while lv <= x_points(end)
    if uv > x_points(end)
        uv = x_points(end);
    end
    ii=ii+1;
    indices = find(x_points>lv & x_points<uv);
    temp_x_points = x_points((indices));
    temp_y_points = y_points((indices));
    if length(temp_x_points) <= 2
        break;
    end
    [b,stats] = robustfit(temp_x_points,temp_y_points);
    reg_data{ii,1} = b(1);
    reg_data{ii,2} = b(2);
    reg_data{ii,3} = length(indices);
    plot(temp_x_points,b(1)+b(2)*temp_x_points,'LineWidth',2);
    lv = lv+x_div;
    uv = lv+x_div;
end

sim_slopes = NaN(length(reg_data),4);
sim_slopes(1,:) = [reg_data{1,1},0,1,1];
idx=1;

for ii=2:length(reg_data)
    coff =sim_slopes(idx,1);
    if abs(reg_data{ii,1}-coff) <= K*coff
        C=zeros(ii-sim_slopes(idx,3)+1,1);
        for kk=sim_slopes(idx,3):ii
            C(kk)=reg_data{kk,1};
        end
        sim_slopes(idx,1)=mean(C);
        sim_slopes(idx,2)=std(C);
        sim_slopes(idx,4)=ii;
    else
        idx = idx + 1;
        sim_slopes(idx,1)=reg_data{ii,1};
        sim_slopes(idx,2)=0;
        sim_slopes(idx,3)=ii;
        sim_slopes(idx,4)=ii;
    end
end
end

很抱歉代码没有得到很好的优化，我对 MATLAB 还是比较陌生。我没有使用导数，因为我的数据相对嘈杂，推导可能会使情况变得更糟。

我设法通过将数据分成名为 x_div_size 的部分然后对每个部分执行 robustfit 来获取代码以找到绘图的最长直线部分，其结果被写入 reg_data。然后代码运行 reg_data 并通过计算图的一部分中的斜率平均值来查找哪些线具有最相似的斜率，由 K 因子确定在 sim_slopes 中记下它。然后它找到最长的间隔 max(sim_slopes(:,4)-sim_slopes(:,3)) 并对其执行回归以给出最终答案。

问题是它只会考虑遇到的第一个直线部分。绘制数据时，它有几个看起来最直的部分:

例如，当我使用 answer2 = {'0.2','0','0.0038','0.3'} 运行脚本时，我得到以下结果，其中黑线是代码找到的最直部分:

我有以下问题:

1. 很明显，从大约 x = 0.04 到 x = 0.2 有一段很长的直线部分，我不确定为什么脚本找不到它。尝试使用不同的值，脚本似乎总是选择第一个最长的直线部分，而忽略后续部分。
2. MATLAB 提示 Warning: Iteration limit reached. 因为要执行 50 多个回归。有没有办法绕过 robustfit 的这个限制？
3. 在生成 sim_slopes 时，可能有部分绘图的斜率与之前斜率的平均值相差太大，因此它被标记为长部分的结尾。但该部分有时夹在两侧的其他几个部分之间，而这些部分具有相似的斜率。如何才能告诉脚本忽略一个任性的部分并继续，就好像它在 K 值允许的公差范围内一样？

Answer 1

看看the Douglas-Peucker algorithm .如果您将 (x,y) 值视为一个(开放的)多边形的顶点，该算法将为您简化它，这样从简化的多边形到原始多边形的最大距离小于您可以选择的某个阈值。简化的多边形将是一组直线。找到相距最远的两个顶点，您就完成了。

MATLAB 在映射工具箱中有一个名为 reducem 的实现.您可能还会在文件交换上找到一个实现(但要小心，那里也有非常糟糕的代码)。或者，您可以推出自己的算法，这是一个非常简单的算法。