function out=Algorithmslib(img,type)

%%%  type:
%%%         pfft2 ：图像fft变换的极坐标表示
%%%         edgehist：灰度图像的边缘方向直方图
%%%         basedcoMatrix：基于灰度共生矩阵的纹理特征
%%%         tamura_feature: 粗糙度; 对比度;方向度


switch type
    case 'pfft2'
        [row,col]=size(img);
        maxW=double(uint16(sqrt(row*row+col*col)));
        W=0:maxW/(row-1):maxW;
        theta=2*pi*(0:90/(col-1):90)/360;        
        fftimg=fft2(img,4*row,4*col);
        U=ceil(W'*cos(theta));
        V=ceil(W'*sin(theta));
        U=reshape(U,row*col,1);
        V=reshape(V,row*col,1);
        UV2=V*2*row+U+1;
        clear W
        clear theta
        clear W
        clear V
%       pfft2out=fftimg(UV2);

%%%     out: row*col
        out=reshape(fftimg(UV2),row,col);
        clear UV2
        clear fftimg
        
        
%     case 'gabor1'
%         %%  S:0,1,2,3,4
%         %%  R:0,1,2,3,4,5
%         Ws=[3/4,3/8,3/16,3/32,3/64];
%         dels=Ws/3/sqrt(2*log(2));
%         
%         thetar=[0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6];
%         delr=pi*pi/36/8/log(2);
%         for i=1:5
%             WO(i,:)=exp(-(W-Ws(i)).^2/2/dels(i)/dels(i));
%         end
%         WTH=zeros(1,300*300);
%         for i=1:6
%         ThetaO=exp(-(theta-thetar(1)).^2/2/delr(1)/delr(1));
%         wwth=[(ThetaO')*WO(1,:);(ThetaO')*WO(2,:);(ThetaO')*WO(3,:);(ThetaO')*WO(4,:);(ThetaO')*WO(5,:)];
%  %       WTH=cat(1,WTH,[ThetaO(i,:).*WO(1,:);ThetaO(i,:).*WO(2,:);ThetaO(i,:).*WO(3,:);ThetaO(i,:).*WO(4,:);ThetaO(i,:).*WO(4,:)]);
%         end
%         
%         %%%%   out :30*row     (W量化数目)
% %       gabor1out=WTH(2:31,:);  
%         out=WTH(2:31,:); 
%         clear WO
%         clear ThetaO
%         clear WTH
%         
    case 'edgehist'
        [row,col]=size(img);
        mm=max(max(img));
        mm=ones(row,col)*double(mm);
        
        gray_img=double(img);
        if((gray_img-mm)==0)
            out=zeros(1,8);
        else
        if(row>800) 
            Gray = blkproc(gray_img,[4 4],@mean);
            Gray = blkproc(Gray',[4 4],@mean)';
            
        else
            Gray=gray_img;
        end
        edge_canny = edge(Gray,'canny');
        gray_img=Gray;
        [row,col]=size(edge_canny);

        edge_num=sum(sum(edge_canny));%边缘总点数

        %计算梯度矢量Gx,Gy
        for cir1=2:row-1
            for cir2=2:col-1
                Gx(cir1,cir2)=sum(gray_img(cir1-1:cir1+1,cir2+1))...
                    -sum(gray_img(cir1-1:cir1+1,cir2-1))...
                    +gray_img(cir1,cir2+1)-gray_img(cir1,cir2-1);
                Gy(cir1,cir2)=sum(gray_img(cir1+1,cir2-1:cir2+1))...
                    -sum(gray_img(cir1-1,cir2-1:cir2+1))...
                    +gray_img(cir1+1,cir2)-gray_img(cir1-1,cir2);
                Gx(cir1,cir2)=Gx(cir1,cir2)+(Gx(cir1,cir2)==0)*1e-6;
                theta(cir1,cir2)=atan(Gy(cir1,cir2)/Gx(cir1,cir2));    %计算边缘方向
            end
        end

        for cir1=1:row
            for cir2=1:col
                if edge_canny(cir1,cir2)==1
                    edge_theta(cir1,cir2)=theta(cir1,cir2);
                else
                    edge_theta(cir1,cir2)=0;
                end
            end
        end

        TH=[tan(pi/16),tan(3*pi/16),tan(5*pi/16),tan(7*pi/16),tan(-1*pi/16),...
            tan(-3*pi/16),tan(-5*pi/16),Inf;tan(-pi/16),tan(1*pi/16),tan(3*pi/16),...
            tan(5*pi/16),tan(-3*pi/16),tan(-5*pi/16),tan(-7*pi/16),tan(7*pi/16)];

        for i = 1:8
            if i~=8
                direct_img = (edge_theta<TH(1,i) & edge_theta>=TH(2,i));
            else
                direct_img = ((abs(edge_theta)<TH(1,i)) & (abs(edge_theta)>=TH(2,i)));
            end
            edgehist(i) = sum(sum(direct_img & edge_canny));
        end

        %归一化
        out = edgehist/edge_num; 
        clear Gray;clear Gx;clear Gy;clear c;clear ca;
        clear direct_img;clear edge_canny;clear edge_theta;
        clear theta;clear gray_img;
        end
%   %----------------------------------------------------------------------------------------%%      
%         
%     %----------------------------
% %计算边缘方向直方图
% %输入：图片序号
% %输出:量化、归一化后的直方图(八个角度的概率统计)
% %----------------------------
% function edgehist=edgedir_his(num)
% 
% %num=1;
% imgdir='C:\MATLAB7\work\';
% imgsuffix='.jpg';
% rgb_img=imread([imgdir,num2str(num),imgsuffix]);
% gray_img=rgb2gray(rgb_img);
% [row,col]=size(gray_img);
% gray_img=double(gray_img);
% 
% %提取canny边缘算子
% edge_canny=edge(gray_img,'canny');
% edge_num=sum(sum(edge_canny));%边缘总点数
% 
% %计算梯度矢量Gx,Gy
% for cir1=2:row-1
%     for cir2=2:col-1
%         Gx(cir1,cir2)=sum(gray_img(cir1-1:cir1+1,cir2+1))...
%             -sum(gray_img(cir1-1:cir1+1,cir2-1))...
%             +gray_img(cir1,cir2+1)-gray_img(cir1,cir2-1);
%         Gy(cir1,cir2)=sum(gray_img(cir1+1,cir2-1:cir2+1))...
%             -sum(gray_img(cir1-1,cir2-1:cir2+1))...
%             +gray_img(cir1+1,cir2)-gray_img(cir1-1,cir2);
%         Gx(cir1,cir2)=Gx(cir1,cir2)+(Gx(cir1,cir2)==0)*1e-6;
%         theta(cir1,cir2)=atan(Gy(cir1,cir2)/Gx(cir1,cir2));    %计算边缘方向
%     end
% end
% 
% %edge_theta=theta&edge_canny;
% for cir1=1:row
%     for cir2=1:col
%         if edge_canny(cir1,cir2)==1
%             edge_theta(cir1,cir2)=theta(cir1,cir2);
%         else
%             edge_theta(cir1,cir2)=0;
%         end
%     end
% end
% 
% TH=[tan(pi/16),tan(3*pi/16),tan(5*pi/16),tan(7*pi/16),tan(-1*pi/16),...
%     tan(-3*pi/16),tan(-5*pi/16),Inf;tan(-pi/16),tan(1*pi/16),tan(3*pi/16),...
%     tan(5*pi/16),tan(-3*pi/16),tan(-5*pi/16),tan(-7*pi/16),tan(7*pi/16)];
% 
% for i = 1:8
%     if i~=8
%         direct_img = (edge_theta<TH(1,i) & edge_theta>=TH(2,i));
%     else
%         direct_img = ((abs(edge_theta)<TH(1,i)) & (abs(edge_theta)>=TH(2,i)));
%     end
%     edgehist(i) = sum(sum(direct_img & edge_canny));
% end
% 
% %归一化
% edgehist = edgehist/edge_num; 
%     
% %平滑直方图
% %for i=2:7
% %    edgehist1(i)=sum(edgehist(i-1:i+1))/3;
% %end
% %edgehist(2:7)=edgehist1(2:7);
% 
%     
%         
%         
%         
%         
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       
    case 'basedcoMatrix'
        [M,N] = size(img);
        Gray=zeros(M,N);
        for n = 1:256/16
            pp=find((n-1)*16<=img&img<n*16); 
            Gray(pp) = n-1;
        end
        P = zeros(16,16,4);
        for ijn = 1:M*N
            i=mod((ijn-1),M)+1;%取余
            j=floor((ijn-1)/M)+1;%取整
            if j<N
                m1=Gray(i,j)+1;n1=Gray(i,j+1)+1;
                P(m1,n1,1) = P(m1,n1,1)+1;
                P(n1,m1,1) = P(m1,n1,1);
            end
            if i>1&j<N
                m2=Gray(i,j)+1;n2=Gray(i-1,j+1)+1;
                P(m2,n2,2) = P(m2,n2,2)+1;
                P(n2,m2,2) = P(m2,n2,2);
            end
            if i<M
                m3=Gray(i,j)+1;n3=Gray(i+1,j)+1;
                P(m3,n3,3) = P(m3,n3,3)+1;
                P(n3,m3,3) = P(m3,n3,3);
            end
            if i<M&j<N
                m4=Gray(i,j)+1;n4=Gray(i+1,j+1)+1;
                P(m4,n4,4) = P(m4,n4,4)+1;
                P(n4,m4,4) = P(m4,n4,4);
            end
        end
        A=ones(1,16);
        B=diag(A);
        [m,n]=find(B~=0);
        P(m,n,:)=P(m,n,:)*2;

        for n = 1:4
            P(:,:,n) = P(:,:,n)/sum(sum(P(:,:,n)));%P是灰度图像Img的灰度共生矩阵
        end

        H = zeros(1,4);
        I = H;
        Ux = H;      Uy = H;
        deltaX= H;  deltaY = H;
        C =H;
        for n = 1:4
            E(n) = sum(sum(P(:,:,n).^2)); %%能量
            for ijn = 1:256
                i=mod(ijn-1,16)+1;
                j=fix((ijn-1)/16)+1;
                if P(i,j,n)~=0
                    H(n) = -P(i,j,n)*log(P(i,j,n))+H(n); %%熵
                end
                I(n) = (i-j)^2*P(i,j,n)+I(n);  %%惯性矩
                Ux(n) = i*P(i,j,n)+Ux(n); %相关性中μx
                Uy(n) = j*P(i,j,n)+Uy(n); %相关性中μy
            end
        end
        for n = 1:4
            for ijn = 1:256
                i=mod(ijn-1,16)+1;
                j=fix((ijn-1)/16)+1;
                deltaX(n) = (i-Ux(n))^2*P(i,j,n)+deltaX(n); %相关性中σx
                deltaY(n) = (j-Uy(n))^2*P(i,j,n)+deltaY(n); %相关性中σy
                C(n) = i*j*P(i,j,n)+C(n);
            end
            C(n) = (C(n)-Ux(n)*Uy(n))/deltaX(n)/deltaY(n); %相关性   
        end
        T(1) = mean(E);   T(2) = sqrt(cov(E));
        T(3) = mean(H);   T(4) = sqrt(cov(H));
        T(5) = mean(I);   T(6) = sqrt(cov(I));
        T(7) = mean(C);   T(8) = sqrt(cov(C));
        out=T;
        clear P;
        clear Gray;
        clear B;
        clear pp;
    case 'tamura_feture'
        
        %计算Tamura纹理特征

        [Nx,Ny] = size(img);
        if(Nx>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
        [Nx,Ny] = size(img);
            
        Ng=256;
        G=double(img);
        %计算粗糙度(coarseness)
        Sbest=zeros(Nx,Ny);
        E0h=zeros(Nx,Ny);E0v=zeros(Nx,Ny);
        E1h=zeros(Nx,Ny);E1v=zeros(Nx,Ny);
        E2h=zeros(Nx,Ny);E2v=zeros(Nx,Ny);
        E3h=zeros(Nx,Ny);E3v=zeros(Nx,Ny);
        E4h=zeros(Nx,Ny);E4v=zeros(Nx,Ny);
        E5h=zeros(Nx,Ny);E5v=zeros(Nx,Ny);
        flag=0;

        for i=1:Nx
            for j=2:Ny
                E0h(i,j)=G(i,j)-G(i,j-1);
            end
        end
        E0h=E0h/2;

        for i=1:Nx-1
            for j=1:Ny
                E0v(i,j)=G(i,j)-G(i+1,j);
            end
        end
        E0v=E0v/2;

        if (Nx<4||Ny<4)
            flag=1;
        end
        if(flag==0)
            for i=1:Nx-1
                for j=3:Ny-1
                    E1h(i,j)=sum(sum(G(i:i+1,j:j+1)))-sum(sum(G(i:i+1,j-2:j-1)));
                end
            end
            for i=2:Nx-2
                for j=2:Ny
                    E1v(i,j)=sum(sum(G(i-1:i,j-1:j)))-sum(sum(G(i+1:i+2,j-1:j)));
                end
            end
            E1h=E1h/4;
            E1v=E1v/4;
        end


        if (Nx<8||Ny<8)
            flag=1;
        end
        if(flag==0)
            for i=2:Nx-2
                for j=5:Ny-3
                    E2h(i,j)=sum(sum(G(i-1:i+2,j:j+3)))-sum(sum(G(i-1:i+2,j-4:j-1)));
                end
            end
            for i=4:Nx-4
                for j=3:Ny-1
                    E2v(i,j)=sum(sum(G(i-3:i,j-2:j+1)))-sum(sum(G(i+1:i+4,j-2:j+1)));
                end
            end
            E2h=E2h/16;
            E2v=E2v/16;
        end


        %图片大小必须大于16*16才能计算E3h、E3v
        if (Nx<16||Ny<16)
            flag=1
        end
        if(flag==0)
            for i=4:Nx-4
                for j=9:Ny-7
                    E3h(i,j)=sum(sum(G(i-3:i+4,j:j+7)))-sum(sum(G(i-3:i+4,j-8:j-1)));
                end
            end
            for i=8:Nx-8
                for j=5:Ny-3
                    E3v(i,j)=sum(sum(G(i-7:i,j-4:j+3)))-sum(sum(G(i+1:i+8,j-4:j+3)));
                end
            end
            E3h=E3h/64;
            E3v=E3v/64;
        end

        %图片大小必须大于32*32才能计算E4h、E4v 
        if (Nx<32||Ny<32)
            flag=1;
        end
        if(flag==0)
            for i=8:Nx-8
                for j=17:Ny-15
                    E4h(i,j)=sum(sum(G(i-7:i+8,j:j+15)))-sum(sum(G(i-7:i+8,j-16:j-1)));
                end
            end
            for i=16:Nx-16
                for j=9:Ny-7
                    E4v(i,j)=sum(sum(G(i-15:i,j-8:j+7)))-sum(sum(G(i+1:i+16,j-8:j+7)));
                end
            end
            E4h=E4h/256;
            E4v=E4v/256;
        end


        %图片大小必须大于64*64才能计算E5h、E5v
        if (Nx<64||Ny<64)
            flag=1;
        end
        if(flag==0)
            for i=16:Nx-16
                for j=33:Ny-31
                    E5h(i,j)=sum(sum(G(i-15:i+16,j:j+31)))-sum(sum(G(i-15:i+16,j-32:j-31)));
                end
            end
            for i=32:Nx-32
                for j=17:Ny-15
                    E5v(i,j)=sum(sum(G(i-31:i,j-16:j+15)))-sum(sum(G(i+1:i+32,j-16:j+15)));
                end
            end
            E5h=E5h/1024;
            E5v=E5v/1024;
        end

        for i=1:Nx
            for j=1:Ny
                [maxv,index]=max([E0h(i,j),E0v(i,j),E1h(i,j),E1v(i,j),E2h(i,j),E2v(i,j),E3h(i,j),E3v(i,j),E4h(i,j),E4v(i,j),E5h(i,j),E5v(i,j)]);
                k=floor((index+1)/2);
                Sbest(i,j)=2.^k;
            end
        end
        Fcoarseness=sum(sum(Sbest))/(Nx*Ny);
        %计算对比度
        [counts,graylevels]=imhist(img);
        PI=counts/(Nx*Ny);
        averagevalue=sum(graylevels.*PI);
        u4=sum((graylevels-repmat(averagevalue,[256,1])).^4.*PI);
        standarddeviation=sum((graylevels-repmat(averagevalue,[256,1])).^2.*PI);
        alpha4=u4/standarddeviation^2;
        Fcontrast=sqrt(standarddeviation)/alpha4.^(1/4);
        %计算方向度
        PrewittH=[-1 0 1;-1 0 1;-1 0 1];
        PrewittV=[1 1 1;0 0 0;-1 -1 -1];
        %计算横向梯度
        deltaH=zeros(Nx,Ny);
        for i=2:Nx-1
            for j=2:Ny-1
                deltaH(i,j)=sum(sum(G(i-1:i+1,j-1:j+1).*PrewittH));
            end
        end
        for j=2:Ny-1
            deltaH(1,j)=G(1,j+1)-G(1,j);
            deltaH(Nx,j)=G(Nx,j+1)-G(Nx,j);   
        end
        for i=1:Nx
            deltaH(i,1)=G(i,2)-G(i,1);
            deltaH(i,Ny)=G(i,Ny)-G(i,Ny-1);   
        end
        %计算竖向梯度
        deltaV=zeros(Nx,Ny);
        for i=2:Nx-1
            for j=2:Ny-1
                deltaV(i,j)=sum(sum(G(i-1:i+1,j-1:j+1).*PrewittV));
            end
        end
        for j=1:Ny
            deltaV(1,j)=G(2,j)-G(1,j);
            deltaV(Nx,j)=G(Nx,j)-G(Nx-1,j);   
        end
        for i=2:Nx-1
            deltaV(i,1)=G(i+1,1)-G(i,1);
            deltaV(i,Ny)=G(i+1,Ny)-G(i,Ny);   
        end
        %梯度向量模
        deltaG=(abs(deltaH)+abs(deltaV))/2;
        %梯度向量方向
        theta=zeros(Nx,Ny);
        for i=1:Nx
            for j=1:Ny
                if (deltaH(i,j)==0)&&(deltaV(i,j)==0)
                elseif deltaH(i,j)==0
                    theta(i,j)=pi;            
                else           
                    theta(i,j)=atan(deltaV(i,j)/deltaH(i,j))+pi/2;
                end
            end
        end
        theta1=reshape(theta,1,[]);
        phai=0:0.0001:pi;
        HD1=hist(theta1,phai);
        HD1=HD1/(Nx*Ny);
        HD2=zeros(size(HD1));
        %定义一个阈值THRESHOLD
        THRESHOLD=0.01;
        for m=1:length(HD2)
            if HD1(m)>=THRESHOLD
                HD2(m)=HD1(m);
            end
        end
        [c,index]=max(HD2);
        phaiP=index*0.0001;
        Fdirection=0;
        for m=1:length(HD2)
            if HD2(m)~=0
                Fdirection=Fdirection+(phai(m)-phaiP)^2*HD2(m);
            end
        end 
        clear E0h;clear E0v;
        clear E1h;clear E1v;
        clear E2h;clear E2v;
        clear E3h;clear E3v;
        clear E4h;clear E4v;
        clear E5h;clear E5v;
        clear G;clear HD1;clear HD2;
        clear PI;clear Sbest;clear counts;
        clear deltaG;clear deltaH;clear deltaV;
        clear phaiP;
        clear theta;clear theta1;
        out=[Fcoarseness,Fcontrast,Fdirection];

        
        


        
        
    otherwise
        error(['Unknown Algorithms .']);  
        
end