blood = imread('blood1.BMP');
[x,y]=size(blood);                % 求出圖象大小
b=double(blood);                  
N =sqrt(100) * randn(x,y);        %  生成方差為10的白噪聲
I=b+N;                            %  噪聲干擾圖象

for i=1:x                         % 實際圖象的灰度為0～255
    for j=1:y
        if (I(i,j)>255)
            I(i,j)=255;
        end
        if (I(i,j)<0)
            I(i,j)=0;
        end
    end
end    


z0=max(max(I));                   % 求出圖象中最大的灰度
z1=min(min(I));                   % 最小的灰度 

T=(z0+z1)/2;                      
TT=0;
S0=0; n0=0;
S1=0; n1=0;
allow=0.5;                       % 新舊閾值的允許接近程度
d=abs(T-TT);
count=0;                         % 記錄幾次循環

while(d>=allow)                 % 迭代最佳閾值分割算法
    count=count+1;
    for i=1:x
        for j=1:y
            if (I(i,j)>=T)
                S0=S0+I(i,j);
                n0=n0+1;
            end
            if (I(i,j)<T)
                S1=S1+I(i,j);
                n1=n1+1;
            end
        end
    end 
    T0=S0/n0;
    T1=S1/n1;
    TT=(T0+T1)/2;
    d=abs(T-TT);
    T=TT;
end

Seg=zeros(x,y);
for i=1:x
    for j=1:y
        if(I(i,j)>=T)
            Seg(i,j)=1;               % 閾值分割的圖象
        end
    end
end

SI=1-Seg;                            % 閾值分割后的圖象求反，便于用腐蝕算法求邊緣
se1=strel('square',3);               % 定義腐蝕算法的結構
SI1=imerode(SI,se1);                 % 腐蝕算法
BW=SI-SI1;                           % 邊緣檢測
%=====傳統的邊緣檢測方法======%
I=uint8(I);
BW1=edge(I,'sobel', 0.11);
BW2=edge(I,'log', 0.015);
%===========圖象顯示==========%
figure(1);
imshow(I);title('Original')                         % 顯示閾值分割的圖象
figure(2);
imshow(BW2);title('Soble')                         % 顯示新算法的邊緣圖象
figure(3)
imshow(BW1);title('Gauss-Laplace') 
figure(4)
imshow(BW);title('New algorithm')