% Chapter  Sparse Shape Representation
%
% (C) Moo K. Chung 2012
% Image registration and template construction was done by S.-G. Kim of SNU
 

% Visualizing amygdala and hippocampus template for MIDUS data set

load midus.mat

amygl=template{1,1}
amygr=template{1,2}
hippol=template{2,1}
hippor=template{2,2}


n=13
dispal= displacement{n,1,1};
dispar= displacement{n,1,2};
disphl= displacement{n,2,1};
disphr= displacement{n,2,2};

figure; %figure_quiver2(amygl, dispal, 10, 'red')
hold on; figure_wire(amygl,'black','red');
hold on; 
hold on; figure_wire(amygr,'black','blue');
hold on; figure_wire(hippol,'black','yellow');
hold on; figure_wire(hippor,'black','green');
legend('Left amygdala', 'Right amygala', 'Left hippocampus', 'Right hippocampus')
legend('boxoff')
hold on; figure_quiver2(amygl, dispal, 10, 'red')
hold on; figure_quiver2(sp_amygr, dispar, 10, 'blue')
hold on; figure_quiver2(sp_hippol, disphl, 10, 'yellow')
hold on; figure_quiver2(sp_hippor, disphr, 10, 'green')
view([-120 15])

print -dpdf -r300 disp13

% Visualizing displacement vector field for the 10th subject

dispal= displacement{10,1,1};
dispar= displacement{10,1,2};
disphl= displacement{10,2,1};
disphr= displacement{10,2,2};

figure; hold on;
figure_patch( sp_amygl, 'red',0.7)
hold on; figure_quiver2(sp_amygl, dispal, 10)
figure_patch( sp_amygr, 'blue',0.7)
hold on; figure_quiver2(sp_amygr, dispar, 10)
figure_patch( sp_hippol, 'yellow',0.7)
hold on; figure_quiver2(sp_hippol, disphl, 10)
figure_patch( sp_hippor, 'green',0.7)
hold on; figure_quiver2(sp_hippor, disphr, 10)
view([-120 15])
axis off
  
%-----------------------------
% Eigenfunctions

[Aal, Cal] =FEM(amygl);
[Val, Dal] = eig(Cal,Aal,'chol');  
Val=Val(:,1:1000);
Dal=diag(Dal(:,1:1000));

[Aar, Car] =FEM(amygr);
[Var, Dar] = eig(Car,Aar, 'chol');  
Var=Var(:,1:1000);
Dar=diag(Dar(:,1:1000));

[Ahl, Chl] =FEM(hippol);
[Vhl, Dhl] = eig(Chl, Ahl, 'chol');
Vhl=Vhl(:,1:1000);
Dhl=diag(Dhl(:,1:1000));

[Ahr, Chr] =FEM(hippor);
[Vhr, Dhr] = eig(Chr, Ahr, 'chol');
Vhr=Vhr(:,1:1000);
Dhr=diag(Dhr(:,1:1000));

%---------------------------------
% Least Squares Estimation of surfaces

LSE_amygl=SPARSEsurf(amygl, Val, 0);
LSE_amygr=SPARSEsurf(amygr, Var, 0);
LSE_hippol=SPARSEsurf(hippol, Vhl, 0);
LSE_hippor=SPARSEsurf(hippor, Vhr, 0);

figure; figure_wire(LSE_amygl,'black','red');
 hold on; figure_wire(LSE_amygr,'black','blue');
hold on; figure_wire(LSE_hippol,'black','yellow');
hold on; figure_wire(LSE_hippor,'black','green');
view([-120 15])
print -dpdf -r300 midus_original_surface
 
%------------------------------
% Least squares estimation 100%
p=amygl.vertices

lambda = 1  
[b1,status]=l1_ls(Val, p(:,1),lambda);
[b2,status]=l1_ls(Val, p(:,2),lambda);
[b3,status]=l1_ls(Val, p(:,3),lambda);

sparsep=Val*[b1 b2 b3]
amyglsp.faces=amygl.faces;
amyglsp.vertices=sparsep;
figure; figure_wire(amyglsp,'black','green');
 
