function t = mri_ttest(images, sequence, smooth, anat_thr, HR_resp_delay);

% function t = mri_ttest(images, sequence, smooth, anat_thr, HR_resp_delay);
%
% Aug 2003, Farshad Moradi
% California Institute of Technology
%
% Examples:
%   seq{1} = [3:10, 23:30, 43:50, 63:70];
%   seq{2} = [13:20, 33:40, 53:60, 73:80];
%   t = mri_ttest(images, seq, 1, 100);
%   mri_display(images(1).data, t);

N = length(images);
if nargin>4,
    for i=1:length(sequence),
        sequence{i} = sequence{i} + HR_resp_delay;
        f = find(sequence{i} <= N);
        sequence{i} = sequence{i}(f);
    end;
end;

n1 = length(sequence{1});
m1 = 0; sd1 = 0;
for i= sequence{1}
    m1 = m1 + images(i).data/n1;
    sd1 = sd1 + (images(i).data.^2)/n1;
end;

n2 = length(sequence{2});
m2 = 0; sd2 = 0;
for i= sequence{2}
    m2 = m2 + images(i).data/n2;
    sd2 = sd2 + (images(i).data.^2)/n2;
end;

if nargin<3, smooth = 0; end;
if nargin<4, anat_thr = -Inf; end;

sd1 = sd1-m1.^2;
sd2 = sd2-m2.^2;

if smooth,
	[xx,yy] = meshgrid(-3:1.5:3); b = exp(-(xx.^2+yy.^2)/2); b = b / sum(b(:));
	disp('smoothing images...');
    for i=1:size(m1,3),
		sd1(:,:,i) = conv2(sd1(:,:,i),b,'same');
		sd2(:,:,i) = conv2(sd2(:,:,i),b,'same');
		m1(:,:,i) = conv2(m1(:,:,i),b,'same');
		m2(:,:,i) = conv2(m2(:,:,i),b,'same');
    end;
end;

a = (n1+n2)/(n1*n2);
b = ((n1-1).*sd1 + (n2-1).*sd2 )/(n1+n2-2);
t = (m1-m2)./sqrt(a.*b);

% Threshold to remove noise in background:

f = find(m1<anat_thr); t(f) = 0;
f = find(m2<anat_thr); t(f) = 0;