function XN = lyngby_prep_detrend(X, varargin)

% lyngby_prep_detrend  - Detrending
%
%       function XN = lyngby_prep_detrend(X, 'PropertyName',
%           'PropertyValue')
%
%       Input:    X      Input data matrix
%
%       Property: Index  The indices (scans/objects) used in the
%                        estimation of the detrending parameters. By
%                        default all data is used. 
%                 Order  [ {1} | integer ] Order of the
%                        detrending. '1' is first order linear
%                        detrending 
%                 Run    Run structure. By default one run across the
%                        entire data set is assumed.
%
%       Output:   XN     Output preprocessed (detrended) data matrix.
%
%       Perform detrending. The detrending is performed independently
%       for each run (provided that the 'Run' property is set) and the
%       estimation of the detrending parameters can be confined to
%       some of the data indicated with the 'Index' property.
%
%       Ref: Goutte et al, "Feature-space clustering for fMRI
%            meta-analysis", Human Brain Mapping, 13(3):165-183.
%
%       See also LYNGBY, LYNGBY_NORMALIZE.  
%
% $Id: lyngby_prep_detrend.m,v 1.3 2003/01/30 13:27:25 fnielsen Exp $


    % Check arguments
    if nargin < 1
      error('Wrong number of arguments')
    end
    if ~isnumeric(X)
      error('Input argument ''X'' should be a matrix')
    end
    
    
    % Default properties
    index = [];
    order = 1;
    run = [];    
    
    % Parse properties
    n = 1;
    while n < nargin-1
      arg = lower(varargin{n}); 

      if strcmp(arg, 'index')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg) & prod(size(arg)) <= size(X,1) & max(arg(:)) <= size(X,1)
	  index = arg(:);
	else
	  error(sprintf([ ...
		'The argument with ''Input'' should be a vector of ' ...
		'integers representing indices in the columns of '...
		'''X''. X: %dx%d, index: %dx%d' ...
		], size(X), size(arg)));
	end

      elseif strcmp(arg, 'order')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg) & length(arg)==1 & arg>=1
	  order = arg;
	else
	  error(sprintf([ ...
		'The argument with ''Order'' should be an integer' ...
		]));
	end
	
      elseif strcmp(arg, 'run')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg) & prod(size(arg)) == size(X,1)
	  run = arg(:);
	else
	  error(sprintf([ ...
		'The argument with ''Run'' should be a vector ' ...
		'of integers with the size of X. X: %dx%d, Run: %dx%d.' ...
		], size(X), size(arg)));
	end
	
      else
	error('Invalid property');
      end
      n = n + 1;
    end

    
    if isempty(index)
      index = (1:size(X,1))';
    end
    if isempty(run)
      run = ones(size(X,1),1);
    end
    
    % Find the run identifiers
    urun = unique(run);

    % Allocated the result matrix
    XN = zeros(size(X));
    
    % Treat each run separately
    for r = urun
      
      % Find the objects that are from run 'r'
      indexRun  = find(run == r);
      indexScan = intersect(indexRun, index);

      % Construct design matrix G
      g = indexScan;
      G  = [ ones(length(g),1) , g ];
      Gall  = [ ones(length(indexRun),1) , indexRun ];
      o = order-1;
      while o,
        G = [ G , G(:,end) .* G(:,2) ];
        Gall = [ Gall , Gall(:,end) .* Gall(:,2) ];
	o = o - 1;
      end

      % Estimate detrending parameters (for this run)
      B = pinv(G) * X(indexScan, :);
      
      % Apply detrend for all scans within the run
      XN(indexRun,:) = X(indexRun,:) - Gall*B;
      
    end
















