function [o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o15,o16,o17,o18,o19,o20]=fisher_opt(x_1,x_2,M_h,N_h,w,r_xx1,r_xx2)
% FISHER_OPT  Design texture filter(s) optimal wrt. the Fisher criterion
%	h=fisher_opt(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=fisher_opt(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]=fisher_opt(...) also return the criterion function, J.
%	
%	[h_1,h_2,h_3,...,h_n,J]=fisher_opt(...) returns the n best filters,
%	with h_1 being the best, and a vector, J, with the corresponding
%	criterion function values.

[M_w,N_w]=size(w);
WaitBarHndl = waitbar(0,'Computing optimal filter');

% Compute the correlation matrices
if exist('r_xx1')~=1,
  r_xx1=xcorrf2(x_1)/(size(x_1,1)*size(x_1,2));
  waitbar(10);
  r_xx2=xcorrf2(x_2)/(size(x_2,1)*size(x_2,2));
end
waitbar(20);
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);
waitbar(30);

% Select the eigenvector yielding the maximum object function
for i=1:M_h*N_h,
  h_=V(:,i);
  h=reshape(h_,N_h,M_h).';
  [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);
  waitbar(30+70*i/(M_h*N_h));
  title(sprintf('Computing optimal filter, J_max=%4.2f, J(%d)=%4.2f',...
    max(J),i,J(i)));
end

% Prepare output
J_sort=-sort(-J);
for i=1:max(1,nargout-1)
  idx=find(J==J_sort(i));
  h_=V(:,idx); 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

close(WaitBarHndl);
