%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                        %
%  A solution to CSE For Your Homework Project 3         %
%  Classified Information: The Data Clustering Problem   %
%  Nargess Memarsadeghi and Dianne P. O'Leary  03/03     %
%  file name = problem3_and_4.m                          %
%                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Load the data and set convergence parameters.         %
%  The data is an rgb image of dimension 500x500,        %
%  so the dimensions of the array are m=500, p=500,      %
%  q = 3.                                                %
%  We will try to cluster the q-vectors into k clusters, %
%  for k=3,4,5.  Then we construct the clustered image   %
%  and compute the relative distance between it and the  %
%  original image.                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global myast myk myell mywhichfunction

load charlie
asdata = double(as);

myell = 2;

[m,p,q] = size(asdata);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Set the range of values to consider in defining       %
%  the clusters.  We restrict these values to columns    %
%  210 and 211 of the image, giving 1000 values in all.  %
%  This is done to shorten the computation time.         %
%  The values are stored in array myast.                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

myast = asdata(:,210:211,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Try function minimization with R and D (Problem 3)    %
%                and k-means (Problem 4)                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for fnct = 1:3,

  if (fnct == 1) 
    mywhichfunction = 'minimizing R';
  elseif (fnct == 2)
    mywhichfunction = 'minimizing D';
  else
    mywhichfunction = 'using k-Means';
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Set up the figure to display results.                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  figure
  subplot(2,2,1)
  image(asdata/255)
  axis image
  axis off
  set(gcf, 'Renderer', 'opengl');
  title('Original Image')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Find the range of values in myast and choose          %
%  candidate cluster centers.                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     mins = min(min(myast,[],1),[],2);
     maxs = max(max(myast,[],1),[],2);

     for qq = 1:q,
       lb(qq,1) = mins(1,1,qq);
       ub(qq,1) = maxs(1,1,qq);
     end

centersav = zeros(q,5);
centersav(q,:) = linspace(lb(q),ub(q),5);
centersav(1,:) = linspace(lb(1),ub(1),5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Try 3, 4, or 5 clusters.                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  for k=3:5
    myk = k;

    disp('')
    disp(sprintf('k=%d %s', ...
                  myk,mywhichfunction))

    ss = sprintf('%d clusters %s',myk,mywhichfunction);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Initialize the cluster centers.                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     centers = centersav(:,1:k);
     
     cs = reshape(centers,myk*q,1);
  
     if (fnct < 3)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Find the  k clusters using minimization.              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

         [cs,fval,exitflag,output] = fminunc(@functr,cs);
         disp(sprintf('Exitflag = %d',exitflag))
         centers = reshape(cs,q,myk);
         [aclus,clustercounts] = map_to_cluster(asdata,centers);
         iters = output.iterations;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Find the  k clusters using k-means.                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     else

         [aclus,centers,clustercounts,iters] = ...
                 mycluster(myast,centers,1,200);
     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Print summary information.                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     clusterradius = computeradius(myast,centers,aclus);
     disp(ss)
     disp(sprintf('Number of iterations = %d',iters))
     disp('')
     disp('    Cluster centers         Counts   Radii')
     disp(sprintf('%7.2f %7.2f %7.2f %10d %7.2f \n', ...
                [centers',clustercounts',clusterradius']'))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Construct the clustered image, measure its error,     %
%  and display the image.                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     [aclus,clustercounts] = map_to_cluster(asdata,centers);
  
     asclus = zeros(m,p,q);
  
     for i=1:m,
        for j=1:p,
          asclus(i,j,:) = centers(:,aclus(i,j));
        end
     end
  
     disp('')
  
     n1 = 0;
     n2 = 0;
     for i=1:q,
       n1 = n1 + norm(asdata(:,:,i)-asclus(:,:,i),'fro')^2;
       n2 = n2 + norm(asdata(:,:,i)              ,'fro')^2;
     end
     disp(sprintf('Norm of relative change in image = %e',sqrt(n1/n2)))
   
     subplot(2,2,myk-1)
     image(asclus/255)
     axis image
     axis off
     set(gcf, 'Renderer', 'opengl');
     title(ss)
     drawnow

    end % for k
  end % for fnct