function [h,J]=vms_opt(x_1,x_2,c,N,N_w,r_1,r_2)
%VMS_OPT	[h,J]=vms_opt(x_1,x_2,c,N,N_w,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 Variance/Mean
%	Sum (VMS). 
%
%	Example use
%	x_1=hipsread('brodatz/D92.hips',0:127,0:127);
%	x_2=hipsread('brodatz/D4.hips',0:127,0:127);
%	h=vms_opt(x_1,x_2,0.5,4,4);

% Compute the correlation matrices
if nargin<6,
  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
J_max=-inf;
for i=1:N*N,
  h=V(:,i);
  my_1=fmean(r_xx,h,N);
  my_2=fmean(r_yy,h,N);
  sigmaSq_1=fvar(r_xx,h,N,N_w);
  sigmaSq_2=fvar(r_yy,h,N,N_w);
  J=c*(my_1-my_2)^2 - (1-c)*(sigmaSq_1 + sigmaSq_2);
  if J>J_max,
    h_max=h;
    J_max=J;
    i_max=i;
  end
end

J=J_max;
h=reshape(h_max,N,N)';

