function [FiltPars,J_min]=weldon(NumFilt,N_tex,data,sigmas,sigma_ps,N_s,p0,p1,p2)
% WELDON  Weldon/Higgins optimal Gabor filter bank design procedure
%	
%	PS! The expression for the variance after post filtering in
%	[weldo94] is probably wrong, and is replaced in this function.
%	See the source code below.
%	
%	Optional parameters are:
%	- 'Verbose'	Run in verbose mode
%	- 'NumFilt',K	Designs K filters (default 1)
%	- 'Postfilter',sigma_p
%			Model a post filter according to Weldon (1994)
%			Default sigma_p value is 8.5

% ---------------------------------------------------------------
% Default values
% ---------------------------------------------------------------
Verbose=0;
nstd=2.0;

% ---------------------------------------------------------------
% Check/parse input
% ---------------------------------------------------------------
[M,N]=size(data);
M=M/N_tex;
if exist('sigmas')~=1,  sigmas=[];          end
if isempty(sigmas),     sigmas=[2 4 8 16];  end
if exist('sigma_ps')~=1,sigma_ps=[];        end
if isempty(sigma_ps),   sigma_ps=2*sigmas;  end
if exist('N_s')~=1,     N_s=[];             end
if isempty(N_s),        N_s=100;            end
ParNo=0;
while exist(sprintf('p%d',ParNo))
  Par=lower(eval(sprintf('p%d',ParNo)));
  if     strcmp(Par,'verbose'),    Verbose=1; 
  else
    error(sprintf('Unsupported parameter "%s"',Par));
  end
  ParNo=ParNo+1;
end



% ---------------------------------------------------------------
% Determine mean and variances for a range of filters
% ---------------------------------------------------------------
% 1) Select a value of sigma (section IV-D): input
for sigma_idx=1:length(sigmas) % 		Loop through list of sigmas
  sigma=sigmas(sigma_idx);
  sigma_p=sigma_ps(sigma_idx);
  
  % Construct the underlying Gaussian
  [dummy,h_G]=gausslp([],sigma,nstd);
  [M_hg,N_hg]=size(h_G);
  
  % 2) Simultaneously apply a large number of Gabor filters
  % 2a) Randomly select a set s of points within a texture
  rand('seed',0);
  s=round((M-M_hg-1) * rand(1,N_s)) + j*round((N-N_hg-1) * rand(1,N_s));
  
  % 2b) Foreach (X,Y) in s compute |F_XY(U,V)|, F_XY is the windowed Fourier
  % transform centered at (X,Y) (section IV-A)
  for Texture=1:N_tex
    fprintf(1,'.');
    i=data((1:M)+(Texture-1)*M,1:N);
    mean_F=0;
    var_F=0;
    for k=1:length(s)
      X=real(s(k));
      Y=imag(s(k));
      i_sub=i(X+(1:M_hg),Y+(1:N_hg));
      F=abs(fftshift(fft2(i_sub .* h_G)))/(M_hg*N_hg);
      % 3a) Compute bar(mu)=mean_F
      mean_F=mean_F + F/length(s);
      eval(sprintf('F_%d=F;',k));
    end
    % 3a) Compute bar(sigma)^2=var_F
    for k=1:length(s)
      F=eval(sprintf('F_%d',k));
      eval(sprintf('clear F_%d',k));
      var_F=var_F + (F - mean_F).^2/length(s);
    end
    % Weldon's model of post filtering
    % var_F=var_F*sigma^2/sigma_p^2; PROBABLY WRONG
    eval(sprintf('mean_F%d_s%d=mean_F;',Texture,sigma_idx));
    eval(sprintf('var_F%d_s%d =var_F ;',Texture,sigma_idx));
  end % 				for Texture
end % 					for sigma
fprintf(1,'\n');

