% Ukazkove prikazy pro analyzu obrazovych dat.
%
% Autor: RNDr. Eva Janousova
% Institut biostatistiky a analyz, Masarykova univerzita, Brno
%
% ---------------------------------------------------------------------- %

% ukazka rozzipovani vsech obrazu v dane slozce:
soubory=ls('c:\Users\janousova\Documents\Neurodata\PacientiAD_kontroly'); % vytvoreni seznamu vsech obrazu
soubory=soubory(3:end,:); % odstraneni prvnich 2 radku
for i=1:size(soubory,1) % postupne rozzipovani jednoho obrazu po druhem
    gunzip(['c:\Users\janousova\Documents\Neurodata\PacientiAD_kontroly\',strtrim(soubory(i,:))]);
end

% ukazka nacteni masky mozku:
brainmask=load_nii('c:\Users\janousova\Documents\Neurodata\brainmask.nii');
brain_inx=find(brainmask.img==1); % vyber voxelu, ktere lezi v mozku 

% ukazka nacteni obrazovych dat pacientu s Alzheimerovou chorobou a kontrolnich subjektu:
X=zeros(100,length(brain_inx)); % alokace matice, kde na radcich budou jednotlive subjekty 
group=[ones(50,1);2*ones(50,1)]; % vektor obsahujici identifikatory skupin (1..pacient, 2..kontrola)
for i=1:50
    img=load_nii(['c:\Users\janousova\Documents\Neurodata\pacient_',num2str(i),'.nii']);
    X(i,:)=img.img(brain_inx);
end
for i=51:100
    img=load_nii(['c:\Users\janousova\Documents\Neurodata\kontrola_',num2str(i),'.nii']);
    X(i,:)=img.img(brain_inx);
end

% ukazka nacteni kovariat:
load 'c:\Users\janousova\Documents\Neurodata\covariates.mat'

% ukazka nastaveni cesty k funkci pro FDR korekci a funkci pro vypocet Fisherova exaktniho testu:
cd 'c:\Users\janousova\Documents\Neurodata\' 

% ukazka vypoctu Mannova-Whitneyova testu pro kazdy voxel zvlast:
pval=zeros(1,size(X,2)); % alokace vektoru s p-hodnotami
for i=1:size(X,2)
    pval(i)=ranksum(X(group==1,i),X(group==2,i));
end
sum(pval<0.05) % pocet voxelu se statisticky vyznamnym vysledkem

% ukazka vizualizace statisticky vyznamnych voxelu:
image=brainmask.img;
image(brain_inx)=(pval<0.05);
imagesc(image(:,:,60)),colormap(gray) % vizualizace 60. rezu

% ukazka provedeni FDR korekce:
[pID,pN] = FDR(pval,0.05);
% pokud je pID prazdne, znamena to, ze po FDR korekci je pocet statisticky
% vyznamnych voxelu roven 0; pokud pID neni prazdne, spocitame:
sum(pval<pID) % pocet statisticky vyznamnych voxelu po FDR korekci

% ukazka provedeni Bonferroniho korekce:
sum(pval*length(pval)<0.05) % pocet statisticky vyznamnych voxelu po Bonferroniho korekci
% nebo: sum(pval<(0.05/length(pval)))

% ukazka vypoctu Spearmanovy korelace hodnot voxelu s vekem u pacientu:
pval=zeros(1,size(X,2));
for i=1:size(X,2)
    [rho,pval(i)]=corr(X(group==1,i),age(group==1),'type','Spearman');
end
% funkce corr umoznuje pocitat p-hodnoty i pro vsechny voxely soucasne:
[rho,pval]=corr(X(group==1,:),age(group==1),'type','Spearman');

% ukazka vypoctu Fisherova exaktniho testu:
pval=zeros(1,size(X,2));
for i=1:size(X,2)
    pval(i)=fexact(X(group==1,i),X(group==2,i));
end