%---------------------------------
% Sparse Representation 5%
sp_amygl=SPARSEsurf(amygl, Val, 1);
sp_amygr=SPARSEsurf(amygr, Var, 1);
sp_hippol=SPARSEsurf(hippol, Vhl, 1);
sp_hippor=SPARSEsurf(hippor, Vhr, 1);

figure; figure_wire(sp_amygl,'black','red');
hold on; figure_wire(sp_amygr,'black','blue');
hold on; figure_wire(sp_hippol,'black','yellow');
hold on; figure_wire(sp_hippor,'black','green');
view([-120 15])



k=6;
figure; figure_surf( sp_amygl , Val(:,k))
hold on;  figure_surf( sp_amygr , Var(:,k))
hold on;  figure_surf( sp_hippol , Vhl(:,k))
hold on;  figure_surf( sp_hippor , Vhr(:,k))
view([-120 15])
caxis([-0.05 0.05])
print -dpdf -r300 midus_basis6

%--------------------------------------------

% Data reformatting
% It is necessary to reformat the data into a 3D array of subjects x
% vertices x coordinates. 


dispal=[];
for i=1:69
    dispal(i,:,:)=displacement{i,1,1} ;
end;

dispar=[];
for i=1:69
    dispar(i,:,:)=displacement{i,1,2};
end;

disphl=[];
for i=1:69
    disphl(i,:,:)=displacement{i,2,1};
end;

disphr=[];
for i=1:69
    disphr(i,:,:)=displacement{i,2,2};
end;

dispal2 = sqrt(sum(dispal.^2,3));
dispar2 = sqrt(sum(dispar.^2,3));
disphl2 = sqrt(sum(disphl.^2,3));
disphr2 = sqrt(sum(disphr.^2,3));

sp_dispal=SPARSEregress(dispal2, Val,1);
sp_dispar=SPARSEregress(dispar2, Var,1);
sp_disphl=SPARSEregress(disphl2, Vhl,1);
sp_disphr=SPARSEregress(disphr2, Vhr,1);

LSE_dispal=SPARSEregress(dispal2, Val,0);
LSE_dispar=SPARSEregress(dispar2, Var,0);
LSE_disphl=SPARSEregress(disphl2, Vhl,0);
LSE_disphr=SPARSEregress(disphr2, Vhr,0);

%--------------------------------------------------
% SPARSE T-test in SurfStat 
%Age & Gender effect

Gender=term(gender);
Age=term(age);
Brain=term(brain);
pvalue= 0.05

LSE_amygl.coord= sp_amygl.vertices';
LSE_amygl.tri=sp_amygl.faces';
slm1 = SurfStatLinMod( dispal2, 1+Brain + Gender+Age, LSE_amygl);
slm = SurfStatT( slm1, gender);
figure; figure_surf( sp_amygl ,slm.t)
threshold=randomfield_threshold(slm, pvalue)

