Как продлить линию 2 PolyFit с любой стороны, чтобы пересечься и получить комбинированную линию посадки

Ниже приведена грубая реализация без набора инструментов для подгонки кривой. Хотя код не требует пояснений, вот схема алгоритма:

1. Мы генерируем некоторые данные.
2. Мы оцениваем точку пересечения, сглаживая данные и находя положение максимального значения.
3. Мы подгоняем линию к каждой стороне предполагаемой точки пересечения.
4. Мы вычисляем пересечение подобранных линий, используя подобранные уравнения.
5. Мы используем mkpp для создания дескриптора функции для "вычисляемого" кусочного полинома.
6. Вывод, ppfunc, представляет собой дескриптор функции 1 переменной, которую вы можете использовать точно так же, как любая обычная функция.

Теперь это решение не является оптимальным ни в каком смысле (например, MMSE, LSQ и т. д.), но, как вы увидите в сравнении с результатом из набора инструментов MATLAB, оно не так уж и плохо!

function ppfunc = q40160257
%% Define the ground truth:
center_x = 6 + randn(1);
center_y = 78.15 + 0.01 * randn(1);
% Define a couple of points for the left section
leftmost_x = 0;
leftmost_y = 78.015 + 0.01 * randn(1);
% Define a couple of points for the right section
rightmost_x = 14.8;
rightmost_y = 78.02 + 0.01 * randn(1);
% Find the line equations:
m1 = (center_y-leftmost_y)/(center_x-leftmost_x);
n1 = getN(leftmost_x,leftmost_y,m1);
m2 = (rightmost_y-center_y)/(rightmost_x-center_x);
n2 = getN(rightmost_x,rightmost_y,m2);
% Print the ground truth:
fprintf(1,'The line equations are: {y1=%f*x+%f} , {y2=%f*x+%f}\n',m1,n1,m2,n2)
%% Generate some data:
NOISE_MAGNITUDE = 0.002;
N_POINTS_PER_SIDE = 1000;
x1 = linspace(leftmost_x,center_x,N_POINTS_PER_SIDE);
y1 = m1*x1+n1+NOISE_MAGNITUDE*randn(1,numel(x1));
x2 = linspace(center_x,rightmost_x,N_POINTS_PER_SIDE);
y2 = m2*x2+n2+NOISE_MAGNITUDE*randn(1,numel(x2));
X = [x1 x2(2:end)]; Y = [y1 y2(2:end)];
%% See what we have:
figure(); plot(X,Y,'.r'); hold on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimating the intersection point:
MOVING_AVERAGE_PERIOD = 10; % Play around with this value.
smoothed_data = conv(Y, ones(1,MOVING_AVERAGE_PERIOD)/MOVING_AVERAGE_PERIOD, 'same');
plot(X, smoothed_data, '-b'); ylim([floor(leftmost_y*10) ceil(center_y*10)]/10);
[~,centerInd] = max(smoothed_data);
fprintf(1,'The real intersection is at index %d, the estimated is at %d.\n',...
           N_POINTS_PER_SIDE, centerInd);
%% Fitting a polynomial to each side:
p1 = polyfit(X(1:centerInd),Y(1:centerInd),1);
p2 = polyfit(X(centerInd+1:end),Y(centerInd+1:end),1);
[x_int,y_int] = getLineIntersection(p1,p2);
plot(x_int,y_int,'sg');

pp = mkpp([X(1) x_int X(end)],[p1; (p2 + [0 x_int*p2(1)])]);
ppfunc = @(x)ppval(pp,x);
plot(X, ppfunc(X),'-k','LineWidth',3)
legend('Original data', 'Smoothed data', 'Computed intersection',...
       'Final piecewise-linear fit');
grid on; grid minor;     

%% Comparison with the curve-fitting toolbox:
if license('test','Curve_Fitting_Toolbox')
  ft = fittype( '(x<=-(n2-n1)/(m2-m1))*(m1*x+n1)+(x>-(n2-n1)/(m2-m1))*(m2*x+n2)',...
                'independent', 'x', 'dependent', 'y' );
  opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
  % Parameter order: m1, m2, n1, n2:
  opts.StartPoint = [0.02 -0.02 78 78];
  fitresult = fit( X(:), Y(:), ft, opts);
  % Comparison with what we did above:
  fprintf(1,[...
    'Our solution:\n'...
    '\tm1 = %-12f\n\tm2 = %-12f\n\tn1 = %-12f\n\tn2 = %-12f\n'...
    'Curve Fitting Toolbox'' solution:\n'...
    '\tm1 = %-12f\n\tm2 = %-12f\n\tn1 = %-12f\n\tn2 = %-12f\n'],...
    m1,m2,n1,n2,fitresult.m1,fitresult.m2,fitresult.n1,fitresult.n2);    
end

%% Helper functions:
function n = getN(x0,y0,m)
% y = m*x+n => n = y0-m*x0;
n = y0-m*x0;

function [x_int,y_int] = getLineIntersection(p1,p2)
% m1*x+n1 = m2*x+n2 => x = -(n2-n1)/(m2-m1)
x_int = -(p2(2)-p1(2))/(p2(1)-p1(1));
y_int = p1(1)*x_int+p1(2);

The result (sample run):

Our solution:
    m1 = 0.022982    
    m2 = -0.011863   
    n1 = 78.012992   
    n2 = 78.208973   
Curve Fitting Toolbox' solution:
    m1 = 0.022974    
    m2 = -0.011882   
    n1 = 78.013022   
    n2 = 78.209127   

Увеличено вокруг перекрестка: