function [Q1]=compute(temp_m,graynum,dist)
%%%%%灰度共生矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%graynum灰度级数
%dist距离

for i=1:graynum
    for j=1:graynum
        matH(i,j)=0;matV(i,j)=0;matRD(i,j)=0;matLD(i,j)=0;
    end
end
[M,N]=size(temp_m);

double total;total=0;

temp_m=double(temp_m);
       for m=1:M
            for n=1:N
                if ceil(temp_m(m,n)/(256/graynum))==0
                    newimage(m,n)=1;
                else
                    newimage(m,n)=ceil(temp_m(m,n)/(256/graynum));%灰度量化
                end
            end
        end
        
        %%% 0
        for m=1:M
            for n=1:N-dist
                a=newimage(m,n);
                b=newimage(m,n+dist);
                matH(a,b)=matH(a,b)+1;
                matH(b,a)=matH(b,a)+1;
            end
        end
%            
%         %%%%%%%%%%%%%%% 90
%         for m=1:16-dist
%             for n=1:16
%                 matV(newimage(m,n),newimage(m+dist,n))=matV(newimage(m,n),newimage(m+dist,n))+1;
%                 matV(newimage(m+dist,n),newimage(m,n))=matV(newimage(m+dist,n),newimage(m,n))+1;
%             end
%         end
%         
%         %%%%%%%%%%%%%%%%%%%   135
%         for m=1:16-dist
%             for n=16-dist
%                 matLD(newimage(m,n),newimage(m+dist,n+dist))= matLD(newimage(m,n),newimage(m+dist,n+dist))+1;
%                 matLD(newimage(m+dist,n+dist),newimage(m,n))= matLD(newimage(m+dist,n+dist),newimage(m,n))+1;
%             end
%         end
%         
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   45
%         for m=1:M-dist
%             for n=1:N-dist
%                  a=newimage(m,n);
%                  b=newimage(m+dist,n+dist);
%                  matH(a,b)=matH(a,b)+1;
%                  matH(b,a)=matH(b,a)+1;
% %                 matRD(newimage(m,n),newimage(m-dist+1,n+dist))=matRD(newimage(m,n),newimage(m-dist+1,n+dist))+1;
% %                 matRD(newimage(m-dist+1,n+dist),newimage(m,n))=matRD(newimage(m-dist+1,n+dist),newimage(m,n))+1;
%             end
%         end
        
          
       %%%%%%%%%%%%%%
            for m=1:graynum
                for n=1:graynum
                    total=total+matH(m,n);
                end
            end
        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%5    
          for m=1:graynum
                for n=1:graynum
                    mat_temp(m,n)=matH(m,n)/total;
                    mat_temp(m,n)=double(mat_temp(m,n));
                end
            end
%%%%%%%%%%计算特征值
Q1=0;Q2=0;
for m=1:graynum
    for n=1:graynum
          Q1=Q1+mat_temp(m,n)*mat_temp(m,n);%%能量
          Q2=Q2+(m-n)*(m-n)*mat_temp(m,n);%%%%对比度
      end
  end
  

A=imread('E:\image\sy24.gif');
  imshow(A);
% % difference(A)
% % laws(a);
% subplot(1,2,1),imshow(A);
[m,n]=size(A);
% A=padarray(A,[21 21],'replicate','post');
t=42;  %窗口大小

for i=1:m-t+1
       for j=1:n-t+1
             A1{i,j}=A(i:i+t-1,j:j+t-1);
             Q1(i,j)=compute(A1{i,j},8,4);

    
       end
end

 [M,N]=size(Q1);
c=reshape(Q1,M*N,1);
[IDC,U]=kmeans(c,2);
       c(IDC==1,1)=0;
       c(IDC==2,1)=0.75;

g=reshape(c,M,N);
g=im2uint8(g);
 figure,imshow(g);
%  %%计算分割正确率
%  [M,N]=size(g);
% g2=zeros(M);
%  for i=1:M/2
%      for j=1:N/2
%          g2(i,j)=191;
%      end
%  end
%  figure,imshow(g2);
%  g3=zeros(M);
%  for i=1:M
%      for j=1:N
%          if g2(i,j)==g(i,j)
%              g3(i,j)=1;
%          end
%      end
%  end
%  g4=sum(g3(:))/(M*N)
%  g5=1-g4
%