LSE_amygr.coord= sp_amygr.vertices';
LSE_amygr.tri=sp_amygr.faces';
slm1 = SurfStatLinMod( sp_dispar', 1+Brain + Gender+Age, LSE_amygr);
slm = SurfStatT( slm1, gender);
hold on; figure_surf( sp_amygr ,slm.t)
threshold=randomfield_threshold(slm, pvalue)

LSE_hippol.coord= sp_hippol.vertices';
LSE_hippol.tri=sp_hippol.faces';
slm1 = SurfStatLinMod( sp_disphl', 1+Brain + Gender+Age, LSE_hippol);
slm = SurfStatT( slm1, gender);
hold on; figure_surf( sp_hippol ,slm.t)
thresholdhl=randomfield_threshold(slm, pvalue)

LSE_hippor.coord= sp_hippor.vertices';
LSE_hippor.tri=sp_hippor.faces';
slm1 = SurfStatLinMod( sp_disphr', 1+Brain + Gender+Age, LSE_hippor);
slm = SurfStatT( slm1, gender);
hold on; figure_surf( sp_hippor ,slm.t)
thresholdhr=randomfield_threshold(slm, pvalue)
colorbar('off')
colorbar('East')
alpha(0.7)

view([-120 15]); camlight
print -dpdf -r300 midus-tstat-gender1.pdf
view([180 45]); camlight
print -dpdf -r300 midus-tstat-gender2.pdf
view([0 -90]); camlight
print -dpdf -r300 midus-tstat-gender3.pdf
view([120 15]); camlight
print -dpdf -r300 midus-tstat-gender4.pdf

pvalue=[0.0001 0.001 0.01 0.05 0.1]
threshold=randomfield_threshold(slm, pvalue)
%5.4876    4.8117    4.0759    3.5001    3.2266 age/gender effect


%--------------------------------------
% LSE T-test 

Gender=term(gender);
Age=term(age);
Brain=term(brain);
pvalue= 0.05

sp_amygl.coord= sp_amygl.vertices';
sp_amygl.tri=sp_amygl.faces';
slm1 = SurfStatLinMod( LSE_dispal, 1+Brain + Gender+Age, sp_amygl);
slm = SurfStatT( slm1, age);
figure; figure_surf( sp_amygl ,slm.t)
threshold=randomfield_threshold(slm, pvalue)

sp_amygr.coord= sp_amygr.vertices';
sp_amygr.tri=sp_amygr.faces';
slm1 = SurfStatLinMod( LSE_dispar, 1+Brain + Gender+Age, sp_amygr);
slm = SurfStatT( slm1, age);
hold on; figure_surf( sp_amygr ,slm.t)
threshold=randomfield_threshold(slm, pvalue)

sp_hippol.coord= sp_hippol.vertices';
sp_hippol.tri=sp_hippol.faces';
slm1 = SurfStatLinMod( LSE_disphl, 1+Brain + Gender+Age, sp_hippol);
slm = SurfStatT( slm1, age);
hold on; figure_surf( sp_hippol ,slm.t)
thresholdhl=randomfield_threshold(slm, pvalue)

sp_hippor.coord= sp_hippor.vertices';
sp_hippor.tri=sp_hippor.faces';
slm1 = SurfStatLinMod( LSE_disphr, 1+Brain + Gender+Age, sp_hippor);
slm = SurfStatT( slm1, age);
hold on; figure_surf( sp_hippor ,slm.t)
thresholdhr=randomfield_threshold(slm, pvalue)

colorbar('off')
colorbar('East')
alpha(0.7)

view([-120 15]); camlight
print -dpdf -r300 midus-tstat-age1.pdf
view([180 45]); camlight
print -dpdf -r300 midus-tstat-age2.pdf
view([0 -90]); camlight
print -dpdf -r300 midus-tstat-age3.pdf
view([120 15]); camlight
print -dpdf -r300 midus-tstat-age4.pdf

pvalue=[0.0001 0.001 0.01 0.05 0.1]
threshold=randomfield_threshold(slm, pvalue)
%5.4876    4.8117    4.0759    3.5001    3.2266 age/gender effect

%---------------------------------------
% rejection region at 5% level. 

pvalue=0.05;

sp_amygl.coord= sp_amygl.vertices';
sp_amygl.tri=sp_amygl.faces';
slm1 = SurfStatLinMod( dispal2, 1+Brain + Gender+Age, sp_amygl);
slm = SurfStatT( slm1, age);
threshold=randomfield_threshold(slm, pvalue)
indal= find(slm.t>=threshold)
n=length(slm.t);
rejectal= zeros(1,n);
rejectal(indal)=slm.t(indal); 
figure; figure_surf( sp_amygl ,rejectal)

sp_amygr.coord= sp_amygr.vertices';
sp_amygr.tri=sp_amygr.faces';
slm1 = SurfStatLinMod( dispar2, 1+Brain + Gender+Age, sp_amygr);
slm = SurfStatT( slm1, age);
threshold=randomfield_threshold(slm, pvalue)
indar= find(slm.t>=threshold)
n=length(slm.t);
rejectar= zeros(1,n);
rejectar(indar)=slm.t(indar); 
hold on; figure_surf( sp_amygr ,rejectar)

sp_hippol.coord= sp_hippol.vertices';
sp_hippol.tri=sp_hippol.faces';
slm1 = SurfStatLinMod( disphl2, 1+Brain + Gender+Age, sp_hippol);
slm = SurfStatT( slm1, age);
thresholdhl=randomfield_threshold(slm, pvalue)
indhl= find(slm.t>=thresholdhl)
n=length(slm.t);
rejecthl= zeros(1,n);
rejecthl(indhl)=slm.t(indhl); 
hold on; figure_surf( sp_hippol ,rejecthl)

sp_hippor.coord= sp_hippor.vertices';
sp_hippor.tri=sp_hippor.faces';
slm1 = SurfStatLinMod( disphr2, 1+Brain + Gender+Age, sp_hippor);
slm = SurfStatT( slm1, age);
thresholdhr=randomfield_threshold(slm, pvalue)
indhr= find(slm.t>=thresholdhr)
n=length(slm.t);
rejecthr= zeros(1,n);
rejecthr(indhr)=slm.t(indhr); 
hold on; figure_surf( sp_hippor ,rejecthr)

colorbar('off')
colorbar('East')
alpha(0.7)

view([-120 15]); camlight
print -dpdf -r300 midus-reject-gender1.pdf
view([120 15]); camlight
print -dpdf -r300 midus-reject-gender2.pdf

view([180 45]); camlight
print -dpdf -r300 midus-tstat-gender3.pdf
view([0 -90]); camlight
print -dpdf -r300 midus-tstat-gender4.pdf
view([120 15]); camlight

%----------------------
% power computation for SPARSE

n_per=5000
Phl=[];
Phr=[];
  
for N=5:5:65
    
    counthl=0;
    counthr=0;
  
    for i=1:n_per
        p=randperm(69);
        pind=p(1:N);
        
        p_brain=brain(pind);
        p_gender=gender(pind);
        p_age=age(pind);
        p_Brain=term(p_brain);
        p_Gender=term(p_gender);
        p_Age=term(p_age);
        
        %     slm1 = SurfStatLinMod( sp_dispal(:,pind)', 1+p_Brain + p_Gender+p_Age, sp_amygl);
        %     slm = SurfStatT( slm1, p_gender);
        %     threshold=randomfield_threshold(slm, pvalue)
        %     indal= find(slm.t<=-threshold)
        %     n=length(slm.t);
        %     rejectal= zeros(1,n);
        %     rejectal(indal)=slm.t(indal);
        %     figure; figure_surf( sp_amygl ,rejectal)
        
        %     slm1 = SurfStatLinMod( sp_dispar', 1+Brain + Gender+Age, sp_amygr);
        %     slm = SurfStatT( slm1, gender);
        %     threshold=randomfield_threshold(slm, pvalue)
        %     indar= find(slm.t<=-threshold)
        %     n=length(slm.t);
        %     rejectar= zeros(1,n);
        %     rejectar(indar)=slm.t(indar);
        %     hold on; figure_surf( sp_amygr ,rejectar)
        
        
        slm1 = SurfStatLinMod( sp_disphl(:,pind)', 1+p_Brain + p_Gender+p_Age, sp_hippol);
        slm = SurfStatT( slm1, p_age);
        
        flag=find(slm.t(indhl)>=thresholdhl);
        if ~isempty(flag)
            counthl=counthl+1;
        end;

        
        slm1 = SurfStatLinMod( sp_disphr(:,pind)', 1+p_Brain + p_Gender+p_Age, sp_hippor);
        slm = SurfStatT( slm1, p_age);
        
        flag=find(slm.t(indhr)>=thresholdhr);
        if ~isempty(flag)
            counthr=counthr+1;
        end;
        
 
    end
    
    Phl(N/5)=counthl/n_per
    Phr(N/5)=counthr/n_per
    
end


%----------------------
% power computation for original data

n_per=5000;
Phl=[];
Phr=[];
  
for N=5:5:65
    
    counthl=0;
    counthr=0;
  
    for i=1:n_per
        p=randperm(69);
        pind=p(1:N);
        
        p_brain=brain(pind);
        p_gender=gender(pind);
        p_age=age(pind);
        p_Brain=term(p_brain);
        p_Gender=term(p_gender);
        p_Age=term(p_age);
        
        %     slm1 = SurfStatLinMod( sp_dispal(:,pind)', 1+p_Brain + p_Gender+p_Age, sp_amygl);
        %     slm = SurfStatT( slm1, p_gender);
        %     threshold=randomfield_threshold(slm, pvalue)
        %     indal= find(slm.t<=-threshold)
        %     n=length(slm.t);
        %     rejectal= zeros(1,n);
        %     rejectal(indal)=slm.t(indal);
        %     figure; figure_surf( sp_amygl ,rejectal)
        
        %     slm1 = SurfStatLinMod( sp_dispar', 1+Brain + Gender+Age, sp_amygr);
        %     slm = SurfStatT( slm1, gender);
        %     threshold=randomfield_threshold(slm, pvalue)
        %     indar= find(slm.t<=-threshold)
        %     n=length(slm.t);
        %     rejectar= zeros(1,n);
        %     rejectar(indar)=slm.t(indar);
        %     hold on; figure_surf( sp_amygr ,rejectar)
        
        
        slm1 = SurfStatLinMod( disphl2(pind,:), 1+p_Brain + p_Gender+p_Age, sp_hippol);
        slm = SurfStatT( slm1, p_age);
        
        flag=find(slm.t(indhl)>=thresholdhl);
        if ~isempty(flag)
            counthl=counthl+1;
        end;

        %threshold = 3.4661
        
        slm1 = SurfStatLinMod( disphr2(pind,:), 1+p_Brain + p_Gender+p_Age, sp_hippor);
        slm = SurfStatT( slm1, p_age);
        
        flag=find(slm.t(indhr)>=thresholdhr);
        if ~isempty(flag)
            counthr=counthr+1;
        end;
        
        %threshold = 3.5001
        
    end
    Phl(N/5)=counthl/n_per
    Phr(N/5)=counthr/n_per
end

% SPARSE result for 5000 iterations 
%Phl
sparsePhl=[
    0.8120
    0.4454
    0.4930
    0.6098
    0.7514
    0.8728
    0.9422
    0.9840
    0.9974
    0.9996
    1.0000
    1.0000
    1.0000]

%Phr
sparsePhr=[
    0.5768
    0.1996
    0.2088
    0.2628
    0.3486
    0.4474
    0.5616
    0.6974
    0.8208
    0.9208
    0.9846
    0.9988
    1.0000]

%ORIGINAL results for 5000 iterations

OPhl=[
    0.8324
    0.4440
    0.4906
    0.5890
    0.7248
    0.8342
    0.9274
    0.9770
    0.9972
    0.9998
    1.0000
    1.0000
    1.0000
    ]

OPhr=[
     0.6304
    0.2278
    0.2264
    0.2662
    0.3248
    0.4146
    0.5116
    0.6240
    0.7298
    0.8422
    0.9350
    0.9888
    1.0000
]

    %LSE result for 5000 iterations
    
    LSEPhl=[
        0.8424
    0.4764
    0.5096
    0.6140
    0.7368
    0.8572
    0.9388
    0.9800
    0.9974
    1.0000
    1.0000
    1.0000
    1.0000
    ]

LSEPhr=[
    0.6388
    0.2432
    0.2472
    0.2970
    0.3536
    0.4558
    0.5530
    0.6648
    0.7822
    0.8790
    0.9676
    0.9954
    1.0000
]

figure; 
subplot(1,2,1);
plot(10:5:65, OPhr(2:end),'--b', 'LineWidth',2)
hold on; plot(10:5:65, sparsePhr(2:end), 'r','LineWidth', 2)
title('Right Hippocampus')
legend('Original','Sparse', 'Location', 'SouthEast')
legend('boxoff')
xlabel('Sample Size')
ylabel('Power')
figure_bg('white')
grid on    
set(gca,'XTick',10:10:70, 'YTick', 0.2:0.1:1)
axis([10 70 0.2 1])
pbaspect([1 1 1])
    
subplot(1,2,2);
plot(10:5:65, OPhl(2:end),'--b', 'LineWidth',2)
hold on; plot(10:5:65, sparsePhl(2:end), 'r','LineWidth', 2)
title('Left Hippocampus')
legend('Original','Sparse', 'Location', 'SouthEast')
legend('boxoff')
xlabel('Sample Size')
ylabel('Power')
figure_bg('white')
grid on
set(gca,'XTick',10:10:70, 'YTick', 0.2:0.1:1)
axis([10 70 0.2 1])
pbaspect([1 1 1])

print -dpdf -r300 power.pdf


