% Skript s ukazkami redukce a klasifikace obrazovych dat mozku.
%
% by Eva Koritakova (Janousova), 11/09/2014, modifikovano 30/10/2018.


%% 0. NACTENI A PREDZPRACOVANI DAT.

% Nastaveni cesty k souborum:
cd 'c:\Users\koritakova\Documents\Vyuka\Vicerozmerky a AKD\Analyza a klasifikace dat\2018-podzim\AKD_predn-07_cviceni\Data_klasifikace\'

% Nacteni souboru se jmeny obrazovych souboru subjektu a zarazenim do skupin:
[num,text]=xlsread('gm_identifikatory.xls','List1'); % num...ciselne hodnoty, text...textove hodnoty
names = text(2:end,1); % jmena souboru s obrazy jednotlivych subjektu
group = num(:,1); % zarazeni subjektu do skupin
inxP = find(group==1); % identifikatory pacientu
inxC = find(group==0); % identifikatory kontrol
clear num text % odstraneni nepotrebnych promennych

% Nacteni masky sede hmoty mozku (GM) - pro odstraneni mimomozkovych a dalsich sktruktur:
maska=load('gm_mask_2D.mat'); % nacteni masky
maska=maska.gm_mask_2D; % vyber matice s maskou ze strukturovaneho pole
zobraz(maska) % zobrazeni masky GM

% Nacteni 2-D obrazu prvniho pacienta na ukazku:
i=1;
img=load([cell2mat(names(inxP(i),1)),'_rez50.mat']); % nacteni obrazu
zobraz(img.rez) % zobrazeni 2-D GM obrazu prvniho pacienta

% Nacteni vsech obrazovych dat (2-D rezy obrazu mozku):
X=zeros(length(inxP)+length(inxC),sum(sum(maska==1))); % alokace datove matice, kde v radcich budou odmaskovane obrazy jednotlivych subjektu
for i=1:length(inxP) % pro pacienty:
    img=load([cell2mat(names(inxP(i),1)),'_rez50.mat']); % nacteni obrazu pacienta
    X(i,:)=img.rez(maska==1); % odmaskovani obrazu a ulozeni do radku matice X
end
for i=1:length(inxC) % pro kontroly:
    img=load([cell2mat(names(inxC(i),1)),'_rez50.mat']); % nacteni obrazu kontrolniho subjektu
    X(i+length(inxP),:)=img.rez(maska==1); % odmaskovani obrazu a ulozeni do radku matice X
end



%% 1. REKDUKCE A KLASIFIKACE DAT - VALIDACE POMOCI RESUBSTITUCE.
% (tzn. trenovaci a testovaci soubor je stejny)

% redukce dat vyberem statisticky vyznamnych pixelu pomoci Mannova-Whitneyova testu:
pvals=zeros(1,size(X,2)); % alokace vektoru s p-hodnotami pro jednotlive pixely
for i=1:size(X,2) % pro vsechny pixely:
    pvals(i)=ranksum(X(group==1,i),X(group==0,i)); % vypocet p-hodnoty pomoci Mannova-Whitneyova testu - srovnani hodnot pacientu a kontrol v danem pixelu
end
sum(pvals<0.05) % pocet statisticky vyznammych pixelu
[pID,pN] = FDR(pvals,0.05); % False Discovery Rate (FDR) korekce
sum(pvals<pID) % pocet statisticky vyznammych pixelu po FDR korekci
X_red=X(:,pvals<pID); % redukce matice X vyberem statisticky vyznamych pixelu po FDR korekci

% zobrazeni statisticky vyznamych pixelu:
vyzn=maska;
vyzn(maska==1)=1+(pvals<pID);
zobraz(vyzn)

