%-----------------------------------
% ZPRACOVÁNÍ A ANALÝZA VIDEOSEKVENCÍ 
%-----------------------------------
% Diplomová práce: Bezkontaktní detekce fyziologických parametrů z obrazových sekvencí.
% Tereza Bršlicová, ID 133958
 
clc;
clear all;
close all;

obrazky = {};          % buňkové pole pro ukládání jednotlivých videosekvencí
fps = 29.1;            % počet snímků za sekundu (frame per second)
                       % tuto hodnotu je nutné změnit před spuštěním
                       % individuálně pro každého dobrovolníka zvlášť

           
%% Cyklus pro načtení jednotlivých snímků

for i=1:291            % počet zachycených snímků
    
   cislo=num2str(i);  
   
   if i<100 && i> 9
       cislo=['0' num2str(i)];
   
   elseif i<10
       cislo=['00' num2str(i)];
   end
   
%    cesta - např. E:\tereza_bršlicová_DP_přílohy\videosekvence
   
%    nazev = ['cesta\1\michalH0' cislo '.bmp'];        % 294
%    nazev = ['cesta\2\gabca0'    cislo '.bmp'];       % 284
%    nazev = ['cesta\3\david_20' cislo '.bmp'];        % 297
%    nazev = ['cesta\4\Mich0' cislo '.bmp'];           % 272
%    nazev = ['cesta\5\Mates0' cislo '.bmp'];          % 268
   nazev = ['C:\Users\terez_000\Desktop\tereza_bršlicová_DP_přílohy\videosekvence\6\Vojtisek0' cislo '.bmp'];       % 291
%    nazev = ['cesta\7\ja0' cislo '.bmp'];             % 297
%    nazev = ['cesta\8\Rene0' cislo '.bmp'];           % 269
%    nazev = ['cesta\9\denca0' cislo '.bmp'];          % 273
%    nazev = ['cesta\10\vasek_20' cislo '.bmp'];       % 278


%   Cyklus načte všchny videosekvence z daného měření do buňkového pole 'obrazky'

    obrazky{i}=imread(nazev);
    
end
   
   
%% Výběr oblasti ROI

imshow(obrazky{1},[])     % zobrazení prvního snímku ze série 
hold on
[a b] = ginput(2)         % dvě souřadnice pro levý horní a pravý dolní roh obdélníku


%% Získání signálů 
%  Z vybrané oblasti ROI je nyní získána průměrná jasová hodnotu pixelů 
%  v celé ploše zvolené oblasti ROI
%  Je proveden rozklad na tři RGB složky

for i=1:length(obrazky)
    
    cc=obrazky{i};
    ccr=double(cc(:,:,1));         % červený kanál  
    
   sigr(i)= mean(mean(ccr(b(1):b(2),a(1):a(2))));
   
   ccg=double(cc(:,:,2));          % zelený kanál
    
   sigg(i)= mean(mean(ccg(b(1):b(2),a(1):a(2))));
   
   ccb=double(cc(:,:,3));          % modrý kanál
    
   sigb(i)= mean(mean(ccb(b(1):b(2),a(1):a(2))));
    
    
end

% Byli vytvořeny 3 vektory sigr (červený kanál), sigb( modrý kanál),
% sigg(zelený kanál)
% Následuje jejich normalizace


close all

%% Normalizace signálů RGB

n_sigr = (sigr-mean(sigr))/std(sigr);   % R
n_sigg = (sigg-mean(sigg))/std(sigg);   % G
n_sigb = (sigb-mean(sigb))/std(sigb);   % B


%% Vykreslení signálů RGB 

figure;
subplot(3,1,1);
plot(linspace(0,10,length(sigr)),(sigr), 'r');
title('Signály RGB složek','FontWeight','bold','FontSize',11);
xlabel('Čas [s]');
ylabel('Průměrný jas')
subplot(3,1,2);
plot(linspace(0,10,length(sigg)),(sigg),'g' );
xlabel('Čas [s]');
ylabel('Průměrný jas')
subplot(3,1,3);
plot(linspace(0,10,length(sigb)),(sigb),'b');
xlabel('Čas [s]');
ylabel('Průměrný jas');


