function [y,t,c,e]=optresh(x,x1,x2,parameter)
%OPTRESH  Optimal threshold
%	y=optresh(x,x1,x2) determine an (sub) optimal threshold for
%	classifying x1 and x2 with thresholding, and apply that threhold to
%	the matrix x. The threshold is found by assuming that the error	is
%	monotonous with respect to the threshold.
%
%	y=optresh(x,x1,x2,'verbose') runs the optimization in verbose mode,
%	displaying intermediate results.  
%
%	[y,t]=optresh(x,x1,x2) also return the threshold value, t. 
%
%	[y,t,c]=optresh(x,x1,x2) return and a value of 1 or 2 in c,
%	indicating which of x1 and x2 will be larger than the threshold.
%
%	[y,t,c,e]=optresh(x,x1,x2) return the error percentage in e.

if nargin<4,
  verbose=0;
else
  if strcmp(lower(parameter),'verbose'),
    verbose=1;
  else
    verbose=0;
  end
end

if verbose,
  fprintf(1,' Threshold      Error\n');
end

% Prepare the inputs
x1=x1(:);
x2=x2(:);
% Make sure x1<x2
if mean(x1) > mean(x2), tmp=x1; x1=x2; x2=tmp; clear tmp; c=1;
else c=2; end

% Find initial threhold
zeros1=zeros(length(x1),1);
ones2 =ones (length(x2),1);
t_cand=[    x1(find(x1>=min(x2)));   x2(find(x2<max(x1)))];
c_cand=[zeros1(find(x1>=min(x2)));ones2(find(x2<max(x1)))];
[t_cand t_candi]=sort(t_cand);
c_cand=c_cand(t_candi);
%error=sum([x1>=t;x2<t]); 
%t=mean([max(x1) min(x2)]);
%[dummy,ti]=min(abs(t_cand-t)); t=t_cand(ti);
max_ti=length(t_cand);
min_ti=1;
if isempty(t_cand)
  continue=0;
  t=mean([min(x2) max(x1)]);
else
  continue=1;
end


while continue,
  ti=round(mean([min_ti max_ti]));
  if      sum(abs((t_cand>=t_cand(round(mean([ti min_ti])))) - c_cand)) ...
	> sum(abs((t_cand>=t_cand(round(mean([ti max_ti])))) - c_cand))
    min_ti=ti;
  else
    max_ti=ti;
  end
  if max_ti-min_ti<=1,
    continue=0;
  end
  t=t_cand(ti);
  e=100*(sum(x1>=t)+sum(x2<t))/length([x1;x2]);
  if verbose,
    fprintf(1,'%11.4e   %6.2f%%\n',t,e);
  end
end
y=x>=t;