function [ClusterA,ClusterB]=txturgroup(Distance,m_x,r_xx,M_h,N_h,w,TestCrit,Separation,IsolateTex,Verbose,Hierarchical,MergeStep)
% TXTURGROUP  Make efficient groupings of textures
%	[ClusterA,ClusterB]=txturgroup(Distance,m_x,r_xx,M_h,N_h,w,...
%	   'MeanSep',1.2,IsolateTex,Verbose,Hierarchical,MergeStep) 
%	clusters the textures by the scheme outlined below. 
%	Verbose=0 means quiet, 1 and 2 means different 
%	levels of user feedback. Distance must be a table with joint distances
%	between textures (by any measure), where Distance(i,j) is the
%	distance between texture i and j. The mean values of each input 
%	texture are given by the vector m_x, and the autocorrelation functions
%	are stacked on top of each other (all of same size) in r_xx.
%	M_h and N_h are the desired filter sizes, while w is the smoothing
%	filter.
%	
%	If IsolateTex is given (and not empty), the texture in it is 
%	kept as one of the groups in all steps.
%	
%	Each line in ClusterA represent a cluster of textures which it
%	should be possible to separate from the corresponding line in
%	ClusterB by an optimal filter.
%	
%	CLUSTERING SCHEME:
%	 1) Start with two most distant textures.
%	 2) For each non-discriminated texture, add it to any of the two 
%	    groups if it increases discrimination and the maximum distance 
%	    within the groups remain less than the minimum distance between.
%	 3) Design filter between the groups according to the Unser criterion.
%	    With the 'MeanSep' option (only option implemented yet), and 
%	    Separation 1.2: If the least feature mean of the group having
%	    largest feature mean is less than 1.2*the maximum mean of the 
%	    other group, remove elements of one of the groups.
%	 4) Continue 3) until satisfying (or only two textures left).
%	 5) For all remaining ungrouped and undiscriminated textures
%	     a) Start with texture most distant from any of the groups.
%	     b) For each non-discriminated texture, add it to this group,
%	        if it increases discrimination and the maximum distance 
%	        within all groups remain less than the minimum distance 
%	        between.
%	 6) Try to add this group to the two others by designing filters
%	    and performing the test of 3). If the test is passed, add this
%	    group to the corresponding group.
%	 7) Repeat steps 5) and 6) until all textures have been examined.
%	 8) Mark off (in a table) the discrimination the resulting grouping
%	    yields, and add the grouping to the two output variables as new 
%	    lines.
%	 9) Go to step 1) until all textures are discriminated.
%	 *) TERMINOLOGY:
%	    Discrimination: if a filter is designed to separate (among others)
%		texture i from texture j, texture i is discriminated from
%		texture j.
%		
%	See also: TXTURSIM

N_T=length(m_x);
[M_r,N_r]=size(r_xx); M_r=M_r/N_T;

if exist('Verbose')~=1,      Verbose=[];      end
if isempty(Verbose),         Verbose=2;       end
if exist('Hierarchical')~=1, Hierarchical=[]; end
if isempty(Hierarchical),    Hierarchical=0;  end
if exist('MergeStep')~=1,    MergeStep=[]; end
if isempty(MergeStep),       MergeStep=1;  end
NotYetSeparated=ones(N_T,N_T)-diag(ones(N_T,1));
DoIsolate=0;
if exist('IsolateTex')==1
  if ~isempty(IsolateTex)
    cau=IsolateTex;
    DoIsolate=1;
    NotYetSeparated=zeros(N_T,N_T);
    NotYetSeparated(IsolateTex,:)=ones(1,N_T);
    NotYetSeparated(IsolateTex,IsolateTex)=0;
  end
end

ClusterA=[]; ClusterB=[];

if Verbose>=1 & ~Hierarchical
  fprintf(1,'Clustering textures\n');