% ---------------------------------------------------------------
% Determine optimal filter combination
% ---------------------------------------------------------------
FiltPars=[]; % 				Parameters for the designed filters
Means=[]; % 				Means of selected filters
Vars=[]; % 				Variances of selected filters
for FiltNo=1:NumFilt
  fprintf(1,'h_%d:',FiltNo);
  % Criterion evaluation
  PP=(1/N_tex)^2;
  J_min=Inf;
  J_mins=Inf * ones(3,length(sigmas));
  for s=1:length(sigmas) % 		Loop through list of sigmas
    sigma=sigmas(s);
    sigma_p=sigma_ps(s);
    [dummy,h_G]=gausslp([],sigma,nstd); % Reconstruct the underlying Gaussian
    [M_hg,N_hg]=size(h_G); % 		Get the size of the filter
    for U_idx=1:M_hg % 			Loop through vertical freq.s
      U=(U_idx-(M_hg+1)/2)/M_hg;
      for V_idx=1:N_hg % 		Loop through horizontal freq.s
	V=(V_idx-(N_hg+1)/2)/N_hg;
	UV=U+sqrt(-1)*V;
	% Check if this or adjacent frequencies is used in already selected
	% filters
	if FiltNo>1
	  NonUsedFreq=min(abs(FiltPars(1:FiltNo-1,1)-UV)) ...
	    > sqrt(2)*abs(1/M_hg + sqrt(-1)*1/N_hg);
	else
	  NonUsedFreq=1;
	end
	if NonUsedFreq
	  % Do not allow same and close frequencies over again
	  J_c=0;
	  J_l=0;
	  % Examine distance for all texture pairs
	  for alpha=1:N_tex-1 % 	Texture 1 in the pair
	    % Determine mu_alpha and C_alpha
	    mu_alpha=eval(sprintf('mean_F%d_s%d(U_idx,V_idx)',alpha,s));
	    C_alpha =eval(sprintf('var_F%d_s%d(U_idx,V_idx)',alpha,s));
	    if size(Means,1)>0
	      mu_alpha=[Means(:,alpha);mu_alpha];
	      C_alpha=[Vars(:,alpha);C_alpha];
	    end
	    C_alpha=diag(C_alpha);
	    for beta=alpha+1:N_tex % 	Texture 2 in the pair
	      % Determine mu_beta and C_beta
	      mu_beta=eval(sprintf('mean_F%d_s%d(U_idx,V_idx)',beta,s));
	      C_beta =eval(sprintf('var_F%d_s%d(U_idx,V_idx)',beta,s));
	      if size(Means,1)>0
		mu_beta=[Means(:,beta);mu_beta];
		C_beta=[Vars(:,beta);C_beta];
	      end
	      C_beta=diag(C_beta);
	      % B-distance
	      B=1/8 * (mu_alpha-mu_beta).' * inv((C_alpha+C_beta)/2) ...
		*     (mu_alpha-mu_beta) ...
		+ 1/2 * log(det((C_alpha+C_beta)/2) ...
		/           sqrt(det(C_alpha)*det(C_beta)));
	      % Composite distance
	      J_c=J_c + exp(-B);
	    end % 			for beta
	  end % 			for alpha
	  J_c = 1/(N_tex-1) * sqrt(PP) * J_c;
	  J_l = 1/FiltNo * 2*N_tex * (sigma^2 + sigma_p^2) / 256^2;
	  for gamma=1:FiltNo-1
	    J_l = J_l + 1/FiltNo * 2*N_tex ...
	      * (FiltPars(gamma,2)^2 + sigma_p^2) / 256^2;
	  end % 			for gamma
	  J = J_c + J_l;
	  % fprintf(1,'J=%4.2f (min=%4.2f) size(mu)=%dx%d size(C)=%dx%d\n',...
	  % J,J_min,size(mu_alpha),size(C_alpha));
	  if J<J_mins(1,s) % 		New best filter for specific sigma
	    J_mins(1:3,s)=[J; J_c; J_l];
	  end
	  if J<J_min % 			New best filter
	    J_min=J;
	    FiltPars(FiltNo,1:2)=[UV sigma];
	    for t=1:N_tex
	      Means(FiltNo,t)=eval(sprintf('mean_F%d_s%d(U_idx,V_idx)',t,s));
	      Vars (FiltNo,t)=eval(sprintf('var_F%d_s%d(U_idx,V_idx)',t,s));
	    end
	  end % 			if J<J_min
	end % 				if min(UV)
      end % 				for V_idx
    end % 				for U_idx
    fprintf(1,'\tJ_min(sigma=%4.1f)=%7.5f+%7.5f=%7.5f\n',...
      sigma,J_mins(2,s),J_mins(3,s),J_mins(1,s));
  end % 				for s
  fprintf(1,' =>\tJ=%7.5f <= U=%6.3f, V=%6.3f, sigma=%4.1f\n',J_min,...
    real(FiltPars(FiltNo,1)),imag(FiltPars(FiltNo,1)),FiltPars(FiltNo,2));
  [SortMeans,SortIdx]=sort(Means(FiltNo,:));
  SortStds=sqrt(Vars(FiltNo,SortIdx));
  fprintf(1,'\tMeans (sorted)=%s\n',sprintf('%5.1e ',SortMeans));
  fprintf(1,'\tStandard dev.s=%s\n',sprintf('%5.1e ',SortStds));
end % 					for FiltNo




