function [I,K,Xdim,Ydim,dl_vec,Kvect,restart]=em_mdl(A,Xdim)

restarts=4;
Kvect=[2 3 4 5];
MDL=30;
%% was 40 for /elib/data/images/blobs/new_blobs/Until980911/ (9/11/98)
%% was 50 until SpatialInit2_MDL40 (9/1/98)

%%%% initialize random number generator
rand('seed',sum(100*clock));

if Xdim==192
   Ydim=128;
elseif Xdim==128
   Ydim=192;
else
   disp(Xdim)
   error('Xdim is not 128 or 192!')
end
Xdim2=Xdim-16;
Ydim2=Ydim-16;

num_attr=size(A,2);
num_ex=size(A,1);
max_num_iter=15;
ll_threshold=0.01;

ll_vec=zeros(1,4);
Ik=cell(length(Kvect));
I=zeros(Ydim2,Xdim2);

for knum=1:length(Kvect)
   K=Kvect(knum);
   Irestarts=cell(restarts);
   
   for r=1:restarts
      [x,y]=meshgrid(1:Xdim,1:Ydim);
      x=trimedge(x,8);
      x=x(:)/Xdim;
      y=trimedge(y,8);
      y=y(:)/Ydim;
      centervect=(x>1/4 & x<3/4 & y>1/4 & y<3/4);
      outervect=~centervect;
      center=find(centervect);

      if K==2
	 means(1,:)=mean(A(center,:),1);
	 means(2,:)=mean(A(find(outervect),:),1);
      elseif K==3
	 means(1,:)=mean(A(center,:),1);
	 if floor(r/2)==(r/2)
	    top= (y<1/2 & outervect);
	    bottom= (~top & outervect);
	    means(2,:)=mean(A(find(top),:),1);
	    means(3,:)=mean(A(find(bottom),:),1);
	 else
	    left= (x<1/2 & outervect);
	    right= (~left & outervect);
	    means(2,:)=mean(A(find(left),:),1);
	    means(3,:)=mean(A(find(center),:),1);
	 end
      elseif K==5
	 means(1,:)=mean(A(center,:),1);
	 top= (y<1/2 & outervect);
	 bottom= (~top & outervect);
	 left= (x<1/2 & outervect);
	 right= (~left & outervect);
	 means(2,:)=mean(A(find(top&left),:),1);
	 means(3,:)=mean(A(find(top&right),:),1);
	 means(4,:)=mean(A(find(bottom&left),:),1);
	 means(5,:)=mean(A(find(bottom&right),:),1);
      elseif (K==4 | K>8)
	 W=sqrt(K);
	 for xi=1:W
	    for yi=1:W
	       use= (x>=(xi-1)/W)&(x<xi/W)&(y>=(yi-1)/W)&(y<yi/W);
	       means((xi-1)*W+yi,:)=mean(A(find(use),:),1);
	    end
	 end
      end
      if (K~=3 & r>1) | (K==3 & r>2)
	 means=means+1.5*randn(size(means)); %% was 1.0 for BigRun2
      end
      
      for k=1:K
	 tmp=0.4*eye(num_attr);
	 if K>2
	    tmp(end-1,end-1)=2;
            tmp(end,end)=2;
         else
            tmp(end-1,end-1)=10;
            tmp(end,end)=10;
         end
	 covars{k}=tmp;
      end
      
      mix_weights=ones(1,K)/K;
      
      % Compute how well we're doing initially using log likelihood
      ccd=class_cond_density_mcm(A,means,covars,mix_weights);
      ll_restarts(r) = sum(log(ccd+eps));
      
      for iter=1:max_num_iter
	 for k=1:K
	    % update mix_weights
	    pmc=prob_mix_comp_mcm(k,A,means,covars,mix_weights,ccd);
	    mix_weights_new = sum(pmc) / num_ex;  %% mean(pmc)
	    
	    % update means
	    m=pmc'*A;
	    denom = mix_weights_new * num_ex;
	    means_new = m(:) / (denom+eps);
	    
	    % update covars
	    A_centered=A'-means_new*ones(1,num_ex);
	    numer=(A_centered.*(ones(num_attr,1)*pmc'))*A_centered';
	    denom = mix_weights_new * num_ex;
	    covars_new = numer / (denom+eps);
	    
	    mix_weights(k) = mix_weights_new;
	    means(k,:) = means_new';
	    covars{k} = covars_new;
	 end
	 
	 % Compute how well we're doing using log likelihood
	 ccd=class_cond_density_mcm(A,means,covars,mix_weights);
	 ll_new = sum(log(ccd+eps));
	 
	 % Decide whether to stop
	 if abs((ll_new - ll_restarts(r)) / ll_restarts(r)) < ll_threshold
	    ll_restarts(r) = ll_new;
	    break
	 end
	 ll_restarts(r) = ll_new;
      end  % for iter
      evaluate_support; % creates Irestarts{r}
   end  % for r (restarts)
   
   [ll,ind]=max(ll_restarts);
   Ik{knum}=Irestarts{ind};
   best(knum)=ind;
   
   ll_vec(knum)=ll;
end

% apply MDL principle
d=size(A,2);
pen = MDL * log(Xdim*Ydim) * ((Kvect-1) + Kvect*d + Kvect*d*(d-1)/2) / 2;
    % hack * log(# data pts)* #parameters for full covariance case
dl_vec=ll_vec-pen;

[junk,k]=max(dl_vec);
K=Kvect(k);
restart=best(k);
I=Ik{k};

%%save(savefile,'I','K','restart','Xdim','Ydim','dl_vec','Kvect');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 1997-1999. The Regents of the University
% of California (Regents). All Rights Reserved.
% 
% Permission to use, copy, modify, and distribute this software and its
% documentation for educational, research, and not-for-profit purposes,
% without fee and without a signed licensing agreement, is hereby
% granted, provided that the above copyright notice, this paragraph and
% the following two paragraphs appear in all copies, modifications, and
% distributions. Contact The Office of Technology Licensing, UC Berkeley,
% 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, (510)
% 643-7201, for commercial licensing opportunities. Created by
% Chad Carson, Serge Belongie, and Jitendra Malik, Department of Electrical
% Engineering and Computer Sciences, University of California, Berkeley.
% For technical information, contact carson@eecs.berkeley.edu.
% 
% IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT,
% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS,
% ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
% REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
% 
% REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
% PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF
% ANY, PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION
% TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
% MODIFICATIONS.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

