% 最大类方差法实现自适应图像阈值分割
% This implements is implemented by WENG, All rights reserved
% Copy Rights (c) WENG
% 2016-3-30 1639clear;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 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~1num1=0; for m=1:threshold-1num1=num1+Histgram_normal(m);endnum2=1-num1;% cal m1, m2m1=0;m2=0;for m=1:gray_levelif m<thresholdm1=m1+(m-1)*Histgram_normal(m);elsem2=m2+(m-1)*Histgram_normal(m);endend% calculate the std betweent the two classesm1=m1/num1;m2=m2/num2;sigma1=num1*(m1-m_overall)^2+num2*(m2-m_overall)^2;% save and updatesigma_B(i)=sigma1;if sigma1>sigma_maxoptical_T=threshold;sigma_max=sigma1;end
endTh_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);