ostu进行遥感图像的分割

% 最大类方差法实现自适应图像阈值分割
% This implements is implemented by WENG, All rights reserved
% Copy Rights (c) WENG
% 2016-3-30 1639


clear;clc;close all;
warning off; %#ok<WNOFF>
IMG = imread('steger_perfect_1_rg.tif');% fdoa1_rs1_rg.tif
IMG_gray=double(rgb2gray(IMG));
[M,N]=size(IMG_gray);

gray_level=256;


% Histgram
Histgram_normal=zeros(256,1);% value span, 1~256
for i_level=1:256,
    Histgram_normal(i_level)=sum(sum(IMG_gray==(i_level+1)));
end
Histgram_normal=Histgram_normal/(M*N); % histgram normalize

m_overall=0;
for m=1:gray_level,m_overall=m_overall+(m-1)*Histgram_normal(m);end

sigma_max=0;
sigma_B=zeros(gray_level,1);
for i=1:gray_level,
    
    threshold=i-1;
    
    % calculate the number of the two class, the number is 0~1
    num1=0; 
    for m=1:threshold-1
        num1=num1+Histgram_normal(m);
    end
    num2=1-num1;
    
    % cal m1, m2
    m1=0;
    m2=0;
    for m=1:gray_level
        if m<threshold
            m1=m1+(m-1)*Histgram_normal(m);
        else
            m2=m2+(m-1)*Histgram_normal(m);
        end
    end
    
    % calculate the std betweent the two classes
    m1=m1/num1;
    m2=m2/num2;
    sigma1=num1*(m1-m_overall)^2+num2*(m2-m_overall)^2;
    
    % save and update
    sigma_B(i)=sigma1;
    if sigma1>sigma_max
        optical_T=threshold;
        sigma_max=sigma1;
    end
end

Th_optical=optical_T;
Th_matalb=graythresh(IMG_gray)*255;%matlab函数求阈值
fprintf(1,'The otsu threshold and Threshold calculate by graythresh are:\n %.1f [best], %.1f\n\n',Th_optical,Th_matalb);

% visual the result
Threshold_Optimal_IMG=zeros(M,N);
Threshold_graythresh_IMG=zeros(M,N);
Threshold_Optimal_IMG(IMG_gray>Th_optical)=1;
Threshold_graythresh_IMG(IMG_gray>Th_matalb)=1;
figure,subplot(2,2,1);imshow(IMG,[]);title('Ori color image');
subplot(2,2,2);imshow(IMG_gray,[]);title('Ori gray image');
subplot(2,2,3);imshow(Threshold_graythresh_IMG,[]);title('Threshold calculated by graythresh');
subplot(2,2,4);imshow(Threshold_Optimal_IMG,[]);title('Threshold calculated by ostu');


figure,subplot(1,2,1);
hold on; 
plot(1:gray_level,sigma_B);title('The two TH marked on \sigma array');
xlabel('Th'); ylabel('\sigma_B');
plot(Th_optical+1,sigma_B(Th_optical+1),'xg','MarkerSize',8,'LineWidth',2); 
text(Th_optical+1,sigma_B(Th_optical+1)-50,'otsu threshold');
plot(Th_matalb+1,sigma_B(Th_matalb+1),'^r','MarkerSize',8,'LineWidth',2); 
text(Th_matalb+1,sigma_B(Th_matalb+1)-50,'graythresh calculated threshold');
hold off;

subplot(1,2,2);plot(1:gray_level,Histgram_normal);title('Gray image Histgram');
hold on;
plot(Th_optical+1,Histgram_normal(Th_optical+1),'xg','MarkerSize',8,'LineWidth',2);
text(Th_optical+1,Histgram_normal(Th_optical+1)-50,'otsu threshold');
plot(Th_matalb+1,Histgram_normal(Th_matalb+1),'^r','MarkerSize',8,'LineWidth',2);
text(Th_matalb+1,Histgram_normal(Th_matalb+1)-50,'graythresh calculated threshold');
hold off;

% Image that can be classified by the threshold evaluate
Sigma_G=std(IMG_gray(:));
eta_ostu=sigma_B(Th_optical+1)/Sigma_G;
eta_matlab=sigma_B(Th_matalb+1)/Sigma_G;
fprintf(1,'The separability of ostu and matlab graythresh func respectively are:\n%.3f, %.3f\n',eta_ostu,eta_matlab);