阈值法进行边缘检测matlab

1、基于全局阈值处理

1）选择初始阈值T

2）用T分割图像，分别是大于T的像素集合和小于T的像素集合。

3）计算两集合的平均灰度m1，m2

4）计算新阈值T=（m1+m2）/2

5）重复上述步骤，直至两次迭代的T的变化小于一定值

6）利用函数im2bw分割图像

例：

clc
clear
f=imread('E:\桌面\数字图像matlab\DIP3E_CH10_Original_Images\DIP3E_Original_Images_CH10\Fig1038(a)(noisy_fingerprint).tif');

count=0;
T=mean2(f);
done=false;
while ~done
    count=count+1;
    g=f>T;
    Tnext=0.5*(mean(f(g))+mean(f(~g)));
    done=abs(T-Tnext)<0.5;
    T=Tnext;
end
g=im2bw(f,T/255);
subplot(121),imshow(f)

subplot(122),imshow(g)
figure,imhist(f)

结果如下：

count=2，迭代两次即收敛。

直方图有两个谷峰，很明显可以分开，程序得到T=125，由图可以看出来是可行的。

2、otsu方法的最佳全局阈值处理

阈值k将灰度分为两组。

首先定义可分性度量，即类间方差与图像灰度方差的比。越大说明可分性越好。

工具箱函数graythresh，[T,SM]=graythresh(f)

例：

clc
clear
f=imread('E:\桌面\数字图像matlab\DIP3E_CH10_Original_Images\DIP3E_Original_Images_CH10\Fig1039(a)(polymersomes).tif');

count=0;
T1=mean2(f);
done=false;
while ~done
    count=count+1;
    g=f>T1;
    Tnext=0.5*(mean(f(g))+mean(f(~g)));
    done=abs(T1-Tnext)<0.5;
    T1=Tnext;
end

g=im2bw(f,T1/255);
subplot(121),imshow(g)

[T2,SM]=graythresh(f);
g1=im2bw(f,T2);
subplot(122),imshow(g1)

 结果如下：

第一张图是利用第一种方法得到的，第二个是otsu方法得到了更好的分割效果。

 

 还有一种基于图像直方图给出阈值的函数otsuthresh。

3、使用图像平滑改进全局阈值

当图像含有噪声时，进行图像平滑（滤波）再进行分割会更好。

例：

clc
clear
f=imread('E:\桌面\数字图像matlab\DIP3E_CH10_Original_Images\DIP3E_Original_Images_CH10\Fig1036(a)(original_septagon).tif');

fn=imnoise(f,'gaussian',0,0.038);%噪声图
subplot(231),imshow(fn)

subplot(232),imhist(fn)
Tn=graythresh(fn);
gn=im2bw(fn,Tn);
subplot(233),imshow(gn)%直接进行otsu

w=fspecial('average',5);
fa=imfilter(fn,w,'replicate');%平均值滤波
subplot(234),imshow(fa)
subplot(235),imhist(fa)
Ta=graythresh(fa);
ga=im2bw(fa,Ta);
subplot(236),imshow(ga)

结果如下：

进行滤波之后，图像可分性变得更强了。

 4、使用边缘改进全局阈值

例：使用基于梯度的边缘信息改进全局阈值

（通过选择梯度较大的像素点，也就是边缘，然后在进行全局阈值处理）

clc
clear

f=imread('E:\桌面\数字图像matlab\DIP3E_CH10_Original_Images\DIP3E_Original_Images_CH10\Fig1041(a)(septagon_small_noisy_mean_0_stdv_10).tif');
subplot(231),imshow(f)
f=tofloat(f);
sx=fspecial('sobel');
sy=sx';
gx=imfilter(f,sx,'replicate');
gy=imfilter(f,sy,'replicate');
grad=sqrt(gx.*gx+gy.*gy);
grad=grad/max(grad(:));

h=imhist(grad);
subplot(232),imhist(f)
Q=percentile2i(h,0.999);%生成阈值，用来找出较大的99.9%的梯度值

markerImage=grad>Q;%梯度选择后的
subplot(233),imshow(markerImage);
fp=f.*markerImage;
subplot(234),imshow(fp)
hp=imhist(fp);

hp(1)=0;%因为大部分是0值，所以把0的统计数去掉
subplot(235),bar(hp)
T=otsuthresh(hp);%以直方图作为输入
T*(numel(hp)-1)

