function out=texturefeature(img,type)

%%%调用了函数gabor2,gaborfilter1,Algorithmslib

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   img 灰度图像
%
% type:
%	a:	Homogeneous Texture Descriptor (HTD)
            
%       textureDescriptors=[Fdc,Fsd,C1,C2,...C30,D1,D2,...,D30]
%
%   b:  Texture Browsing Descriptor	 
%           TBD=[V1,V2,V3,V4,V5];
%%          V1 :  纹理图像的规则程度｛1，2，3，4｝
%           V2,V4 : 纹理图像主方向性｛0，1，2，3，4，5，6｝
%%          V3,V5 : 纹理图像主要纹理大小｛1，2，3，4｝
%
%   c:   Edge Histogram Descriptor 
%
%   d:  coMatrix：基于灰度共生矩阵的纹理特征
%
%   e:  autocorrelation
%    
%   f:  edge frequency 边缘频率
%
%   g: tamura_feature out=[粗糙度,对比度,方向度]
%
%   h: Extract the Primitive Length vector from an image：基元长度（行程）
%
%   k: Gabor 纹理
%
%   l:  Fourier 纹理
%
%   m:  Simultaneous Autoregression texture model  自回归纹理模型
%
%   n:  分形纹理描述符



switch type
    case 'a'
    %%% 同质纹理 Homogeneous Texture Descriptor (HTD)
        [row,col]=size(img);
        img1=reshape(img,1,row*col);
        DC=mean(img1);  %1:均值Fdc
        SD=std(img1);   %标准差Fsd
        clear img1
        pfft2img=Algorithmslib(img,'pfft2');
        Ws=[3/4,3/8,3/16,3/32,3/64];         
        thetar=[0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6];
        
        for ijn=0:29
            j=mod(ijn,6)+1;
            i=floor(ijn/6)+1;                
            Hij = pfft2img.*(gabor2(img,Ws(i),thetar(j)));
            bb=ijn+1;
            p(bb)=sum(sum(Hij.*Hij));%平均能量
            q(bb)=sqrt(sum(sum(Hij.*Hij-p(bb))));%能量标准差
            ee(bb)=log(1+p(bb));
            dd(bb)=log(1+q(bb));
        end
        clear pfft2img
        clear Hij
        out=[DC,SD,ee,dd];
        
    case 'b'
        [row,col]=size(img);
        if(row>800)            
            img=double(img);
            [c,s]=wavedec2(img,2,'db2');
            ca=appcoef2(c,s,'db2',1);
            Max=max(max(ca));
            Min=min(min(ca));
            img=255*(ca-Min)/(Max-Min);
        end
        [row,col]=size(img);
    %%%   纹理浏览描述符  Texture Browsing Descriptor
        %%%%%%%%
        
       %   TBD=[V1,V2,V3,V4,V5]; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     求： V2  V4   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %   例子： f:0,2,4,8
        %   theta = 0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6
        %   [G,gabout] = gaborfilter1(img,2,4,16,pi/3);
        %   figure,imshow(uint8(gabout));
        
        Wsk=ones(row*4,col*6);
        %Wskseg=mat2cell(Wsk,[row,row,row,row],[col,col,col,col,col,col]);
        SWsk=zeros(row*4,col);
        SWskseg=mat2cell(SWsk,[row,row,row,row],[col]);
        S=[0,2,4,8];
        K=[0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6];

         freq=pi/1.4;
