Update 13/5/2019

Code/CardiacPhase.m

@@ -1,8 +1,13 @@
 classdef CardiacPhase < Phase
         properties
         FullPhase = Phase;
+        IMF = uint32(0);
         SystolePositions = [];
         DiastolePositions = [];
+        
+        end
+        
+        properties % (GetAccess = private) should become private
         SystoleTreshold = double(0);
         DiastoleTreshold = double(0);
         end
@@ -12,32 +17,30 @@ classdef CardiacPhase < Phase
             if nargin==0;return;end
             % gets properties of superclass
             obj.FullPhase = Phase.ModePhases;
-            obj.WorkLoad = Phase.WorkLoad;
             obj.Video = Phase.Video;
         end
-%     
+
-%     methods %set methods
-% %         function set.ToBeSetByUser(obj,val)
-% %             if ~isboolean(val);error('provide a boolean value please');end
-% %             obj.ToBeSetByUser=val;
-% %         end
-%     end
-%
         % Choose the correct IMF depending on the work load
         function SelectPhase(obj)
-            if strcmp(obj.WorkLoad,'Rest')
+            prompt = {'Choose the IMF representing cardiac phase:','Combine with additional IMF:'};
-                imf_number = 3;
+            dlgtitle = 'Cardiac Phase';
-            else if strcmp(obj.WorkLoad,'Exercise')
+            dims = [1 60];
-                    imf_number = 2;
+            definput = {'3',''};
-                else if strcmp(obj.WorkLoad,'Continuous Acquisition')
+            opts.Interpreter = 'tex';
-                        imf_number = 3;
+            answer = inputdlg(prompt,dlgtitle,dims,definput,opts);
-                    end
+            obj.IMF = str2num(answer{1});
-                end
+            
+            if  length(answer{2}) > 0.
+                 new_imf = obj.FullPhase(obj.IMF,:) + obj.FullPhase(str2num(answer{2}),:); % if the IMFs are unorthogonal combine them into one IMF
+                 obj.FullPhase(obj.IMF,:) = new_imf;
+                 obj.FullPhase(str2num(answer{2}),:) = [];
+                 
+                 Complot(obj.FullPhase);
             end
             % SHOULD TRY TO USE GETPHASE FUNCTION
 %             obj.IMF = imf_number;
 %             GetPhase@Phase(obj.FullPhase,obj.IMF);
-            obj.Values = obj.FullPhase(imf_number,:); % Phase corresponding to the cardiac phase
+            obj.Values = obj.FullPhase(obj.IMF,:); %Phase corresponding to respiratory phase
         end
         
         % Detects the frame's positions of diastolic events

Code/Complot.m

