% Statistical and Computational Methods in Brain Image Analysis
% Chapter 3. General Linear Models
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu

% data is in the path /CCdensity
path(path, './CCdensity');

% Subject identifiers in the directory
c1=[01 02 05 08 09]'
c2=[10 11 12 13 14 15 16]'

a1=[01 02 03 04 05 06 07 08 09]'
a2=[10 11 13 14 16 17 18]'

% Image names for control subjects
file_c1=strcat('CCcontrol0',int2str(c1),'.txt')
file_c2=strcat('CCcontrol',int2str(c2),'.txt')
file_c=[file_c1; file_c2]

% Image names for autistic subjects
file_a1=strcat('CCautism0',int2str(a1),'.txt')
file_a2=strcat('CCautism',int2str(a2),'.txt')
file_a=[file_a1; file_a2]

% Subject numbers
% 12 controls
n_c=size(file_c,1)

% 16 autistic
n_a=size(file_a,1)

% Image dimension 95 x 68
d1=95; d2=68;

density_a=zeros(n_c,d1,d2);
density_a=zeros(n_a,d1,d2);

for i=1:n_c
    temp=load(file_c(i,:));
    temp=reshape(temp,d1,d2);
    density_c(i,:,:)=temp;
end;

for i=1:n_a
    temp=load(file_a(i,:));
    temp=reshape(temp,d1,d2);
    density_a(i,:,:)=temp;
end;

img=[];
for i=1:n_c 
    img = [img, squeeze(density_c(i,:,:))];
end;
figure; imagesc(img);colormap('bone');colorbar;

img=[];
for i=1:n_a
    img = [img, squeeze(density_a(i,:,:))];
end;
figure; imagesc(img);colormap('bone');colorbar

% general linar model 1:
% density = lambda0 + lambda1*age
% we test for the significance of the coefficient beta1

age_c = [15 18 18 16 15 13 18 15 21 17 16 23]'
age_a = [15 20 17 13 12 15 25 14 15 14 24 18 10 12 22 12]'
age=[age_c;age_a];

const=ones(n_c+n_a,1);
X=[const, age];

p_c=reshape(density_c,12,d1*d2);
p_a=reshape(density_a,16,d1*d2);
p=[p_c;p_a];

lambda=zeros(2,d1*d2);
for i=1:(d1*d2)
    lambda(:,i)=pinv(X)*p(:,i);
end
lambda0=reshape(lambda(1,:),d1,d2);
lambda1=reshape(lambda(2,:),d1,d2);

figure;imagesc(lambda0);colorbar; colormap('bone')
axis off
set(gcf,'Color','w')
figure;imagesc(lambda1);colorbar; colormap('bone')
axis off
set(gcf,'Color','w')

SSE0 = sum((p - (X*lambda)).^2);
SSE0=reshape(SSE0,d1,d2);

%figure;imagesc(SSE0);colorbar;
% density = lambda0 + lambda1*age + lambda2*group

group = [zeros(12,1);ones(16,1)]
lambda=zeros(3,d1*d2);
X=[const,age,group];

for i=1:(d1*d2)
    lambda(:,i)=pinv(X)*p(:,i);
end

SSE1 = sum((p - (X*lambda)).^2);
SSE1=reshape(SSE1,d1,d2);

F=25*(SSE0-SSE1)./SSE0;
imagesc(F);colorbar; colormap('bone')
axis off

pvalue=1 - fcdf(F,1,25);
figure;imagesc(-log10(pvalue));colorbar; colormap('bone')



% [beta,rho,pvalue]=regress_image(control,age_control);
% figure; imagesc(squeeze(beta(:,:,1))); colorbar;
% figure; imagesc(squeeze(beta(:,:,2))); colorbar;
% figure; imagesc(rho); colorbar;
% figure; imagesc(pvalue); colorbar;
% 
% [p1,p2]=image_threshold(pvalue,0.1);
% figure; imagesc(p1); colorbar;
% 
% 
% %autism
% [beta_a,rho_a,pvalue_a]=regress_image(autism,age_autism);
% figure; imagesc(squeeze(beta_a(:,:,1))); colorbar;
% figure; imagesc(squeeze(beta_a(:,:,2))); colorbar;
% figure; imagesc(rho_a); colorbar;


