%============================================s3p.m=============================================
%Step3：對二值圖像進行霍夫變換
%Mmax和Nmax分別為離散化的分割數，這里設置為1000
Mmax = 500;
Nmax = 500;

H = zeros(Mmax,Nmax);               %初始化霍夫數組
L = size(P);                        %獲得二值圖像的大小

%霍夫變換，對每個值為1的像素點進行遍歷，并累加到霍夫數組中
for i = 3:L(1)-2,
    for j = 3:L(2)-2,
        if P(i,j) == 1,
            for n = 1:Nmax,
                th = pi*(n - 1)/(Nmax - 1)/2;   %得到n表示的θ值
                r = j*cos(th) + i*sin(th);      %得到θ對應的ρ值
                m = (Mmax - 1)*r/3000 + 1;      %得到ρ值對應的離散數
                m = round(m);
                H(m,n) = H(m,n) + 1;            %霍夫數組加一
            end
        end
    end
end

%求出霍夫數組中最大的和第三大的數據
i1=1;
j1=1;
i2=1;
j2=1;
L = size(H);
for i = 1:L(1),
    for j = 1:L(2),
        if H(i,j)>H(i1,j1),
            j3 = j2;
            i3 = i2;
            i2 = i1;
            j2 = j1;
            i1 = i;
            j1 = j;
        end
    end
end

L = size(I);                                    %獲得原圖像的大小
x = 1:L(2)+2;                                   %x的取值范圍是1~L(2)+2

%求出兩個邊緣直線的參數
th1 = pi*(j1 - 1)/(Nmax - 1)/2;
th2 = pi*(j3 - 1)/(Nmax - 1)/2;
r1 = 3000*(i1 -1)/(Mmax - 1);
r2 = 3000*(i3 -1)/(Mmax - 1);
y1 = round((r1 - x*cos(th1))/sin(th1));
y2 = round((r2 - x*cos(th2))/sin(th2));

%為了使直線在圖像中比較明顯，使用紅色，并且厚度取3個像素
for i = 2:L(2)+1,
    I(y1(i)-1,i-1,1) = 255;
    I(y1(i)-1,i-1,2) = 0;
    I(y1(i)-1,i-1,3) = 0;
    I(y1(i),i-1,1) = 255;
    I(y1(i),i-1,2) = 0;
    I(y1(i),i-1,3) = 0;
    I(y1(i)+1,i-1,1) = 255;
    I(y1(i)+1,i-1,2) = 0;
    I(y1(i)+1,i-1,3) = 0;
    I(y2(i)-1,i-1,1) = 255;
    I(y2(i)-1,i-1,2) = 0;
    I(y2(i)-1,i-1,3) = 0;
    I(y2(i),i-1,1) = 255;
    I(y2(i),i-1,2) = 0;
    I(y2(i),i-1,3) = 0;
    I(y2(i)+1,i-1,1) = 255;
    I(y2(i)+1,i-1,2) = 0;
    I(y2(i)+1,i-1,3) = 0;
end
imwrite(I,'pipe_out.jpg','jpg');                %將處理出來的圖像輸出到pipe_out.jpg中
%=========================================s3p.m(end)==========================================