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

一、孔洞填充

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

算法实现步骤：
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 flagT=imdilate(T,SE);T=T&AC;               %将T与背景相与，注意这里T和AC用1表示需要填充的空洞即黑色局域 if T==T_last          %如果操作前后没有变化则认为填充完成flag=0;endT_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-5temp1=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):256temp2=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:mfor j=1:nif A1(i,j)>T0A1(i,j)=255;elseA1(i,j)=0;endend
end
str1=['迭代法 阈值为',num2str(T0)];
subplot(121);
imshow(A1);title(str1);

2.2 Otsu方法, 大津法

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

代码实现：

%otsu方法
for i=1:256P1=sum(h(1:i));P2=1-P1;m1=0;m2=0;for j=1:256if j<=im1=m1+h(j)*(j-1)/P1;elsem2=m2+h(j)*(j-1)/P2;endendop(i)=P1*P2*(m1-m2)^2;
end
t=find(op==max(op))-1;
A2=A;
for i=1:mfor j=1:nif A2(i,j)>tA2(i,j)=255;elseA2(i,j)=0;endend
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)>1t1=255*graythresh(block);elset1=0;endt=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:m1for j=1:n1if A(i,j)>T2_1(i,j)A(i,j)=255;elseA(i,j)=0;endend
endfigure;
%进行开操作
mask=[0,1,0;1,1,1;0,1,0]; %膨胀操作的结构元
SE=strel('arbitrary',mask);
A=imerode(A,SE);
A=imdilate(A,SE)
imshow(A);

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