function [XN, meanImage, stdImage, meanSeq,stdSeq] = ...
    lyngby_normalize(X, arg1, arg2, arg3, arg4, arg5, arg6, arg7, ...
    arg8, arg9, arg10);

% lyngby_normalize     - Normalize (Preprocess) a datamatrix
%
%       function [XN, meanImage, stdImage, meanSeq, stdSeq] =     
%                    lyngby_normalize(X, PropertyName, PropertyValue);
%
%       Input:    X   Datamatrix. 
%         
%       Property: Centering       [ 0 | {1} ] Extract the sequence
%                                 mean. This will give a mean image
%                 RunCentering    [ {0} | 1 ] Extract the mean within
%                                 each run
%                 Run             Run specification
%                 ImageCentering  [ {0} | 1 ] Extract the image
%                                 mean. This will give the mean sequence
%                 Normalization   [ {0} | 1 ] Divide by the standard
%                                 deviation 
%	
%       Output:   XN         Output matrix
%                 meanImage  Mean values in each voxel
%                 stdImage   Standard deviation in each voxel
%                 meanSeq    Mean sequence
%                 stdSeq     Standard deviation in each image
%
%       The function performs centering (extract the mean) for a
%       datamatrix. This can either be done on the whole data set
%       ('Centering') or for each run separately ('RunCentering').
%       The mean of each image can also be extracted with
%       'ImageCentering'.
%
%       The data can be adjusted to common range by dividing by the
%       standard deviation. The standard diviation (stdImage and
%       stdSeq) is computed by normalizing with N rather than N-1.
%
%       See also: LYNGBY, LYNGBY_UI_PREPROC_NORMALIZE.
%
% $Id: lyngby_normalize.m,v 1.24 2002/04/22 09:54:32 fnielsen Exp $


    % Sizes
    [rX, cX] = size(X);
    
    % Default values
    bCentering = 1;
    bRunCentering = 0;
    bNormalization = 0;
    bImageCentering = 0;
    bWriteToGlobal = 0;
    R = [];
    
    % Parse arguments
    n = 1;
    while n < nargin
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'centering')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if arg == 0
	  bCentering = 0;
	else
	  bCentering = 1;
	end
      elseif strcmp(arg, 'runcentering')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if arg == 0
	  bRunCentering = 0;
	else
	  bRunCentering = 1;
	end
      elseif strcmp(arg, 'run')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if length(arg) == rX | isempty(arg)
	  R = arg;
	else
	  error(sprintf(['The run specification and the datamatrix ' ...
	      'do not have the same length: Run: %dx%d, datamatrix: ' ...
	      '%dx%d'], size(arg), rX, cX));
	end
      elseif strcmp(arg, 'imagecentering')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if arg == 0
	  bImageCentering = 0;
	else
	  bImageCentering = 1;
	end
      elseif strcmp(arg, 'normalization')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if arg == 0
	  bNormalization = 0;
	else
	  bNormalization = 1;
	end
      else
	error('Invalid property');
      end
      n = n + 1;
    end

    if bRunCentering
      if isempty(R)
	error('Run centering is on, but there is no run structure');
      end
    end

    [rX, cX] = size(X);
    XN = zeros(rX, cX);

    lyngby_log('Computing mean image ...');
    meanImage = mean(X,1); 
    lyngby_log('Computing std image ...');
    if 0
      % This is a memory consuming operation
      stdImage = std(X,0,1); 
    else
      % Less memory
      stdImage = zeros(1,cX);
      blocks = ceil(cX*rX / 10^5);
      indices = unique(round((0:1/blocks:1) * cX));
      for n = 1:length(indices)-1
	index = indices(n)+1:indices(n+1);
	xc = X(:, index) - repmat(meanImage(index), rX, 1);
	stdImage(index) = sqrt(sum(xc.*xc, 1) / rX);
	lyngby_log(sprintf('Computing std image (%.0d%%)', ...
	    round(indices(n+1)/cX*100)));
      end
    end
    lyngby_log('Computing mean sequence ...');
    meanSeq = mean(X,2);
    lyngby_log('Computing std sequence ...');
    if 0 
      % This is a memory consuming operation
      stdSeq = std(X,0,2);
    else
      % Less memory
      stdSeq = zeros(rX,1);
      blocks = ceil(cX*rX / 10^5);
      indices = unique(round((0:1/blocks:1) * rX));
      for n = 1:length(indices)-1
	index = indices(n)+1:indices(n+1);
	xc = X(index, :) - repmat(meanSeq(index), 1, cX);
	stdSeq(index) = sqrt(sum(xc.*xc, 2) / cX);
	lyngby_log(sprintf('Computing std squence (%.0d%%)', ...
	    round(indices(n+1)/rX*100)));
      end
    end
    
    if bCentering & ~bRunCentering
      lyngby_log('Centering');
      if 0
	% This is a memory consuming operation
	XN = X - ones(rX, 1) * meanImage;
      else
	% Less memory
	XN = zeros(rX,cX);
	blocks = ceil(cX*rX / 10^5);
	indices = unique(round((0:1/blocks:1) * cX));
	for n = 1:length(indices)-1
	  index = indices(n)+1:indices(n+1);
	  XN(:,index) = X(:, index) - repmat(meanImage(index), rX, 1);
	  lyngby_log(sprintf('Centering (%.0d%%)', ...
	      round(indices(n+1)/cX*100)));
	end
      end
	
    elseif bRunCentering
      Rcut = R;
      n = 1;
      while length(Rcut) > 0
	lyngby_log(sprintf('Run centering: %d. run', n));
	run = Rcut(1);
	Rcut = Rcut(find(Rcut~=run));
	indices = find(R == run);
	runMean(n,:) = mean(X(indices,:));
	XN(indices,:) = X(indices,:) - ...
	    ones(length(indices),1) * runMean(n,:);
	n = n + 1;
      end
      
    else
      XN = X;
      
    end

    if bImageCentering
      lyngby_log('Image centering ...');
      XN = XN - meanSeq * ones(1, cX);
    end

    if bNormalization
      lyngby_log('Normalization ...');
      XN = XN ./ repmat(stdImage + realmin, rX, 1);
    end







