diff --git a/Code/CardiacPhase.m b/Code/CardiacPhase.m index e7acd93..714a111 100644 --- a/Code/CardiacPhase.m +++ b/Code/CardiacPhase.m @@ -1,16 +1,17 @@ classdef CardiacPhase < Phase properties FullPhase = Phase; - SelectedPhase = []; SystolePositions = []; DiastolePositions = []; + SystoleTreshold = double(0); + DiastoleTreshold = double(0); end methods %constructor function obj=CardiacPhase(Phase) if nargin==0;return;end % gets properties of superclass - obj.FullPhase = Phase.Values; + obj.FullPhase = Phase.ModePhases; obj.WorkLoad = Phase.WorkLoad; obj.Video = Phase.Video; end @@ -33,27 +34,30 @@ classdef CardiacPhase < Phase end end end - obj.SelectedPhase = obj.FullPhase(imf_number,:); % Phase corresponding to the cardiac phase + % 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 end % Detects the frame's positions of diastolic events function Diastole(obj) - max_heart_phase = max(obj.SelectedPhase); % cardiac phase maxima - treshold_diast = 0.8*max_heart_phase; % diastole treshold - index1 = obj.SelectedPhase > treshold_diast; % frames with phase above the diastolic treshold + max_heart_phase = max(obj.Values); % cardiac phase maxima + obj.DiastoleTreshold = 0.8*max_heart_phase; % diastole treshold + index1 = obj.Values > obj.DiastoleTreshold; % frames with phase above the diastolic treshold indexNumeric1 = find(index1); - values1 = obj.SelectedPhase(indexNumeric1); + values1 = obj.Values(indexNumeric1); [ordered_values1 oredered_index1] = sort(values1); obj.DiastolePositions = indexNumeric1(oredered_index1); % sorts the frames in ascending cardiac phase end % Detects the frame's positions of systolic events function Systole(obj) - min_heart_phase = min(obj.SelectedPhase); % cardiac phase minima - treshold_syst = 0.8*min_heart_phase; % systole treshold - index2 = obj.SelectedPhase < treshold_syst; % frames with phase below the systolic treshold + min_heart_phase = min(obj.Values); % cardiac phase minima + obj.SystoleTreshold = 0.8*min_heart_phase; % systole treshold + index2 = obj.Values < obj.SystoleTreshold; % frames with phase below the systolic treshold indexNumeric2 = find(index2); - values2 = obj.SelectedPhase(indexNumeric2); + values2 = obj.Values(indexNumeric2); [ordered_values2 oredered_index2] = sort(values2); obj.SystolePositions = indexNumeric2(oredered_index2); % sorts the frames in ascending cardiac phase end diff --git a/Code/Phase.m b/Code/Phase.m index 78ff25d..13efff1 100644 --- a/Code/Phase.m +++ b/Code/Phase.m @@ -3,9 +3,10 @@ classdef Phase < handle properties MeanIntensity = []; Video = []; - IMFs = []; - Values = []; WorkLoad = ''; + IMFs = []; + ModePhases = []; + Values = []; end properties % (GetAccess = private) should become private @@ -76,15 +77,19 @@ classdef Phase < handle % Apply Hilber Transform to all the IMFs and extract the % instantneous phase from each one of them - function GetPhase(obj) + function GetModes(obj) [number_imf lenght] = size(obj.IMFs); analytic_EMD = zeros(number_imf,lenght); for i=1:1:number_imf analytic_EMD(i,:) = hilbert(obj.IMFs(i,:)); %Apply Hilbert Transform - obj.Values(i,:) = angle(analytic_EMD(i,:))/pi; % Get instataneous phase + obj.ModePhases(i,:) = angle(analytic_EMD(i,:))/pi; % Get instataneous phase end end + function GetPhase(obj,imf) + obj.Values = obj.ModePhases(imf,:); + end + % Shows the frames of the original image sequence specified in the % array 'positions' function ShowFrames(obj,positions,title) @@ -93,6 +98,10 @@ classdef Phase < handle video = implay(triggered,3); set(video.Parent, 'Name', title); end + + function Display(obj) + plot(obj.Values,'-s','LineWidth',0.8,'MarkerSize',4); + end end end \ No newline at end of file diff --git a/Code/RespiratoryPhase.m b/Code/RespiratoryPhase.m index dfe6030..460dfa8 100644 --- a/Code/RespiratoryPhase.m +++ b/Code/RespiratoryPhase.m @@ -1,16 +1,23 @@ classdef RespiratoryPhase < Phase properties FullPhase = Phase; - SelectedPhase = []; InspirationPositions = []; ExpirationPositions = []; end + properties % (GetAccess = private) should become private + IMF = uint32(0); + InspirationTreshold = double(0); + ExpirationTreshold = double(0); + MaxRespPhase = double(0); + MinRespPhase = double(0); + end + methods %constructor function obj=RespiratoryPhase(Phase) if nargin==0;return;end % gets properties of superclass - obj.FullPhase = Phase.Values; + obj.FullPhase = Phase.ModePhases; obj.WorkLoad = Phase.WorkLoad; obj.Video = Phase.Video; end @@ -26,30 +33,33 @@ classdef RespiratoryPhase < Phase end end end - obj.SelectedPhase = obj.FullPhase(imf_number,:); %Phase corresponding to respiratory phase + % 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 end % Detects the frame's positions of inspiration moments function Inspiration(obj) - max_resp_phase = max(obj.SelectedPhase); % respiratory phase maxima - treshold_insp = 0.8*max_resp_phase; % inspiration treshold - index1 = obj.SelectedPhase > treshold_insp; % frames above inpiration treshold + obj.MaxRespPhase = max(obj.Values); % respiratory phase maxima + obj.InspirationTreshold = 0.8*obj.MaxRespPhase; % inspiration treshold + index1 = obj.Values > obj.InspirationTreshold; % frames above inpiration treshold indexNumeric1 = find(index1); - values1 = obj.SelectedPhase(indexNumeric1); + values1 = obj.Values(indexNumeric1); [ordered_values1 oredered_index1] = sort(values1); obj.InspirationPositions = indexNumeric1(oredered_index1); % sorts the frames in ascending respiratory phase end % Detects the frame's positions of expiration moments function Expiration(obj) - min_resp_phase = min(obj.SelectedPhase); % respiratory phase minima - treshold_exp = 0.8*min_resp_phase; % systole treshold - index2 = obj.SelectedPhase < treshold_exp; % frames below the systole treshold + obj.MinRespPhase = min(obj.Values); % respiratory phase minima + obj.ExpirationTreshold = 0.8*obj.MinRespPhase; % systole treshold + index2 = obj.Values < obj.ExpirationTreshold; % frames below the systole treshold indexNumeric2 = find(index2); - values2 = obj.SelectedPhase(indexNumeric2); + values2 = obj.Values(indexNumeric2); [ordered_values2 oredered_index2] = sort(values2); obj.ExpirationPositions = indexNumeric2(oredered_index2); % sorts the frames in ascending respiratory phase end - + end end \ No newline at end of file diff --git a/Code/test_OOP.m b/Code/test_OOP.m index 85923b9..bddb828 100644 --- a/Code/test_OOP.m +++ b/Code/test_OOP.m @@ -1,16 +1,23 @@ % Test Script % TO DO LIST: -% - Display phases % - Heart beat reconstruction with steady respiratory phase % - Multislice analysis +% - Choose IMF based on Instantaneous Frequency % - Error messages and robustness of the methods clear classes +%% Extracts data phase = Phase; % creates a Phase object phase.GetData; % extracts data +%% Time Information +time = 0:1:phase.nFrames-1; +repetitionTime = (1/15); % 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'}; @@ -18,17 +25,47 @@ list = {'Rest','Exercise','Continuous Acquisition'}; str = list{indx}; phase.WorkLoad = str; -phase.GetSignal % gets mean intensity level of the images -phase.GetEMD % decompose the siganl using EMD -phase.GetPhase % turn the modes into phase information +%% Instantaneous Phase Extraction (cardiac ad respiratory) +phase.GetSignal; % gets mean intensity level of the images +phase.GetEMD; % decompose the siganl using EMD +phase.GetModes; % turn the modes into phase information -cardiac_phase = CardiacPhase(phase) % creates a CradiacPhase object -cardiac_phase.SelectPhase % choose the right imf besed on the work load +cardiac_phase = CardiacPhase(phase); % creates a CradiacPhase object +cardiac_phase.SelectPhase; % choose the right imf besed on the work load IMPROVE + +respiratory_phase = RespiratoryPhase(phase); % creates a RespiratoryPhase object +respiratory_phase.SelectPhase; % choose the right imf based on the work load IMPROVE + +cardiac_phase.Display +hold on +respiratory_phase.Display +% 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)'); + +%% Cardiac and Respiratory events cardiac_phase.Diastole % detects diastolic events cardiac_phase.Systole % detects systolic events - -respiratory_phase = RespiratoryPhase(phase) % creates a RespiratoryPhase object -respiratory_phase.SelectPhase % choose the right imf besed on the work load respiratory_phase.Expiration % detects expirtion moments respiratory_phase.Inspiration % detects inspiration moments @@ -43,6 +80,27 @@ phase.ShowFrames(SE,'End-systolic frames in expiration'); phase.ShowFrames(DI,'End-diastolic frames in inpiration'); phase.ShowFrames(DE,'End-diastolic frames in expiration'); +%% 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; +phase.ShowFrames(new_index{selected_range},'Heart Beat Reconstruction'); + % c = ImageSequence % c.ToBeSetByUser=15 % c.ToBeSetByUser=true