differentialFeature( img, features, stepSize )

function [ fVal, idx ] = differentialFeature( img, features, stepSize )
% differentialFeature method providing a 3D differential feature.
%   
%   Differential feature is computed in LBP manner without computing
%   decimal reduction.
%   
%   Input:
%       img = 3 dimensional array containing intensity values. (volume)
%       features = vector containing the displacement vectors used as
%       neighborhood.
%       stepSize = step size used for sampling
%
%   Output:
%       fVal = binary feature matrix.
%       idx = linear index of the binary feature response.
%
%   Example:
%       
%       img = load('...');
%       feature = randn(3, 100) * 10; % 100 neighbours normal distributed
%       stepSize = 2;
%
%       result = differentialFeature(img, feature, stepSize);

%% sanity check
assert ( ndims( img ) == 3 )

%% initialize
[X, Y, Z] = size(img);

%% compute

% center voxel intensity matrix
C = img(1:stepSize:X, 1:stepSize:Y, 1:stepSize:Z);
fVal = false(size(features, 2), numel(C));

for n = 1:size(features, 2)
   
    N = zeros(size(C));
    
    x = max(1, features(1,n));
    xx = min(X, X - features(1,n));
    
    y = max(1, features(2,n));
    yy = min(Y, Y - features(2,n));
    
    z = max(1, features(3,n));
    zz = min(Z, Z - features(3,n));
    
    % neighboor n voxel intensity matrix
    A = img(x:stepSize:xx, y:stepSize:yy, z:stepSize:zz);
    
    [aX, aY, aZ] = size(A);
         
    startX =  numel(1:stepSize:x);
    endX = aX + startX - 1;
     
    startY =  numel(1:stepSize:y);
    endY = aY + startY - 1;
     
    startZ =  numel(1:stepSize:z);
    endZ = aZ + startZ - 1;

    N(startX:endX, startY:endY, startZ:endZ) = A;
    
    % compare values
    B = N > C;
    fVal(n,:) = B(:);
    
    if (n == (size(features, 2) / 2))
        disp('50%')
    end
end

%% add linear index
[P, Q, R] = ndgrid(1:stepSize:X, 1:stepSize:Y, 1:stepSize:Z);
idx = sub2ind(size(img), P, Q, R);
idx = single(idx(:));
end