218 lines
7.4 KiB
Matlab
218 lines
7.4 KiB
Matlab
% 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);
|