function [CC,prior] = MergeCluster(R,prior,option)

% if option==1, skip the first class
if option==1
  Kstart = 2;
else
  Kstart =1;
end

% If number of clusters>1, merge two clusters
K = length(prior);
if (((K>1)&(option~=1))|((K>2)&(option==1)))

  % Determine closest clusters
  maxinvDist = 0;
  for k1=Kstart:K
  for k2=(k1+1):K
    R1 = squeeze(R(k1,:,:));
    R2 = squeeze(R(k2,:,:));
    pi1 = prior(k1);
    pi2 = prior(k2);
    invDist = 1/(ClusterDistance(R1,R2,pi1,pi2)+1E-30);
  
    if invDist>maxinvDist
      maxinvDist = invDist;
      k1min = k1;
      k2min = k2;
    end
  end
  end

  % compute parameters of new class 
  R(k1min,:,:) = R(k1min,:,:) + R(k2min,:,:);
  prior(k1min) = prior(k1min) + prior(k2min);
  
  % Delete class k2min using the fact that k2min>k1min>=1
  if k2min==K
    sample = 1:(k2min-1);
  else
    sample = [1:(k2min-1),(k2min+1):K];
  end
  R = R(sample,:,:);  
  prior = prior(sample);  

end

% Compute clustered component for each class
for k = 1:length(prior)
  [U,S,V] = svd(squeeze(R(k,:,:)));
  
  % Use principal component
  CC(:,k) = U(:,1);
end

% Check for option 1
if option==1
  CC(:,1) = 0*CC(:,1);
end



