function [newimages, driftmap] = mri_removeLinearDrift(images, seq);

% function [newimages, driftmap] = mri_removeLinearDrift(images, seq);
%
% Aug 2003, Farshad Moradi
% California Institute of Technology
%
% Fits a line for MRI signal change over time in each condition (for each
% voxel). Then subtracts the line form the data to remove the baseline change
%

if nargin<2,
    seq{1} = 1:length(images);
end;

n = 0; a = 0;
for i=1: length(seq);
    m = mean(seq{i});
    xy = 0; 
    for j=seq{i};
        xy = xy + (j-m)*images(j).data;
    end;
    x2 = sum((seq{i}-m).^2);
    a = a + length(seq{i})* xy/x2;
    n = n+length(seq{i});
end;

a = a / n; 

n = length(images);
bias = a*(n-1);
f = find(bias<0); bias(f) = 0;
bias = max(bias(:));

if nargout>=2,
    driftmap = a;
end;

for i=1: n,
    newimages(i) = images(i);
    newimages(i).data = images(i).data - a*(i-1) + bias;
end;