end
while sum(NotYetSeparated(:))
  if Verbose>=2
    fprintf(1,'Constructing cluster pair %d: ',size(ClusterA,1)+1);
  elseif Verbose==1
    fprintf(1,'.');
  end
  [I,J]=find(Distance==max(max(Distance(find(NotYetSeparated)))));
  if ~DoIsolate
    cau=I(1);
  end
  if J(1)==cau
    cbu=J(2);
  else
    cbu=J(1);
  end
  
  % Add textures while distance within groups is less than distance between
  Prospects=notin(find(sum(NotYetSeparated([cau cbu],:))),[cau cbu]);
  Rejected=[];
  for p=Prospects
    da2a=nonzero(Distance(cau,cau)); if isempty(da2a), da2a=0; end
    db2b=nonzero(Distance(cbu,cbu)); if isempty(db2b), db2b=0; end
    dap2b=nonzero(Distance([cau p],cbu));
    dbp2a=nonzero(Distance([cbu p],cau));
    dap2a=nonzero(Distance([cau p],cau));
    dbp2b=nonzero(Distance([cbu p],cbu));
    if (max(dap2a) < min(dap2b)) & (max(db2b) < min(dap2b)) ...
	& sum(NotYetSeparated(cbu,p))>0 & ~DoIsolate
      cau=[cau p];
    elseif (max(dbp2b) < min(dbp2a)) & (max(da2a) < min(dbp2a)) ...
	& sum(NotYetSeparated(cau,p))>0
      cbu=[cbu p];
    else
      Rejected=[Rejected p];
    end
  end
  
  % Test if cluster pairs satisfy the condition for successful grouping
  h=moptxf(m_x,r_xx,M_h,N_h,'Unser',w,cau,cbu,[],[],M_h,N_h,N_T);
  if strcmp(lower(TestCrit),'meansep')
    m_v=[];
    for i=1:N_T
      m_v(i)=fstats(m_x(i),r_xx((1:M_r)+(i-1)*M_r,:),h,w);
    end
    while (min(m_v(cau)) < Separation*max(m_v(cbu))) ...
	& (min(m_v(cbu)) < Separation*max(m_v(cau))) ...
	& (length(cau)>1 | length(cbu)>1)
      if Verbose>=2
	fprintf(2,'  WARNING: First stage of clustering did not yield %s\n',...
	  sprintf('well separable clusters.'));
	fprintf(2,'  Problem:%s vs.%s yeilds feature means%s and%s.\n',...
	  sprintf(' %d',cau),sprintf(' %d',cbu),...
	  sprintf(' %4.2f',m_v(cau)),sprintf(' %4.2f',m_v(cbu)));
      end
      if length(cau)>1 | length(cbu)>1
	if Verbose>=2
	  fprintf(2,'  Removing problematic textures by distance.\n');
	end
	% Identify what pair of textures from the different clusters
	% are closest to each other
	[aIdx,bIdx]=find(Distance(cau,cbu)==min(nonzero(Distance(cau,cbu))));
	% Remove texture of the closest pair in the largest cluster
	if length(cau)>length(cbu) & ~DoIsolate
	  Rejected=[Rejected cau(aIdx)];
	  cau=notin(cau,cau(aIdx));
	else
	  Rejected=[Rejected cbu(bIdx)];
	  cbu=notin(cbu,cbu(bIdx));
	end
	h=moptxf(m_x,r_xx,M_h,N_h,'Unser',w,cau,cbu,[],[],M_h,N_h,N_T);
	for i=1:N_T
	  m_v(i)=fstats(m_x(i),r_xx((1:M_r)+(i-1)*M_r,:),h,w);
	end
      end
    end
  else
    error(sprintf('Test criterion "%s" not supported',TestCrit));
  end
  
  if MergeStep
    % Group rejected textures into clusters where within distance is less
    % than distance to cau and cbu. Test each new cluster for inclusion to
    % cau or cbu. 
    while length(Rejected)>0
      if Verbose>=2
	fprintf(1,'.');
      end
      % Note that all cluster memebers of each new formed cluster are
      % never components of new Rejected.
      Prospects=Rejected;
      Rejected=[];
      % Select most distant texture (seed)
      MaxDist=-Inf;
      for p=Prospects
	dp2a=min(nonzero(Distance(p,cau)));
	dp2b=min(nonzero(Distance(p,cbu)));
	if mean([dp2a dp2b])>MaxDist
	  ccu=p;
	  MaxDist=mean([dp2a dp2b]);
	end
      end
      Prospects=notin(Prospects,ccu);
      % Grow from seed as long as possible
      for p=Prospects
	da2a=nonzero(Distance(cau,cau));
	db2b=nonzero(Distance(cbu,cbu));
	dcp2a=nonzero(Distance(p,cau));
	dcp2b=nonzero(Distance(p,cbu));
	dcp2c=nonzero(Distance([ccu p],ccu));
	if (max(dcp2c) < min(dcp2a)) & (max(da2a) < min(dcp2a)) ...
	    & (max(dcp2c) < min(dcp2b)) & (max(db2b) < min(dcp2b))
	  ccu=[ccu p];
	else
	  Rejected=[Rejected p];
	end
      end
      % Test new cluster for inclusion into both cau and cbu by filter design
      % and mean analysis.
      if strcmp(lower(TestCrit),'meansep')
	% Test new cluster for inclusion into cbu
	h_B=moptxf(m_x,r_xx,M_h,N_h,'Unser',w,cau,[cbu ccu],[],[],M_h,N_h,N_T);
	m_v=[];
	for i=1:N_T, m_v(i)=fstats(m_x(i),r_xx((1:M_r)+(i-1)*M_r,:),h,w); end
	ccu2cbuOK=(min(m_v([cbu ccu])) > Separation*max(m_v(cau))) ...
	  | (min(m_v(cau)) > Separation*max(m_v([cbu ccu])));
	% Test new cluster for inclusion into cau
	h=moptxf(m_x,r_xx,M_h,N_h,'Unser',w,[cau ccu],cbu,[],[],M_h,N_h,N_T);
	m_v=[];
	for i=1:N_T, m_v(i)=fstats(m_x(i),r_xx((1:M_r)+(i-1)*M_r,:),h,w); end
	ccu2cauOK = ((min(m_v([cau ccu])) > Separation*max(m_v(cbu))) ...
	  | (min(m_v(cbu)) > Separation*max(m_v([cau ccu])))) ...
	  & ~DoIsolate;
	% Include new cluster if possible
	if ccu2cauOK & (~ccu2cbuOK | ccu2cbuOK & (length(cau)<length(cbu)))
	  cau=[cau ccu];
	elseif ccu2cbuOK 
	  cbu=[cbu ccu];
	else
	  % New cluster not possible to include to neither cau nor cbu
	end
      else
	error(sprintf('Test criterion "%s" not supported',TestCrit));
      end
    end
  end % 				if MergeSte
  if Verbose>=2
    fprintf(1,'%d vs. %d textures\n',length(cau),length(cbu));
  end
  
  % Set elements that will now be separated to zero in NotYetSeparated
  NotYetSeparated(cau,cbu)=zeros(length(cau),length(cbu));
  NotYetSeparated(cbu,cau)=zeros(length(cbu),length(cau)); 
  if Hierarchical
    NotYetSeparated(cau,cau)=zeros(length(cau),length(cau));
    NotYetSeparated(cbu,cbu)=zeros(length(cbu),length(cbu)); 
  end
  % Get sizes
  [M_CA,N_CA]=size(ClusterA); [M_cau,N_cau]=size(cau); 
  [M_CB,N_CB]=size(ClusterB); [M_cbu,N_cbu]=size(cbu);
  % Update cluster databases
  ClusterA=[[ClusterA zeros(M_CA,N_cau-N_CA)];[cau zeros(M_cau,N_CA-N_cau)]];
  ClusterB=[[ClusterB zeros(M_CB,N_cbu-N_CB)];[cbu zeros(M_cbu,N_CB-N_cbu)]];
end
if Verbose==1
  fprintf(1,'\n');
end
anc=sum(sum([ClusterA ClusterB]>0))/N_T;

% Hierarchical clustering of clusters
if Hierarchical
  [M_CA,N_CA]=size(ClusterA);
  [M_CB,N_CB]=size(ClusterB);
  Clusters=[];
  Clusters(1:M_CA,1:N_CA)=ClusterA;
  Clusters(M_CA+(1:M_CB),1:N_CB)=ClusterB;
  Clusters=sortrows(sort(Clusters')');
  [M_C,N_C]=size(Clusters);
  for ClusterNo=1:M_C
    ThisCluster=sort(nonzero(Clusters(ClusterNo,:)));
    if length(ThisCluster)<=1
      DoCluster=0;
    else % 				length(ThisCluster)<=1
      if ClusterNo==1
	DoCluster=1;
      else % 				k~=1
	if Clusters(ClusterNo,:)==Clusters(ClusterNo-1,:)
	  DoCluster=0;
	else
	  DoCluster=1;
	end
      end % 				if k==1
    end % 				if length(ThisCluster)

    if DoCluster
      m_x_sub=[];
      r_xx_sub=[];
      for t=ThisCluster
	m_x_sub=[m_x_sub;m_x(t)];
	r_xx_sub=[r_xx_sub;r_xx((1:M_r)+(t-1)*M_r,:)];
      end
      [cau,cbu]=txturgroup(Distance(ThisCluster,ThisCluster),...
	m_x_sub,r_xx_sub,M_h,N_h,w,TestCrit,...
	Separation,IsolateTex,Verbose,Hierarchical,MergeStep);
      for l=1:size(cau,1)
	[M_CA,N_CA]=size(ClusterA);
	[M_CB,N_CB]=size(ClusterB); 
	ClusterA(M_CA+1,1:length(nonzero(cau(l,:)))) ...
	  =ThisCluster(nonzero(cau(l,:)));
	ClusterB(M_CB+1,1:length(nonzero(cbu(l,:)))) ...
	  =ThisCluster(nonzero(cbu(l,:)));
      end
    end % 				if DoCluster
  end % 				for k
end % 					if Hierarchical

% Prune duplicate texture group pairs
[M_CA,N_CA]=size(ClusterA);
[M_CB,N_CB]=size(ClusterB);
if M_CA>1 & N_CA>1, ClusterA=sort(ClusterA')'; end
if M_CB>1 & N_CB>1, ClusterB=sort(ClusterB')'; end
[dummy,idx]=sortrows([ClusterA ClusterB]);
ClusterA=ClusterA(idx,:);
ClusterB=ClusterB(idx,:);
[M_CA,N_CA]=size(ClusterA);
idx=1;
for k=2:M_CA
  if ~( min(ClusterA(k,:)==ClusterA(k-1,:))>0 ...
      & min(ClusterB(k,:)==ClusterB(k-1,:))>0)
    idx=[idx k];
  end
end
ClusterA=ClusterA(idx,:);
ClusterB=ClusterB(idx,:);
