全局阈值处理的改进

前言

在之前的学习中，我们采用的全局阈值处理的原图像往往本身是十分理想的。然而在实际的操作过程中几乎不可能存在恰好能使用全局阈值处理的原图，例如一幅图像中总会存在一些噪声，而噪声的存在会对图像的全局阈值处理带来毁灭性的打击。因此如何实现全局阈值处理的改进就是我们本次学习的重点内容。

 

一、使用图像平滑改进全局阈值处理

正如前言中所说，我们拿到的一手的图像其中总会存在着许多干扰噪声，而噪声的存在会把原本简单的阈值处理问题变得完全无法解决。当我们不能在源头降低噪声时，如果要进行阈值处理则必须先对噪声进行滤除，也就是说要先进行图像平滑来改进原图的质量，再进行全局阈值处理。

图像平滑的常用方法在空间域滤波的学习过程中我已经详细的记录过，在这里不再赘述。本质上就是先用滤波器滤除噪声，再进行图像全局阈值处理。

二、代码示例

给定一幅人为加入均值为0，方差为0.038的高斯噪声的图像，先进行全局阈值处理观察图像情况；再利用5*5的均值滤波器滤波后，再进行全局阈值处理观察图像情况。最后进行二者的对比。

I1=imread('septagon.tif');
I2=imnoise(I1,'gaussian',0,0.038);%加入噪声
T1=graythresh(I2);%对图像进行Otsu的全局阈值处理
I3=im2bw(I2,T1);
w=fspecial('average',5);%引入均值滤波器模板
I4=imfilter(I2,w,'replicate');%进行图像平滑滤除噪声
T2=graythresh(I4);%再进行Otsu的全局阈值处理
I5=im2bw(I4,T2);
subplot(231),imshow(I2),title('噪声污染图');
subplot(232),imhist(I2),title('噪声污染图直方图');
subplot(233),imshow(I3),title('基于Otsu方法的阈值处理图像');
subplot(234),imshow(I4),title('滤波后图像');
subplot(235),imhist(I4),title('滤波后图像直方图');
subplot(236),imshow(I5),title('基于Otsu方法的阈值处理图像');

三、结果展示及总结

加入噪声后我们可以通过观察图像的直方图发现该图像的直方图分布并不具有适合图像分割的明显特征，对这样的图像直接进行阈值处理往往会失败，事实也确实如此，在233图像中，我们采用了Otsu方法进行全局阈值处理，结果得到了许多亮暗点（也就是阈值处理过程中产生的误差），因此这次阈值处理是完全不成功的。

采用图像平滑后，观察235的直方图可知，针对这种直方图的阈值处理一定可以得到非常完美的分割图像，结果也正如236展示的那样，图像的分割非常完美。

四、使用拉普拉斯算子边缘信息改进全局阈值处理

因为在图像分割的时候最重要的是判断该像素点是否在图像的边缘上，而判断其是否在边缘上最常用的方法就是通过计算其梯度或拉普拉斯算子的绝对值来获得的（拉普拉斯算子有正有负），因此我们使用拉普拉斯算子预先对图像进行处理，得到一幅具有强边缘的图像，再利用Otsu方法进行图像分割往往可以得到非常好的图像分割效果。

具体步骤为：

1.使用拉普拉斯算子由f计算一幅边缘图像。边缘图像我们设为拉普拉斯算子的绝对值。

2.指定一个阈值

3.使用来自步骤2的阈值对来自步骤1的图像进行阈值处理，产生二值图像g。这幅图像在步骤4中用作一幅标记图像，以便从f中选取对应于强边缘像素的像素。

4.仅使用f中对应于g中1值像素的位置的像素计算直方图

5.使用来自步骤4的直方图采用Otsu方法进行全局分割

通常在指定阈值时我们会使用自定义函数percentile2i，这是一个可以给定对应百分位灰度值的函数，例如：

Q=percentile2i(f,0.999)指的就是给定输入图像f中百分位为99.9%的灰度值，赋值给Q。

对于该自定义函数，其代码如下：

function I=percentile2i(h,P)
%PERCENTILE2I Computes an intensity value given a percentile.
%I=PERCENTILE2I(H,P) Given a percentile,P,and a histogram,
%H,this function computes an intensity,I,representing the
%Pth percentile and returns the value in I.P must be in the
%range [0, 1] and I is returned as a value in the range [0, 1]
%also.

%Check value of P.
if P<0||P>1
erro('The percentile must be in the range [0, 1]');
end
%Normalized the histogram to unit area.If it is already normalized
%the following computation has no effect.
h=h/sum(h);
%Cummulative distribution.
C=cumsum(h);
%Calculations.
idx=find(C>=P,1,'first');
%Subtract 1 from idx because indexing starts at 1,but intensities
%start at 0.Also,normalize to the range [0, 1].
I=(idx-1)/(numel(h)-1);

五、代码示例

现在我们对一幅酵母细胞的图像进行全局阈值处理。我们会对比直接使用Otsu方法和先采用拉普拉斯算子处理再使用Otsu方法得到的两幅图像来进行对比。

I1=imread('cell.tif');
I1=tofloat(I1);%改变图像类型
hf=imhist(I1);
[Tf SMf]=graythresh(I1);%先对图像直接使用Otsu方法
I2=im2bw(I1,Tf);
w=[-1 -1 -1;-1 8 -1;-1 -1 -1];%引入拉普拉斯模板
I3=abs(imfilter(I1,w,'replicate'));%对图像用拉普拉斯模板进行处理
I3=I3/max(I3(:));%将I3的值包含在[0,1]内
h=imhist(I3);
Q=percentile2i(h,0.995);%取h中0.995为阈值，对应于灰度值的99.5%位，以便边缘图像中仅有很少的像素用于阈值的计算
markerimage=I3>Q;%取标记图像为大于阈值的部分
I4=I1.*markerimage;%做相乘图像得到强边缘像素图像
hp=imhist(I4);
hp(1)=0;%去除0对直方图的贡献，把hp第一个元素排除在外
T=otsuthresh(hp);%利用Otsu方法再次进行图像分割
I5=im2bw(I1,T);
subplot(231),imshow(I1),title('原图');
subplot(232),imhist(I1),title('原图直方图');
subplot(233),imshow(I2),title('使用Otsu方法直接进行分割');
subplot(234),imshow(I4),title('图像乘积');
subplot(235),bar(hp),title('棒状图');
subplot(236),imshow(I5),title('处理后图像');

六、结果展示及分析

根据图像展示我们可以看到，直接使用Otsu方法确实可以提取出酵母细胞的轮廓，但是许多细胞都互相粘连在一起而没能完全分割出来，另外我们可以发现在原图中酵母细胞的细胞核明显比细胞质更加明亮，如果仅仅提取细胞核则显然直接使用Otsu方法无法做到，它只能大致实现细胞质的提取。而采用拉普拉斯算子进行强处理之后，我们可以看到图像乘积中的亮点基本已经聚集在细胞核周围（这是因为我们设置了阈值为99.5%的缘故，如果仅仅想提取细胞质则应该适当降低阈值），利用该相乘图像的直方图再次进行Otsu方法的全局阈值处理，最终我们得到了酵母细胞的细胞核的提取图。提取的非常成功！