function [Y] = lyngby_histeq(X, varargin)

% lyngby_histeq	       - Perform histogram equalization
%
%       function [Y] = lyngby_histeq(X, 'PropertyName', PropertyValue)
%
%       Input:	  X         Column vector to be histogram equalized
%
%       Property: Bins      {256} Number of bins 
%                 Reference { not defined } Vector which histogram
%                           should be used as a reference for the
%                           equalization 
%
%       Output:	  Y         The equalized vector
%
%       Without the 'Reference' vector the histogram equalization will be
%       uniform.
%       'Tied ranks' cannot be handle proberly: if two or values are
%       the same the sorting will be ambigouos.
%
%       See also LYNGBY.
%
% $Id: lyngby_histeq.m,v 1.4 2002/03/12 02:01:02 fnielsen Exp $


    % Default values
    bins = 256; 
    R = [];       

    n = 1;
    while n < nargin-1
      arg = lower(varargin{n});
      
      if strcmp(arg, 'bins')
	n = n + 1;
	arg = varargin{n};
	if isnumeric(arg) 
	  if length(arg) == 1 & arg > 0  
	    bins = arg;
	  else
	    error(sprintf(['Property value to ''Bins'' should be an ' ...
		  'integer. Length: %d'], length(arg)));
	  end
	else
	  error(sprintf(['Property value to ''Bins'' should be a ' ...
		'value. Is a string: ''%s'''], arg));
	end
	
      elseif strcmp(arg, 'reference')
	n = n + 1;
	arg = varargin{n};
	if isnumeric(arg) 
	  if size(X,2) == size(arg,2)
	    R = arg;
	  else
	    error(sprintf([ 'Property value to ''Reference'' should be ' ...
		  'of the same size as X. X: %dx%d, Ref: %dx%d'], ...
		size(X), size(arg)));
	  end
	else
	  error(sprintf(['Property value to ''Reference'' should be a ' ...
		'value. Is a string: ''%s'''], arg));
	end
	
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end
    
    if isempty(R)
      % Uniform 
      
      Xn = round((X - min(X)) * (bins-1) / (max(X) - min(X))) + 1;
      N = cumsum(hist(Xn, bins),2)';
      Nn = round((N - min(N)) * (bins-1) / (max(N) - min(N)));
      Y = Nn(Xn);
      
    else
      % Target histogram
      
      histR = hist(R,bins);
      [dummy,I] = sort(X,2);

      Y = X;
      for n = 1:size(X,1)
	i1 = 1;
	for j = 1:bins
	  i2 = histR(j) + i1 - 1;
	  Y(n,I(n,i1:i2)) = ones(size(I(n,i1:i2))) * j;
	  i1 = i2 + 1;
	end
      end
    end


    
    
    














