% Test Script

% TO DO LIST:
% - Error messages and robustness of the methods

clear classes

%% Extracts data
phase = Phase; % creates a Phase object
phase.GetData; % extracts data

%% Visualize data
phase.ShowFrames(:,'cine Images',13);

%% Time Information
time = 0:1:phase.nFrames-1;
repetitionTime = (1/13); % time between samples
time = repetitionTime*time; % sampling moments in time

%% Instantaneous Phase Extraction (cardiac and respiratory)
phase.GetSignal; % gets mean intensity level of the images
phase.DisplayIntensity;
title('Mean Intensity Value');
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values
xlabel('Time (sec)');

phase.GetEMD; % decompose the siganl using EMD

Complot(phase.IMFs);
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values
title('Empirical Mode Decomposition');
xlabel('Time (sec)');
legend('IMF1','IMF2','IMF3','IMF4','IMF5','IMF6','IMF7');

phase.GetModes; % turn the modes into phase information

figure
Complot(phase.ModePhases)
title('IMFs Phase');
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values
xlabel('Time (sec)');

cardiac_phase = CardiacPhase(phase); % creates a CradiacPhase object
cardiac_phase.SelectPhase; % choose the cardiac imf

respiratory_phase = RespiratoryPhase(phase); % creates a RespiratoryPhase object
respiratory_phase.SelectPhase; % choose the respiratory imf

figure
cardiac_phase.DisplayValues
hold on
respiratory_phase.DisplayValues
title('Cardiac and Respiratory Phase');
legend('Cardiac Phase','Respiratory Phase');
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values

%% File with Cardiac and Respiratory Phase
file_matrix = [cardiac_phase.Values; respiratory_phase.Values];
csvwrite('phase_file.csv',file_matrix);

%% Instantaneous Frequency
cardiac_phase.GetFrequency(time);
respiratory_phase.GetFrequency(time);

%% Cardiac and Respiratory events
cardiac_phase.Diastole % detects diastolic events
cardiac_phase.Systole % detects systolic events
respiratory_phase.Expiration % detects expirtion moments
respiratory_phase.Inspiration % detects inspiration moments

SI = intersect(cardiac_phase.SystolePositions,respiratory_phase.InspirationPositions); % systole in inspiration
SE = intersect(cardiac_phase.SystolePositions,respiratory_phase.ExpirationPositions); % systole in expiration
DI = intersect(cardiac_phase.DiastolePositions,respiratory_phase.InspirationPositions); % diastole in inspiration
DE = intersect(cardiac_phase.DiastolePositions,respiratory_phase.ExpirationPositions); % diastole in expiration

% Show frames of the different cardiac/respiratory events
% Systole in Inpiration
phase.ShowFrames(SI,'End-systolic frames in inpiration',3);
directory = 'Systole in Inspiration'
mkdir(directory);
for i=1:length(SI)
   baseFileName = sprintf('SI%d.dcm', i);
   fullFileName = fullfile(directory, baseFileName); 
   dicomwrite(phase.Video(:,:,1,SI(i)), fullFileName);
end

figure
cardiac_phase.DisplayValues
hold on
respiratory_phase.DisplayValues
hold on
hline = refline([0 cardiac_phase.SystoleTreshold]); % Reference line for systole treshold phase
hline.Color = 'b';
hline.LineStyle = ':';
hold on
hline = refline([0 respiratory_phase.InspirationTreshold]); % Reference line for systole treshold phase
hline.Color = 'r';
hline.LineStyle = ':';
title('Systole in Inspiration');
legend('Cardiac Phase','Respiratory Phase');
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values

%Systole in Expiration
phase.ShowFrames(SE,'End-systolic frames in expiration',3);
directory = 'Systole in Expiration'
mkdir(directory);
for i=1:length(SE)
   baseFileName = sprintf('SE%d.dcm', i);
   fullFileName = fullfile(directory, baseFileName); 
   dicomwrite(phase.Video(:,:,1,SE(i)), fullFileName);
