I am struggling to understand Matlab implementation of LBP algorithm found here. I am trying to find how it calculates the binaries for every pixel? It just calculates where the neighbor pixel are greater than the actual center pixel size. I want to calculate the binaries for every pixel in order to use local histograms to calculate the features of 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 depends
% 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.
C = image(origy:origy+dy,origx:origx+dx);
d_C = double(C);

bins = 2^neighbors;

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

%Compute the LBP code image
% the whole process here
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);
    D = N >= C; 
  else
  % Interpolation needed, use double type images 
  ty = y - fy;
  tx = x - fx;

  % Calculate the interpolation weights.
  w1 = roundn((1 - tx) * (1 - ty),-6);
  w2 = roundn(tx * (1 - ty),-6);
  w3 = roundn((1 - tx) * ty,-6) ;
  % w4 = roundn(tx * ty,-6) ;
  w4 = roundn(1 - w1 - w2 - w3, -6);

  % 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);
  N = roundn(N,-4);
  D = N >= d_C; 
 end  
   % Update the result matrix.
  v = 2^(i-1);
  result = result + v*D;
end

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

 if (strcmp(mode,'h') || strcmp(mode,'hist') || strcmp(mode,'nh'))
  % Return with LBP histogram if mode equals 'hist'.
  result=hist(result(:),0:(bins-1));
  if (strcmp(mode,'nh'))
    result=result/sum(result);
  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
 size(result)
end

Iteratively it adds some value in results for all 8 neighbors of every pixel. But how it is correlated with LBP binaries? How is it correlate with the following code for the following c++ LBP approach: 

uchar lbp(const Mat_ & img, int x, int y)
 {
  // this is pretty much the same what you already got..
  uchar v = 0;
  uchar c = img(y,x);
  v += (img(y-1,x  ) > c) << 0;
  v += (img(y-1,x+1) > c) << 1;
  v += (img(y  ,x+1) > c) << 2;
  v += (img(y+1,x+1) > c) << 3;
  v += (img(y+1,x  ) > c) << 4;
  v += (img(y+1,x-1) > c) << 5;
  v += (img(y  ,x-1) > c) << 6;
  v += (img(y-1,x-1) > c) << 7;
  return v;

It's a vectorized implementation of LBP, rather well-suited for Matlab.

After the initialization instructions, let's look the main loop, beginning at the line "for i = 1:neighbors". The loop is pretty clear: it computes the comparison of one neighbor with the center pixel, the loop iterates over all neighbors. You've got this point, so now enter deep into the loop to understand how it accumulates all results.

The core of the loop is in fact over complicated because it takes into account the real circle instead of an approximate integer circle. So the purpose of the major part of the instructions is to compute the interpolated intensity of the neighbor pixel. Here it differs from the C++ code you have as reference, where it takes only the integer, 1-pixel-wide-radius circle. Remember that with the lbp.m code you can -- theoretically, I will discuss that later -- compute the LBP along a circle of radius R with N sampling points, so the C++ would correspond to a circle of radius 1 and with 8 sampling points, if only there was no interpolation. But there is an interpolation when the neighbor does not fit the pixel grid of the image, when (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6) is false).

If (abs(x - rx) < 1e-6) && (abs(y - ry) < 1e-6) is true, there is no interpolation, so the computation of all comparisons between the central pixel and the current neighbor is stored directly into D. Else, it computes a bilinear interpolation of the intensity at the sampling neighbor point, over the entire image: 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);.

And finally, turn to the update part: v = 2^(i-1); result = result + v*D;. v is the equivalent of the shift: for the ith neighbor, you shift the value of the comparison by i-1 to the left, or equivalently multiplying be 2^(i-1). Then you sum with result. So at the end of the loop, the computation is really equivalent to your C++ code, except that it's done over the entire image instead of one single pixel. And the C++ code can be seen as a unrolled version of the matlab loop with the neighbor circle of radius 1 and 8 sampling points. At this point, the LBP map is computed, the following blocks are additional processing of the LBP map (remap through a mapping table, and optionally computing the histogram of the LBP image instead of the LBP image itself).

Now, a little discussion about the whole script. There is a flaw here that is hidden at the end of the script. In fact, through the code, you are limited to 32 neighbors, no more, because at the end the LBP image is cast to int32. The flaw is that the variable result is allocated as a double matrix and not an integer matrix, so I really hope that there is no approximation problem when updating result and later when casting into integer, leading to changing bits in the LBP. Normally there should not be as there is at least 52 precision bits (according to wikipedia for IEEE 754 spec). I think it's risky here ... and on the contrary I am not aware of a matlab type for long fixed-sized, efficient bit vector. I would use int64 instead of int32, but the limit will be there at 64 sampling neighbors.

EDIT

Now, if your wish is to commpute some local binary patterns restricted on the 3*3 neighborhood, this Matlab function is way too generic for you, and the best thing is to unroll the loop for this neighborhood, and thus be really close to the C++ code. Here is a piece of code for that (I use bitwise or instead of addition, but it's equivalent): 

result = uint8(ysize, xsize);
result = (image(1:end-2,2:end-1) > image(2:end-1,2:end-1));                                 % <=> v += (img(y-1,x  ) > c) << 0;
result = result|bitshift((image(1:end-2,3:end) > image(2:end-1,2:end-1)), 1, 'uint8');      % <=> v += (img(y-1,x+1) > c) << 1;
result = result|bitshift((image(2:end-1,3:end) > image(2:end-1,2:end-1)), 2, 'uint8');      % <=> v += (img(y  ,x+1) > c) << 2;
result = result|bitshift((image(3:end,3:end) > image(2:end-1,2:end-1)), 3, 'uint8');        % <=> v += (img(y+1,x+1) > c) << 3;
result = result|bitshift((image(3:end,2:end-1) > image(2:end-1,2:end-1)), 4, 'uint8');      % <=> v += (img(y+1,x  ) > c) << 4;
result = result|bitshift((image(3:end,1:end-2) > image(2:end-1,2:end-1)), 5, 'uint8');      % <=> v += (img(y+1,x-1) > c) << 5;
result = result|bitshift((image(2:end-1,3:end) > image(2:end-1,2:end-1)), 6, 'uint8');      % <=> v += (img(y  ,x-1) > c) << 6;
result = result|bitshift((image(1:end-2,1:end-2) > image(2:end-1,2:end-1)), 7, 'uint8');    % <=> v += (img(y-1,x-1) > c) << 7;

It's the exact translation of the C code to a Matlab script, using the powerful vectorization. With this in hand, it's pretty simple to change for another order or different tests in this neighborhood. I also mention this point because there is an error in the Matlab script for this case, line 53 there is a wrong sign: neighobrhood is better asspoints=[-1 -1; -1 0; -1 1; 0 -1; 0 -1; 1 -1; 1 0; 1 1]; instead of spoints=[-1 -1; -1 0; -1 1; 0 -1; -0 1; 1 -1; 1 0; 1 1];.