function [ln, fourPoints] = findLines(E);

%Use matlab's radon function to compute the hough transform:
theta = (0 : 179)';
[R, xp] = radon(E, theta);

peaks = 70;
i = find(R > peaks);
[foo,ind] = sort(-R(i));
k = i(ind(1 : size(i)));  
[y, x] = ind2sub(size(R), k); 
t = - theta(x) * pi / 180;
r = xp(y);
lines = [cos(t) sin(t) -r];  
cx = size(E, 2) / 2 - 1;
cy = size(E, 1) / 2 - 1;
lines(:, 3) = lines(:, 3) - lines(:, 1) * cx - lines(:, 2) * cy;  

t_r = [lines(:, 3), t];
    
for k = 1 : 4
    
    [min_r ind_r] = min(t_r(:, 1));
    t_t = t_r(ind_r, 2);
    
    for i = 1 : length(lines)
        if abs(t_r(i, 2) - t_t) <= 0.2
            t_r(i, :) = Inf;
        end
    end
    
    ln(k, :) = [cos(t_t) sin(t_t) min_r];    
    
end

for i = 1 : 4
    foundLines(i, 1) = - ln(i, 1) / ln(i, 2);
    foundLines(i, 2) = - ln(i, 3) / ln(i, 2);
end

points = zeros(6, 2);
k = 0;

for i = 1 : 4
    for j = i + 1 : 4
        k = k + 1;  
        points(k, 1) = - (foundLines(i, 2) - foundLines(j, 2)) / (foundLines(i, 1) - foundLines(j, 1));
        points(k, 2) = foundLines(i, 1) * points(k, 1) + foundLines(i, 2);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Get the "center of mass" of the points
mX = mean(points(:,1));
mY = mean(points(:,2));

cPoints = zeros(6, 3);
% Translate all points to the centeroid coord system
cPoints(:, 1) = points(:, 1) - mX;
cPoints(:, 2) = points(:, 2) - mY;

% Compute distance of points to their "center of mass"
cPoints(:, 3) = sqrt(points(:, 1) .^ 2 + points(:, 2) .^ 2);

% Get the 4th biggest distance
sorted = sort(cPoints(:, 3));
thres = sorted(4);

fourPoints = zeros(4, 2);

k = 0;
for i = 1 : 6
    
    if cPoints(i, 3) <= thres
        k = k + 1;
        fourPoints(k, :) = points(i, :);
    end
    
end

% Swap centers so that they go clock - wise
bubble = fourPoints(3, :);
fourPoints(3, :) = fourPoints(4, :);
fourPoints(4, :) = bubble;
