% **************************************************************************************************
% File: vyska_jamky.m
% 
% Brief: subor obsahuje algoritmus spracovania signalu z experimentalnych merani 
%        a vypocet vysky jamky
% 
% Version: 1.0  Revision: 0.0   Name: peter novotny
% 
% Date: 2011-05-21
%
% Matlab version: R2009a
% **************************************************************************************************

clear all; close all; clc;

%% *************************************************************************************************
% Nacitanie experimentalne nameranych dat z textoveho suboru a odstrananie redundantnych dat.
% !! Tato cast kodu je kvoli velkosti .txt suborov vykonavana len raz a nasledne sa data nacitaju uz
% len z ulozenych .mat suborov a tato cast sa zakomentuje.!!
% z - je relativna hodnota posunutia zrkadla v meracej vetve (ziskana z mikroposuvu PZT)
% I_960xN - matica obsahujuca v riadkoch zaznamenane hodnoty relativnej intenzity pre vsetky pixely
%           s trojnasobnou redundanciou
% I_320xN - matica obsahujuca v riadkoch zaznamenane hodnoty relativnej intenzity pre vsetky pixely
%           po odstraneni trojnasobnej redundancie
% **************************************************************************************************

% % nacitanie textoveho suboru
% data = load('namerane_data\data_130511_1.txt');
% 
% % vyber hodnot odpovedajucich posunutiu posuvneho zrkadla a zistenie postu zaznamenanych hodnot N
% z = data(1:end, 1)';    
% N = length(z);
% 
% % vyber hodnot odpovedajucich zaznamenanym intenzitam pre vsetky pixely
% I_960xN = data(1:end, 3:end)';
% 
% % nadefinovanie rozmerov snimaca
% rows = 16; colls = 20;
% 
% % odstranenie redundantnych dat - pretoze su vzdy 3 susedne stlpce rovnake vyberie sa len kazdy treti.
% I_320xN = zeros(rows*colls, N);
% k = 1;
% for i=1:rows*colls*3
%     if ~mod(i,3)
%         I_320xN(k,:) = I_960xN(i,:);
%         k = k + 1;
%     end
% end
% 
% % vycistenie workspace od dalej nadbytocnych premennych
% clear data I_960xN;
% 
% % ulozenie pozadovanyh premennych do suboru s priponou .mat, ktory je podstatne mensi ako textovy
% % subor a bude sluzit pre nasledovne spracovanie dat
% save('namerane_data\data_130511_1.mat','I_320xN','z');

%% *************************************************************************************************
% Nacitanie experimentalne nameranych dat z vopred pripraveneho suboru .mat, ktory uz neobsahuje
% redundantne data. 
% rozmery snimaca a dlzka dat sa definuju opat, pretoze predosli kod sa vykona len raz a potom sa 
% zakomentuje kvoli velkosti textovych dat
% **************************************************************************************************

subor = 'namerane_data\data_130511_2.mat';
load(subor);

% nadefinovanie rozmerov snimaca
rows = 16; colls = 20;

% dlzka zaznamenanych dat (odpoveda poctu zaznamenanych snimkov z kamery)
N = length(z);

%% *************************************************************************************************
% Zobrazenie priebehu intenzit pri vertikalnom skenovani v zvolenom riadku a stlpci
% !! k zobrazeniu nutne odkomentovat 1x celu sekciu!
% **************************************************************************************************

% % vyber zobrazovaneho riadku a stlpca
% riadok = 8;
% stlpec = 10;
% 
% % jeden riadok snimaca bude mat dlzku podla poctu stlpcov
% I_v_riadku(:,:) = I_320xN(colls*(riadok - 1) + 1:colls*riadok, :)';
% 
% % jeden stlpec bude mat dlzku podla poctu riadkov
% I_v_stlpci = zeros(rows, N);
% % prejdu sa vsetky riadky z I_320xN a vyberu sa len data pre odpovedajuci stlpec
% 
% for i=1:rows 
%     I_v_stlpci(i,:) = I_320xN((colls*i) - colls + stlpec, :)';
% end
% 
% % figure; imshow(imresize(I_320xN, [640 1000]),[]);
% % zobrazenie 
% figure; imshow(imresize(I_v_riadku, [1200 640]),[]);
% title(strcat(num2str(riadok),'. riadok (',num2str(colls),' hodnot)'));
% xlabel('stlpce -->'); ylabel('cislo zaznamu -->'); 
% figure; imshow(imresize(I_v_stlpci, [600 1200]),[]);
% title(strcat(num2str(stlpec),'. stlpec (',num2str(rows),' hodnot)'));
% xlabel('cislo zaznamu -->'); ylabel('riadky -->'); 

%% *************************************************************************************************
% aproximacia hodnot posunutia "z" metodou najmensich stvorcov 
% **************************************************************************************************

% n urcuje rad aproximacneho polynomu
n = 2; 
% vypocet koeficientov aproximacneho polynomu
koef = polyfit(1:N, z, n); 
% vypocet hodnot aproximacneho polynomu z koeficientov
x = polyval(koef, 1:N);
% posunutie polynomu na zaciatocnu hodnotu 0 
x = x - x(1); 

% % graficke znazornenie aproximacie
% figure; plot(z-z(1),'b','linewidth', 2); grid on; hold on;
% plot(x,'r','linewidth', 2);

% krok merania (posuvania PZT) v [um] sa urci ako priemerny rozdiel medzi hodnotami posunutia
d = diff(x);
step = mean(d);

%% *************************************************************************************************
% Vyber algoritmu spracovania dat a jeh aplikacia postupne na kazdy riadok matice I_320xN
% oznacenie algoritmov pre vyber je nasledujuce:
% 'wc' = vahovanie stredu
% 'cc' = koherentna korelacia
% 'fda' = analyza vo frekvencnej oblasti
% Vysledna matica "stredy_pruzkov" (320x1) obsahuje informacie o vypocitanom strede pruzkov zvolenym
% algoritmom
% **************************************************************************************************

