function y=mll(m,beta,K,neigh)
%         mll(n,beta)
% Generates a sample a multi level logistic  distribution 
% Uses 1st or 2nd order neighborhoods
%
% Method: Gibbs sampler
% Input parameters ---------------------------------------
% m -> image size
% b -> MLL parameter
% K -> number of labels
% neigh -> neighborhood {1,2}

% Author: Jose M. Bioucas-Dias, 2005

 n= m+2;
 y=zeros(n);
 % histogram of the neighbors
 h = zeros(K,1);
  if neigh == 1
    % 1 st order neighborhood
    index=[0 1; -1 0; 0 -1; 1 0];
 else
    % 2nd order neighborhood
    index=[0 1; -1 1; -1 0; -1 -1; 0 -1; 1 -1; 1 0; 1 1];
 end

 iter = 100;
 % start with a uniform distribution and set boundary
 % to zeros
 y(2:n-1,2:n-1) = floor(rand(m)*K)+1;
 % image sweep randonmly
 for t=1:iter/2
     for i=2:n-1
        for j=2:n-1
           i= floor(rand()*m)+2;
           j= floor(rand()*m)+2;
           % comute neighbor histogram 
           h=0*h;
           for v =1:size(index,1)
               label = y(i+index(v,1),j+index(v,2));
               % if not on the boudary
               if label ~= 0
                  h(label) = h(label)+1;
               end
           end
           % compute conditional probabilities
           p = exp(beta*h);
           p=p/sum(p);
           % generate label
           cp = cumsum(p);
           cp = [0; cp];
           w = rand();
           for u=1:K
              if (cp(u)<= w) &  (cp(u+1)>= w)
                  y(i,j) = u;
                  break;
              end
           end
               
        end
     end
 end
 % sweep image linewise
  for t=1:iter/2
     for i=2:n-1
        for j=2:n-1
           % comute neighbor histogram 
           h=0*h;
           for v =1:size(index,1)
               label = y(i+index(v,1),j+index(v,2));
               % if not on the boudary
               if label ~= 0
                  h(label) = h(label)+1;
               end
           end
           % compute conditional probabilities
           p = exp(beta*h);
           p=p/sum(p);
           % generate label
           cp = cumsum(p);
           cp = [0; cp];
           w = rand();
           for u=1:K
              if (cp(u)<= w) &  (cp(u+1)>= w)
                  y(i,j) = u;
                  break;
              end
           end
               
        end
     end
 end

 
 y=y(2:n-1,2:n-1);
