Update 14/05/2'19

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

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

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

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);