From a6ad9b0c7499bb381ddabee21d03486d680110b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jo=C3=A3o=20Narciso?= Date: Tue, 14 May 2019 14:54:03 +0200 Subject: [PATCH] Update 14/05/2'19 --- Code/CardiacPhase.m | 15 ++++- Code/Phase.m | 12 ++-- Code/RespiratoryPhase.m | 13 ++++ Code/test_OOP.m | 139 +++++++++++++++++++++++++++++----------- 4 files changed, 136 insertions(+), 43 deletions(-) diff --git a/Code/CardiacPhase.m b/Code/CardiacPhase.m index bbb5e8b..2bf7c01 100644 --- a/Code/CardiacPhase.m +++ b/Code/CardiacPhase.m @@ -4,7 +4,6 @@ classdef CardiacPhase < Phase IMF = uint32(0); SystolePositions = []; DiastolePositions = []; - end properties % (GetAccess = private) should become private @@ -18,6 +17,7 @@ classdef CardiacPhase < Phase % gets properties of superclass obj.FullPhase = Phase.ModePhases; obj.Video = Phase.Video; + obj.InstantaneousFreq = Phase.InstantaneousFreq; end % Choose the correct IMF depending on the work load @@ -43,6 +43,19 @@ classdef CardiacPhase < Phase obj.Values = obj.FullPhase(obj.IMF,:); %Phase corresponding to respiratory phase end + function GetFrequency(obj,time) + heart_inst_freq = (abs(diff(obj.Values))./diff(time))/(2*pi()); + obj.InstantaneousFreq = medfilt1(heart_inst_freq,3); + figure + scatter(time,[obj.InstantaneousFreq,0],8); + mean_freq1 = mean(obj.InstantaneousFreq); + hline = refline([0 mean_freq1]); % Reference line for maximum phase + hline.Color = 'm'; + title('Cardiac Instantaneous Frequency'); + xlabel('Time (sec)') + ylabel('Frequency (Hz)'); + end + % Detects the frame's positions of diastolic events function Diastole(obj) max_heart_phase = max(obj.Values); % cardiac phase maxima diff --git a/Code/Phase.m b/Code/Phase.m index 99793f3..ed4b700 100644 --- a/Code/Phase.m +++ b/Code/Phase.m @@ -6,6 +6,7 @@ classdef Phase < handle IMFs = []; ModePhases = []; Values = []; + InstantaneousFreq = []; NumberOfSlices = uint32(0); CurrentSlice = uint32(0); end @@ -92,6 +93,7 @@ classdef Phase < handle end function DisplayIntensity(obj) + figure plot(obj.MeanIntensity,'-s','LineWidth',0.8,'MarkerSize',4); end @@ -119,15 +121,17 @@ classdef Phase < handle % Shows the frames of the original image sequence specified in the % array 'positions' - function ShowFrames(obj,positions,title) + function ShowFrames(obj,positions,title,fps) triggered = []; - triggered = cat(2,triggered,obj.Video(:,:,positions)); - video = implay(triggered,3); - set(video.Parent, 'Name', title); + triggered = cat(2,triggered,obj.Video(:,:,positions)); + video = implay(triggered,fps); + set(video.Parent, 'Name', title); end function DisplayValues(obj) plot(obj.Values,'-s','LineWidth',0.8,'MarkerSize',4); + xlabel('Time (sec)'); + ylabel('Phase(rad)'); end end diff --git a/Code/RespiratoryPhase.m b/Code/RespiratoryPhase.m index c52f4be..e9142cd 100644 --- a/Code/RespiratoryPhase.m +++ b/Code/RespiratoryPhase.m @@ -37,6 +37,19 @@ classdef RespiratoryPhase < Phase obj.Values = obj.FullPhase(obj.IMF,:); %Phase corresponding to respiratory phase end + function GetFrequency(obj,time) + resp_inst_freq = (abs(diff(obj.Values))./diff(time))/(2*pi()); + obj.InstantaneousFreq = medfilt1(resp_inst_freq,3); + figure + scatter(time,[obj.InstantaneousFreq,0],8); + mean_freq1 = mean(obj.InstantaneousFreq); + hline = refline([0 mean_freq1]); % Reference line for maximum phase + hline.Color = 'm'; + title('Respiratory Instantaneous Frequency'); + xlabel('Time (sec)') + ylabel('Frequency (Hz)'); + end + % Detects the frame's positions of inspiration moments function Inspiration(obj) obj.MaxRespPhase = max(obj.Values); % respiratory phase maxima diff --git a/Code/test_OOP.m b/Code/test_OOP.m index d023605..d5a17a5 100644 --- a/Code/test_OOP.m +++ b/Code/test_OOP.m @@ -10,7 +10,7 @@ phase = Phase; % creates a Phase object phase.GetData; % extracts data %% Visualize data -phase.ShowFrames(:,'cine Images'); +phase.ShowFrames(:,'cine Images',13); %% Time Information time = 0:1:phase.nFrames-1; @@ -19,28 +19,28 @@ 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.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'); +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('IMF Phases'); +title('IMFs Phase'); te=repetitionTime; xt=get(gca,'xtick'); set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time values @@ -56,34 +56,20 @@ 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 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)'); %% 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 @@ -96,7 +82,8 @@ DI = intersect(cardiac_phase.DiastolePositions,respiratory_phase.InspirationPosi 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'); +% Systole in Inpiration +phase.ShowFrames(SI,'End-systolic frames in inpiration',3); directory = 'Systole in Inspiration' mkdir(directory); for i=1:length(SI) @@ -105,7 +92,26 @@ for i=1:length(SI) dicomwrite(phase.Video(:,:,1,SI(i)), fullFileName); end -phase.ShowFrames(SE,'End-systolic frames in expiration'); +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) @@ -114,7 +120,26 @@ for i=1:length(SE) dicomwrite(phase.Video(:,:,1,SE(i)), fullFileName); end -phase.ShowFrames(DI,'End-diastolic frames in inpiration'); +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) @@ -123,7 +148,26 @@ for i=1:length(DI) dicomwrite(phase.Video(:,:,1,DI(i)), fullFileName); end -phase.ShowFrames(DE,'End-diastolic frames in expiration'); +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) @@ -132,6 +176,25 @@ for i=1:length(DE) 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 @@ -151,4 +214,4 @@ for i=1:1:n_ranges new_index{i} = k(j); % ordered indexes end selected_range = n_ranges-2; % optional range -phase.ShowFrames(new_index{selected_range},'Heart Beat Reconstruction'); +phase.ShowFrames(new_index{selected_range},'Heart Beat Reconstruction',3);