g=im2bw(f,T);%大于T的为1，小于为0，得到二值图像
subplot(236),imshow(g)

结果如下：

       进行梯度选择后的图像，可分性变强。

       这个例子中中间白色物体的大小比背景小的多，对直方图的贡献很小，所以如果运用上面的平均值滤波，得不到想要的结果，只能考虑边缘信息。

例：使用拉普拉斯边缘信息改进全局阈值处理

得到原图酵母细胞的亮点区域，图三是直接用于graythresh函数

clc
clear

f=imread('E:\桌面\数字图像matlab\DIP3E_CH10_Original_Images\DIP3E_Original_Images_CH10\Fig1043(a)(yeast_USC).tif');
subplot(231),imshow(f),title('原图')
f=tofloat(f);
subplot(232),imhist(f)
[Tf,SMf]=graythresh(f);
gf=im2bw(f,Tf);
subplot(233),imshow(gf),title('直接otsu')%直接otsu

w=[-1 -1 -1;-1 8 -1;-1 -1 -1];
lap=abs(imfilter(f,w,'replicate'));
lap=lap/max(lap(:));
h=imhist(lap);
Q=percentile2i(h,0.995);
markerImage=lap>Q;
fp=f.*markerImage;
subplot(234),imshow(fp);
hp=imhist(fp);
hp(1)=0;
subplot(235),bar(hp)
T=otsuthresh(hp);%原理也是让类间方差最大
g=im2bw(f,T);
subplot(236),imshow(g)

 其实上面的基于梯度的改进与拉普拉斯类似，也能得到很好的结果。

5、使用移动平均的图像阈值处理

对图像像素进行‘己’字扫描，书上称为zigzag模式。在每个扫描点处计算平均灰度

例：

对一副由斑点灰度模式遮蔽了的手写文本图像，进行分割

clc
clear

f=imread('E:\桌面\数字图像matlab\DIP3E_CH10_Original_Images\DIP3E_Original_Images_CH10\Fig1049(a)(spot_shaded_text_image).tif');

subplot(131),imshow(f)
T=graythresh(f);
g1=im2bw(f,T);
subplot(132),imshow(g1)%直接应用otsu

g2=movingthresh(f,20,0.5);%移动平均
subplot(133),imshow(g2)

function g = movingthresh(f, n, K)
%MOVINGTHRESH Image segmentation using a moving average threshold.
% G = MOVINGTHRESH(F, n, K) segments image F by thresholding its
% intensities based on the moving average of the intensities along
% individual rows of the image. The average at pixel k is formed
% by averaging the intensities of that pixel and its n − 1
% preceding neighbors. To reduce shading bias, the scanning is
% done in a zig-zag manner, treating the pixels as if they were a
% 1-D, continuous stream. If the value of the image at a point
% exceeds K percent of the value of the running average at that
% point, a 1 is output in that location in G. Otherwise a 0 is
% output. At the end of the procedure, G is thus the thresholded
% (segmented) image. K must be a scalar in the range [0, 1].
% Preliminaries.
f = tofloat(f);
[M, N] = size(f);
if (n < 1) || (rem(n, 1) ~= 0)
error('n must be an integer >= 1.')
end
if K < 0 || K > 1
error('K must be a fraction in the range [0, 1].')
end
% Flip every other row of f to produce the equivalent of a zig-zag
% scanning pattern. Convert image to a vector.
%按己字扫描
f(2:2:end, :) = fliplr(f(2:2:end, :));%列翻转，假设有三列第一列与第三列互换
f = f'; 	% Still a matrix.
f = f(:)'; % Convert to row vector for use in function filter.
% Compute the moving average.
maf = ones(1, n)/n; 	% The 1-D moving average filter.
ma = filter(maf, 1, f); % Computation of moving average.
%https://blog.csdn.net/qq_38559814/article/details/86521602?utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromMachineLearnPai2%7Edefault-4.control&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7EBlogCommendFromMachineLearnPai2%7Edefault-4.control
% Perform thresholding.
g = f > K * ma;
% Go back to image format (indexed subscripts).
g = reshape(g, N, M)';
% Flip alternate rows back.
g(2:2:end, :) = fliplr(g(2:2:end, :));

图1是原图，图二是直接用otsu全局阈值处理，第三个是移动平均