function y=optexatr(x,x_1,x_2,N,N_w_mod,N_w,DesignMethod,Smoother)
%OPTEXATR	y=optexatr(x,x_1,x_2,N,N_w_mod,N_w,DesignMethod,Smoother)
%	Perform texture attribute extraction with optimal an filter and
%	a local energy function on the image x, with the texture regions 
%	(which the filter is optimized with respect to) defined by
%	m1, n1, m2 and n2.
%	DesignMethod may be 'er', 'mnsed', 'vsnsed', 'vsnsedg'
%	Smoother may be 'Gauss', 'Rectangular', 'Maximum filter', 'Minimum filter'
%	NaN are replaced by 0 (zero) prior to the filtering operation.
%
%	Example:
%	v=optexatr(z,z(351:450,331:430),z(501:600,1051:1150));


WaitbarHndl=waitbar(5,'Optimal texture feature extraction',100);
% Get the size of the input
[xM,xN]=size(x);

if strcmp(DesignMethod,'er') | strcmp(DesignMethod,'mnsed'),
  MaxLag=N;
else
  MaxLag=N+N_w_mod;
end

% Check the size of the two training regions
if (size(x_1,1) < MaxLag) | (size(x_1,2) < MaxLag),
  close(WaitbarHndl);
  error('Region 1 too small');
end
if (size(x_2,1) < MaxLag) | (size(x_2,2) < MaxLag),
  close(WaitbarHndl);
  error('Region 2 too small');
end

% Replace all NaNs by zeros, and retain a matrix with zeros where x have
% values, and NaNs where x have NaNs
NaNs=x-x; 
x=NaN2Val(x,0);

% Determine the autocorrelation functions
waitbar(10);
r_1=corrcoeff(x_1,x_1,MaxLag);
waitbar(30);
r_2=corrcoeff(x_2,x_2,MaxLag);

% Optimal filter design
waitbar(50);
if strcmp(DesignMethod,'er'),
  [h,J]=er_opt(x_1,x_2,N,r_1,r_2);
elseif strcmp(DesignMethod,'mnsed'),
  [h,J]=mnsed_opt(x_1,x_2,N,r_1,r_2);
elseif strcmp(DesignMethod,'vsnsed'),
  [h,J]=vsnsed_opt(x_1,x_2,N,N_w_mod,r_1,r_2);
elseif strcmp(DesignMethod,'vsnsedg'),
  %  [h,J]=vsnsed_opt(x_1,x_2,N,N_w_mod,r_1,r_2);
  %  [h,J]=vsnsed_gopt(h,x_1,x_2,N,N_w_mod,r_1,r_2);
  close(WaitbarHndl);
  error('Gradient search not supported yet');
end
waitbar(70);

% Symmetric extension
top=N+N_w; bot=N+N_w; lft=N+N_w; rgt=N+N_w;
x=symmext(x,top,bot,lft,rgt);
waitbar(72);

% Feature extraction
y=filter2(h,x);
waitbar(85);

y=y.^2;
if strcmp(Smoother,'Gauss'),
  y=gaussmth(y);
elseif strcmp(Smoother,'Rectangular'),
  y=avgsmth(y,N_w);
elseif strcmp(Smoother,'Maximum filter'),
  y=maxfilt2(y,[N_w N_w]);
elseif strcmp(Smoother,'Minimum filter'),
  y=minfilt2(y,[N_w N_w]);
else
  close(WaitbarHndl);
  error('Unsupported smoothing operator requested');
end
waitbar(100);

% Remove the symmetric extension
y=y(top+1:top+xM,lft+1:lft+xN);

% Put NaN's where x had NaN's
y=NaNs + y;
close(WaitbarHndl);