%%
function resultTot = slbp(varargin) 


narginchk(1,7);

image = varargin{1};
d_image = double(image);

if nargin==1
    spoints = [-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];
    neighbours = 8;
    mapping = 0;
    threshold = 3;
    steps = threshold*2+1;
    mode = 'h';
end

if (nargin == 2) && (length(varargin{2}) == 1)
    error('Input arguments');
end

if (nargin >= 5) && (length(varargin{2}) == 1)
    radius = varargin{2};
    neighbors = varargin{3};
    
    steps = varargin{4};
    threshold = varargin{5};
     
    spoints = zeros(neighbors,2);
    
    % Angle step.
    a = 2*pi/neighbors;
    
    for i = 1:neighbors
        spoints(i,1) = -radius*sin((i-1)*a);
        spoints(i,2) = radius*cos((i-1)*a);
    end
    
    if(nargin >= 6)
        mapping=varargin{6};
        if(isstruct(mapping) && mapping.samples ~= neighbors)
            error('Incompatible mapping');
        end
    else
        mapping=0;
    end
    
    if(nargin >= 7)
        mode=varargin{7};
    else
        mode='h';
    end
end

% Determine the dimensions of the input image.
[ysize xsize] = size(image);

miny=min(spoints(:,1));
maxy=max(spoints(:,1));
minx=min(spoints(:,2));
maxx=max(spoints(:,2));

% Block size, each LBP code is computed within a block of size bsizey*bsizex
bsizey=ceil(max(maxy,0))-floor(min(miny,0))+1;
bsizex=ceil(max(maxx,0))-floor(min(minx,0))+1;

% Coordinates of origin (0,0) in the block
origy=1-floor(min(miny,0));
origx=1-floor(min(minx,0));

% Minimum allowed size for the input image 2depends
% on the radius of the used LBP operator.
if(xsize < bsizex || ysize < bsizey)
    error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');
end

% Calculate dx and dy;
dx = xsize - bsizex;
dy = ysize - bsizey;

% Fill the center pixel matrix C.
d_C = double(image(origy:origy+dy,origx:origx+dx));

% compute interval based on the threshold value and number of steps, steps
interval = round(linspace(-threshold,threshold,steps));

% Initialize the result matrix with zeros.
result = zeros(length(interval),dy+1,dx+1);

bins = 2^neighbors;

%     result(k,:,:)=zeros(dy+1,dx+1);


%Compute the LBP code image
for i = 1:neighbors
    y = spoints(i,1)+origy;
    x = spoints(i,2)+origx;
    % Calculate floors, ceils and rounds for the x and y.
    fy = floor(y); cy = ceil(y); ry = round(y);
    fx = floor(x); cx = ceil(x); rx = round(x);
    % Check if interpolation is needed.
    if (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6)
        % Interpolation is not needed, use original datatypes
        N = image(ry:ry+dy,rx:rx+dx);
        %% EDIT: displacement added
        %             D = N >= d_C;
        %             D = (N >= (C + interval(k)));
    else
        % Interpolation needed, use double type images
        ty = y - fy;
        tx = x - fx;
        
        % Calculate the interpolation weights.
        w1 = (1 - tx) * (1 - ty);
        w2 =      tx  * (1 - ty);
        w3 = (1 - tx) *      ty ;
        w4 =      tx  *      ty ;
        % Compute interpolated pixel values
        N = w1*d_image(fy:fy+dy,fx:fx+dx) + w2*d_image(fy:fy+dy,cx:cx+dx) + ...
            w3*d_image(cy:cy+dy,fx:fx+dx) + w4*d_image(cy:cy+dy,cx:cx+dx);
        
    end
    % Update the result matrix.
    v = 2^(i-1);
    for k = 1:length(interval)
        D = N >= d_C+interval(k);
        result(k,:,:) = squeeze(result(k,:,:)) + v*D;
    end
end


%% Apply mapping if it is defined
if isstruct(mapping)
    for n = 1:length(interval)
        bins = mapping.num;
        for i = 1:size(result,2)
            for j = 1:size(result,3)
                result(n,i,j) = mapping.table(result(n,i,j)+1);
            end
        end
    end
end

if (strcmp(mode,'h') || strcmp(mode,'hist') || strcmp(mode,'nh'))
    % Return with LBP histogram if mode equals 'hist'.
    
    resultTot = zeros(1,bins);
    for k=1:size(result,1)
        res = result(k,:,:);
        resultTot = resultTot + hist(res(:),0:(bins-1));
    end
    resultTot = resultTot./(length(interval));
    
    if (strcmp(mode,'nh'))
        resultTot=resultTot./sum(resultTot);
    end
else
    %Otherwise return a matrix of unsigned integers
    if ((bins-1)<=intmax('uint8'))
        result=uint8(result);
    elseif ((bins-1)<=intmax('uint16'))
        result=uint16(result);
    else
        result=uint32(result);
    end
end

end




