数字图像处理（八）——孔洞填充+全局阈值+自适应阈值

一、孔洞填充

基本思想：
基于形态学算法，膨胀后与上取反的原图

算法实现步骤：
1.首先找出所有孔洞的位置，只需知道洞中的一个点的坐标即可，下面直接以改点代替该洞
2.新建一张全零图，用0表示背景，1表示前景，大小与原图相同
（1）取出一个洞的坐标，在新图中该位置表1
（2）对该新图用一个结构元进行膨胀，然后再与原图的反（孔洞的地方应全为1）求与
（3）如果检测到一次操作完的结果与操作前相同，则结束迭代
（4）取出下一个洞的坐标，返回到（2）操作，直到所有洞都补完
3.此时新图中为1的位置即为需要填补的孔洞，将新图与原图相或即得到孔洞填充完成的图像
注：
边缘检测：图-腐蚀后的图

代码实现：

clc;clear;
A=imread('Experiment4_tak1.png');
BW=im2bw(A);
AC=~BW;
mask=[0,1,0;1,1,1;0,1,0]; %膨胀操作的结构元
SE=strel('arbitrary',mask);
[m,n]=size(BW);
T=zeros(m,n);

%膨胀操作 某位置只要与结构元有交际，则该位置就为集合
%x0找出了所有的孔洞位置
x0={[29,24],[92,20],[187,17],[48,79],[2,91],[132,72],[92,102],[236,70],[213,117],[53,153],[122,155],[197,189],[262,181],[29,230],[119,233],[204,240],[260,250]};
for i=1:length(x0)
flag=1;
T_last=zeros(m,n);
col_temp=x0{i}(1);
row_temp=x0{i}(2);

T(row_temp,col_temp)=1; 	  %令该位置为1 然后利用模板进行膨胀
    while flag

        T=imdilate(T,SE);
        T=T∾               %将T与背景相与，注意这里T和AC用1表示需要填充的空洞即黑色局域 

        if T==T_last          %如果操作前后没有变化则认为填充完成
            flag=0;
        end
        T_last=T;
    end
end
New=BW|T;

edge=New-imerode(New,SE);

subplot(131);
imshow(BW);title('处理前');
subplot(132);
imshow(New);title('处理后');
subplot(133);
imshow(edge);title('边缘');

结果：

二、全局阈值

2.1 迭代阈值

基本思想：
Ti1=(小于Ti的平均灰度+大于Ti的平均灰度)/2，一直迭代到Ti1与Ti相同

注：
matlab默认操作的数据类型是double，所以转int时一定要慎重，一个表达式中有一个int整个表达式就是int型了！

代码实现：

clc;clear;
A=imread('Experiment4_task2.jpg');
A=rgb2gray(A);

H=imhist(A);
h=H./sum(H);
h=h';
T0=120;
error=100;

%循环或用矩阵相乘做
%这里注意表达式中存在int整个式子就默认为int型
while (error)>=1e-5
    temp1=0;
    sum1=sum(h(1:(T0+1)));
    temp2=0;
    sum2=sum(h((T0+2):256));
    for i=1:double(T0+1)
        temp1=temp1+h(i)*(i-1)/sum1;
    end
     
%     x=double(0:T0);
%     temp1=x*h(1:T0+1)'/sum(h(1:T0+1));
    
    for j=double(T0+2):256
        temp2=temp2+h(j)*(j-1)/sum2;
    end
%     y=double(T0+1:255);
%     temp2=y*h(T0+2:256)'/sum(h(T0+2:256));
    
    T1=(1/2*(temp1+temp2));
    error=abs(T1-T0);
    T0=int32(T1);
end

% while error>1
%     temp1=(find(A>T0));
%     temp2=(find(A<T0));
%     T1=1/2*(mean(A(temp1))+mean(A(temp2)));
%     error=abs(T1-T0);
%     T0=T1;
% end

[m,n]=size(A);
A1=A;
for i=1:m
    for j=1:n
        if A1(i,j)>T0
            A1(i,j)=255;
        else
            A1(i,j)=0;
        end
    end
end
str1=['迭代法 阈值为',num2str(T0)];
subplot(121);
imshow(A1);title(str1);

2.2 Otsu方法, 大津法

基本思想：
阈值从0取到255，计算每一阈值对应类间方差，取类间方差取max时的阈值

代码实现：

%otsu方法
for i=1:256
    P1=sum(h(1:i));
    P2=1-P1;
    m1=0;
    m2=0;
    for j=1:256
        if j<=i
            m1=m1+h(j)*(j-1)/P1;
        else
            m2=m2+h(j)*(j-1)/P2;
        end
    end

    op(i)=P1*P2*(m1-m2)^2;
end
t=find(op==max(op))-1;
A2=A;
for i=1:m
    for j=1:n
        if A2(i,j)>t
            A2(i,j)=255;
        else
            A2(i,j)=0;
        end
    end
end
str2=['大津法 阈值为',num2str(t)];
subplot(122);
imshow(A2);title(str2);

结果：

三、自适应阈值

基本思想：
如果图像光照不均匀等因素，则不能使用一个规定的全局阈值，选择的阈值时坐标的函数。

算法实现步骤：
1.将整幅图像分成一系列互相有50%重叠的子图像；
2.得到每个子图像的直方图；
3.检测各子图像是否为双峰，若是则采用最优阈值法，否则不进行处理；
4.根据对直方图为双峰的子图像得到的阈值，通过插值得到所有子图像的插值；
5.根据各子图像的阈值，然后对图像进行分割；
注：
1.这里判断双峰采用了findpeaks()函数，并规定而间距过小的峰值不算，检测到的极大值的个数超过1个则认为是双峰；
2.插值利用的是imresize()函数，直接将得到的一系列阈值模板插值得到与原图等大的阈值矩阵；
3.这里处理的具体方法为：
（1）先不管峰的个数，直接利用大津法直接得到阈值模板后插值得到阈值矩阵T；
（2）对峰的个数进行判断，如果仅有1个的话该处阈值设为0，得到阈值模板后进行插值得到阈值矩阵T1；
（3）最终的阈值模板应为T2=(T1&T).*T

代码实现：

clc;clear;
A=imread('Experiment4_task3.jpg');
A=rgb2gray(A);
[m1,n1]=size(A);
row_num=floor(m1/5);
col_num=floor(n1/5);
for i=0:(row_num-2)
    for j=0:(col_num-2)
        block=A(1+5*i:1+5*i+9,1+5*j:1+5*j+9);
        B=imhist(block);
        pks=findpeaks(B,'minpeakdistance',35);
        if length(pks)>1
            t1=255*graythresh(block);
        else
            t1=0;
        end
        t=255*graythresh(block);
        T(i+1,j+1)=t;
        T1(i+1,j+1)=t1;
    end
end
T_1=imresize(T,[m1,n1],'bilinear');
T1_1=imresize(T1,[m1,n1],'bilinear');
T2=T_1&T1_1;
T2_1=T2.*T_1;
for i=1:m1
    for j=1:n1
        if A(i,j)>T2_1(i,j)
            A(i,j)=255;
        else
            A(i,j)=0;
        end
    end
end

figure;
%进行开操作
mask=[0,1,0;1,1,1;0,1,0]; %膨胀操作的结构元
SE=strel('arbitrary',mask);
A=imerode(A,SE);
A=imdilate(A,SE)
imshow(A);

结果：
效果还行，如果把窗减得更小估计能更好，但就是跑得慢。。。