% klasifikace redukovanych obrazovych dat:
SVMStruct = fitcsvm(X_red,group);
group_klasif = predict(SVMStruct, X_red); % klasifikace dat pomoci SVM klasifikatoru
group_klasif=group_klasif(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 89.6%; 87.5%; 91.7%

% klasifikace redukovanych obrazovych dat pri C=100:
SVMStruct = fitcsvm(X_red,group,'BoxConstraint',100); % BoxConstraint odpovida parametru C (vychozi je C=1)
group_klasif = predict(SVMStruct, X_red); % klasifikace dat pomoci SVM klasifikatoru
group_klasif=group_klasif(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 100%; 100%; 100%

% klasifikace redukovanych obrazovych dat pri C=0.1:
SVMStruct = fitcsvm(X_red,group,'BoxConstraint',0.1); % BoxConstraint odpovida parametru C (vychozi je C=1)
group_klasif = predict(SVMStruct, X_red); % klasifikace dat pomoci SVM klasifikatoru
group_klasif=group_klasif(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 85.4%; 79.2%; 91.7%

% ve starsich verzich Matlabu pouziti funkci svmtrain a svmclassify misto fitcsvm a predict:
%SVMStruct = svmtrain(X_red,group); % natrenovani SVM klasifikatoru (ve starsich verzich Matlabu)
%SVMStruct = svmtrain(X_red,group,'boxconstraint', C); %alternativne lze menit hodnotu parametru C (boxconstraint), vychozi je C=1
%group_klasif = svmclassify(SVMStruct, X_red); % klasifikace dat pomoci SVM klasifikatoru
%Poznamka: pri pouziti funkce svmtrain vyjdou pro C=1 jine vysledky nez pri pouziti
%          fitcsvm, protoze funkce svmtrain automaticky data standardizovala.



%% 2. REDUKCE A KLASIFIKACE DAT - VALIDACE POMOCI METODY LEAVE-ONE-OUT (NESPRAVNA VARIANTA)
% (nespravna varianta proto, ze krosvalidace je delana jen pro klasifikaci
% a ne i pro redukci dat)

% redukce dat vyberem statisticky vyznamnych pixelu pomoci Mannova-Whitneyova testu:
pvals=zeros(1,size(X,2)); % alokace vektoru s p-hodnotami pro jednotlive pixely
for i=1:size(X,2) % pro vsechny pixely:
    pvals(i)=ranksum(X(group==1,i),X(group==0,i)); % vypocet p-hodnoty pomoci Mannova-Whitneyova testu - srovnani hodnot pacientu a kontrol v danem pixelu
end
sum(pvals<0.05) % pocet statisticky vyznammych pixelu
[pID,pN] = FDR(pvals,0.05); % False Discovery Rate (FDR) korekce
sum(pvals<pID) % pocet statisticky vyznammych pixelu po FDR korekci
X_red=X(:,pvals<pID); % redukce matice X vyberem statisticky vyznamych pixelu po FDR korekci

% klasifikace redukovanych obrazovych dat:
group_klasif2=zeros(1,length(group)); % alokace vektoru s vysledky klasifikace (tzn. indexy skupin urcenymi klasifikatorem)
for i=1:length(group)
    i
    x_test_red=X_red(i,:); % vyber testovaciho subjektu ze zredukovanych dat
    X_train_red=X_red(~ismember(1:length(group),i),:); % vyber trenovacich subjektu ze zredukovanych dat
    group_train=group(~ismember(1:length(group),i)); % vyber identifikatoru trenovacich subjektu
    SVMStruct = fitcsvm(X_train_red,group_train); % natrenovani SVM klasifikatoru na trenovacich datech (bez testovaciho subjektu)
    group_klasif2(i)=predict(SVMStruct, x_test_red); % klasifikace dat pomoci SVM klasifikatoru
end
group_klasif2=group_klasif2(:); % prevod vektoru s indexy skupin na sloupcovy vektor
spravnost = sum(group_klasif2 == group) / length(group) % vypocet spravnosti klasifikace
senzitivita = sum(group_klasif2(group==1) == group(group==1)) / length(group(group==1)) % vypocet senzitivity klasifikace
specificita = sum(group_klasif2(group==0) == group(group==0)) / length(group(group==0)) % vypocet specificity klasifikace
% 79.2%; 79.2%; 79.2%



%% 3. REDUKCE A KLASIFIKACE DAT - VALIDACE POMOCI METODY LEAVE-ONE-OUT (SPRAVNA VARIANTA)

% UKOL - vypocitejte spravnost, senzitivitu a specificitu pri pouziti
%        spravne krizove validace
%      - zamyslete se nad tim, jakym zpusobem by bylo mozne vizualizovat 
%        statisticky vyznamne pixely pri pouziti spravne krizove validace