%         width=5;
        for s=1:4       %s表示Gabor滤波器的尺度数
            for k=1:6       %k表示Gabor滤波器的方向数
                sigma=pi/2.^(s-1);
                theta=pi*(k-1)/8;
                [G,Wskseg{s,k}]=gaborfilter1(img,sigma,freq,theta,5);
                SWskseg{s,1}=SWskseg{s,1}+Wskseg{s,k};
            end
             Uu(s)=mean(mean(SWskseg{s,1}/6));
             sigu(s)=std(reshape(SWskseg{s,1}/6,1,row*col));
             Tu(s)=Uu(s)+sigu(s);           
        end
        Wsk=cell2mat(Wskseg);
        clear SWsk;
        clear SWskseg;
        
        NumK=zeros(4,1);
        for skn=0:23
            k=mod(skn,6)+1;
            s=floor(skn/6)+1;
            Num=size(find(Wskseg{s,k}>Tu(s)));
            Nsk(skn+1)=Num(1)*Num(2);
            NumK(s)=NumK(s)+Nsk(skn+1);
        end
        clear Wskseg;

        Hsk=[Nsk(1:6)/NumK(1),Nsk(7:12)/NumK(2),Nsk(13:18)/NumK(3),Nsk(19:24)/NumK(4)];
        Hsk=reshape(Hsk,6,4); 
        Hsk=Hsk';
        Hsk00=[[0;0;0;0],Hsk,[0;0;0;0]];
        for skn=0:23
            s=mod(skn,4)+1;
            k=floor(skn/4)+1;
            CC0(s,k)=(2*Hsk00(s,k+1)-Hsk00(s,k)-Hsk00(s,k+2))/2;
        end
        [SignS0,SignK0]=find(CC0==max(max(CC0)));
        
       if(size(SignS0)<2)            
            CC01=[min(min(CC0))*ones(1,6);CC0;min(min(CC0))*ones(1,6)];
            maxS=max(CC01'); %  每行最大值       
            if(CC01(SignS0,SignK0)==maxS(SignS0)|CC01(SignS0+2,SignK0)==maxS(SignS0+2))%最大值的上方值或下方值仍是所在行的最大值
                V2=SignK0;
                Hsk11=Hsk00;
                Hsk11(:,V2)=[];
                %%再次找最佳方向
                for skn=0:19
                    s=mod(skn,4)+1;
                    k=floor(skn/4)+1;
                    CC1(s,k)=(2*Hsk11(s,k+1)-Hsk11(s,k)-Hsk11(s,k+2))/2;
                end

                    [SignS1,SignK1]=find(CC1==max(max(CC1)));
                    
                    if(size(SignS1)<2)
                        CC11=[min(min(CC1))*ones(1,5);CC1;min(min(CC1))*ones(1,5)];
                        maxS1=max(CC11');
                        if(CC11(SignS1,SignK1)==maxS1(SignS1)|CC11(SignS1+2,SignK1)==maxS1(SignS1+2))
                            if(SignK1<V2)
                                V4=SignK1;
                            else
                                V4=SignK1+1;
                            end
                        else
                            V4=0;
                        end
                    else
                        if(SignK1<V2)
                            V4=SignK1(1);
                        else
                            V4=SignK1(1)+1;
                        end
                    end
           else
                V2=0;V4=0;
            end
      else
          V2=SignK0(1);
          V4=SignK0(2);                               
      end
  

        clear CC0;clear CC01;clear CC1;clear CC11
        clear Hsk;clear Hsk00;clear Hsk11
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     求： V3  V5   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
        [Phsk,Pvsk] = PVHsk(Wsk,V2,V4);
        
        Phskseg=mat2cell(Phsk,[1,1,1,1],[256,256,256,256,256,256]);
        Pvskseg=mat2cell(Pvsk,[1,1,1,1],[256,256,256,256,256,256]);  

        for skn=0:23
            s=mod(skn,4)+1;
            k=floor(skn/4)+1;
            aa=1./sqrt(sum(Phskseg{s,k}.^2));
            bb=ones(256,1)*Phskseg{s,k};
            bb=triu(bb)';
            aa=aa./sqrt(sum(bb.^2));
            cc=conv(Phskseg{s,k},fliplr(Phskseg{s,k}));
            cc=cc(1:256);cc=triu(cc);
            abc=cc.*aa;            
            NACsk{s,k}=abc;
            aa1=1./sqrt(sum(Pvskseg{s,k}.^2));
            bb1=ones(256,1)*Pvskseg{s,k};
            bb1=triu(bb1)';
            aa1=aa1./sqrt(sum(bb1.^2));
            cc1=conv(Pvskseg{s,k},fliplr(Pvskseg{s,k}));
            cc1=cc1(1:256);cc1=triu(cc1);
            abc1=cc1.*aa1;
            NACsk1{s,k}=abc1;            
        end
        
        NACkh=cell2mat(NACsk);NACkv=cell2mat(NACsk1);
        %clear NACsk ;clear NACsk1; 
        clear aa ; clear bb ; clear ab ; clear cc ; clear abc ;
        clear aa1; clear bb1; clear ab1; clear cc1; clear abc1;
        %%%%%%%%%%%%%%%%        找峰值和谷值      %%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        for skn=0:23
            s=mod(skn,4)+1;
            k=floor(skn/4)+1;
            NACsk01=NACsk{s,k};                            NACsk11=NACsk1{s,k};
            NACsk00=[NACsk01(1),NACsk01(1:255)];           NACsk10=[NACsk11(1),NACsk11(1:255)];
            NACsk02=[NACsk01(2:256),NACsk01(256)];         NACsk12=[NACsk11(2:256),NACsk11(256)];
            sigP=find(NACsk01>=NACsk00&NACsk01>=NACsk02);  sigP1=find(NACsk11>=NACsk10&NACsk11>=NACsk12);
            avP=mean(NACsk01(sigP));                       avP1=mean(NACsk11(sigP1));
            sigV=find(NACsk01<=NACsk00&NACsk01<=NACsk02);  sigV1=find(NACsk11<=NACsk10&NACsk11<=NACsk12);
            avV=mean(NACsk01(sigV));                       avV1=mean(NACsk11(sigV1));
            cc(skn+1)=avP-avV;                             cc1(skn+1)=avP1-avV1;
            Gph(skn+1)=mean(sigP);               Gpv(skn+1)=mean(sigP1);
            %/std(sigP) /std(sigP1)
        end
        clear NACsk ;clear NACsk1;
        clear NACsk01 ;clear NACsk00 ; clear NACsk02 ;
        clear NACsk11 ;clear NACsk10 ; clear NACsk12 ;
        V3sk=find(max(cc));
        V3=mod(V3sk-1,4)+1;
        V5sk=find(max(cc1));
        V5=mod(V5sk-1,4)+1;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%    求： V1        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Gphposi=find(Gph<0.33);
        Gh=zeros(1,24);
        Gh(Gphposi)=1;
        Gpvposi=find(Gpv<0.33);
        Gv=zeros(1,24);
        Gv(Gpvposi)=1;
        Mh=0;Mv=0;
        if(Gh~=0)
            Gh=reshape(Gh,4,6);
            Gh00=[zeros(6,1),[zeros(1,6);Gh;zeros(1,6)],zeros(6,1)];
            jj=size(Gphposi);
            for ij=1:jj
                skn=Gphposi(ij);
                s=mod(skn-1,4)+1;
                k=floor(skn/4)+1;
                if(Gh00(s+1,k)~=0|Gh00(s+1,k+2)~=0|Gh00(s,k+1)~=0|Gh00(s+2,k+1)~=0)
                    Mh=Mh+1;
                else
                    Gh00(s+1,k+1)=0;
                    if (Gh00(s+1,:)~=0|Gh00(:,k+1)~=0)
                        Mh=Mh+0.3;
                    else
                        Mh=Mh+0;
                    end
                end
            end
        else
            Mh=0;
        end
        clear Gh00;clear Gh;clear Gph;
        
        if(Gv~=0)
            Gv=reshape(Gv,4,6);Gv=reshape(Gv,4,6);
            Gv00=[zeros(6,1),[zeros(1,6);Gv;zeros(1,6)],zeros(6,1)];
            jj=size(Gpvposi);
            for ij=1:jj
                skn=Gpvposi(ij);
                s=mod(skn-1,4)+1;
                k=floor(skn/4)+1;
                if(Gv00(s+1,k)~=0|Gv00(s+1,k+2)~=0|Gv00(s,k+1)~=0|Gv00(s+2,k+1)~=0)
                    Mv=Mv+1;
                else
                    Gv00(s+1,k+1)=0;
                    if (Gv00(s+1,:)~=0|Gv00(:,k+1)~=0)
                        Mv=Mv+0.3;
                    else
                        Mv=Mv+0;
                    end
                end
            end
        else
            Mv=0;
        end
        clear Gv00;clear Gv;clear Gpv;
        Mimg=Mh+Mv;
        if(Mimg<3)
            V1=0;
        elseif(Mimg<6)
            V1=1;
        elseif(Mimg<15)
            V1=2;
        else
            V1=3;
        end
        out=[V1,V2,V3,V4,V5];       
     
    case 'c'
    %%% 边缘方向直方图  Edge Histogram Descriptor
        img=double(img);
        out=Algorithmslib(img,'edgehist');
        
    case 'd'
    %%   基于共生矩阵的纹理特征
        out=Algorithmslib(img,'basedcoMatrix');
        
    case 'e'
    %%   自相关
        %out=corrcoef(img);
        img=double(img);
        [M,N] = size(img);
        if(M>800)
            [c,s]=wavedec2(img,2,'db2');
            ca=appcoef2(c,s,'db2',1);
            Max=max(max(ca));
            Min=min(min(ca));
            img=255*(ca-Min)/(Max-Min);
        end
        
        [M,N] = size(img);
        img=double(img);
        for p=1:15-1
            for q=1:15-1
                img1=img(1:M-p,1:N-q);
                img2=img(1+p:M,1+q:N);
                s1=img1.*img2;
                s1=sum(sum(s1));
                s2=img.*img;
                s2=sum(sum(s2));
                Acor(p,q)=M*N*s1/(M-p)/(N-q)/s2;
            end
        end
        out=Acor;
        out=reshape(out,1,[]);
        clear img1;clear img2;clear Acor;
        
    case 'f'
    %%% 边缘频率
        img=double(img);
        [M,N,R]=size(img);
        if(M>800)
            [c,s]=wavedec2(img,2,'db2');
            ca=appcoef2(c,s,'db2',1);
            Max=max(max(ca));
            Min=min(min(ca));
            img=255*(ca-Min)/(Max-Min);
            [c,s]=wavedec2(img,2,'db2');
            ca=appcoef2(c,s,'db2',1);
            Max=max(max(ca));
            Min=min(min(ca));
            img=255*(ca-Min)/(Max-Min);
        end
        
        [M,N]=size(img);
        dd=max(M,N)-1;
        for d=1:dd
            if(d<M)
                img1=[img(d+1:M,:,:);img(M+1-d:M,:,:)];
                img1=img-img1;
                upd=sum(sum(abs(img1)));
                clear img1;
            else
                upd=0;
            end
            if(d<M)
                img2=[img(1:d,:,:);img(1:M-d,:,:)];
                img2=img-img2;
                downd=sum(sum(abs(img2))); 
                clear img2;
            else
                upd=0;
                
            end
            if(d<N)
                img3=[img(:,d+1:N,:) img(:,N+1-d:N,:)];
                img3=img-img3;
                leftd=sum(sum(abs(img3))); 
                clear img3;
            else
                leftd=0;
            end
            if(d<N)
                img4=[img(:,1:d,:) img(:,1:N-d,:)];
                img4=img-img4;
                rightd=sum(sum(abs(img4)));
                clear img4;
            else
                rightd=0;
            end
            out(d,:,:)=upd+downd+leftd+rightd;     
        end        
        
    case 'g'
    %%%     
         out=Algorithmslib(img,'tamura_feture');   
        
    case 'h'
    %%%    基元长度
        [M,N]=size(img);
        img=uint8(img);
        img=uint16(img)+1;
        B=zeros(256,N);
        kj=1;
        for i=1:M
            kj=1;
            while(kj<=N)
                n=1;
                while (img(i,kj)==img(i,kj+1)) %统计第i行，相邻像素点相同灰度级的个数n
                    n=n+1;
                    kj=kj+1;
                    if(kj==N)break;
                    end
                end
                B(img(i,kj),n)=B(img(i,kj),n)+1;
                kj=kj+1;
                if(kj==N)
                    B(img(i,kj),1)=B(img(i,kj),1)+1; %B矩阵256*1,256代表灰度级，第1列存储的是所有行里对应的相邻的灰度级的个数
                    break;
                end                
            end                
        end
        maxc1=1;
        for i=1:256
            if(B(i,:)==0)
                cc1=0;
            else
                cc1=max(find(B(i,:)~=0));
            end
            maxc1=max(maxc1,cc1);
        end
        B=double(B(:,1:maxc1));
        
        %%K为行程的总数
        K=sum(sum(B));
        %%% 短基元加重 out(1)
        RR=repmat([1:maxc1],256,1);
        RR=RR.*RR;
        B1=B./RR;
        out(1)=sum(sum(B1))/K;
        clear B1;
        %%  长基元加重 out(2)
        B2=B.*RR;
        out(2)=sum(sum(B2))/K;
        clear B2;
        %%  灰度级的均匀性 out(3)
        out(3)=(sum(sum(B').^2))/K;
        %%  基元长度的均匀性  out(4)
        out(4)=(sum(sum(B).^2))/K;
        %% 基元百分比    out(5)
        out(5)=K/M/N;
        clear B;
        
    case 'k'
    %%%%  Gabor 纹理
        [row,col]=size(img);
%        Wsk=ones(row*4,col*6);

         freq=pi/1.4;
        for s=1:4       %s表示Gabor滤波器的尺度数
            for k=1:6       %k表示Gabor滤波器的方向数
                sigma=pi/2.^(s-1);
                theta=pi*(k-1)/8;
                [G,Wskseg{s,k}]=gaborfilter1(img,sigma,freq,theta,5);
                Wskseg{s,k}=Wskseg{s,k}.^2;
                sumG(s,k)=sum(sum(Wskseg{s,k}));
            end
        end
        clear Wskseg;clear G;
        out=sumG/sum(sum(sumG));
        out=reshape(out,1,[]);
        clear sumG;
       
    case 'l1'
    %%%Fourier 纹理
        img=double(img);
        [M,N] = size(img);
        if(M>800)
            [c,s]=wavedec2(img,2,'db2');
            ca=appcoef2(c,s,'db2',1);
            Max=max(max(ca));
            Min=min(min(ca));
            img=255*(ca-Min)/(Max-Min);
        end
        [M,N] = size(img);
        imgfft2=fft2(img,16,16);
        Rimgfft=real(imgfft2);  %   实部
        Iimgfft=imag(imgfft2);  %   虚部
        clear imgfft2;
        %%  幅度
        Apimgf=sqrt(Rimgfft.^2+Iimgfft.^2);      
        clear Rimgfft;clear Iimgfft;
        
        out=Apimgf;out=reshape(out,1,[]);
        clear Apimgf;
          
    case 'l2'
    %%%Fourier 纹理
        img=double(img);
        [M,N] = size(img);
        if(M>800)
            [c,s]=wavedec2(img,2,'db2');
            ca=appcoef2(c,s,'db2',1);
            Max=max(max(ca));
            Min=min(min(ca));
            img=255*(ca-Min)/(Max-Min);
        end
        imgfft2=fft2(img,16,16); %频谱矩阵16*16
        Rimgfft=real(imgfft2);  %   实部
        Iimgfft=imag(imgfft2);  %   虚部
        clear imgfft2;
         
        %% 相角
        Agimgf=atan2(Rimgfft,Iimgfft)/pi*180;
        
        clear Rimgfft;clear Iimgfft;
        out=Agimgf;out=reshape(out,1,[]);
        clear Agimgf;
        
%     case m
%         %%自回归模型纹理特征
%         
%         
%         
%     case n
%         %%%分形纹理描述符
        
        
        
    otherwise
        error(['Unknown texture feature .']);     
        
end