function [ F ] = calcHaralick( glcm )

% HARALICK Fast Calculation of Haralick Features

%   IN:   glcm = Co-Occurrence Matrix     

%   OUT:  F = Feature Vector   

%

%   Stefan Winzeck 2012   

%   winzeck@hm.edu

% 

%   Feature Calculation according to:

%   [1] R. Haralick: 'Textural Feature for Image Classification' (1979)

%   [2] E. Miyamoto: 'Fast Calculation of Haralick Texture Features' 

% 

% MISSING:   f14  [1]

%% ALLOCATION

S=size(glcm,1);

f_2=zeros(S);

f_3=zeros(S);

f_4=zeros(S);

f_5=zeros(S);

f_6=zeros(1,2*S);

f_7=zeros(1,2*S);

f_8=zeros(1,2*S);

f_9=zeros(S);

f_11=zeros(1,S);

pxy=zeros(1,2*S);

px_y=zeros(1,S);

HX_=zeros(1,S);

HY_=zeros(1,S);

HXY_1=zeros(S);

HXY_2=zeros(S);

%% CALCULATION

% Normalization

M = glcm/sum(glcm(:));

% Energy

f_1 = M.^2;

f1 = sum(f_1(:));

%-------------------------------------------------------------------------%

u = mean2(M);

py = sum(M,1);

px = sum(M,2);

%OPTIMIZATION NEEDED: change the indicies to j = i and multiple the 

% end summations by 2.

for i=1:S

    for j=1:S

        Mij = M(i,j);

        f_3(i,j) = i*j*Mij;

        f_4(i,j) = (i-u)^2*Mij;

        f_5(i,j) = Mij/(1+(i-j)^2);

        %OPTIMIZATION NEEDED: Use of log tables

        f_9(i,j) = Mij*log(Mij+eps);

        %Equation (5) from [2]

        pxy(i+j) = pxy(i+j)+Mij;

        %Equation (6) from [2]

        px_y(abs(i-j)+1) = px_y(abs(i-j)+1)+Mij;

        %Different than the paper

        %OPTIMIZATION NEEDED: Use of log tables

        %Related to Equation (20) from [2]

        HX_(i)= px(i)*log(px(i)+eps);

        %Related to Equation (20) from [2]

        HY_(j)= py(j)*log(py(j)+eps);

        %Equation (7) from [2]

        HXY_1(i,j) = Mij*log(px(i)*py(j)+eps);

        %Equation (8) from [2]

        HXY_2(i,j) = px(i)*py(j)*log(px(i)*py(j)+eps);

    end

end

% Correlation

ux = mean(px); sx=std(px);

uy = mean(py); sy=std(py);

f3 =(sum(f_3(:))-(ux*uy))/(sx*sy);

% Sum of Variances

f4 = sum(f_4(:));

% Inverse Difference Moment

f5 = sum(f_5(:));

% Entropy

f9 = -sum(f_9(:));

% Information Measures of Correlation 1&2

HX = -sum(HX_);

HY = -sum(HY_);

HXY = f9;

HXY1 = -sum(HXY_1(:));

HXY2 = -sum(HXY_2(:));

f12 = (HXY-HXY1)/max([HX, HY]);

f13 = (1 - exp((-2)*(HXY2 - HXY)))^0.5;

%-------------------------------------------------------------------------%

for i=2:2*S

    f_6(i) = i*pxy(i);

    %OPTIMIZATION NEEDED: Use of log tables

    f_8(i) = pxy(i)*log(pxy(i)+eps);

end

% Sum Average

%f_6(1) = [];       % not necessary f_6(1) is zero anyway

f6 = sum(f_6);

% Sum Entropy

%f_8(1)=[];         % not necessary f_8(1) is zero anyway

f8 = -sum(f_8);

%-------------------------------------------------------------------------%

%Different than the paper, f7 should be calculated along side f6 and f8

for i=2:2*S

    f_7(i)=(i-f8)^2*pxy(i);

end

% Sum Variance

%f_7(1) = [];       % not necessary f_7(1) is zero anyway

f7 = sum(f_7);

% Difference Variance

f10 = var(px_y);

%-------------------------------------------------------------------------%

for k=1:S

    f_2(k) = (k-1)^2*px_y(k);

    %OPTIMIZATION NEEDED: Use of log tables

    f_11(k) = px_y(k)*log(px_y(k)+eps);

end

% Contrast

f2 = sum(f_2(:));

% Difference Entropy

f11 = -sum(f_11);

%-------------------------------------------------------------------------%

F = [f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11;f12;f13];

end