From cb821f88486325da261ca60021c61a768c037a56 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jo=C3=A3o=20Narciso?= Date: Mon, 13 May 2019 10:15:01 +0200 Subject: [PATCH] Update 13/5/2019 --- Code/CardiacPhase.m | 39 +++++++++-------- Code/Complot.m | 41 ++++++++++++++++++ Code/Phase.m | 75 +++++++++++++++++++++----------- Code/RespiratoryPhase.m | 22 +++++----- Code/rParabEmd__L.m | 3 ++ Code/test_OOP.m | 94 +++++++++++++++++++++++++++++++---------- 6 files changed, 197 insertions(+), 77 deletions(-) create mode 100644 Code/Complot.m diff --git a/Code/CardiacPhase.m b/Code/CardiacPhase.m index 714a111..bbb5e8b 100644 --- a/Code/CardiacPhase.m +++ b/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') - imf_number = 3; - else if strcmp(obj.WorkLoad,'Exercise') - imf_number = 2; - else if strcmp(obj.WorkLoad,'Continuous Acquisition') - imf_number = 3; - end - end + prompt = {'Choose the IMF representing cardiac phase:','Combine with additional IMF:'}; + dlgtitle = 'Cardiac Phase'; + dims = [1 60]; + definput = {'3',''}; + opts.Interpreter = 'tex'; + answer = inputdlg(prompt,dlgtitle,dims,definput,opts); + obj.IMF = str2num(answer{1}); + + 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 diff --git a/Code/Complot.m b/Code/Complot.m new file mode 100644 index 0000000..942f290 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/Code/Phase.m b/Code/Phase.m index 13efff1..99793f3 100644 --- a/Code/Phase.m +++ b/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 - 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 + + prompt = {'Number of slices:'}; + dlgtitle = 'Slice'; + dims = [1 60]; + definput = {'1'}; + opts.Interpreter = 'tex'; + answer = inputdlg(prompt,dlgtitle,dims,definput,opts); + obj.NumberOfSlices = str2num(answer{1}); + + 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 diff --git a/Code/RespiratoryPhase.m b/Code/RespiratoryPhase.m index 460dfa8..c52f4be 100644 --- a/Code/RespiratoryPhase.m +++ b/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') - imf_number = 5; - else if strcmp(obj.WorkLoad,'Exercise') - imf_number = 4; - else if strcmp(obj.WorkLoad,'Continuous Acquisition') - imf_number = 5; - end - end - end + prompt = {'Choose the IMF representing respiratory phase:'}; + dlgtitle = 'Respiratory Phase'; + dims = [1 60]; + definput = {'5'}; + opts.Interpreter = 'tex'; + answer = inputdlg(prompt,dlgtitle,dims,definput,opts); + obj.IMF = str2num(answer{1}) % 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 diff --git a/Code/rParabEmd__L.m b/Code/rParabEmd__L.m index 1fe6041..2b9e286 100644 --- a/Code/rParabEmd__L.m +++ b/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 diff --git a/Code/test_OOP.m b/Code/test_OOP.m index bddb828..d023605 100644 --- a/Code/test_OOP.m +++ b/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 -% 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) +%% 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('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