[I0,map]=imread('lines','bmp');
colormap(map);
I0=double(I0); %increasing of calculation precision
[row col deep]=size(I0); %acquisition of the input image size

%gray-level image creation
IG0=I0;
%conversion to grayscale
IG0(:,:,1)=(I0(:,:,1)+I0(:,:,2)+I0(:,:,3))/3;
IG0(:,:,2)=IG0(:,:,1);
IG0(:,:,3)=IG0(:,:,1);

%Accumulator
max_ro=floor((sqrt(row*row+col*col))/2);
min_ro=-max_ro;
rozmiar_aku=max_ro-min_ro;
max_alfa=180;
A=zeros(rozmiar_aku,max_alfa);

%Hough transform
for x=1:col
    for y=1:row
        if IG0(row-y+1,x,1)==0
            for ro=min_ro:max_ro
                for alfa=1:max_alfa
                    alf=pi*alfa/180.0;
                    %if ro==round(x.*cos(alf)+y.*sin(alf))
                    if ro==round((x+(-col/2)).*cos(alf)+(y+(-row/2)).*sin(alf))
                        A(max_ro-ro+1,alfa)= A(max_ro-ro+1,alfa)+1;
                    end;
                end;
            end;
        end; %if
    end;
end;


%Accumulator
IA=double(zeros(rozmiar_aku,max_alfa,3));
IA(:,:,1)=A;
minIA=min(min(IA(:,:,1)));
maxIA=max(max(IA(:,:,1)));
IA(:,:,1)=255*(IA(:,:,1)-minIA)/(maxIA-minIA);
IA(:,:,2)=IA(:,:,1);
IA(:,:,3)=IA(:,:,1);


%Lines
IG1=IG0;
for ilosc_prostych=1:5
max_lokalne=max(max(A));
    for x=1:col
        for y=1:row
            if I0(row-y+1,x,1)==0
                for ro=min_ro:max_ro
                    for alfa=1:max_alfa
                    alf=pi*alfa/180.0;
                    if A(max_ro+1-round((x+(-col/2)).*cos(alf)+(y+(-row/2)).*sin(alf)),alfa)== max_lokalne
                        IG1(y,x,1)=(50*ilosc_prostych);
                        IG1(y,x,2)=(10*ilosc_prostych);
                        IG1(y,x,3)=(2*ilosc_prostych);
                    %A(max_ro+1-round((x+(-col/2)).*cos(alf)+(y+(-row/2)).*sin(alf)),alfa)=0;
                    end;
                    end;
                end;
            end; %if
        end;
    end;
end;


% OUTPUT FIGURES

%input grayscale image
IG0=uint8(IG0);
figure(1);
image(IG0);
title('INPUT IMAGE');

%Accumulator
IA=uint8(IA);
figure(2);
image(IA);
title('ACCUMULATOR');

%Accumulator - 3D
figure(3);
surf(A);
title('ACCUMULATOR - 3D');

%input grayscale image with color lines
IG1=uint8(IG1);
figure(4);
image(IG1);
title('DETECTED LINES');