function bestLines = getLines(I,Iedge)

% set 10 degree variance between diagonals, horizontals and verticals
thetaThresh = 10*(pi/180);

Isobel = edge(I,'sobel');

[row, col] = size(Iedge);

[H,theta,rho] = hough(Iedge);

P = houghpeaks(H,10,'threshold',ceil(0.3*max(H(:))));
x = theta(P(:,2)); 
y = rho(P(:,1));

lines = houghlines(Iedge,theta,rho,P,'FillGap',5,'MinLength',7);

thetarho = [lines.theta; lines.rho]';
[finalPolar,m,n] = unique(thetarho,'rows');

[fRow,fCol] = size(finalPolar);
finalLines = struct('start',zeros(fRow,2),'stop',zeros(fRow,2),'slope',zeros(fRow,1),'yint',zeros(fRow,1),'s1',zeros(fRow,1),'s2',zeros(fRow,1),'s3',zeros(fRow,1),'type',zeros(fRow,1));

figure, imshow(I), hold on

for k =1:size(finalPolar,1),
    xy = [lines(m(k)).point1 lines(m(k)).point2];
    tf = ismember(thetarho,[lines(m(k)).theta lines(m(k)).rho],'rows');
    indi = find(tf);
    
    for j = 1:numel(indi)
        p1(j,:) = lines(indi(j)).point1;
        p2(j,:) = lines(indi(j)).point2;
    end
    
    slopes = (p1(:,2)-p2(:,2))./(p1(:,1)-p2(:,1));
    slope = mean(slopes);
    
    yint = xy(2) - slope*xy(1);
    
    lengths = sqrt((p1(:,2)-p2(:,2)).^2 + (p1(:,1)-p2(:,1)).^2);
    s1 = sum(lengths);
    
    allpts = [p1; p2];
    maxX   = max(allpts(:,1));
    minX   = min(allpts(:,1));
    maxY   = max(allpts(:,2));
    minY   = min(allpts(:,2));
    s2     = sqrt( (maxX-minX)^2 + (maxY-minY)^2 );
    
    
    
    if (abs(slope) == Inf),
        start = [xy(1) 1];
        stop = [xy(1) row];
        
        sobelS = Isobel(:,xy(1));
    else
        start = [1 yint];
        stop  = [col (slope.*col + yint)];
        
        x = [minX:maxX]';
        y = x.*slope + yint;
        sobelS = zeros(numel(x),1);
        for c = 1:numel(x),
            sobelS(c) = Isobel(min(max(round(y(c)),1),row),min(max(round(x(c)),1),col));
        end
    end
    
    % Length of all the line segments combined
    finalLines(k).s1 = s1;
    % Proportional to the length of L
    finalLines(k).s2 = s2;
    % Median value of the norm of the gradient
    finalLines(k).s3 = mean(sobelS);
    
    finalLines(k).slope = slope;
    finalLines(k).yint = yint;
    finalLines(k).start = start;
    finalLines(k).stop = stop;
    clear p1 p2 s1 s2 s3;
end
%% Normalize s1, s2, s3

maxS1 = max([finalLines.s1],[],2);
maxS2 = max([finalLines.s2],[],2);
maxS3 = max([finalLines.s3],[],2);
%%

for k = 1:size(finalLines,2),
    finalLines(k).s1    = finalLines(k).s1/maxS1;
    finalLines(k).s2    = finalLines(k).s2/maxS2;
    finalLines(k).s3    = finalLines(k).s3/maxS3;
    
    finalLines(k).salience = (finalLines(k).s1 + finalLines(k).s2 + finalLines(k).s3)/3;     
end
%%