@@ -0,0 +1,41 @@
+function [] = Complot(varargin)
+
+error(nargchk(1,2,nargin));
+[msg,x,S] = Cplotchk(varargin);
+error(msg)
+
+% Get some spacing between the various signals
+alpha = mean(max(S')-min(S'));
+% save alpha alpha
+% load alpha alpha
+for i = 1:size(S,1)
+  SS(i,:) = S(i,:) - alpha*(i-1)*ones(size(S(1,:)));
+end
+
+% Plotting per se
+if nargin == 1
+   plot(SS')
+   else
+   plot(x,SS')
+end
+axis([x(1) x(length(SS)) min(min(SS)) max(max(SS)) ])
+% axis([0 x(length(SS)) min(min(SS)) max(max(SS)) ])
+set(gca,'YTick',-size(S,1)*alpha:alpha:0);
+set(gca,'YTicklabel',size(S,1)+1:-1:1);
+
+% Reading routine
+function [msg,x,S] = Cplotchk(P);
+msg = [];
+
+if length(P) > 2
+   msg = 'Too many inputs.';
+else
+   if length(P) == 1
+      S = P{1};
+      x = 1:length(S);
+   end
+   if length(P) == 2
+      x = P{1};
+      S = P{2};
+   end
+end

Code/Phase.m

@@ -3,10 +3,11 @@ classdef Phase < handle
     properties 
         MeanIntensity = [];
         Video = [];
-        WorkLoad = '';
         IMFs = [];
         ModePhases = [];
         Values = [];
+        NumberOfSlices = uint32(0);
+        CurrentSlice = uint32(0);
     end
     
     properties % (GetAccess = private) should become private
@@ -25,37 +26,54 @@ classdef Phase < handle
         end
     end
     
-    methods %set methods
-        function set.WorkLoad(obj,str)
-            obj.WorkLoad = str;
-%             if indx == 1
-%                 obj.WorkLoad = 'Rest';
-%             else if indx == 2
-%                     obj.WorkLoad = 'Exercise';
-%                 else
-%                     obj.WorkLoad = 'Continuous Acquisition';
-%             end
-        end
-    end
-    
     methods
         % Extracts the image sequence from the folder
         function GetData(obj)
+            
             directory = uigetdir(obj.InitialDirectory); % Choose the cineMRI folder 
             dicomFiles = dir(fullfile(directory, '*.dcm' ));
             obj.nFrames = length(dicomFiles); % Number of frames of the cineMRI
             info = dicominfo(fullfile(directory,dicomFiles(1).name));
             obj.xImageSize = info.Width; % Image width
             obj.yImageSize = info.Height; % Image height
-            obj.Video = zeros(obj.xImageSize, obj.yImageSize, 1, obj.nFrames, 'uint8');
+            
-            for p=1:obj.nFrames
+            prompt = {'Number of slices:'};
-               filename = fullfile(directory, dicomFiles(p).name );
+            dlgtitle = 'Slice';
-               Info=dicominfo(filename);
+            dims = [1 60];
-               if Info.Width~=obj.xImageSize || Info.Height~=obj.yImageSize
+            definput = {'1'};
-                   error('invalid image dimensions')
+            opts.Interpreter = 'tex';
-               end
+            answer = inputdlg(prompt,dlgtitle,dims,definput,opts);
-               obj.Video(:,:,1,p) = dicomread(filename); % Read the DICOM files in the folder
+            obj.NumberOfSlices = str2num(answer{1});
-            end 
+            
+            if obj.NumberOfSlices == 1
+                obj.Video = zeros(obj.xImageSize, obj.yImageSize, 1, obj.nFrames, 'uint8');
+                for p=1:obj.nFrames
+                   filename = fullfile(directory, dicomFiles(p).name );
+                   Info=dicominfo(filename);
+                   if Info.Width~=obj.xImageSize || Info.Height~=obj.yImageSize
+                       error('invalid image dimensions')
+                   end
+                   obj.Video(:,:,1,p) = dicomread(filename); % Read the DICOM files in the folder
+                end
+            else
+                list_array = 2:1:obj.NumberOfSlices;
+                list = {'2','3','4','5','6'};
+                [indx,tf] = listdlg('ListString',list);
+                obj.CurrentSlice = indx;
+                
+                obj.nFrames = obj.nFrames/obj.NumberOfSlices;
+                
+                obj.Video = zeros(obj.xImageSize, obj.yImageSize, 1, obj.nFrames, 'uint8');
+                for p=1:obj.nFrames
+                   shift = (obj.CurrentSlice-1)*obj.nFrames;
+                   filename = fullfile(directory, dicomFiles(shift+p).name );
+                   Info=dicominfo(filename);
+                   if Info.Width~=obj.xImageSize || Info.Height~=obj.yImageSize
+                       error('invalid image dimensions')
+                   end
+                   obj.Video(:,:,1,p) = dicomread(filename); % Read the DICOM files in the folder
+                end
+            end
         end
 
         % Gets the mean intensity value of the image sequence usin the
@@ -63,12 +81,21 @@ classdef Phase < handle
         function GetSignal(obj)
              d2IFT = ifft2(obj.Video); % performs Inverse Fourier Transform
              d2IFT_shifted = fftshift(d2IFT);
+             
+%              video = implay(abs(d2IFT_shifted));
+%              set(video.Parent, 'Name', '2D IFT of the Images Sequence along time');
+             
              for i=1:obj.nFrames
                 obj.MeanIntensity(:,i) = d2IFT_shifted(floor(obj.yImageSize/2),floor(obj.xImageSize/2),i); % extracts intensity of the central point in k-space for every frame
              end
              obj.MeanIntensity = abs(obj.MeanIntensity);
         end
         
+        function DisplayIntensity(obj)
+            plot(obj.MeanIntensity,'-s','LineWidth',0.8,'MarkerSize',4);
+        end
+        
+        
         % Perfomrs Empirical Mode Decomposition
         function GetEMD(obj)
             obj.IMFs = rParabEmd__L(obj.MeanIntensity,40,40,1);
@@ -99,7 +126,7 @@ classdef Phase < handle
                 set(video.Parent, 'Name', title);
         end
         
-        function Display(obj)
+        function DisplayValues(obj)
             plot(obj.Values,'-s','LineWidth',0.8,'MarkerSize',4);
         end
     

Code/RespiratoryPhase.m

@@ -1,12 +1,12 @@
 classdef RespiratoryPhase < Phase
     properties
         FullPhase = Phase;
+        IMF = uint32(0);
         InspirationPositions = [];
         ExpirationPositions = [];
     end
     
     properties % (GetAccess = private) should become private
-        IMF = uint32(0);
         InspirationTreshold = double(0);
         ExpirationTreshold = double(0);
         MaxRespPhase = double(0);
@@ -18,25 +18,23 @@ classdef RespiratoryPhase < Phase
             if nargin==0;return;end
             % gets properties of superclass
             obj.FullPhase = Phase.ModePhases;
-            obj.WorkLoad = Phase.WorkLoad;
+%             obj.WorkLoad = Phase.WorkLoad;
             obj.Video = Phase.Video;
        end
    
     % Choose the correct IMF depending on the work load
     function SelectPhase(obj)
-            if strcmp(obj.WorkLoad,'Rest')
+            prompt = {'Choose the IMF representing respiratory phase:'};
-                imf_number = 5;
+            dlgtitle = 'Respiratory Phase';
-            else if strcmp(obj.WorkLoad,'Exercise')
+            dims = [1 60];
-                    imf_number = 4;
+            definput = {'5'};
-                else if strcmp(obj.WorkLoad,'Continuous Acquisition')
+            opts.Interpreter = 'tex';
-                        imf_number = 5;
+            answer = inputdlg(prompt,dlgtitle,dims,definput,opts);
-                    end
+            obj.IMF = str2num(answer{1})
-                end
-            end
             % SHOULD TRY TO USE GETPHASE FUNCTION
 %             obj.IMF = imf_number;
 %             GetPhase@Phase(obj.FullPhase,obj.IMF);
-            obj.Values = obj.FullPhase(imf_number,:); %Phase corresponding to respiratory phase
+            obj.Values = obj.FullPhase(obj.IMF,:); %Phase corresponding to respiratory phase
     end
         
         % Detects the frame's positions of inspiration moments

Code/rParabEmd__L.m

@@ -1,3 +1,6 @@
+% This function was extracted from https://nl.mathworks.com/matlabcentral/fileexchange/21409-empirical-mode-decomposition
+%  Authors: Raul Rato (rtr@uninova.DOT.pt) and Manuel Ortigueira (mdortigueira@uninova.pt or mdo@fct.unl.pt)
+
 function rParabEmd = rParabEmd__L (x, qResol, qResid, qAlfa)
 dbstop if warning
 if(nargin~=4), error('rParabEmd__L: Use with 4 inputs.'), end

Code/test_OOP.m

@@ -1,9 +1,6 @@
 % 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
@@ -12,33 +9,53 @@ clear classes
 phase = Phase; % creates a Phase object
 phase.GetData; % extracts data
 
+%% Visualize data
+phase.ShowFrames(:,'cine Images');
+
 %% Time Information
 time = 0:1:phase.nFrames-1;
-repetitionTime = (1/15); % time between samples
+repetitionTime = (1/13); % time between samples
 time = repetitionTime*time; % sampling moments in time
 
-%% Workload information
+%% Instantaneous Phase Extraction (cardiac and respiratory)
-% 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.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('IMF Phases');
+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 right imf besed on the work load IMPROVE
+cardiac_phase.SelectPhase; % choose the cardiac imf
 
 respiratory_phase = RespiratoryPhase(phase); % creates a RespiratoryPhase object
-respiratory_phase.SelectPhase; % choose the right imf based on the work load IMPROVE
+respiratory_phase.SelectPhase; % choose the respiratory imf
 
-cardiac_phase.Display
+figure
+cardiac_phase.DisplayValues
 hold on
-respiratory_phase.Display
+respiratory_phase.DisplayValues
 % hold on
 % hline = refline([0 cardiac_phase.SystoleTreshold]); % Reference line for systole treshold phase
 % hline.Color = 'b';
@@ -63,6 +80,10 @@ set(gca,'xticklabel',arrayfun(@num2str,xt*te,'un',0)); % adjust x axis to time v
 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);
+
 %% Cardiac and Respiratory events
 cardiac_phase.Diastole % detects diastolic events
 cardiac_phase.Systole % detects systolic events
@@ -75,10 +96,41 @@ 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'); 
+phase.ShowFrames(SI,'End-systolic frames in inpiration');
+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
+
 phase.ShowFrames(SE,'End-systolic frames in expiration');
+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
+
 phase.ShowFrames(DI,'End-diastolic frames in inpiration');
+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
+
 phase.ShowFrames(DE,'End-diastolic frames in expiration');
+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
 
 %% Heart Beat Reconstruction
 n_ranges = 30; % number of phase ranges
@@ -98,9 +150,5 @@ for i=1:1:n_ranges
     j = oredered_index{i};
     new_index{i} = k(j); % ordered indexes
 end
-selected_range = n_ranges-2;
+selected_range = n_ranges-2; % optional range
 phase.ShowFrames(new_index{selected_range},'Heart Beat Reconstruction');
-
-% c = ImageSequence
-% c.ToBeSetByUser=15
-% c.ToBeSetByUser=true