function [x1,y1,x2,y2,weights1] = points_2d_em(pts);

% Initialization 

minx = min(pts(1,:));
miny = min(pts(2,:));
maxx = max(pts(1,:));
maxy = max(pts(2,:));

x1 = minx + (maxx-minx)*rand;
y1 = miny + (maxy-miny)*rand;
x2 = minx + (maxx-minx)*rand;
y2 = miny + (maxy-miny)*rand;

sigma_cur = (maxx-minx)/4;

dists1 = sqrt((x1-pts(1,:)).^2 + (y1-pts(2,:)).^2);
dists2 = sqrt((x2-pts(1,:)).^2 + (y2-pts(2,:)).^2);

for i = 1:5
    
    % Expectation step
    
    prob1 = exp(- (dists1.^2 / sigma_cur^2));
    prob2 = exp(- (dists2.^2 / sigma_cur^2));
    weights1 = prob1./(prob1+prob2);
    weights2 = prob2./(prob1+prob2);
    
   if mod(i,1) == 0
    figure(i)
    plot(pts(1,:), pts(2,:), 'ko');
    hold on
    plot([x1,x2], [y1, y2], 'ro');
    hold off
    end
    
    % Maximization step

    x1 = sum(pts(1,:).*weights1)/sum(weights1);
    y1 = sum(pts(2,:).*weights1)/sum(weights1);  
    x2 = sum(pts(1,:).*weights2)/sum(weights2);
    y2 = sum(pts(2,:).*weights2)/sum(weights2);      
    
    dists1 = sqrt((x1-pts(1,:)).^2 + (y1-pts(2,:)).^2);
    dists2 = sqrt((x2-pts(1,:)).^2 + (y2-pts(2,:)).^2);
    
     sigma_cur = max(eps, sqrt( (sum((dists1.^2).*weights1) ...
         + sum((dists2.^2).*weights2)) / (length(pts)-1)));
end