A tentative answer for stackoverflow.com/questions/33423107 written in Octave.

so_33423107.m

# https://stackoverflow.com/questions/33423107/how-to-find-2-lines-given-only-points

clear;
pkg load image
img = rgb2gray(imread("VGUM6.png"));
img = img < 128;
rp = regionprops(bwlabel(img));

X = zeros(1,size(rp)(1));
Y = zeros(1,size(rp)(1));
for ( i =1:size(rp)(1))
  X(1,i)=rp(i).Centroid(1);
  Y(1,i)=rp(i).Centroid(2);
 endfor
 

%%% THE EM ALGORITHM
NumModel = 2;
NumData = 64;
Sigma = 0.1; % IN E-STEP
Noise = 0.1; % IN DATA
NIterate = 10;
%%%% MAKE SYNTHETIC DATA
%A = 2*rand(1,NumModel) - 1;
%B = 2*rand(1,NumModel) - 1;
%X = [];
%Y = [];
%for i = 1 : NumModel
%x = 2*rand(1,NumData) - 1;
%y = A(i) * x + B(i) + Noise*(rand(1,NumData)-0.5);
%X = [X x];
%Y = [Y y];
%end

size(X)
size(Y)

%%% INITIALIZE MODEL
a = 2*rand(1,NumModel) - 1;
b = 2*rand(1,NumModel) - 1;
for j = 1 : NIterate
%%% E-STEP
for i = 1 : NumModel
r(i,:) = a(i)*X + b(i) - Y;
end
den = sum( exp( -(r.^2)/Sigma^2 ) );
for i = 1 : NumModel
w(i,:) = exp( (-(r(i,:).^2)/Sigma^2) ) ./ (den+eps);
end
%%% M-STEP
for i = 1 : NumModel
M = [ sum(w(i,:).*X.^2) sum(w(i,:).*X) ; sum(w(i,:).*X) sum(w(i,:)) ];
V = [ sum(w(i,:).*X.*Y) ; sum(w(i,:).*Y) ];
U = pinv(M) * V;
a(i) = U(1);
b(i) = U(2);
end
%%% SHOW MODEL FIT AND ASSIGNMENT
xc = [-1 : 0.1 : 1];
subplot(2,1,1); cla; hold on;
plot( X, Y, 'bo' );
for i = 1 : NumModel
yc = a(i)*xc + b(i);
plot( xc, yc, 'r' );
end
hold off;
subplot(2,1,2); cla; imagesc(w); colormap gray; pause(1);
end