% 最大类方差法实现自适应图像阈值分割% This implements is implemented by WENG, All rights reserved% Copy Rights (c) WENG% 2016-3-30 1639clear;clc;close all;warning off; %#ok
      
       IMG = imread('steger_perfect_1_rg.tif');% fdoa1_rs1_rg.tifIMG_gray=double(rgb2gray(IMG));[M,N]=size(IMG_gray);gray_level=256;% HistgramHistgram_normal=zeros(256,1);% value span, 1~256for i_level=1:256,    Histgram_normal(i_level)=sum(sum(IMG_gray==(i_level+1)));endHistgram_normal=Histgram_normal/(M*N); % histgram normalizem_overall=0;for m=1:gray_level,m_overall=m_overall+(m-1)*Histgram_normal(m);endsigma_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
       
        sigma_max        optical_T=threshold;        sigma_max=sigma1;    endendTh_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 resultThreshold_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 evaluateSigma_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);