%-----------------------------------
%computing two sample t-statistic

% Group 1
control1=[
    
/home/chung/Jasmine/control/snI001S3_flip_seg1.img
/home/chung/Jasmine/control/snI002S3_flip_seg1.img
/home/chung/Jasmine/control/snI005S3_flip_seg1.img
/home/chung/Jasmine/control/snI008S3_flip_seg1.img
/home/chung/Jasmine/control/snI009S3_flip.img
/home/chung/Jasmine/control/snI010S3_flip.img
/home/chung/Jasmine/control/snI011S3_flip.img
/home/chung/Jasmine/control/snI012S3_flip.img
/home/chung/Jasmine/control/snI013S3_flip.img
/home/chung/Jasmine/control/snI014S3_flip.img
/home/chung/Jasmine/control/snI015S3_flip.img
/home/chung/Jasmine/control/snI016S3_flip.img];


%----------------------------------
% Group 2

autism=[
    /home/chung/Jasmine/autism/I001S3_flip.img
/home/chung/Jasmine/autism/I002S3_flip.img
/home/chung/Jasmine/autism/I003S3_flip.img
/home/chung/Jasmine/autism/I004S3_flip.img
/home/chung/Jasmine/autism/I005S3_flip.img
/home/chung/Jasmine/autism/I006S3_flip.img
/home/chung/Jasmine/autism/I007S3_flip.img
/home/chung/Jasmine/autism/I008S3_flip.img
/home/chung/Jasmine/autism/I009S3_flip.img
/home/chung/Jasmine/autism/I010S3_flip.img
/home/chung/Jasmine/autism/I011S3_flip.img
% Subject 12 will not be inlcuded as the primary diagnosis is fetal alcohol syndrome.
/home/chung/Jasmine/autism/I013S3_flip.img
/home/chung/Jasmine/autism/I014S3_flip.img
/home/chung/Jasmine/autism/I016S3_flip.img
/home/chung/Jasmine/autism/I017S3_flip.img
/home/chung/Jasmine/autism/I018S3_flip.img]


%seg1 = gray, seg2=white

directory='/study/morsub/5mmFWHM/';
id_control1=strcat(directory,'snE',int2str(control1),'S3_seg2');
id_control2=strcat(directory,'snE',int2str(control2),'S3_seg2');

id_def1=strcat(directory, 'snE', int2str(def1), 'S3_seg2');
id_def2=strcat(directory, 'snE', int2str(def2), 'S3_seg2');

i=read_analyze('/study/morsub/5mmFWHM/snE10002S3_seg2');
i=read_analyze(id_control1(1,:));
n_control=size(control1,1) +size(control2,1);
n_def=size(def1,1) + size(def2,1);
dim=i.dim;

%imagesc(i.data(:,:,40));

%computing mean
mean_control= zeros(dim(1),dim(2),dim(3));
for i=1:size(control1,1)
    i
    image=read_analyze(id_control1(i,:));
    mean_control=mean_control + image.data;
end;
for i=1:size(control2,1)
    i
    image=read_analyze(id_control2(i,:));
    mean_control=mean_control + image.data;
end;

mean_control=mean_control/(n_control);
%imagesc(mean_control(:,:,40))

mean_def=zeros(dim(1),dim(2),dim(3));
for i=1:size(def1,1)
    i
    image=read_analyze(id_def1(i,:));
    mean_def=mean_def + image.data;
end;

for i=1:size(def2,1)
    i
    image=read_analyze(id_def2(i,:));
    mean_def=mean_def + image.data;
end;

mean_def=mean_def/n_def;

%computing variance
var_control=zeros(dim(1),dim(2),dim(3));

for i=1:size(control1,1)
    i
    image=read_analyze(id_control1(i,:));
    var_control=var_control + (image.data-mean_control).^2;
end;

for i=1:size(control2,1)
    i
    image=read_analyze(id_control2(i,:));
    var_control= var_control + (image.data-mean_control).^2;
end;

var_control=var_control/n_control;

var_def=zeros(dim(1),dim(2),dim(3));

for i=1:size(control1,1)
    i
    image=read_analyze(id_control1(i,:));
    var_def=var_def + (image.data-mean_control).^2;
end;

for i=1:size(control2,1)
    i
    image=read_analyze(id_control2(i,:));
    var_def= var_def + (image.data-mean_control).^2;
end;

var_def=var_def/n_control;

% group t test
den=(var_def/(n_def-1)+var_control/(n_control-1)).^(1/2);
den2=squeeze(reshape(den,dim(1)*dim(2)*dim(3),1));
for i=1:size(den2,1)
    if den2(i)==0,
        den2(i)=10000;
    end;
end;

den3=squeeze(reshape(den2,dim(1),dim(2),dim(3)));
t= (mean_def-mean_control)./den3;
df= n_def+ n_control-2;


%-----------final results--------------

save /study/morsub/matlab/white_mc_t_5mm.mat t

%load /study/morsub/matlab/white_mc_t.mat;

tem=read_analyze('/home/chung/projects/morsub/avg305T1');
%tem=read_analyze('/study/morsub/MNItemplate/navg305T1_flip_ap');
imagesc(squeeze(tem.data(:,:,34)))
colormap(gray);



imagesc(squeeze(t(:,:,43)));
colorbar
