Update 7/5/2019

This commit is contained in:
João Narciso 2019-05-07 16:36:02 +02:00
parent 431ff5f7d4
commit cd47edb1fd
4 changed files with 117 additions and 36 deletions

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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