contrast_line = diff(image);
imshow(contrast_line); % shows the line
thresholded_image=mean(contrast_line,3) > 240; % set here a value between 0 and 255 that would fit with your contrast line. It could be found perhaps with max(max(max(contrast_line)))
imshow(thresholded_image) % show the thresholded difference, in case you need to extract only the points above a certain threshold (or below, in case you set a < in the line above)
%% the other lines below require this function (see link here below)
%% based on this post (https://de.mathworks.com/matlabcentral/answers/512870-how-to-perform-robust-line-fitting-in-a-binary-image)
[y,x]=find(flipud(thresholded_image));
p_best=fit_line(x,y);
p_poly=polyfit(x,y,1);%polyfit for reference
x_bounds=[1 size(thresholded_image,2)];
figure(1),clf(1)
subplot(1,2,1)
imshow(thresholded_image)
subplot(1,2,2)
plot(x,y,'.','DisplayName','data'),hold on
plot(x_bounds,polyval(p_best ,x_bounds),'k--',...
'DisplayName',sprintf('R^2=%.2f (fit_line)',get_r_square(p_best ,x,y)))
plot(x_bounds,polyval(p_poly ,x_bounds),'r--',...
'DisplayName',sprintf('R^2=%.2f (polyfit)',get_r_square(p_poly ,x,y)))
hold off
axis equal
axis([1 size(thresholded_image,2) 1 size(thresholded_image,1)])
legend('Location','northeastoutside')
function p=fit_line(x,y)
%Fit a line to the data. Use the RMS of the orthogonal distance as a cost
%function instead of MSE_y, as polyfit probably does.
%
%This works with fminsearch, which is sensitive to initial values that are
%far from the optimum, sometimes returning local optima.
x=x(:);y=y(:);
pt=[x y];
%root mean square of the distance between the line defined by the two
%points in v and the points defined by x and y.
RMS=@(x) sqrt(mean(x.^2));
cost_fun=@(v) RMS(point_to_line_distance(pt, v([1 3]), v([2 4])));
%convert [x1 y1 phi2] to [x1 y1 x2 y2]
vr_2_v=@(vr) [vr(1:2) vr(1:2)+[cos(vr(3)) sin(vr(3))]];
%wrap the cost function and the converter
fun=@(vr) cost_fun(vr_2_v(vr));
%initialize to around the center of the data
init=[mean(x) mean(y) 0*pi];
%execute fit
opts = optimset('MaxFunEvals',50000, 'MaxIter',10000);
fit_val = fminsearch(fun, init, opts);
tmp=vr_2_v(fit_val);
p=polyfit(tmp(1:2),tmp(3:4),1);
end
function r2=get_r_square(p,x,y_real)
y_fit=polyval(p,x);
err=y_real-y_fit;
SSres=sum(err.^2);
SStot=sum((y_real-mean(y_real)).^2);
r2=1-(SSres./SStot);
end
