% Test Script

% TO DO LIST:
% - Heart beat reconstruction with steady respiratory phase
% - Multislice analysis
% - Choose IMF based on Instantaneous Frequency
% - Error messages and robustness of the methods

clear classes

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

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

%% Workload information
% Should become a set method
% User sets the workload, important to choose the correct IMF
list = {'Rest','Exercise','Continuous Acquisition'};
[indx,tf] = listdlg('ListString',list);
str = list{indx};
phase.WorkLoad = str;

%% Instantaneous Phase Extraction (cardiac ad respiratory)
phase.GetSignal; % gets mean intensity level of the images
phase.GetEMD; % decompose the siganl using EMD
phase.GetModes; % turn the modes into phase information

cardiac_phase = CardiacPhase(phase); % creates a CradiacPhase object
cardiac_phase.SelectPhase; % choose the right imf besed on the work load IMPROVE

respiratory_phase = RespiratoryPhase(phase); % creates a RespiratoryPhase object
respiratory_phase.SelectPhase; % choose the right imf based on the work load IMPROVE

cardiac_phase.Display
hold on
respiratory_phase.Display
% hold on
% hline = refline([0 cardiac_phase.SystoleTreshold]); % Reference line for systole treshold phase
% hline.Color = 'b';
% hline.LineStyle = ':';
% 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 = ':';
% hold on
% hline = refline([0 respiratory_phase.ExpirationTreshold]); % Reference line for diastole treshold phase
% hline.Color = 'r';
% hline.LineStyle = ':';
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
xlabel('Time (sec)');
ylabel('Phase(rad)');

%% 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
phase.ShowFrames(SI,'End-systolic frames in inpiration'); 
phase.ShowFrames(SE,'End-systolic frames in expiration');
phase.ShowFrames(DI,'End-diastolic frames in inpiration');
phase.ShowFrames(DE,'End-diastolic frames in expiration');

%% 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;
phase.ShowFrames(new_index{selected_range},'Heart Beat Reconstruction');

% c = ImageSequence
% c.ToBeSetByUser=15
% c.ToBeSetByUser=true