图像分割——采用梯度算法改善全局阈值设定

clc;
clear all;
close all;
%用边缘信息改善全局阈值设定
%Using Edges to Improve Global Thresholding

%%
%====================采用梯度算法改善全局阈值设定==========================
I=imread('D:\Gray Files\10-41.tif');
[M,N]=size(I);
%是否统计灰度值为0的像素点
g_gradient=SobelFilter(im2double(I));
% 梯度阈值
th=0.997;
g_T=zeros(M,N);
g_T=uint8(g_T);
ind=find(g_gradient>th);
g_T(ind)=1;
% imshow(g_T)
%将梯度二值图像与原图相乘
g_edge=I.*g_T;
% imshow(g_edge)
%通过Otsu算法寻找门限
L_start=1;
T=Otsu_FindThreshold(g_edge,L_start);
g=zeros(M,N);
ind=find(I>T);
g(ind)=1;
g=uint8(g).*I;
imshow(g)

Sobel滤波函数，SobelFilter如下：

%利用梯度来增强图像的边缘算子，一阶导数，锐化图像
function g=SobelFilter(f)
len=1;
%对原始图像进行扩展，此处采用了镜像扩展，目的是解决边缘计算的问题
f_pad=padarray(f,[len,len]);
[M,N]=size(f_pad);
%sobel算子
Lx=[-1 -2 -1;
    0 0 0;
    1 2 1];
Ly=[-1 0 1;
     -2 0 2;
     -1 0 1];     
%M(x,y)=sqrt(gx^2+gy^2)
for i=1+len:M-len
    for j=1+len:N-len
        %从扩展图像中，取出局部图像
        Block=f_pad(i-len:i+len,j-len:j+len);
        gx=sum(sum(Block.*Lx));
        gy=sum(sum(Block.*Ly));
        %取这组数的中值,赋值给输出图像        
%         g(i-len,j-len)=f(i-len,j-len)+ sqrt(gx^2+gy^2);
        g(i-len,j-len)=sqrt(gx^2+gy^2);
    end
end

Otsu函数，Otsu_FindThreshold如下：

%%%%%%%灰度图像二值分割门限
%I为输入灰度图像
function T=Otsu_FindThreshold(I,L_start)
    if isempty(find(I>1))
        Display('Please input image with integer type!');
        T=0;
        return;
    end
    %统计直方图
    [Pi,Pk]=Histeq_My(I,L_start);
    L=L_start:255;
    %计算全局的灰度平均值
    m_G=sum(L.*Pi);
    k_optimum=0;
    delta_max=0;
    %保存多个最大值
    ks=[];
    for k=L_start:255
        %计算第k平均值
        i=L_start:k;    
        Pi_k=Pi(1,1:k+1-L_start);
        m_k=sum(i.*Pi_k);
        %如果被除数为0或1，其值直接为0。详细请参见公式
        if Pk(k+1-L_start)==0 ||Pk(k+1-L_start)==1
            delta=0;
        else
           delta=(m_G*Pk(k+1-L_start)-m_k)^2/((Pk(k+1-L_start))*(1-Pk(k+1-L_start)));
        end   
        %寻找最大值
        if delta>delta_max
            delta_max=delta;
            k_optimum=k;
            ks=[];            
        elseif delta==delta_max %如果有多个最大值，将其记录
            ks=cat(2,ks,k);
        end
        %求多个最大值下标的平均值
        if ~isempty(ks)
            k_optimum=mean(ks);
        end
    end
   T=k_optimum;
end

Histeq_My函数如下：

%%%直方图计算
%L_start取值为0和1
function [Pi,Pk]=Histeq_My(I,L_start)
    [M,N]=size(I);
    if L_start
        hist_I=zeros(1,255);
        %遍历图像，对每一个灰度分量进行数量统计
        for i=1:M*N
            if I(i)~=0
                hist_I(I(i))=hist_I(I(i))+1;
            end
        end      
        %找出图像中的非零元素的个数
        c=length(find(I>0));
        Pi=hist_I/c;
        Pk=zeros(1,255);
        for k=1:255
            Pi_k=Pi(1:k);
            Pk(k)=sum(Pi_k);
        end            
    else
        hist_I=zeros(1,256);
        %每一个值加1，消除零值，以与矩阵的下标对应
        I=I+1;
        %遍历图像，对每一个灰度分量进行数量统计
        for i=1:M*N
            hist_I(I(i))=hist_I(I(i))+1;
        end
        Pi=hist_I/(M*N);
        Pk=zeros(1,256);
        for k=1:256
            Pi_k=Pi(1:k);
            Pk(k)=sum(Pi_k);
        end
    end
    %显示直方图统计信息
%     for x=L_start:255
%         x_ps=[x,x];
%         y_ps=[Pi(x+1-L_start),0];
%         plot(x_ps,y_ps);
%         hold on
%     end
%     hold off    
end