classdef Phase < handle
    
    properties 
        MeanIntensity = [];
        Video = [];
        IMFs = [];
        ModePhases = [];
        Values = [];
        InstantaneousFreq = [];
        NumberOfSlices = uint32(0);
        CurrentSlice = uint32(0);
    end
    
    properties % (GetAccess = private) should become private
        nFrames = uint32(0);
        xImageSize = uint16(0);
        yImageSize = uint16(0);
    end
    
    properties(Constant)
        InitialDirectory = 'C:\Uss\Jero�o\Desktop\Tese\Data';
    end
  
    methods %constructor
        function obj=Phase(varargin)
            if nargin==0;return;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
            
            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
        % central point in k-space for all the images
        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)
            figure
            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);
            obj.IMFs = obj.IMFs.';
        end
        
        % Apply Hilber Transform to all the IMFs and extract the
        % instantneous phase from each one of them
        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.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,fps)
            triggered = [];
            triggered = cat(2,triggered,obj.Video(:,:,positions));
            video = implay(triggered,fps);
            set(video.Parent, 'Name', title);
        end
        
        function DisplayValues(obj)
            plot(obj.Values,'-s','LineWidth',0.8,'MarkerSize',4);
            xlabel('Time (sec)');
            ylabel('Phase(rad)');
        end
    
    end
end