How to extend the line of 2 PolyFit from either side to intersect and get a combined fit line

Question

I am trying to get combined fit line made from two linear polyfit from either side (should intersect), here is the picture of fit lines:

I am trying to make the two fit (blue) lines intersect and produce a combined fit line as shown in the picture below:

Note that the crest can happen anywhere so I cannot assume to be in the center.

Here is the code that creates the first plot:

xdatPart1 = R;
zdatPart1 = z;

n = 3000;
ln = length(R);

[sX,In] = sort(R,1);

sZ = z(In);       

xdatP1 = sX(1:n,1);
zdatP1 = sZ(1:n,1);

n2 = ln - 3000;

xdatP2 = sX(n2:ln,1);
zdatP2 = sZ(n2:ln,1);

pp1 = polyfit(xdatP1,zdatP1,1);
pp2 = polyfit(xdatP2,zdatP2,1);

ff1 = polyval(pp1,xdatP1);
ff2 = polyval(pp2,xdatP2);

xDat = [xdatPart1];
zDat = [zdatPart1];

axes(handles.axes2);
cla(handles.axes2);
plot(xdatPart1,zdatPart1,'.r');
hold on
plot(xdatP1,ff1,'.b');
plot(xdatP2,ff2,'.b');
xlabel(['R ',units]);
ylabel(['Z ', units]);
grid on
hold off

Answer 1

Below's a rough implementation with no curve fitting toolbox. Although the code should be self-explanatory, here's an outline of the algorithm:

1. We generate some data.
2. We estimate the intersection point by smoothing the data and finding the location of the maximum value.
3. We fit a line to each side of the estimated intersection point.
4. We compute the intersection of the fitted lines using the fitted equations.
5. We use mkpp to construct a function handle to an "evaluateable" piecewise polynomial.
6. The output, ppfunc, is a function handle of 1 variable, that you can use just like any regular function.

Now, this solution is not optimal in any sense (such as MMSE, LSQ, etc.) but as you will see in the comparison with the result from MATLAB's toolbox, it's not that bad!

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   

Zoomed in around the intersection: