function [o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o15,o16,o17,o18,o19,o20]=optxf(Criterion,x_1,x_2,hh,w,r_xx1,r_xx2)
% OPTXF	Design optimal texture filter(s)
%	h=optxf('Criterion',x_1,x_2,[M_h N_h],w) designs an optized 
%	FIR filter of size M_h by N_h and returns the filter coefficients h.
%	The autocorrelation function for x_1 and x_2 are inherently estimated.
%	
%	h=optxf('Criterion',x_1,x_2,[M_h N_h],w,r_xx1,r_xx2) uses the given 
%	estimates for the correlation functions, r_xx1 and r_xx2.
%	
%	[h,J]=optxf(...) also return the criterion function, J.
%	
%	...=optxf('Criterion',x_1,x_2,h_init,...) for any non-closed form
%	solution supplies the initial filter h_init for the iterative
%	optimization. Currently only supported with 'FisherG' criterion.
%	
%	[h_1,h_2,h_3,...,h_n,J]=optxf(...) returns the n best filters,
%	with h_1 being the best, and a vector, J, with the corresponding
%	criterion function values.
%	
%	Supported criterion functions are
%	- 'Fisher'	Closed form solution wrt. the Fisher criterion
%	- 'FisherG'	Gradient search solution wrt. the Fisher criterion
%	- 'Singh'	Closed form solution wrt. the Singh criterion
%	- 'Unser'	Closed form solution wrt. the Unser criterion
%	- 'Histogram'	Closed form solution wrt. the generalized criterion,
%			select filter yielding minimum feature histogram 
%			overlap, or better, maximum distance between hist.'s
%	
%	See also: FISHER_GOPT, FSTATS, DUNN

if exist('w')
  [M_w,N_w]=size(w);
end

if strcmp(lower(Criterion),'fisher') | strcmp(lower(Criterion),'histogram')
  Verbose=1;
else
  Verbose=0;
end

% Compute the correlation matrices
if exist('r_xx1')~=1,
  r_xx1=xcorrf2(x_1)/(size(x_1,1)*size(x_1,2));
  r_xx2=xcorrf2(x_2)/(size(x_2,1)*size(x_2,2));
end

if strcmp(lower(Criterion),'fisherg')
  % Gradient search solution
  options=[];		% Clear default options
  options(1)=1;		% Display intermediate gradient search results
  options(14)=1000; 	% Max number of iterations in gradient search
  % options(100)=1; 	% Constrained search
  %%% [o1,o2]=fisher_gopt(hh,x_1,x_2,w,options,r_xx1,r_xx2);
  error('Gradient search temporarily unavailable');
else
  % Closed form solution
  M_h=hh(1);
  N_h=hh(2);

  % Trim correlation function to only contain necessary samples (minimize
  % parameter passing)
  K=M_h+M_w-1; % Max vertical correlation-lag required
  L=N_h+N_w-1; % Max horizontal correlation-lag required
  K0_x1=(size(r_xx1,1)+1)/2; % Vertical zero lag position, r_xx1
  L0_x1=(size(r_xx1,2)+1)/2; % Horizontal zero lag position, r_xx1
  K0_x2=(size(r_xx2,1)+1)/2; % Vertical zero lag position, r_xx2
  L0_x2=(size(r_xx2,2)+1)/2; % Horizontal zero lag position, r_xx2
  r_xx1=r_xx1((-K:K)+K0_x1,(-L:L)+L0_x1);
  r_xx2=r_xx2((-K:K)+K0_x2,(-L:L)+L0_x2);

  R_xx1=r2R(r_xx1,M_h,N_h);
  R_xx2=r2R(r_xx2,M_h,N_h);
  % Solve the optimization equation
  [V,D]=eig(inv(R_xx2)*R_xx1);
  D=diag(D);

  % Select the eigenvector yielding the maximum object function
  if Verbose
    fprintf(1,'Optimal filter design');
  end
  for i=1:M_h*N_h,
    if Verbose
      fprintf(1,'.');
    end
    h_=V(:,i);
    h=reshape(h_,N_h,M_h).';
    if strcmp(lower(Criterion),'fisher')
      [m_1,var_1]=fstats(mean(x_1(:)),r_xx1,h,w);
      [m_2,var_2]=fstats(mean(x_2(:)),r_xx2,h,w);
      J(i)=(m_1-m_2)^2 / (var_1 + var_2);
    elseif strcmp(lower(Criterion),'singh')
      J(i)=D(i);
    elseif strcmp(lower(Criterion),'unser')
      J(i)=D(i)+D(i)^(-1)-2;
    elseif strcmp(lower(Criterion),'histogram')
      x=[x_1 x_2];
      [M_x,N_x]=size(x);
      % Feature extraction, use separable smoothing if possible
      w_h=w(1,:);
      w_v=w(:,1);
      w_s=conv2(w_h,w_v);
      Scale=mean(mean(w./w_s));
      w_h=sqrt(Scale)*w_h;
      w_v=sqrt(Scale)*w_v;
      w_s=conv2(w_h,w_v);
      dw=w-w_s;
      if std(dw(:))/std(w(:))<1e-10
	% The smoothing filter is separable
	v=filter2(w_h,filter2(w_v,abs(filter2(h,x)).^2));
      else	
	v=filter2(w,abs(filter2(h,x)).^2); % Feature extraction
      end
      [M_v,N_v]=size(v);
      v=v(M_h+M_w-1:M_v-M_h-M_w+2,N_h+N_w-1:N_v-N_h-N_w+2); % Rem. edge eff.
      [M_v,N_v]=size(v);
      v_1=v(:,1:N_v/2);     v_1=v_1(:)';
      v_2=v(:,N_v/2+1:N_v); v_2=v_2(:)';
      
      % Determine the interesting histogram range
      % [n_1,j_1]=hist(v_1,128); [n_2,j_2]=hist(v_2,128);
      % j=(0:511) .* (max([j_1 j_2])-min([j_1 j_2]))/511 + min([j_1 j_2]);
      j_lo=max([min([max(v_1) min(v_2)]) min([max(v_2) min(v_1)])]);
      j_hi=min([max([max(v_1) min(v_2)]) max([max(v_2) min(v_1)])]);
      j=(0:0.005:1) * (j_hi-j_lo) + j_lo;
      dj=j(2)-j(1);
      j=[j(1)-dj j j(length(j))+dj];
      
      % Determine the relative histogram overlap
      [n_1,j]=hist(v_1,j); [n_2,j]=hist(v_2,j);
      J(i)=-sum(min([n_1;n_2]))/sum([n_1 n_2]); % Classification error
      if J(i)==0
	% No overlap, let criterion be distance between tails
	J(i)=(max(min(v_1),min(v_2)) - min(max(v_1),max(v_2)))...
	  /(max([v_1 v_2]) - min([v_1 v_2]));
      end
    else
      error(['Not supported criterion function: ' Criterion]);
    end
  end
  if Verbose
    fprintf(1,'\n');
  end

  % Prepare output
  J_sort=-sort(-J);
  for i=1:max(1,nargout-1)
    idx=find(J==J_sort(i));
    h_=V(:,idx(1)); h=reshape(h_,N_h,M_h).';
    eval(sprintf('o%d=h;',i));
  end
  if nargout>1
    eval(sprintf('o%d=J_sort(1:nargout-1);',nargout));
  end
end