图像处理：自动阈值选取算法

      阈值分割是图像处理的基本问题，在图像分析和识别中起着重要作用。有关图像阈值自动选取方法的已多达数十种。

以下是 循环分割算法 和 最大类间方差法 的matlab代码实现自动阈值选取~

1. 循环分割算法

       适用情况：

  循环分割算法的优点在于实现起来十分简单,容易提取灰度值差异十分明显的变化像元 , 忽略灰度值相对较小变化像元。

  以下是对一副NDVI影像的自动阈值选取的matlab代码，目的是分割植被与非植被区域：

% 读取 tif 数据
[N2a,~]=geotiffread('NDVI.tif');
% 读取 tif 最值
[N2min,N2max]=deal(min(min(N2a)),max(max(N2a)));
% 循环分割算法，参考论文《基于遥感影像的变化检测方法在滑坡体提取中的应用-陈》第3.2.1节描述
% 灰度级
T1=(N2min+N2max)/2;
u1=sum(sum(N2a(N2a>=T1)))/sum(sum(N2a>=T1));
u2=sum(sum(N2a(N2a<T1)))/sum(sum(N2a<T1));
T2=(u1+u2)/2;
while T2-T1>=0.01   % 0.01为最终迭代的差，该值可依据需要设定
    T1=T2;
    u1=sum(sum(N2a(N2a>=T1)))/sum(sum(N2a>=T1));
    u2=sum(sum(N2a(N2a<T1)))/sum(sum(N2a<T1));
    T2=(u1+u2)/2;
end
Tg=T2;  % Tg 即为分割阈值

2. 最大类间方差法（又称Otsu或大津法）

      适用情况：

      这类算法是面向均匀性的算法 , 所得到二值图像比较均匀且目标性状较好。该方法具有抗亮度变化和对比度变化的特点，因此非常适合于实时图像处理。

     以下是对一副山体阴影影像的自动阈值选取的matlab代码，目的是提取阴影分割阈值：

% 读取 tif 数据
[N2a,~]=geotiffread('hillshade.tif');
% 读取 tif 最值
[N2min,N2max]=deal(min(min(N2a)),max(max(N2a)));
% 最大类间方差法，参考论文《一些新的图像阈值选取方法-付》第3节描述，文中频数f和频率p说反了
% 将tif数据直方图归一化
N2t=(255-0)/(N2max-N2min)*(N2a-N2min)+0;
% 改为1维，方便使用tabulate函数（若用循环实现频数、频率计算，运行速度慢的很）
N2t_row=reshape(N2t,1,size(N2t,1)*size(N2t,2));
N2t_count=tabulate(N2t_row); % 三列分别为灰度值、频数、频率
N2t_ave_mat=N2t_count(:,1).*N2t_count(:,3); % 均值分子矩阵

G=zeros(1,254);
GG=zeros(1,254);
for t=1:254
    % 求取目标个数
    num_pre=sum(N2t_ave_mat<=t);
    % 求取背景个数
    num_post=sum(N2t_ave_mat>t); 
    % 目标均值
    delta_pre=sum(N2t_ave_mat(1:num_pre,1))/sum(N2t_count(1:num_pre,3));
    % 背景均值
    delta_post=sum(N2t_ave_mat(num_pre+1:end,1))/sum(N2t_count(num_pre+1:end,3));
    % 目标部分比例
    w_pre=sum(N2t_count(1:num_pre,2))/(size(N2t,1)*size(N2t,2));
    % 背景部分比例
    w_post=sum(N2t_count(num_pre+1:end,2))/(size(N2t,1)*size(N2t,2));
    % 全图总均值
    delta=w_pre*delta_pre+w_post*delta_post;
    G(t)=(delta_post-delta)*(delta-delta_pre);
    GG(t)=(delta_post-delta_pre)^2*w_pre*w_post;
end
g=find(G==max(G)); % G最大值时的t值
gg=find(GG==max(GG)); 
% 反算到g, 按比例计算(255-0)/(N2max-N2min)=(255-gg)/(N2max-g)
tre=N2max-(N2max-N2min)/(255-0)*(255-g);
tre2=N2max-(N2max-N2min)/(255-0)*(255-gg);

       最大类间方差法代码计算了论文《一些新的图像阈值选取方法-付》第5节结尾中的两个公式的阈值分割结果，分别为tre，tre2。