How to find SGLD matirx?

I got the following code from net

% getting the size of the input image
im_final=floor(double(image)/(2^output_bits));
[im_final_x im_final_y]=size(im_final);

% setting the size of the co-occurence matrices depending on the grey level depth
CO_size=2^input_bits/(2^output_bits);
SGLD=zeros(CO_size,CO_size);

switch theta
   case {0}
     for i=1:im_final_x
       for j=1:(im_final_y-d)
          SGLD(im_final(i,j)+1,im_final(i,j+d)+1)=SGLD(im_final(i,j)+1,im_final(i,j+d)+1)+1;
       end
     end

% make the SGLD matrix symmetric by adding it's transpose to it
SGLD=SGLD+SGLD';

% normalize the SGLD matrix to values between 0 and 1
SGLD=SGLD/sum(sum(SGLD));

But what is the logic of this program,ie how it actually works?