Calculating first-order histogram based features

I'm trying to calculate first order histogram features in matlab. It looks like I'm missing something, because current implementation of FOF doesn't work correctly - it should classify correctly like 80% of textures, but current best result is 20%.

Database - Brodatz, 576 images (180x180); each texture has 6 images; total number of classes - 96.

I believe that classify algorithm is working fine, because more complex algorithms are returning expected results:

• Co-occurrence matrix ~ 94%,
• LBP - 99,1%.

Could you please take a look at my implementation? First order features are explained here, pages: 7-9.

Matlab Implementation:

    function [ result ] = GetFOF( img1 )
    %GetFOF img1 must be two-dimensional.

    [pixelCounts, GLs] = imhist(img1);
    NM = sum(pixelCounts); % number of pixels

    % Calculate the mean
    meanGL = sum(GLs .* (pixelCounts / NM));

    varresult = 0;  % variance temp var
    skewresult = 0; % skewness temp var
    kurtresult = 0; % kurtosis temp var

    for i=0:1:length(pixelCounts)-1
        varresult = varresult + (i-meanGL)^2 * (pixelCounts(i+1)/NM);
        skewresult = skewresult + (i-meanGL)^3 * (pixelCounts(i+1)/NM);
        kurtresult = kurtresult + (i-meanGL)^4 * (pixelCounts(i+1)/NM)-3;
    end

    skewresult = skewresult * varresult ^-3; % skewness
    kurtresult = kurtresult * varresult ^-4; % kurtosis

    %energy
    energy = sum((pixelCounts / NM) .^ 2);

    %entropy 
    pI = pixelCounts / NM;
    entropy1 = -sum(pI(pI~=0) .* log2(pI(pI~=0)));

    % returns the same result as above
    %entropy2 = -sum(pI .* log2(pI+eps))
    %entropy3 = entropy(img1)

    result = [meanGL, varresult, skewresult, kurtresult, energy, entropy1];

    % calculate features using the matlab-optimized way; currently broken 
    % (different results)

    % Subtract the mean from the gray levels
    % m0 = pixelCounts - meanGL;
    % diffSquared = m0 .^ 2;
    % % variance
    % variance = sum(diffSquared .* (pixelCounts / NM))
    % 
    % %skewness
    % diffThird = m0 .^ 3;
    % skewness = variance ^ -3 * sum(diffThird .* (pixelCounts / NM))
    % 
    % %kurtosis
    %  diffFour = m0 .^ 4;
    %  kurtosis = variance ^ -4 * sum(diffFour .* ((pixelCounts / NM) - 3))
    end

I have tried to calculate features manually for sample image:

    im = [5 0 1; 4 2 3; 0 5 2] (gray levels = 0..5 = 6)

And I have received almost the same values, with the exception of kurtosis. (If we change loop to 0..6, kurtosis will be the same too).

Edit

Classification details: For each texture there are 6 images. With external script (not matlab code) I have created text file with following structure (Image path,class Id):

    ...\group000_1.bmp,1
    ...\group001_1.bmp,2
    ...\group002_1.bmp,3
    ...\group003_1.bmp,4
      (...)
    ...\group000_2.bmp,1
    ...\group001_2.bmp,2
    ...\group002_2.bmp,3
    ...\group003_2.bmp,4
      (...)

Following structure and order should ensure that data will be still valid after splitting it, in into two parts. By 'valid' I mean: 3 images from each texture goes into set 1, and the rest 3 goes into set 2. (Yes,the total number of images is even)

After this, load all images, process them with specified algorithm. Algorithm must return vector of features/data. Split data & group information into two sets. Implementation:

    [paths, groups] = LoadPathGroupPairFile();

    % temp buffer for all feature vectors
    temp = zeros(length(paths), 256);

    % size of each set
    maxlength = floor(length(temp(:,1))/2);

    validSampleClasses = groups((1:maxlength), :);
    validTestClasses = groups((maxlength+1:length(groups)), :);

    % check if grouping script outputed valid data (each set has all classes)
    if(length(validSampleClasses) ~= length(validTestClasses) || ...
       isempty(setdiff(validSampleClasses, validTestClasses)) == 0 || ... 
       isempty(setdiff(validTestClasses, validSampleClasses)) == 0)
        error('invalid classes');
    end

    % load images, process them with specified algorithm; 
    % put feature vector into temp(i,:) = ...
    % (...)
    % f.e.: temp(i,:) => feature vector for 'i' image.

    % split feature vector buffer into two sets
    sampleValues = temp((1:maxlength), :);
    testValues = temp((maxlength+1:length(temp(:,1))), :);

    k = 1;
    resultknn = knnclassify(sampleValues, testValues, validTestClasses, k);
      
    invalidvalues = 0;
    for i=1:1:length(resultknn)
        if(resultknn(i) ~= validSampleClasses(i))
        invalidvalues = invalidvalues + 1;
      end
    end

    classfresult = ((maxlength-invalidvalues)*100/maxlength);
    display(sprintf('result %0.2f %%', classfresult));

[dpb]

Seems to implement the paper definition altho seems a little peculiar. Depending on the image and the resulting distribution, two things can happen--

a) If, as the paper mentions, the histogram is not unimodal, the central moments don't help a lot as the moosh everything together and don't reflect the various modes at all, and

b) Depending on the granularity the weighted estimates don't necessarily reproduce the values computed directly from the observations terribly closely.

Those may be issues in your classification; would need to look in the details to ascertain.

One thing in the computation; you can replace all the loops

varresult = 0;  % variance temp var
skewresult = 0; % skewness temp var
kurtresult = 0; % kurtosis temp var

for i=0:1:length(pixelCounts)-1
  varresult = varresult + (i-meanGL)^2 * (pixelCounts(i+1)/NM);
  skewresult = skewresult + (i-meanGL)^3 * (pixelCounts(i+1)/NM);
  kurtresult = kurtresult + (i-meanGL)^4 * (pixelCounts(i+1)/NM)-3;
end

with vectorized versions as

    L=length(pixelCounts);
    ii=0:L-1;
    vr=sum((ii-meanGL).^2.*pixelCounts/NM);
    sk=sum((ii-meanGL).^3.*pixelCounts/NM);
    ku=sum((ii-meanGL).^4.*pixelCounts/NM)-3*L;

[Pawel]

@dpb

Thanks for reply and vectorized version of loops. Added classification details to question body.