二维最大类间方差(otsu)matlab实现源码

今天整理电脑资料，没想到翻到以前上学时候写的小论文，电脑里居然还存着以前用matlab写得最大类间方差的二维实现源码。当时持续写了好几个星期才搞定，忘了哪个是最终版本了，只好按照日期把最后的那个源码版本贴出来，不知道里面有没有bug。

%%传统二维Otsu方法function [idx idy]=otsuHistgram2(I)%I=imread('rice.png');noise=imnoise(I,'gaussian',0,0.03);I=noise;figure,imhist(I);figure,imshow(I);%% 求出图像的邻域和中值h=fspecial('average',[3 3]);[m,n]=size(I);totalnum=m*n;flt=filter2(h,I);%邻域均值flt=uint8(flt);I=uint8(I);%duality 二元性（灰度，邻域）%% 求（灰度，邻域）出现的频率num_bins = 256;f(1:num_bins,1:num_bins)=0;for r=1:m    for l=1:n        gry=I(r,l);  %灰度gray        adj=flt(r,l);%邻域adjacent        f(gry+1,adj+1)=f(gry+1,adj+1)+1;    endendPij=f./totalnum;%% 画出二维直方图[X,Y] = meshgrid(1:num_bins);Z=f; plot3(X,Y,Z);figure;mesh(X,Y,Z);colormap(cool); %% 二维直方图上总的均值矢量Ut=(uti,utj)'r=1:num_bins;muti=sum((r-1)*Pij(r,:));mutj=sum((r-1)*Pij(:,r)');%% 两类出现的概率 W0，W1=1-W0 w0=cumsum(Pij); w0=cumsum(w0,2);  w0=w0+(w0==0)*eps; %% 两类对应的均值矢量 Ui,Uj%%将[0:255]分别乘以Pij的每行和每列，然后再求(s,t)内的和aa=zeros(1,num_bins);cc=zeros(num_bins,1);for i=1:num_bins    aa=[aa;(i-1)*Pij(i,:)];    cc=[cc (i-1)*Pij(:,i)];endbb=aa(2:end,:);mui=cumsum(bb);mui=cumsum(mui,2); dd=cc(:,2:end);muj=cumsum(dd);muj=cumsum(muj,2); %% 求出使类间方差最大的阈值sigma_b_squared=((w0*muti-mui).^2+(w0*mutj-muj).^2)./(w0.*(1 - w0));maxval = max(sigma_b_squared(:));isfinite_maxval = isfinite(maxval);if isfinite_maxval    [idx idy]=find(sigma_b_squared==maxval);     idx=mean(idx)-1;    idy=mean(idy)-1; else    idx=0;    idy=0;end%% 根据阈值对图像进行分割out1=im2bw(I,graythresh(I));figure,imshow(out1);out=((I<=idx)&(flt<=idy));out=1-out;figure,imshow(out);