function [h,J]=pdf_opt(x_1,x_2,N,r_1,r_2)
%PDF_OPT	[h,J]=pdf_opt(x_1,x_2,N,r_1,r_2)
%	Computes an optized FIR filter of size NxN and returns the filter 
%	coefficients h.r_1 and r_2 are the 2D autocorrelation matrices
%	of x_1 and x_2, and if supplied, the x_1 and x_2 inputs may be [].
%	The filter is optimal with respect to the true Probability
%	Density Function (PDF) for the local energy.
%
%	Example use
%	x_1=hipsread('brodatz/D92.hips',0:127,0:127);
%	x_2=hipsread('brodatz/D4.hips',0:127,0:127);
%	h=pdf_opt(x_1,x_2,6);


% Compute the correlation matrices
if nargin<4,
  r_1=corrcoeff(x_1,x_1,N);
  r_2=corrcoeff(x_2,x_2,N);
end
R_1=r2R(r_1,N);
R_2=r2R(r_2,N);

% Solve the optimization equation
[V,D]=eig(inv(R_2)*R_1);

% Select the eigenvector yielding the maximum object function
i_max=-inf;
J_max=-inf;
J_2_max=-inf;
for i=1:N*N,
  fprintf(1,'%03d:%03d',i,N*N);
  % Select filter from candidate filters
  h=reshape(V(:,i),N,N)';
  
  % Compute the local energies in each image
  [L_1,M_1]=size(x_1);
  [L_2,M_2]=size(x_2);
  z_1=smooth(filter2(h,x_1).^2);
  z_2=smooth(filter2(h,x_2).^2);
  z_1=z_1(20:L_1-20,20:M_1-20);
  z_2=z_2(20:L_2-20,20:M_2-20);
  [L_1,M_1]=size(z_1);
  [L_2,M_2]=size(z_2);
  z_1=reshape(z_1,1,L_1*M_1);
  z_2=reshape(z_2,1,L_2*M_2);
  
  % Determine the interesting histogram range
  [n_1,j_1]=hist(z_1,128);
  [n_2,j_2]=hist(z_2,128);
  j=(0:511) .* (max([j_1 j_2])-min([j_1 j_2]))/511 + min([j_1 j_2]);
  
  % Determine the histogram overlap
  [n_1,j]=hist(z_1,j);
  [n_2,j]=hist(z_2,j);
  J=-sum(min([n_1;n_2]));
  J_2=(max(min(z_1),min(z_2)) - min(max(z_1),max(z_2)))...
      /(max([z_1 z_2]) - min([z_1 z_2]));

  % Terminal feedback (information to the user)
  % clf; hold on;
  % plot(j,n_1,'y-');
  % plot(j,n_2,'c:');
  % title(sprintf('%d:%d  -  J=%d,  J_2=%5.3f (max: i=%d,  J=%d,  J_2=%5.3f)',...
  %     i,N*N,J,J_2,i_max,J_max,J_2_max));
  % hold off;
  % pause(2)
  
  if J==0  & J_2>J_2_max,
    i_max=i;
    h_max=h;
    J_2_max=J_2;
    J_max=J;
  elseif (J>J_max),
    i_max=i;
    h_max=h;
    J_max=J;
  end
  fprintf(1,'\b\b\b\b\b\b\b');
end

J=J_max;
h=h_max;