end

figure
cardiac_phase.DisplayValues
hold on
respiratory_phase.DisplayValues
hold on
hline = refline([0 cardiac_phase.SystoleTreshold]); % Reference line for systole treshold phase
hline.Color = 'b';
hline.LineStyle = ':';
hold on
hline = refline([0 respiratory_phase.ExpirationTreshold]); % Reference line for diastole treshold phase
hline.Color = 'r';
hline.LineStyle = ':';
title('Systole in Expiration');
legend('Cardiac Phase','Respiratory Phase');
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values

% Disatole in Inspiration
phase.ShowFrames(DI,'End-diastolic frames in inpiration',3);
directory = 'Diastole in Inspiration'
mkdir(directory);
for i=1:length(DI)
   baseFileName = sprintf('DI%d.dcm', i);
   fullFileName = fullfile(directory, baseFileName); 
   dicomwrite(phase.Video(:,:,1,DI(i)), fullFileName);
end

figure
cardiac_phase.DisplayValues
hold on
respiratory_phase.DisplayValues
hold on
hline = refline([0 cardiac_phase.DiastoleTreshold]); % Reference line for diastole treshold phase
hline.Color = 'b';
hline.LineStyle = ':'
hold on
hline = refline([0 respiratory_phase.InspirationTreshold]); % Reference line for systole treshold phase
hline.Color = 'r';
hline.LineStyle = ':';
title('Disatole in Inspiration');
legend('Cardiac Phase','Respiratory Phase');
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values

% Distole in Expiration
phase.ShowFrames(DE,'End-diastolic frames in expiration',3);
directory = 'Diastole in Expiration'
mkdir(directory);
for i=1:length(DE)
   baseFileName = sprintf('DE%d.dcm', i);
   fullFileName = fullfile(directory, baseFileName); 
   dicomwrite(phase.Video(:,:,1,DE(i)), fullFileName);
end

figure
cardiac_phase.DisplayValues
hold on
respiratory_phase.DisplayValues
hold on
hline = refline([0 cardiac_phase.DiastoleTreshold]); % Reference line for diastole treshold phase
hline.Color = 'b';
hline.LineStyle = ':'
hold on
hline = refline([0 respiratory_phase.ExpirationTreshold]); % Reference line for diastole treshold phase
hline.Color = 'r';
hline.LineStyle = ':';
title('Diastole in Expiration');
legend('Cardiac Phase','Respiratory Phase');
te=repetitionTime;
xt=get(gca,'xtick');
set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values


%% Heart Beat Reconstruction
n_ranges = 30; % number of phase ranges
ranges = linspace(respiratory_phase.MinRespPhase,respiratory_phase.MaxRespPhase,n_ranges+1); % splits the respiratory phase in n_ranges equal parts
for i=1:1:n_ranges
    Ltreshold = ranges(i);
    Utreshold = ranges(i+1);
    index(i,:) = respiratory_phase.Values < Utreshold & respiratory_phase.Values > Ltreshold; % the respiratory phase is within the upper and lower treshold of this range
    indexNumeric{i} = find(index(i,:)); % indexes of the frames satifying the previous condition
    if iscolumn(indexNumeric{i})
        indexNumeric{i} = indexNumeric{i}.'; % puts the indexes in a vector
    end
    values{i} = cardiac_phase.Values(indexNumeric{i}); % cardiac phase values of frames within this respiratory range
    % Sorts cardiac phase values and their indexes by ascending cardiac phase
    [ordered_values{i}, oredered_index{i}] = sort(values{i});
    k = indexNumeric{i};
    j = oredered_index{i};
    new_index{i} = k(j); % ordered indexes
end
selected_range = n_ranges-2; % optional range
phase.ShowFrames(new_index{selected_range},'Heart Beat Reconstruction',3);