function project()
%% load image
% load image
img0 = imread( 'A41.bmp.bmp' );

% convert RGB format to grayscale
img = rgb2gray( img0 );
figure;
subplot(2,2,1), imshow( img );

% crop image
img = img( 5:end-60, 10:end-5 );

% remove background
background = imopen(img, strel('disk',64));
img = imsubtract(img, background);

%% detect lines
%  Create a Binary Version and Skeletonize image
bw = edge(img, 'canny', [0.0375 0.115]);
subplot(2,2,2), imshow( bw );

% Hough transform
[H,theta,rho] = hough(bw);
subplot(2,2,3);
imshow(H, [], 'XData', theta, 'YData', rho, ...
    'InitialMagnification', 'fit');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;

P = houghpeaks(H,200,'threshold',ceil(0.3*max(H(:))));

x = theta(P(:,2));
y = rho(P(:,1));
plot(x,y,'s','color','red');

lines = houghlines(bw, theta, rho, P, ...
    'FillGap', 10, 'MinLength', 38);
nLines = length( lines );

% plot all detected lines
subplot(2,2,4);
imshow(img), hold on
for k = 1:nLines;
    % Plot Line
    xy = [lines(k).point1; lines(k).point2];
    plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
    % Plot beginnings and ends of lines
    plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
    plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
end

%% section II compare lines with each other
% sort and compile lines
figure;
imshow(img), hold on
mainLines = lines;
nMainL = 0;
for i = 1:nLines
    % replot lines in new figure
    xy = [lines(i).point1; lines(i).point2];
    plot(xy(:,1),xy(:,2),'LineWidth',1,'Color','green');
    % Plot beginnings and ends of lines
    plot(xy(1,1),xy(1,2),'x','LineWidth',1,'Color','white');
    plot(xy(2,1),xy(2,2),'x','LineWidth',1,'Color','white');
    
    % similar to whom...
    similarTo = 0;
    distance = inf;
    % try to find similar lines
    % detect line position
    isLeft = left_right(img, lines(i));
    % compare with lines detected till now
    for j = 1:nMainL
        % have a same slope..?
        if abs( lines(i).theta - mainLines(j).theta ) <= 10
            % check two line distance from each other
            theta1 = lines(i).theta * pi/180;
            theta2 = mainLines(j).theta * pi/180;
            theta = (theta1 + theta2) / 2;
            x = ( lines(i).point1(1) + lines(i).point2(1) ) / 2;
            y1 = ( lines(i).rho - x*cos(theta1) ) / sin(theta1);
            y2 = ( mainLines(j).rho - x*cos(theta2) ) / sin(theta2);
            if isLeft
                z = (y1-y2) * cos(pi/2-abs(theta));
            else
                z = (y2-y1) * cos(pi/2-abs(theta));
            end
            z = sign(theta) * z;
            if z<-20 || z>40
                continue;
            elseif  z<15 || xor(isLeft, mainLines(j).isLeft)
                if z < distance
                    similarTo = j;
                    distance = z;
                end
            end
        end
    end
    % new line
    if similarTo == 0
        nMainL = nMainL + 1;
        mainLines(nMainL).theta = lines(i).theta;
        mainLines(nMainL).rho = lines(i).rho;
        mainLines(nMainL).point1 = lines(i).point1;
        mainLines(nMainL).point2 = lines(i).point2;
        mainLines(nMainL).repeat = 1;
        mainLines(nMainL).isLeft = left_right(img, mainLines(nMainL));
        text(lines(i).point1(1), lines(i).point1(2), sprintf('%d,%d', nMainL, i), 'color', 'red');
    else
        % compile two lines into one
        n = mainLines(similarTo).repeat + 1;
        mainLines(similarTo).theta = ( mainLines(similarTo).repeat*mainLines(similarTo).theta ...
            + lines(i).theta ) / n;
        mainLines(similarTo).rho = ( mainLines(similarTo).repeat*mainLines(similarTo).rho ...
            + lines(i).rho ) / n;
        % count compiled lines
        mainLines(similarTo).repeat = n;
        text(lines(i).point1(1), lines(i).point1(2), sprintf('%d,%d', similarTo, i), 'color', 'yellow');
        % move line boundry
        %  check lower x
        mainLines(similarTo).point1(1) = min( [ lines(i).point1(1) mainLines(similarTo).point1(1) ] );
        mainLines(similarTo).point1(2) =  ( mainLines(similarTo).rho - ...
            mainLines(similarTo).point1(1) * cos(  mainLines(similarTo).theta*pi/180 ) ...
            ) / sin(  mainLines(similarTo).theta*pi/180 );
        %  check higher x
        mainLines(similarTo).point2(1) = max( [ lines(i).point2(1) mainLines(similarTo).point2(1) ] );
        mainLines(similarTo).point2(2) =  ( mainLines(similarTo).rho - ...
            mainLines(similarTo).point2(1) * cos(  mainLines(similarTo).theta*pi/180 ) ...
            ) / sin(  mainLines(similarTo).theta*pi/180 );
    end
end

%% results...
% plot detected lines
for k = 1:nMainL;
    % Plot Line
    xy = [mainLines(k).point1; mainLines(k).point2];
    plot(xy(:,1),xy(:,2), '--','LineWidth',1,'Color','blue');
end

% report resultes
clc;
disp('  #     deg.         no.');
disp('------------------------');
for j = 1:nMainL
    theta = sign(mainLines(j).theta)*(90 - abs(mainLines(j).theta));
    disp( sprintf(' %2d %8.2f    %8d', j, theta, mainLines(j).repeat) );
end

%% assisting function
function b = left_right(img, line)
x = ( line.point1(1) + line.point2(1) ) / 2;
y = ( line.rho - x*cos(line.theta*pi/180) ) / sin(line.theta*pi/180);
A = 5 * cos(line.theta*pi/180);
B = 5 * sin(line.theta*pi/180);
b = img(round(y-B), round(x-A)) > img(round(y+B), round(x+A));