They key is in how you define what a good fit is. The code below uses the root-mean-square of the orthogonal distance between all points and the line as the cost function. The code below needs my point_to_line_distance function, which you can find here on the FEX.
im=imread('image.png');
im=mean(im,3)>240;
[y,x]=find(flipud(im));
p_best=fit_line(x,y);
p_poly=polyfit(x,y,1);%polyfit for reference
x_bounds=[1 size(im,2)];
figure(1),clf(1)
subplot(1,2,1)
imshow(im)
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(im,2) 1 size(im,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
