function [h,J]=fisher_gopt(h,x_1,x_2,w,IntermSaveProb,options,r_xx1,r_xx2,w_1,o_1,w_2,o_2,w_3,o_3,w_4,o_4,w_5,o_5)
% FISHER_GOPT  Gradient search texure filter optimization, Fisher criterion
%	[h,J]=fisher_gopt(h,x_1,x_2,w,IntermSaveProb,options,r_xx1,r_xx2)
%	starts gradient search with the filter h as initial and returns the
%	optimized filter h, and the Fischer criterion value J. The parameter
%	w is the smoothing filter to be assumed in the model. The 
%	autocorrelation functions of x_1 and x_2 are inherently 
%	estimated if they are empty or not given. Intermediate results
%	are saved to the file fisher_gopt_h.mat with probability IntermSaveProb
%	at each criterion function evaluation.
%	
%	[h,J]=fisher_gopt(h,x_1,x_2,w,IntermSaveProb,options,r_xx1,r_xx2,...
%	                  w_1,o_1,w_2,o_2,...) 
%	performs the optimization iteratively, first with smoothing filter
%	w_1, then w_2,... and finally with w. Uses options o_1 in first pass,
%	options o_2 in second etc. Useful for speeding up
%	the search, doing fast searches with small smoothing filters
%	initially.
%	
%	For information on the options parameter, see help foptions
%	Additional: options(100)=1 => use constrained search.

% Estimate the autocorrelation functions, if not supplied
if exist('r_xx2')~=1, r_xx2=[]; end
if exist('r_xx1')~=1, r_xx1=[]; end
if isempty(r_xx1), r_xx1=xcorrf2(x_1); end
if isempty(r_xx2), r_xx2=xcorrf2(x_2); end

% Initialize
[M_h,N_h]=size(h);
[M_w,N_w]=size(w);
m_x1=mean(x_1(:));
m_x2=mean(x_2(:));

% Iterate through list of smoothing filters
i=1;
while exist(sprintf('w_%d',i)) & exist(sprintf('o_%d',i))
  TheseOpts=eval(sprintf('o_%d',i));
  ThisW=eval(sprintf('w_%d',i));
  h_=reshape(h.',M_h*N_h,1);
  J=-fisher_gopt_J(h_,w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h);
  fprintf(1,'Starting initial gradient search %d: J=%6.4f\n',i,J);
  h=fisher_gopt(h,x_1,x_2,ThisW,IntermSaveProb,TheseOpts,r_xx1,r_xx2);
  i=i+1;
end
h_=reshape(h.',M_h*N_h,1);
J=-fisher_gopt_J(h_,w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h);
fprintf(1,'Starting gradient search, J=%6.4f (max %d iterations)\n',...
  J,options(14));

% Parse input arguments
if exist('options')~=1
  options=[];
end
if length(options)>=100,
  constrained=options(100);
  options=trimend(options(1:99));
else
  constrained=0;
end

% 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);

% Perform the gradient search
if options(14)>0
  if constrained,
    [h_,options]=constr('fisher_gopt_J',h_,options,[],[],'fisher_gopt_DJ',...
      w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h,IntermSaveProb);
  else
    [h_,options]=fminu('fisher_gopt_J',h_,options,'fisher_gopt_DJ',...
      w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h,IntermSaveProb);
  end
end

% Compute the final criterion function value
if nargout>1
  J=-fisher_gopt_J(h_,w,m_x1,m_x2,r_xx1,r_xx2,M_h,N_h);
end

% Transform filter back from vector representation to 2D
h=reshape(h_,M_h,N_h).';


