p147685/CardiacPhase
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Initial class construction

Browse Source
This commit is contained in:
João Narciso
2019-05-06 16:34:28 +02:00
parent 67f2d57e03
commit 431ff5f7d4
5813 changed files with 1622108 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

62
Code/CardiacPhase.m Normal file
View File

@ -0,0 +1,62 @@
classdef CardiacPhase < Phase
properties
FullPhase = Phase;
SelectedPhase = [];
SystolePositions = [];
DiastolePositions = [];
end
methods %constructor
function obj=CardiacPhase(Phase)
if nargin==0;return;end
% gets properties of superclass
obj.FullPhase = Phase.Values;
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
end
obj.SelectedPhase = 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
indexNumeric1 = find(index1);
values1 = obj.SelectedPhase(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
indexNumeric2 = find(index2);
values2 = obj.SelectedPhase(indexNumeric2);
[ordered_values2 oredered_index2] = sort(values2);
obj.SystolePositions = indexNumeric2(oredered_index2); % sorts the frames in ascending cardiac phase
end
end
end

98
Code/Phase.m Normal file
View File

@ -0,0 +1,98 @@
classdef Phase < handle
properties
MeanIntensity = [];
Video = [];
IMFs = [];
Values = [];
WorkLoad = '';
end
properties % (GetAccess = private) should become private
nFrames = uint32(0);
xImageSize = uint16(0);
yImageSize = uint16(0);
end
properties(Constant)
InitialDirectory = 'C:\Uss\Jero<EFBFBD>o\Desktop\Tese\Data';
end
methods %constructor
function obj=Phase(varargin)
if nargin==0;return;end
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
end
% Gets the mean intensity value of the image sequence usin the
% central point in k-space for all the images
function GetSignal(obj)
d2IFT = ifft2(obj.Video); % performs Inverse Fourier Transform
d2IFT_shifted = fftshift(d2IFT);
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
% Perfomrs Empirical Mode Decomposition
function GetEMD(obj)
obj.IMFs = rParabEmd__L(obj.MeanIntensity,40,40,1);
obj.IMFs = obj.IMFs.';
end
% Apply Hilber Transform to all the IMFs and extract the
% instantneous phase from each one of them
function GetPhase(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
end
end
% Shows the frames of the original image sequence specified in the
% array 'positions'
function ShowFrames(obj,positions,title)
triggered = [];
triggered = cat(2,triggered,obj.Video(:,:,positions));
video = implay(triggered,3);
set(video.Parent, 'Name', title);
end
end
end

55
Code/RespiratoryPhase.m Normal file
View File

@ -0,0 +1,55 @@
classdef RespiratoryPhase < Phase
properties
FullPhase = Phase;
SelectedPhase = [];
InspirationPositions = [];
ExpirationPositions = [];
end
methods %constructor
function obj=RespiratoryPhase(Phase)
if nargin==0;return;end
% gets properties of superclass
obj.FullPhase = Phase.Values;
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
obj.SelectedPhase = 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
indexNumeric1 = find(index1);
values1 = obj.SelectedPhase(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
indexNumeric2 = find(index2);
values2 = obj.SelectedPhase(indexNumeric2);
[ordered_values2 oredered_index2] = sort(values2);
obj.ExpirationPositions = indexNumeric2(oredered_index2); % sorts the frames in ascending respiratory phase
end
end
end

306
Code/rParabEmd__L.m Normal file
View File

@ -0,0 +1,306 @@
function rParabEmd = rParabEmd__L (x, qResol, qResid, qAlfa)
dbstop if warning
if(nargin~=4), error('rParabEmd__L: Use with 4 inputs.'), end
if(nargout>1), error('rParabEmd__L: Use with just one output.'), end
ArgCheck_s(x, qResol, qResid, qAlfa)
% Actual computation -------------------------------------
kc = x(:); % ket copy of the input signal
Wx= kc'*kc; % Original signal energy
quntN = length(kc); % Signal length
% loop to decompose the input signal into successive IMFs
rParabEmd= []; % Matrix which will contain the successive IMFs, and the residue
rParabEmdCnt= 0;
qDbResid= 0; %Equal energies at start
quntOscCnt= quntNOsc_s(kc);
while ((qDbResid<qResid) && (quntOscCnt>2) ) % c has some energy and oscilates
kImf = kc; % at the beginning of the sifting process, kImf is the signal
rPMOri= rGetPMaxs_s(kImf); % rPM= [xM(M), yM(M)];
rPmOri= rGetPMins_s(kImf); % rPm= [xm(m), ym(m)];
rPM= rPMaxExtrapol_s(rPMOri, rPmOri, quntN);
rPm= rPMinExtrapol_s(rPMOri, rPmOri, quntN);
quntLM= length(rPM); quntLm= length(rPm);
% if (abs(quntLM-quntLm)>2), disp('Debug: Max-Min count mismatch.'),keyboard,end;
if (abs(quntLM-quntLm)>2), disp('Debug: Max-Min count mismatch.'),end;
if(sum(abs(diff(sign(rPM(1:min(quntLM,quntLm),1)- rPm(1:min(quntLM,quntLm),1)))))>0)
% disp('Debug: Max-Min sequence mismatch.'),keyboard;
disp('Debug: Max-Min sequence mismatch.');
end
if(sum(abs(diff(sign(rPm(1:min(quntLM,quntLm),1)- rPM(1:min(quntLM,quntLm),1)))))>0)
% disp('Debug: Max-Min reverse sequence mismatch.'),keyboard;
disp('Debug: Max-Min reverse sequence mismatch.');
end
bTenv= spline(rPM(:,1), rPM(:,2), 1:quntN); % Top envelop: bTenv[n];
bDenv= spline(rPm(:,1), rPm(:,2), 1:quntN); % Down envelop: bDenv[n];
bBias= (bTenv+bDenv)/2; % first bias estimate
while true(1) % inner loop to find each IMF
WImf= kImf'*kImf; %current IMF energy
WBias= bBias*bBias'; %bias energy
if WBias*WImf<0 , warning('rParabEmd__L: Ooops, negative energy detected.'), end
if WBias> 0, DbqResol= 10*log10(WImf/WBias); else DbqResol= Inf; end
if (DbqResol>qResol), break, end %Resolution reached
%Resolution not reached. More work is needed
kImf = kImf- qAlfa*bBias'; % subtract qAlfa bias from kImf
rPMOri= rGetPMaxs_s(kImf); % rPM= [xM(M), yM(M)];
rPmOri= rGetPMins_s(kImf); % rPm= [xm(m), ym(m)];
rPM= rPMaxExtrapol_s(rPMOri, rPmOri, quntN);
rPm= rPMinExtrapol_s(rPMOri, rPmOri, quntN);
bTenv= spline(rPM(:,1), rPM(:,2), 1:quntN); % Top envelop: bTenv[n];
bDenv= spline(rPm(:,1), rPm(:,2), 1:quntN); % Down envelop: bDenv[n];
bBias= (bTenv+bDenv)/2; % new bias estimate
end % Wend true
%
rParabEmd = [rParabEmd; kImf']; % store the extracted rParabEmd in the matrix rParabEmd
kc = kc - kImf; % subtract the extracted rParabEmd from the signal
quntOscCnt= quntNOsc_s(kc);
rParabEmdCnt=rParabEmdCnt+1;
if (kc'*kc)>0
qDbResid= 10*log10(Wx/(kc'*kc));
else
qDbResid = Inf
end
%
end % Wend ((DbR... ))
if ((kc'*kc)/Wx)>(10^-12)
rParabEmd=[rParabEmd; kc']; %The residual is the last IMF
rParabEmdCnt=rParabEmdCnt+1;
NumOscqResiduais= quntNOsc_s(kc);
end
rParabEmd= rParabEmd';
end %main function
%SubFunctions ------------------------------------------------------------
%-------------------------------------------------------------------------
function ArgCheck_s(x, qResol, qResid, qAlfa)
[qL, qC] = size(x);
if ((qL*qC)~= max(qL,qC)), error('rParabEmd__L: Input signal must be a one dim vector.'), end
if ((qL*qC)<= 1), error('rParabEmd__L: Input signal must be a vector.'), end
[qL,qC] = size(qResol);
if ( ~((qL==1)&(qC==1)) ), error('rParabEmd__L: Input resolution must be a scalar.'), end
if ( qResol<=0 ), error('rParabEmd__L: Input resolution must strictly positive.'), end
[qL,qC] = size(qResid);
if ( ~((qL==1)&(qC==1)) ), error('rParabEmd__L: Input residual must be a scalar.'), end
if ( qResid<=0 ), error('rParabEmd__L: Input residual must strictly positive.'), end
[qL,qC] = size(qAlfa);
if ( ~((qL==1)&(qC==1)) ), error('rParabEmd__L: qAlfa step must be a scalar.'), end
if ( qAlfa<=0 ), error('rParabEmd__L: qAlfa step must be strictly positive.'), end
end
%--------------------------------------------------------------------------
%---------- make at 17-Jul-07 10:16:59.44
% quntNOsc_s v1.01
% build 20070409001
% Returns the oscilation count, no steps
function quntNOsc = quntNOsc_s (x)
y=0; qisTop= false; qisDown= false;
for i=2:(length(x)-1)
if( ((x(i-1)) < (x(i))) && ((x(i+1))< (x(i))) ) %Max /-\
y=y+1;
end
if( ((x(i-1)) > (x(i))) && ((x(i+1))> (x(i))) ) %min \_/
y=y+1;
end
%Top
if( ((x(i-1)) < (x(i))) && ((x(i+1))== (x(i))) ) %StepL /-
qisTop= true; qisDown= false;
end
if( ((x(i-1)) == (x(i))) && ((x(i+1))< (x(i))) ) %stepR -\
if qisTop; y=y+1; end;
qisTop= false;
end
%Downs
if( ((x(i-1)) > (x(i))) && ((x(i+1))== (x(i))) ) %stepL \_
qisTop= false; qisDown= true;
end
if( ((x(i-1)) == (x(i))) && ((x(i+1))> (x(i))) ) %StepR _/
if qisDown; y=y+1; end
qisDown=false;
end
end % for i=2:(length(x)-1)
quntNOsc= y;
end % function y = quntNOsc_s (x)
%---------- make at 17-Jul-07 10:16:59.44
function rPMaxExtrapol= rPMaxExtrapol_s(rPM, rPm, quntL)
%rPMaxExtrapol_s V1.00
% build 2007407001
% Time-mirrored top extrema (Parabolic Maxs) extrapolation
%Init ------------------------------------
rPM= sortrows(rPM); %assumes nothing on rPM sort order
rPm= sortrows(rPm); %assumes nothing on rPm sort order
kTopTim1= rPM(:,1); kTopVal= rPM(:,2);
kDwnTim1= rPm(:,1); kDwnVal= rPm(:,2);
%Start extrapolation ---------------------
if ( (kTopTim1(1)== 1) && (kDwnTim1(1)== 1) )
disp ('rPMaxExtrapol_s: Poliextrema at signal''s start');
elseif ( (kTopTim1(1)<1) || (kDwnTim1(1)< 1) )
disp ('rPMaxExtrapol_s: Invalid extrema at signal''s start');
else
kTopTim1=[2-kDwnTim1(1); kTopTim1]; % New first Top at the (one based) specular Min
kTopVal=[kTopVal(1); kTopVal]; % Same Val as old first Top
end
% End extrapolation -----------------------
if ( (kTopTim1(end)== quntL) && (kDwnTim1(end)== quntL) )
disp ('rPMaxExtrapol_s: Poliextrema at signal''s end');
elseif ( (kTopTim1(end)> quntL) || (kDwnTim1(end)> quntL) )
disp ('rPMaxExtrapol_s: Invalid extrema at signal''s end');
else
kTopTim1=[kTopTim1; (2*quntL - kDwnTim1(end))]; % New last Top at the specular Min
kTopVal=[ kTopVal; kTopVal(end)]; % Same Val as old last Top
end
% return value ------------------------
rPMaxExtrapol= sortrows([kTopTim1, kTopVal]);
end
%-------------------------------------------------------------------------
%---------- make at 17-Jul-07 10:16:59.44
function rPMinExtrapol= rPMinExtrapol_s(rPM, rPm, quntL)
%rPMinExtrapol_s V1.00
% build 2007407001
% Time-mirrored down extrema (Parabolic Mins) extrapolation
%Init ------------------------------------
rPM= sortrows(rPM); %assumes nothing on rPM sort order
rPm= sortrows(rPm); %assumes nothing on rPm sort order
kTopTim1= rPM(:,1); kTopVal= rPM(:,2);
kDwnTim1= rPm(:,1); kDwnVal= rPm(:,2);
%Start extrapolation ---------------------
if ( (kTopTim1(1)== 1) && (kDwnTim1(1)== 1) )
disp ('rPMinExtrapol_s: Poliextrema at signal''s start');
elseif ( (kTopTim1(1)<1) || (kDwnTim1(1)< 1) )
disp ('rPMinExtrapol_s: Invalid extrema at signal''s start');
else
kDwnTim1=[2-kTopTim1(1); kDwnTim1]; % New first Dwn at the (one based) specular Max
kDwnVal=[kDwnVal(1); kDwnVal]; % Same Val as old first Dwn
end
% End extrapolation -----------------------
if ( (kTopTim1(end)== quntL) && (kDwnTim1(end)== quntL) )
disp ('rPMinExtrapol_s: Poliextrema at signal''s end');
elseif ( (kTopTim1(end)> quntL) || (kDwnTim1(end)> quntL) )
disp ('rPMinExtrapol_s: Invalid extrema at signal''s end');
else
kDwnTim1=[kDwnTim1; (2*quntL - kTopTim1(end))]; % New last Dwn at the specular Max
kDwnVal=[ kDwnVal; kDwnVal(end)]; % Same Val as old last Dwn
end
% return value ------------------------
rPMinExtrapol= sortrows([kDwnTim1, kDwnVal]);
end
%-------------------------------------------------------------------------
%---------- make at 17-Jul-07 10:16:59.44
function rPMax= rGetPMaxs_s(aS) %Get Parabolic Maxs, plateaus out
% build 20070612001
kS= aS(:);
quntLenS=length(kS);
quntMaxCnt=0;
kSMNdx1= []; kSMVal=[]; %signal S Maxima indices and values
kSPMTim1= []; kSPMVal=[]; %signal S Parabolic Maxima times and values
if (quntLenS>2) %if signal has enough length
for Cnt=2:(quntLenS-1) %search the Maxs
if ( ((kS(Cnt) > kS(Cnt+1))) && ((kS(Cnt) >= kS(Cnt-1))) || ((kS(Cnt) >= kS(Cnt+1))) && ((kS(Cnt) > kS(Cnt-1))) )
quntMaxCnt=quntMaxCnt+1;
kSMNdx1= [kSMNdx1; Cnt]; kSMVal=[kSMVal; kS(Cnt)];
end
end
end
% Now we have the Maxs, lets get the Parabolic Maxs
oldxv= -Inf; oldyv= -Inf;
intGapMax= max(kS)-min(kS);
for jj=1:quntMaxCnt %for all Maxs
%xa= -1; xb= 0; xc= 1;
ya= kS(kSMNdx1(jj)-1); % Sample point before
yb= kS(kSMNdx1(jj)); % Sample point, == kSMVal(jj)
yc= kS(kSMNdx1(jj)+1); % Sample point after
D= (-4*yb+2*ya+2*yc);
if (D==0), xv= kSMNdx1(jj);
else xv= kSMNdx1(jj)+(ya-yc)/D; end; % Vertix abscissa
D= (-16*yb+ 8*ya+ 8*yc);
if (D==0), yv= yb;
else yv= yb+ (2*yc*ya- ya*ya- yc*yc)/D; end;
% Lets check for double maxima
if ( (xv==oldxv)||(abs(yv-oldyv)/abs(xv-oldxv))> (2*intGapMax) )
xv= (xv+ oldxv)/2; yv= max(yv,oldyv); %Double found
kSPMTim1(length(kSPMTim1))= xv; kSPMVal(length(kSPMVal))= yv;
else
kSPMTim1= [kSPMTim1; xv]; kSPMVal=[kSPMVal; yv];
end
oldxv= xv; oldyv= yv;
end % for jj=1:quntMaxCnt
if quntMaxCnt>0
if ( kS(1) >= kSPMVal(1) )
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal ]; %Start must be included as a Max
end
if ( kS(end) >= kSPMVal(end))
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Max
end
end
if quntMaxCnt==0
if ( kS(1) > kS(2) )
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal ]; %Start must be included as a Max
end
if ( kS(end) > kS(end-1))
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Max
end
end
if quntMaxCnt<0
error('rGetPMaxs_s: Invalid MaxCnt value');
end
rPMax= sortrows([kSPMTim1, kSPMVal]);
end
%---------- make at 17-Jul-07 10:16:59.44
function rPMin= rGetPMins_s(aS) %Get Parabolic Mins, plateaus out
% build 20070612001
kS= aS(:);
quntLenS=length(kS);
quntMinCnt=0;
kSMNdx1= []; kSMVal=[]; %signal S Minima indices and values
kSPMTim1= []; kSPMVal=[]; %signal S Parabolic Minima times and values
if (quntLenS>2) %if signal has enough length
for Cnt=2:(quntLenS-1) %search the Mins
if ( ((kS(Cnt) < kS(Cnt+1))) && ((kS(Cnt) <= kS(Cnt-1))) || ((kS(Cnt) <= kS(Cnt+1))) && ((kS(Cnt) < kS(Cnt-1))) )
quntMinCnt=quntMinCnt+1;
kSMNdx1= [kSMNdx1; Cnt]; kSMVal=[kSMVal; kS(Cnt)];
end
end
end
% Now we have the Mins, lets get the Parabolic Mins
oldxv= -Inf; oldyv= -Inf;
intGapMax= max(kS)-min(kS);
for jj=1:quntMinCnt %for all Mins
%xa= -1; xb= 0; xc= 1;
ya= kS(kSMNdx1(jj)-1); % Sample point before
yb= kS(kSMNdx1(jj)); % Sample point, == kSMVal(jj)
yc= kS(kSMNdx1(jj)+1); % Sample point after
D= (-4*yb+2*ya+2*yc);
if (D==0), xv= kSMNdx1(jj);
else xv= kSMNdx1(jj)+(ya-yc)/D; end; % Vertix abscissa
D= (-16*yb+ 8*ya+ 8*yc);
if (D==0), yv= yb;
else yv= yb+ (2*yc*ya- ya*ya- yc*yc)/D; end;
% Lets check for double minima
if ( (xv==oldxv)||(abs(yv-oldyv)/abs(xv-oldxv))> (2*intGapMax) )
xv= (xv+ oldxv)/2; yv= min(yv,oldyv); %Double found
kSPMTim1(length(kSPMTim1))= xv; kSPMVal(length(kSPMVal))= yv;
else
kSPMTim1= [kSPMTim1; xv]; kSPMVal=[kSPMVal; yv];
end
oldxv= xv; oldyv= yv;
end % for jj=1:quntMinCnt
if quntMinCnt>0
if ( kS(1) <= kSPMVal(1) )
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal ]; %Start must be included as a Min
end
if ( kS(end) <= kSPMVal(end))
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Min
end
end
if quntMinCnt==0
if ( kS(1) < kS(2) )
kSPMTim1= [1; kSPMTim1]; kSPMVal=[kS(1); kSPMVal]; %Start must be included as a Min
end
if ( kS(end) < kS(end-1))
kSPMTim1= [kSPMTim1; quntLenS]; kSPMVal=[kSPMVal; kS(end)]; %End must be included as a Min
end
end
if quntMinCnt<0
error('rGetPMins_s: Invalid MinCnt value');
end
rPMin= sortrows([kSPMTim1, kSPMVal]);
end
%---------- make at 17-Jul-07 10:16:59.44

48
Code/test_OOP.m Normal file
View File

@ -0,0 +1,48 @@
% Test Script
% TO DO LIST:
% - Display phases
% - Heart beat reconstruction with steady respiratory phase
% - Multislice analysis
% - Error messages and robustness of the methods
clear classes
phase = Phase; % creates a Phase object
phase.GetData; % extracts data
% 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;
phase.GetSignal % gets mean intensity level of the images
phase.GetEMD % decompose the siganl using EMD
phase.GetPhase % 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.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
SI = intersect(cardiac_phase.SystolePositions,respiratory_phase.InspirationPositions); % systole in inspiration
SE = intersect(cardiac_phase.SystolePositions,respiratory_phase.ExpirationPositions); % systole in expiration
DI = intersect(cardiac_phase.DiastolePositions,respiratory_phase.InspirationPositions); % diastole in inspiration
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(SE,'End-systolic frames in expiration');
phase.ShowFrames(DI,'End-diastolic frames in inpiration');
phase.ShowFrames(DE,'End-diastolic frames in expiration');
% c = ImageSequence
% c.ToBeSetByUser=15
% c.ToBeSetByUser=true