[sorted,sInd] = sortrows([finalLines.salience]',-1);

sortedLines = finalLines(sInd);

%% remove lines around border of image
threshb     = 2;
removeLine  = false(size(sortedLines,2));
for i = 1:size(sortedLines,2),
    if (sortedLines(i).slope == 0) && ( ((sortedLines(i).start(2)>=1) && (sortedLines(i).start(2)<=1+threshb)) || ((sortedLines(i).start(2)<=row) && (sortedLines(i).start(2)>=row-threshb)) ),
        removeLine(i) = true;
    elseif (abs(sortedLines(i).slope) == Inf) && ( ((sortedLines(i).start(1)>=1) && (sortedLines(i).start(1)<=1+threshb)) || ((sortedLines(i).start(2)<=col) && (sortedLines(i).start(2)>=col-threshb)) ),
        removeLine(i) = true;
    end
end

sortedLines(removeLine) = [];

%% Print unoptomized Lines
figure; imshow(I);
numLines = size(sortedLines,2);
for k = 1:numLines,
        coords = [sortedLines(k).start; sortedLines(k).stop];
        hold on
        plot(coords(:,1),coords(:,2),'LineWidth',2,'Color',[sortedLines(k).salience (numLines-(k-1))/numLines 0 ]);
        clear coords;
end

%% Now, make the vertical and horizontal lines the most important ones...

numImp  = 1;
bestLines = sortedLines;

%% put horizons before everything else
hThresh         = abs(sin(0+thetaThresh)/cos(0+thetaThresh));
removeLine      = false(numLines,1);
horizLines      = false(numLines,1);

for k=1:numLines,
    if (abs(sortedLines(k).slope) <= hThresh),
       % this is a horizon line, super important
       horizLines(numImp)       = true; 
       bestLines(numImp)        = sortedLines(k);
       removeLine(k)            = true;
       bestLines(numImp).type   = 1;
       numImp = numImp + 1;
    end
end

%% Normalize the Salience of the horizontal lines

maxHorizSalience = max([bestLines(horizLines).salience]);
horizInd = find(horizLines);

for k = 1:numel(horizInd),
    bestLines(horizInd(k)).salience = bestLines(horizInd(k)).salience/maxHorizSalience;
end

%% put verticals after horizons

vThresh     = abs(sin(90-thetaThresh)/cos(90-thetaThresh));
vertLines   = false(numLines,1);

for k=1:numLines,
    if (abs(sortedLines(k).slope) >= vThresh),
       % this is a vertical size, slightly less super important
       vertLines(numImp)        = true; 
       bestLines(numImp)        = sortedLines(k);
       removeLine(k)            = true;
       bestLines(numImp).type   = 2;
       numImp = numImp + 1;
    end
end

sortedLines(removeLine) = [];

%% Normalize vertical lines

maxVertSalience = max([bestLines(vertLines).salience]);
vertInd = find(vertLines);

for k = 1:numel(vertInd),
    bestLines(vertInd(k)).salience = .75*(bestLines(vertInd(k)).salience/maxVertSalience);
end

%% Designate Diagonals
op = (row-1);
ad = (col-1);
hy = sqrt(op^2+ad^2);

% sin(THETA) = op/hy;
diagTheta   = asin(op/hy);

smSlope     = sin(diagTheta-thetaThresh)/cos(diagTheta-thetaThresh);
lgSlope     = sin(diagTheta+thetaThresh)/cos(diagTheta+thetaThresh);

%% Now, the rest of the lines
[sortS,locS] = sort([sortedLines.salience]','descend');
maxSalience = max([sortedLines.salience]);
for k = 1:numel(locS)
    bestLines(numImp) = sortedLines(locS(k));
    bestLines(numImp).salience = .5*(bestLines(numImp).salience/maxSalience);   
    if (abs(bestLines(numImp).slope) <= lgSlope)&&(abs(bestLines(numImp).slope) >= smSlope),
        bestLines(numImp).type = 3;
    else
        bestLines(numImp).type = 0;
    end
    numImp = numImp +1;
end

%% Print Best Lines
figure; imshow(I);
numLines = size(bestLines,2);
for k = 1:4,
        coords = [bestLines(k).start; bestLines(k).stop];
        hold on
        plot(coords(:,1),coords(:,2),'LineWidth',2,'Color',[bestLines(k).salience (numLines-(k-1))/numLines 0 ]);
        clear coords;
end
