plotting two lines intersecting at a certain point

Hello friends, I have got some data and i want to draw two lines of best fit, one from the left and the other from the right that intersect at certain point. I need to read that point for further analysis. But am not sure how to go about that. Here is my data and some ploting

>>Temp=[25 26 26 26.5 27 27 27.5 28 28 28.5 29 29 29 29 29 29 29 29 29 29 29 29 29 28.5 28.5 28]';
>> vol=[0:25]';
>> scatter(vol,Temp);
>>a=gca;ax=axis;ax(3:4)=[25 32];a.YLim=ax(3:4)

[Rik]

This data has two pivots: one at 10 and another at 22. Which one do you want?

My first try would probably be to loop over all possible pivots and see which pivot will result in the best fit for the two lines. I suspect polyfit will already be enough, your data looks clean enough.

[Samuel Katongole]

Thanks guys for that.

But I am looking for two lines: one should best fit the points before the first pivot, which is 10, and the other to best fit the points after the second pivot, that is, 22. The intersection of these two lines is what I am interested in

Answer 1

[Image Analyst]

You can scan the data, fitting each portion (left and right) to a line. Then find where the difference in slopes of the lines is greatest -- that is the crossing point. See demo below, and attached.

% Find optimal pair of lines to fit noisy data, one line on left side and one line on right side.  Separation/crossing x value is identified.
clc;	% Clear command window.
clear;	% Delete all variables.
close all;	% Close all figure windows except those created by imtool.
workspace;	% Make sure the workspace panel is showing.
fontSize = 18;
markerSize = 20;
%============================================================================================================
% FIRST CREATE X-Y DATA. 
Temp=[25 26 26 26.5 27 27 27.5 28 28 28.5 29 29 29 29 29 29 29 29 29 29 29 29 29 28.5 28.5 28]';
vol=[0:25]';
% Rename to x and yData
x = vol;
yData = Temp;
% Done defining sample data.
%============================================================================================================
% Plot lines.
subplot(2, 2, 1);
plot(x, yData, 'b.', 'MarkerSize', markerSize);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
title('Noisy Y vs. X data', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
%================================================================================================================================================
% NOW SCAN ACROSS DATA FITTING LEFT AND RIGHT PORTION TO LINES.
% Assume the crossing point will be somewhere in the middle half of the points.
% If you go much more outward than that there are too few points to get a good line and the slopes of the lines will vary tremendously.
% Fit a line through the right and left parts and get the slopes.
% Keep the point where the slope difference is greatest.
numPoints = length(x);
index1 = round(0.25 * numPoints); % 25% of the way through.
index2 = round(0.75 * numPoints); % 75% of the way through.
% In other words, assume that we need at least 25 percent of the points to make a good estimate of the line.
% Obviously if we took only 2 or 3 points, then the slope could vary quite dramatically,
% so let's use at least 25% of the points to make sure we don't get crazy slopes.
% Initialize structure array
for k = 1 : numPoints
    lineData(k).slopeDifferences = 0;
    lineData(k).line1 = [0,0];
    lineData(k).line2 = [0,0];
end
for k = index1 : index2
    % Get data in left side.
    x1 = x(1:k);
    y1 = yData(1:k);
    % Fit a line through the left side.
    coefficients1 = polyfit(x1, y1, 1); % The slope is coefficients1(1).
    % Get data in right side.
    x2 = x(k+1:end);
    y2 = yData(k+1:end);
    % Fit a line through the left side.
    coefficients2 = polyfit(x2, y2, 1); % The slope is coefficients2(1).
    
    % Compute difference in slopes, and store in structure array along with line equation coefficients.
    lineData(k).slopeDifferences = abs(coefficients1(1) - coefficients2(1));
    lineData(k).line1 = coefficients1;
    lineData(k).line2 = coefficients2;
end
% Find index for which slope difference is greatest.
slopeDifferences = [lineData.slopeDifferences]; % Extract from structure array into double vector of slope differences only
% slope1s = struct2table(lineData.line1); % Extract from structure array into double vector of slopes only
% slope2s = [lineData.line2(1)]; % Extract from structure array into double vector of slopes only
[maxSlopeDiff, indexOfMaxSlopeDiff] = max(slopeDifferences)
% Plot slope differences.
subplot(2, 2, 2);
plot(slopeDifferences, 'b.', 'MarkerSize', markerSize);
xlabel('Index', 'FontSize', fontSize);
ylabel('Slope', 'FontSize', fontSize);
grid on;
caption = sprintf('Slope Differences Maximum at Index = %d, x value = %.2f', indexOfMaxSlopeDiff, x(indexOfMaxSlopeDiff));
title(caption, 'FontSize', fontSize);
% Mark it with a red line.
line([indexOfMaxSlopeDiff, indexOfMaxSlopeDiff], [0, maxSlopeDiff], 'Color', 'r', 'LineWidth', 2);
% Show everything together all on one plot.
% Plot lines.
subplot(2, 2, 3:4);
plot(x, yData, 'b.', 'MarkerSize', markerSize);
grid on;
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
hold on;
% Use the equation of line1 to get fitted/regressed y1 values.
slope1 = lineData(indexOfMaxSlopeDiff).line1(1);
intercept1 = lineData(indexOfMaxSlopeDiff).line1(2);
y1Fitted = slope1 * x + intercept1;
% Plot line 1 over/through data.
plot(x, y1Fitted, 'r-', 'LineWidth', 2);
% Use the equation of line2 to get fitted/regressed y2 values.
slope2 = lineData(indexOfMaxSlopeDiff).line2(1);
intercept2 = lineData(indexOfMaxSlopeDiff).line2(2);
y2Fitted = slope2 * x + intercept2;
% Plot line 2 over/through data.
plot(x, y2Fitted, 'r-', 'LineWidth', 2);
%================================================================================================================================================
% FIND THE CROSSING POINT.  IT IS WHERE THE Y VALUES ARE EQUAL.
% Just set the equations equal to each other and solve for x.  
% So y1Fitted = y2Fitted, which means (slope1 * xc + intercept1) = (slope2 * xc + intercept2).  Solving for xc gives:
xc = (intercept2 - intercept1) / (slope1 - slope2);
y1c = slope1 * xc + intercept1; % Will be the same as y2c.
y2c = slope2 * xc + intercept2; % Will be the same as y1c.
% Mark crossing with a magenta line.
% Draw a line up from the x axis to the crossing point.
ylim([24, 30]);
yl = ylim(); % Get yrange.
line([xc, xc], [yl(1), y1c], 'Color', 'm', 'LineWidth', 2);
% Draw a line over from the y axis to the crossing point.
line([0, xc], [y1c, y1c], 'Color', 'm', 'LineWidth', 1);
caption = sprintf('Data with left and right lines overlaid.  Lines cross at (x,y) = (%.4f, %.4f)', xc, y1c);
title(caption, 'FontSize', fontSize);
message1 = sprintf('Left  Equation: y = %.3f * x + %.3f', slope1, intercept1);
message2 = sprintf('Right Equation: y = %.3f * x + %.3f', slope2, intercept2);
message = sprintf('%s\n%s', message1, message2);
fprintf('%s\n', message);
text(12, 26, message, 'Color', 'r', 'FontSize', 15, 'FontWeight', 'bold', 'VerticalAlignment', 'middle');
uiwait(helpdlg(message));

Left Equation: y = 0.359 * x + 25.341

Right Equation: y = -0.045 * x + 29.670

Crossing point at (10.7232, 29.1915)

Answer 2

[Matt J]

This is equivalent to a first order free-knot spline fit. This FEX submission might be useful,

https://www.mathworks.com/matlabcentral/fileexchange/25872-free-knot-spline-approximation?s_tid=srchtitle