两次利用大津法对图像进行分割，将前景与背景，白色区与背景分割开，并二值化。

close all; clear; clc; I=imread('dan11.jpg'); g=rgb2gray(I); %g imshow(g),title('原始灰度图') %figure(2),imhist(g),title('原始灰度图的直方图') %使用Matlab函数计算阀值 %level=graythresh(g); %level %bw=im2bw(g,level);%阀值分割 %figure(2),subplot(132),imshow(bw),title('graythresh计算阀值'); %disp(strcat('graythresh计算灰度阀值',num2str(uint8(level*255)))); %Matlab实现简化大津法计算阀值 zmax=max(max(g));%取出最大灰度值，max只能找到每一列的最大值。 zmin=min(min(g));%取出最小灰度值 t1=double(zmin:zmax);%求初始阀值 %t1 gsize=size(g);%图像大小 %gsize muxsize=gsize(1)*gsize(2);%图像像素总个数 for k=1:length(t1) tk=t1(1,k);%从最小灰度值到最大值的计算方法 %定义前景和背景像素数 iforeground=0; ibackground=0; %定义前景和背景的灰度总和 foreground=0; background=0; for i=1:gsize(1) for j=1:gsize(2) tmp=g(i,j); if(tmp>=tk) %前景灰度值 iforeground=iforeground+1; foreground=foreground+double(tmp); else %背景灰度值 ibackground=ibackground+1; background=background+double(tmp); end end end w0=iforeground/muxsize; w1=ibackground/muxsize; u0=foreground/iforeground; u1=background/ibackground; t1(2,k)=w0*w1*(u0-u1)*(u0-u1); end %遍历后寻找g的第二行的最大值 omax=max(t1(2,:));%第二行方差的最大值，忽略NaN idx=find(t1(2,:)>=omax);%方差最大值所对应的列号，find只能检查一行中的数，或一维数组 t1=uint8(t1(1,idx));%从第一行取出灰度值作为阀值 %t1 disp(strcat('一次简化大津法计算灰度阀值t1=',num2str(t1))); s1=1; for i=1:gsize(1) for j=1:gsize(2) if (g(i,j) < t1) g(i,j)=0; else bw(1,s1)=g(i,j); s1=s1+1; end end end s1=s1-1; disp(strcat('米粒图像中米粒像素个数s=',num2str(s1)));%米粒图像中米粒像素个数 figure(2),imshow(g),title('简化大津法计算灰度阀值'); bwsize=size(bw); %bwsize muxsize1=bwsize(1)*bwsize(2);%第一次大津法后的图像像素总个数 zmax=max(bw);%取出最大灰度值 %max只能找到每一列的最大值。 zmin=min(bw);%取出最小灰度值 t2=double(zmin:zmax);%求初始阀值 for k=1:length(t2) tk=t2(1,k);%从最小灰度值到最大值的计算方法 %定义前景和背景像素数 iforeground=0; ibackground=0; %定义前景和背景的灰度总和 foreground=0; background=0; for i=1:bwsize(1) for j=1:bwsize(2) tmp=bw(i,j); if(tmp>=tk) %前景灰度值 iforeground=iforeground+1; foreground=foreground+double(tmp); else %背景灰度值 ibackground=ibackground+1; background=background+double(tmp); end end end w0=iforeground/muxsize1; w1=ibackground/muxsize1; u0=foreground/iforeground; u1=background/ibackground; t2(2,k)=w0*w1*(u0-u1)*(u0-u1); end %遍历后寻找g的第二行的最大值 omax=max(t2(2,:));%第二行方差的最大值，忽略NaN idx=find(t2(2,:)>=omax);%方差最大值所对应的列号，find只能检查一行中的数，或一维数组 t2=uint8(t2(1,idx));%从第一行取出灰度值作为阀值 %t2(1,end) %可能大津法计算的阀值有多个，去第一个 disp(strcat('第二次简化大津法计算灰度阀值t2=',num2str(t2(1,end)))); s2=1; for i=1:gsize(1) for j=1:gsize(2) if (g(i,j) < t2(1,end)) g(i,j)=0; else s2=s2+1; end end end s2=s2-1; disp(strcat('米粒图像中垩白像素个数s2=',num2str(s2)));%单米粒图像中米粒像素个数 figure(3),imshow(g),title('第二次简化大津法计算灰度阀值');