%一個小波邊緣檢測的MATLAB程序
%load woman;
X=imread('Miss.bmp'); 
clear coef1;
clear coef2;
clear coef3;
clear form;
x=X;
x=double(x);%200 line,300 column; 
line=512;%250;
column=512;%250;
J=8;
Tu=7;
smooth_area_num=20;
edge_area_num=40;
for i=1:line
    coef1(i,:)=cwt(x(i,:),1,'bior4.4');
    coef2(i,:)=cwt(x(i,:),2,'bior4.4');
    coef3(i,:)=cwt(x(i,:),J,'bior4.4');
end
%正值1，負值0的矩陣temp表示(突變點左位置）
max=0;
max_loc=0;
max_num=13;%15;
for i=1:line
    j=1;
    while j<column-1
        
        if (coef3(i,j)>0 & coef3(i,j+1)<0) | (coef3(i,j)<0 & coef3(i,j+1)>0 )%找到左位置j
            for k=j+1:column-1
                
                if coef3(i,k)*coef3(i,j)>=0%找到右位置k
                    for m=j:k
                        if abs(coef3(i,m))>max
                            max=abs(coef3(i,m));
                            max_loc=m;
                        end
                    end
                    if max>max_num
                        form(i,max_loc)=1;
                    end
                    j=k;
                    max=0;
                    break;
                end
            end%for k
            if k==column-1
                j=column-1;
            end
        else
            j=j+1;
        end%if
    end%while
end
%for i=1:line
%    for j=1:column
%        if abs(coef1(i,j))>=2*abs(coef2(i,j)) %是否加絕對值？
%           form(i,j)=0;%第一種方法去噪聲
%        end
%    end
%end
%對列作小波變換
x=x';
for i=1:line
    coef1(i,:)=cwt(x(i,:),1,'bior4.4');
    coef2(i,:)=cwt(x(i,:),2,'bior4.4');
    coef3(i,:)=cwt(x(i,:),J,'bior4.4');
end
%正值1，負值0的矩陣temp表示(突變點左位置）
max=0;
max_loc=0;
for i=1:line
    j=1;
    while j<column-1
        
        if (coef3(i,j)>0 & coef3(i,j+1)<0) | (coef3(i,j)<0 & coef3(i,j+1)>0 )%找到左位置j
            for k=j+1:column-1
                
                if coef3(i,k)*coef3(i,j)>=0%找到右位置k
                    for m=j:k
                        if abs(coef3(i,m))>max
                            max=abs(coef3(i,m));
                            max_loc=m;
                        end
                    end
                    if max>max_num
                        form(max_loc,i)=1;
                    end
                    j=k;
                    max=0;
                    break;
                end
            end%for k
            if k==column-1
                j=column-1;
            end
        else
            j=j+1;
        end%if
    end%while
end
%包括噪聲在內的所有極值點①
for i=1:line
    for j=1:column
        if abs(coef1(i,j))>=2*abs(coef2(i,j)) %是否加絕對值？!
            form(i,j)=0;%第一種方法去噪聲
        end
    end
end
%此處添加用第二種方法去噪!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!② 結束
for i=2:line-2
    for j=2:column-2
        %鏈長為1
        if form(i,j-1)==0 & form(i,j+1)==0 & form(i-1,j-1)==0 & form(i-1,j)==0 & form(i-1,j+1)==0 & form(i+1,j-1)==0 & form(i+1,j)==0 & form(i+1,j+1)==0
            form(i,j)=0;%去除孤立點
        end
        
        %****************************************
        %鏈長為2
        %0000
        %0110
        %0000
        if form(i,j)==1 & form(i,j+1)==1 & sum(sum(form(i-1:i+1,j-1:j+2)))==2
            form(i,j)=0;
            form(i,j+1)=0;%行去除兩點
        end
        %000
        %010
        %010
        %000
        if form(i,j)==1 & form(i+1,j)==1 & sum(sum(form(i-1:i+2,j-1:j+1)))==2
            form(i,j)=0;
            form(i+1,j)=0;%列去除兩點 
        end
        %0000
        %0100
        %0010
        %0000
        if form(i,j)==1 & form(i+1,j+1)==1 & sum(sum(form(i:i+1,j:j+1)))==2
            if sum(sum(form(i-1:i+2,j-1:j+2)))==3 & ( form(i-1,j+2)==1 | form(i+2,j-1)==1)
                form(i,j)=0;
                form(i+1,j+1)=0;
            end
            if sum(sum(form(i-1:i+2,j-1:j+2)))==4 &  form(i-1,j+2)==1 & form(i+2,j-1)==1)
                form(i,j)=0;
                form(i+1,j+1)=0;
            end
        end
        %0000
        %0010
        %0100
        %0000
        if form(i,j+1)==1 & form(i+1,j)==1 & sum(sum(form(i:i+1,j:j+1)))==2
            if sum(sum(form(i-1:i+2,j-1:j+2)))==3 & ( form(i-1,j-1)==1 | form(i+2,j+2)==1)
                form(i,j+1)=0;
                form(i+1,j)=0;
            end
            if sum(sum(form(i-1:i+2,j-1:j+2)))==4 &  form(i-1,j-1)==1 & form(i+2,j+2)==1
                form(i,j+1)=0;
                form(i+1,j)=0;
            end
        end
    end
end
%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!!
max=0;






for i=3:line-3
    for j=3:column-3
        s=sum(sum(form(i-2:i+2,j-2:j+2)));
        for m=i-2:i+2
            for n=j-2:j+2
                if s>=7 & abs(coef3(m,n))<=s*3
                    form(m,n)=0;
                end
            end
        end
    end
end
end


%顯示form即為邊緣檢測圖像