function [h,U,V,J_min,z]=dunn(i_1,i_2,sigma,N_s,p0,p1,p2)
% DUNN	Dunn/Higgins optimal Gabor filter design procedure
%	h=dunn(x_1,x_2) designs a Gabor filter, h, optimal for separating
%	textures x_1 and x_2. Uses sigma=8 and samples with 20 spatial 
%	samples, s.
%	
%	h=dunn(x_1,x_2,sigma,N_s) specifies a different sigma and number of 
%	spatial samples s, N_s.
%	
%	[h,U,V,J,z]=dunn(x_1,x_2,...) returns also the center 
%	frequency specifications U and V in the range 0 to 0.5, with 0.5 
%	correspoding to half the sampling frequency, the criterion J, and the
%	filtered composite image, z. 
%	
%	Optional parameters are:
%	- 'Verbose'	Run in verbose mode
%	- 'Error'	Use the minimum error criterion (defualt)
%	- 'ModError'	Use the modified minimum error criterion
%	- 'Fisher'	Use the Fisher criterion
%	- 'Postfilter',sigma_p
%			Model a post filter according to Weldon (1994)
%			Default sigma_p value is 8.5
%	
%	Reference:
%	[1] D. Dunn and W.E. Higgins, "Optimal Gabor Filters for Texture
%	    Segmentation", IEEE Trans. Image Proc., Vol. 4, No. 7,
%	    July 1995, pp. 947-964.

Verbose=0;
Criterion='error';
sigma_p=8.5;
UsePostFilter=0;

% Check input
if sum(abs(size(i_1)-size(i_2)))>0
  error('Currently, only equally sized textures are supported.');
end
[M,N]=size(i_1);
if exist('sigma')~=1, sigma=[]; end
if isempty(sigma),    sigma=8;  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; 
  elseif strcmp(Par,'error'),      Criterion=Par;
  elseif strcmp(Par,'moderror'),   Criterion=Par;
  elseif strcmp(Par,'fisher'),     Criterion=Par;
  elseif strcmp(Par,'postfilter'), 
    UsePostFilter=1;
    if exist(sprintf('p%d',ParNo+1))
      if ~isstr(sprintf('p%d',ParNo+1))
	sigma_p=eval(sprintf('p%d',ParNo+1));
	ParNo=ParNo+1;
      end
    end
  else
    error(sprintf('Unsupported parameter "%s"',Par));
  end
  ParNo=ParNo+1;
end


% 1) Select a value of sigma (section IV-D): input
nstd=2.5;

% 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:2
  i=eval(sprintf('i_%d',Texture));
  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)));
    % 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
  eval(sprintf('mean_F%d=mean_F;',Texture));
  eval(sprintf('var_F%d =var_F ;',Texture));
end

% Weldon's model of post filtering
if UsePostFilter
  var_F1=var_F1*sigma^2/sigma_p^2;
  var_F2=var_F2*sigma^2/sigma_p^2;
end

% 3b) Compute P_E(U,V)=error probability
% Compute integration limits, x_1 and x_2
a=1./var_F1 - 1./var_F2;
b=-2*mean_F1./var_F1 + 2*mean_F2./var_F2;
c=mean_F1.^2./var_F1 - mean_F2.^2./var_F2 + 2*log(sqrt(var_F1./var_F2));
x_1=(-b + sqrt(b.^2 - 4*a.*c)) ./ (2*a);
x_2=(-b - sqrt(b.^2 - 4*a.*c)) ./ (2*a);

% Make mean_F1 and var_F1 correspond to the distribution most likely
% (midway) between x_1 and x_2
x_m=(x_1+x_2)/2;
p_m1=1./(sqrt(2*pi*var_F1)) * exp(-0.5*(x_m-mean_F1).^2./var_F1);
p_m2=1./(sqrt(2*pi*var_F2)) * exp(-0.5*(x_m-mean_F2).^2./var_F2);
SwapMask=p_m2 > p_m1;
[mean_F1,mean_F2]=maskswap(mean_F1,mean_F2,SwapMask);
[var_F1 ,var_F2 ]=maskswap(var_F1 ,var_F2 ,SwapMask);
[x_1    ,x_2    ]=maskswap(x_1    ,x_2    ,SwapMask);

% Criterion evaluation
if strcmp(Criterion,'error')
  % Integrate the error under the Gaussian distributions
  %        x_1                  x_2                  Inf         
  %        /-                   /-                   /-          
  % P_e =  |  P_2 p_2(x) dx  +  |  P_1 p_1(x) dx  +  |  P_2 p_2(x) dx
  %       -/                   -/                   -/           
  %       -Inf                 x_1                   x_2    
  % under the assumption P_1=P_2=1/2
  P_1=1/2; P_2=1/2;
  J=  P_2*gaussint(mean_F2,sqrt(var_F2),-Inf,x_1) ...
    + P_1*gaussint(mean_F1,sqrt(var_F1),x_1,x_2) ...
    + P_2*gaussint(mean_F2,sqrt(var_F2),x_2,Inf);
elseif strcmp(Criterion,'moderror')
  % Only integrate with one divisor
  % Expected to yield better features after smoothing
  P_1=1/2; P_2=1/2;
  P_ea=P_2*gaussint(mean_F2,sqrt(var_F2),-Inf,x_1) ...
    +  P_1*gaussint(mean_F1,sqrt(var_F1),x_1,Inf);
  P_eb=P_2*gaussint(mean_F1,sqrt(var_F1),-Inf,x_2) ...
    +  P_1*gaussint(mean_F2,sqrt(var_F2),x_2,Inf);
  J=min(P_ea,P_eb);
elseif strcmp(Criterion,'fisher')
  J=-(mean_F1-mean_F2).^2 ./ (var_F1+var_F2);
else
  error('Internal error, illegal state');
end

% 4) Rank J(U,V) to indicate most effective filters
[J_min,U_idx]=min(min(J.')); U=(U_idx-(M_hg+1)/2)/M_hg;
[J_min,V_idx]=min(min(J));   V=(V_idx-(N_hg+1)/2)/N_hg;
UV=U+sqrt(-1)*V;

% User feedback
if Verbose
  m1=mean_F1(U_idx,V_idx);  std1=sqrt(var_F1(U_idx,V_idx));
  m2=mean_F2(U_idx,V_idx);  std2=sqrt(var_F2(U_idx,V_idx));
  x1=x_1(U_idx,V_idx);
  x2=x_2(U_idx,V_idx);
  fprintf(1,'min(J)=%4.2f <= m=%1.0f,%1.0f, sigma=%1.0f,%1.0f. U=%f, V=%f\n',...
    J_min,m1,m2,std1,std2,U,V);
  [Hist,Bins] ...
    =hist([m1+randn(1,30000)*std1; m2+randn(1,30000)*std2]',100);
  fprintf(1,'Histogram overlap %4.2f%%, theoretical %4.2f%%\n',...
    100*sum(min(Hist'))/sum(Hist(:)),...
    100* (0.5*gaussint(m2,std2,-Inf,x1) + 0.5*gaussint(m1,std1,x1,x2) ...
    + 0.5*gaussint(m2,std2,x2,Inf)));
  plot(Bins,Hist); title(Criterion); drawnow;
end

% Construct filter
[dummy,h]=gabor([],abs(UV),angle(UV),sigma,sigma,nstd);
if nargout>1
  [M_h,N_h]=size(h);
  H=zeros(size([i_1 i_2])); H(1:M_h,1:N_h)=h; H=fft2(H);
  z=ifft2(H .* fft2([i_1 i_2]));
end