% vyber algoritmu
algoritmus = 'cc';

% predpriprava matice stredy_pruzkov (320x1) pre informacie o strede pruzkov
stredy_pruzkov = zeros(rows*colls,1);

% prechadzanie matice I_320xN po riadkoch a aplikacia vybraneho algoritmu 
for i=1:rows*colls
    % vyber jedneho riadku matice I_320xN (odpoveda jednemu pixelu snimaca kamery)
    I_1D = I_320xN(i,:);
    switch algoritmus
        case 'wc'
            run weight_center; 
        case 'cc'
            run coherent_correlation; 
        case 'fda'
            run frequency_domain_analysis;
        otherwise
            error('Nespravne zadany alebo neznamy algoritmus!');
    end
    stredy_pruzkov(i) = fringe_center;
end

%% *************************************************************************************************
% Prevedenie matice "stredy_pruzkov" z rozmermi 320x1 na maticu "profil_povrchu", ktora odpoveda 
% fyzickym rozmerom snimaca 16x20 
% **************************************************************************************************

profil_povrchu = zeros(rows,colls);

for i=1:rows
    % Naukladanie riadkov (s dlzkou odpovedajucou poctu stlpcov) postupne pod seba
    profil_povrchu(i,:) = stredy_pruzkov(colls*i-(colls-1):colls*i)';
end

% zobrazenie ziskaneho profilu povrchu 
% figure; surf(profil_povrchu(:,:));
% title('Vvyskovy profil casti povrchu jamky');
% xlabel('stlpce [-]'); ylabel('riadky [-]'); zlabel('vyska [um]');
% shading INTERP; % SHADING FLAT; SHADING FACETED;

%% *************************************************************************************************
% Nakolko nebola osvetlena rovnomerne cela plocha snimaca vyrata sa vyska len z jeho vybranej casti
% kde sa dalo osvetlenie povazovat za rovnomerne
% **************************************************************************************************

% vyber casti snimaca s dostatosne rovnomernym osvetlenim (na zaklade zobrazeneho grafu povrchu)
ubrat_v_riadkoch = 3;
ubrat_v_stlpcoch = 2;
vybrana_cast = profil_povrchu(1 + ubrat_v_riadkoch:end - ubrat_v_riadkoch + 2 ...
                             ,1 + ubrat_v_stlpcoch:end - ubrat_v_stlpcoch); 

% zobrazenie vybranej dobre osvetlenej casti 
figure; surf(vybrana_cast);
title('Vyskovy profil jamky z dobre osvetlenej casti snimaca');
xlabel('stlpce [-]'); ylabel('riadky [-]'); zlabel('Vyska [um]');
% shading INTERP;% FLAT; FACETED - default

%% *************************************************************************************************
% Zrovnanie vybranej plochy len podla riadkov nakolko stlpce obsahuju zlom ukazalo sa,
% ze rovnanie nema vyraznejsi vplyv na presnost
% **************************************************************************************************

% Aproximacia a zrovnanie cez riadky
pomocna_os_riadky = 1:size(vybrana_cast,2);

% Vyber riadku pomocou ktoreho sa bude zrovnavat
rovnany_riadok = vybrana_cast(1,:);

% vypocet aproximacnych koeficientov 
koef_zrovnanie_riadky = polyfit(pomocna_os_riadky,rovnany_riadok,1);
% a = polyval(koef_zrovnanie_riadky, pomocna_os_riadky);
b = koef_zrovnanie_riadky(1)*pomocna_os_riadky;

% zrovnanie odpocitanim prislusnych hodnot 
zrovnany_povrch = vybrana_cast - meshgrid(b, 1:size(vybrana_cast,1));

% % zobrazenie zrovnaneho povrchu
% figure; surf(zrovnany_povrch);
% title('zrovnany povrch'); xlabel('stlpce [-]'); ylabel('riadky [-]'); zlabel('vyska [um]');

%% *************************************************************************************************
% rozdelenie vybranej casti na vrchnu a spodnu plochu podla priemernej hodnoty
povrch_pre_vypocet = vybrana_cast;
vrch = povrch_pre_vypocet; 
spodok = povrch_pre_vypocet;
stred = mean(mean(povrch_pre_vypocet));
vrch(vrch < stred) = 0; 
spodok(spodok > stred) = 0;

% odstranenie hodnot moc vzdialenych medianu plochy
max_presah = 1; % [um]
vrch_median = median(median(vrch(end-1:end,:)));
spodok_median = median(median(spodok(1:2,:)));
vrch(vrch < (vrch_median - max_presah)) = 0; 
vrch(vrch > (vrch_median + max_presah)) = 0; 
spodok(spodok < (spodok_median - max_presah)) = 0; 
spodok(spodok > (spodok_median + max_presah)) = 0; 

% Vypocet vysky jamky ako rozdiel priemernej vysky spodnej a vrchnej casti 
% priemer sa zrata len z hodnot ktore niesu nulove
priemer_vrch = sum(sum(vrch))/(numel(vrch) - numel(find(vrch == 0)));
priemer_spodok = sum(sum(spodok))/((numel(spodok) - numel(find(spodok == 0))));

% Vysledna vyska jamky
hole_height = priemer_vrch - priemer_spodok;

vypis = strcat('Vyska jamky pre meranie c.',subor(end-4),' zistena algoritmom "',algoritmus, '" je rovna:');
disp(vypis); 
disp(hole_height);

% **************************************************************************************************
%  END OF FILE: vyska_jamky.m
% **************************************************************************************************