%% Vykreslení normalizovaných signálů RGB 
figure;
subplot(3,1,1);
plot(linspace(0,10,length(n_sigr)),(n_sigr), 'r');
xlabel('Čas [s]');
title('Signály normalizovaných RGB složek','FontWeight','bold','FontSize',11);
ylabel('Průměrný jas')
subplot(3,1,2);
plot(linspace(0,10,length(n_sigg)),(n_sigg),'g' );
xlabel('Čas [s]');
ylabel('Průměrný jas')
subplot(3,1,3);
plot(linspace(0,10,length(n_sigb)),(n_sigb),'b');
xlabel('Čas [s]');
ylabel('Průměrný jas');


%% ICA (Independent Component Analysis)

% funkce JADE algoritmu (volána z jade.m)
[~,Out1] = jade([n_sigr;n_sigg;n_sigb]); 

% vstup funkce jade - normalizované signály
% výspup funkce jade -  matice Out1 obsahující 3 sloupce, z nichž 
%                       každý reprezentuje jednu komponentu 

%% Vykreslení ICA komponent v závislosti na čase

figure;
subplot(3,1,1);
plot(linspace(0,10,length(Out1(1,:))),(Out1(1,:)),'k');
title('Komponenty ICA analýzy','FontWeight','bold','FontSize',11);
xlabel('Čas [s]');
ylabel('IC1','FontWeight','bold','FontSize',11);
subplot(3,1,2);
plot(linspace(0,10,length(Out1(2,:))),(Out1(2,:)),'k');
% title('IC2');
xlabel('Čas [s]');
ylabel('IC2','FontWeight','bold','FontSize',11);
subplot(3,1,3);
plot(linspace(0,10,length(Out1(3,:))),(Out1(3,:)),'k');
% title('IC3');
xlabel('Čas [s]');
ylabel('IC3','FontWeight','bold','FontSize',11);


%% Spektrální analýza pomocí FFT (rychlá Fourierova transformace)

presnost = 20;  % přesnost umožní vypočítat více spektrálních koeficientů

% FFT
spect1 = fft(Out1(1,:),length(sigr)*presnost); 
spect2 = fft(Out1(2,:),length(sigg)*presnost); 
spect3 = fft(Out1(3,:),length(sigb)*presnost); 

% Vykreslení výkonového spektra FFT
figure;
subplot(3,1,1);
plot(linspace(0,fps,length(spect1)),abs(spect1),'k');
xlabel('Frekvence [Hz]');
ylabel('IC1','FontWeight','bold','FontSize',11);
title('FFT','FontWeight','bold','FontSize',11);
subplot(3,1,2);
plot(linspace(0,fps,length(spect2)),abs(spect2),'k');
xlabel('Frekvence [Hz]');
ylabel('IC2','FontWeight','bold','FontSize',11);
subplot(3,1,3);
plot(linspace(0,fps,length(spect3)),abs(spect3),'k');
xlabel('Frekvence [Hz]');
ylabel('IC3','FontWeight','bold','FontSize',11);

% FFT s Hannovým oknem
w = hann(length(Out1(1,:)))';  % Hannovo váhovací okno

wspect1 = fft(w.*Out1(1,:),length(sigr)*presnost); 
wspect2 = fft(w.*Out1(2,:),length(sigg)*presnost); 
wspect3 = fft(w.*Out1(3,:),length(sigb)*presnost);

% Vykreslení výkonového spektra FFT s Hannovým oknem
figure;
subplot(3,1,1);
plot(linspace(0,fps,length(wspect1)),abs(wspect1),'k');
xlabel('Frekvence [Hz]');
ylabel('IC1','FontWeight','bold','FontSize',12);
title('FFT s Hannovým oknem','FontWeight','bold','FontSize',11);
subplot(3,1,2);
plot(linspace(0,fps,length(wspect2)),abs(wspect2),'k');
xlabel('Frekvence [Hz]');
ylabel('IC2','FontWeight','bold','FontSize',12);
subplot(3,1,3);
plot(linspace(0,fps,length(wspect3)),abs(wspect3),'k');
xlabel('Frekvence [Hz]');
ylabel('IC3','FontWeight','bold','FontSize',12);
 