PWV_2DPC.m

function varargout = PWV_2DPC(varargin)
    % PWV_2DPC MATLAB code for PWV_2DPC.fig
    %      PWV_2DPC, by itself, creates a new PWV_2DPC or raises the existing
    %      singleton*.
    %
    %      H = PWV_2DPC returns the handle to a new PWV_2DPC or the handle to
    %      the existing singleton*.
    %
    %      PWV_2DPC('CALLBACK',hObject,eventData,handles,...) calls the local
    %      function named CALLBACK in PWV_2DPC.M with the given input arguments.
    %
    %      PWV_2DPC('Property','Value',...) creates a new PWV_2DPC or raises the
    %      existing singleton*.  Starting from the left, property value pairs are
    %      applied to the GUI before PWV_2DPC_OpeningFcn gets called.  An
    %      unrecognized property name or invalid value makes property application
    %      stop.  All inputs are passed to PWV_2DPC_OpeningFcn via varargin.
    %
    %      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
    %      instance to run (singleton)".
    %
    % See also: GUIDE, GUIDATA, GUIHANDLES

    % Edit the above text to modify the response to help PWV_2DPC

    % Last Modified by GUIDE v2.5 15-Nov-2022 10:16:55
    % Developed by Grant S Roberts, University of Wisconsin-Madison, 2019
    
    
    % Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @PWV_2DPC_OpeningFcn, ...
                       'gui_OutputFcn',  @PWV_2DPC_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end

    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % End initialization code - DO NOT EDIT


function PWV_2DPC_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to PWV_2DPC (see VARARGIN)

    % Choose default command line output for PWV_2DPC
    handles.output = hObject;

    % Update handles structure
    guidata(hObject, handles);

    % Get anatomical, pc, and bSSFP data from LoadPWV GUI
    handles.pcDatasets(1).ROI = []; %initialize ROI field
    handles.pcDatasets(1).Interp = []; %throgh-plane flow curve data
    handles.pcDatasets(1).Gaussian = []; %flow smoothed with Gauss.
    handles.pcDatasets(1).Shifted = [];
    
    handles.global.dataType = 'Cartesian';
    handles.global.interpType = 'None'; %interpolation type (i.e. Gaussian)
    handles.global.startAnalyzing = 0; %flag to begin PWV calculations
    handles.global.totalROIs = 0;
    handles.global.pgShift = 0;
    handles.global.pcIter = 1;
    handles.global.analysisType = 'Absolute_Time'; % (or 'Relative_Tme') use timestamps for each acq separately
    handles.global.minTime = NaN;
    %handles.global.analysisType = 'Relative_Time';
    
    set(handles.load2DPCbutton,'Enable','off');
    set(handles.pcPlanePopup,'Enable','off');
    set(handles.pcDatasetPopup,'Enable','off');
    set(handles.drawROIbutton,'Enable','off');
    set(handles.loadROIbutton,'Enable','off');
    set(handles.pcSlider,'Enable','off');
    set(handles.interpolatePopup,'String',{'None','Gaussian','Shifted'}); %set all possible interpolation types
    set(handles.interpolatePopup,'Enable','off'); %initialize radios and buttons
    set(handles.errorBarRadio,'Enable','off');
    set(handles.pgShift,'Enable','off');
    set(handles.ttpointRadio,'Value',1); 
    set(handles.ttpointRadio,'Enable','off');
    set(handles.ttuRadio,'Value',1);
    set(handles.ttuRadio,'Enable','off');
    set(handles.ttfRadio,'Value',1);
    set(handles.ttfRadio,'Enable','off');
    set(handles.xcorrRadio,'Value',1);
    set(handles.xcorrRadio,'Enable','off');
    set(handles.exportAnalysisButton,'Enable','off');
    
    guidata(hObject, handles);

% --- Outputs from this function are returned to the command line.
function varargout = PWV_2DPC_OutputFcn(hObject, eventdata, handles) 
    varargout{1} = handles.output;


    
%%%%%%%%%%%% LOAD 2DPC PLANE %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% --- PLANE PLOT - CREATE FUNCTION
function pcPlanePlot_CreateFcn(hObject, eventdata, handles)


% --- LOAD CENTERLINE DATA - CALLBACK
function loadCLpush_Callback(hObject, eventdata, handles)
    [clFile, clDir] = uigetfile({'*.mat;','Useable Files (*.mat)';
       '*.mat','MAT-files (*.mat)'; ...
       '*.*',  'All Files (*.*)'}, 'Select the centerline dataset (anatCLdataset.mat)');
    load([clDir clFile]);
    handles.centerline = anatCLdataset;
    
    cd(clDir);
    handles.global.homeDir = clDir;
    set(handles.load2DPCbutton,'Enable','on');
    set(handles.planeLoadedText,'String','Centerline Loaded');
    guidata(hObject, handles);
    
    
% --- LOAD 2DPC DATASETS - CALLBACK
function load2DPCbutton_Callback(hObject, eventdata, handles)
    [pcFile, pcDir] = uigetfile({'*.dcm;*.dat;*.mat;*.h5','Useable Files (*.dcm,*.dat,*.mat,*.h5)';
       '*.dcm',  'DICOM files (*.dcm)'; ...
       '*.dat',  'DAT-files (*.dat)'; ...
       '*.mat',  'MAT-files (*.mat)'; ...
       '*.h5',   'HDF5-files (*.h5)'; ...
       '*.*',    'All Files (*.*)'}, 'Select ONE 2DPC file in the dataset');
    pcIter = handles.global.pcIter;
    [~,~,extension] = fileparts(pcFile);
    dirInfo = dir(fullfile(pcDir,['*' extension]));
    if isequal(extension,'.dcm') %if our extension is a dicom file
        handles.global.dataType = 'Cartesian';
        handles.pcDatasets(pcIter).Info = dicominfo(fullfile(pcDir,dirInfo(1).name)); %get dicom metadata (from 1st dicom)
        for i=1:length(dirInfo)
            hold(:,:,i) = single(dicomread(fullfile(pcDir,dirInfo(i).name))); %read dicoms and cast to single
        end  
        mag = hold(:,:,floor(length(dirInfo)/2)+1:end); %magnitude is last half
        v = hold(:,:,1:floor(length(dirInfo)/2)); %velocity is first half of images
        %mag = circshift(mag,10,3); %account for PG gating
        %v = circshift(v,10,3);
        MAG = mean(mag,3); %time-averaged magnitude
        VMEAN = mean(v,3); %time-averaged velocity
        CD = MAG.*sin( pi/2*abs(VMEAN)/max(VMEAN(:)) );
        
        handles.pcDatasets(pcIter).Images.MAG = MAG;
        handles.pcDatasets(pcIter).Images.CD = CD;
        handles.pcDatasets(pcIter).Images.V = VMEAN;
        handles.pcDatasets(pcIter).Images.mag = mag;
        handles.pcDatasets(pcIter).Images.v = v; 
        handles.pcDatasets(pcIter).Names = ['Plane' num2str(pcIter)];
    elseif isequal(extension,'.dat')
        fid = fopen([pcDir filesep 'pcvipr_header.txt'], 'r'); %open header
        dataArray = textscan(fid,'%s%s%[^\n\r]','Delimiter',' ', ...
            'MultipleDelimsAsOne',true,'ReturnOnError',false); %parse header info
        fclose(fid);
        dataArray{1,2} = cellfun(@str2num,dataArray{1,2}(:),'UniformOutput',false);
        pcviprHeader = cell2struct(dataArray{1,2}(:),dataArray{1,1}(:),1); %turn to structure
        handles.pcDatasets(pcIter).Info = pcviprHeader; %add pcvipr header to handles
        resx = pcviprHeader.matrixx; %resolution in x
        resy = pcviprHeader.matrixy; %resolution in y
        nframes = pcviprHeader.frames; %number of cardiac frames
        if nframes<80
            if contains(pcDir,'expiration')
                handles.global.dataType = 'Radial_Expiration';
            else
                handles.global.dataType = 'Radial_LowRes';
            end 
        elseif nframes>=80
            handles.global.dataType = 'Radial_HighRes';
        end 
        MAG = load_dat(fullfile(pcDir,'MAG.dat'),[resx resy]); %Average magnitude
        CD = load_dat(fullfile(pcDir,'CD.dat'),[resx resy]); %Average complex difference
        VMEAN = load_dat(fullfile(pcDir,'comp_vd_3.dat'),[resx resy]); %Average velocity

        % Initialize data time-resolved data arrays
        mag = zeros(resx,resy,nframes); %Time-resolved magnitude
        cd = zeros(resx,resy,nframes); %Time-resolved complex difference
        v = zeros(resx,resy,nframes); %Time-resolved velocity 
        for j = 1:nframes  %velocity is placed in v3 for 2D (through-plane)
            mag(:,:,j) = load_dat(fullfile(pcDir,[filesep 'ph_' num2str(j-1,'%03i') '_mag.dat']),[resx resy]);
            cd(:,:,j) = load_dat(fullfile(pcDir,[filesep 'ph_' num2str(j-1,'%03i') '_cd.dat']),[resx resy]);
            v(:,:,j) = load_dat(fullfile(pcDir,[filesep 'ph_' num2str(j-1,'%03i') '_vd_3.dat']),[resx resy]);
        end
        
        if resy==640
            MAG = MAG(:,161:480);
            CD = CD(:,161:480);
            VMEAN = VMEAN(:,161:480);
            mag = mag(:,161:480,:);
            cd = cd(:,161:480,:);
            v = v(:,161:480,:);
        end 
        handles.pcDatasets(pcIter).Images.MAG = flipud(MAG);
        handles.pcDatasets(pcIter).Images.CD = flipud(CD);
        handles.pcDatasets(pcIter).Images.V = flipud(VMEAN);
        handles.pcDatasets(pcIter).Images.mag = flipud(mag);
        handles.pcDatasets(pcIter).Images.cd = flipud(cd);
        handles.pcDatasets(pcIter).Images.v = flipud(v); 
        handles.pcDatasets(pcIter).Names = ['Plane' num2str(pcIter)];
    elseif isequal(extension,'.h5')
        fid = fopen([pcDir filesep 'pcvipr_header.txt'], 'r'); %open header
        dataArray = textscan(fid,'%s%s%[^\n\r]','Delimiter',' ', ...
            'MultipleDelimsAsOne',true,'ReturnOnError',false); %parse header info
        fclose(fid);
        dataArray{1,2} = cellfun(@str2num,dataArray{1,2}(:),'UniformOutput',false);
        pcviprHeader = cell2struct(dataArray{1,2}(:),dataArray{1,1}(:),1); %turn to structure
        handles.pcDatasets(pcIter).Info = pcviprHeader; %add pcvipr header to handles
        handles.global.dataType = 'Radial_SMS';

        % Initialize data time-resolved data arrays
        mag = h5read(fullfile(pcDir,pcFile),'/MAG'); %Time-resolved magnitude
        cd = h5read(fullfile(pcDir,pcFile),'/CD'); %Time-resolved complex difference
        v = h5read(fullfile(pcDir,pcFile),'/VZ'); %Time-resolved velocity 
        MAG = mean(mag,3); %Average magnitude
        CD = mean(cd,3); %Average complex difference
        VMEAN = mean(v,3); %Average velocity
%         aliasBIN = (v > 1000);
%         nonaliased = ~aliasBIN.*v;
%         aliased = aliasBIN.*v;
%         fixed = (aliased - 3000).*aliasBIN;
%         v = fixed + nonaliased;
        
        handles.pcDatasets(pcIter).Images.MAG = flipud(MAG);
        handles.pcDatasets(pcIter).Images.CD = flipud(CD);
        handles.pcDatasets(pcIter).Images.V = flipud(VMEAN);
        handles.pcDatasets(pcIter).Images.mag = flipud(mag);
        handles.pcDatasets(pcIter).Images.cd = flipud(cd);
        handles.pcDatasets(pcIter).Images.v = flipud(v); 
        handles.pcDatasets(pcIter).Names = ['Plane' num2str(pcIter)];
    else %if a single matlab file (with all images)
        handles.global.isRadial = 1;
        hold = load([pcDir pcFile]);
        planeName = fieldnames(hold);
        planeName = planeName{1};
        handles.pcDatasets(pcIter).Info = hold.(planeName).Info;
        images = hold.(planeName).Images;
        images = single(images);
        
        handles.pcDatasets(pcIter).Images.MAG = mean(images(:,:,:,1),3);
        handles.pcDatasets(pcIter).Images.CD = mean(images(:,:,:,2),3);
        handles.pcDatasets(pcIter).Images.V = mean(images(:,:,:,3),3);
        handles.pcDatasets(pcIter).Images.mag = images(:,:,:,1);
        handles.pcDatasets(pcIter).Images.cd = images(:,:,:,2);
        handles.pcDatasets(pcIter).Images.v = images(:,:,:,3);
        handles.pcDatasets(pcIter).Names = planeName;
    end
    set(handles.planeLoadedText,'String',['Plane #' num2str(handles.global.pcIter) ' Loaded']);
    handles.global.pcIter = handles.global.pcIter + 1;
    
    set(handles.pcPlanePopup,'Enable','on');
    set(handles.pcDatasetPopup,'Enable','on');
    set(handles.drawROIbutton,'Enable','on');
    set(handles.loadROIbutton,'Enable','on');
    set(handles.pcSlider,'Enable','on');
    set(handles.pcPlanePopup,'String',{handles.pcDatasets.Names}); %list of all planes (AAo, AbdAo, etc.)
    set(handles.pcDatasetPopup,'String',fieldnames(handles.pcDatasets(pcIter).Images)); %list of all datasets (CD, MAG, v, etc.)
    guidata(hObject, handles);
    updatePCImages(handles);
    
    
% --- PLANE DROPDOWN - CALLBACK
function pcPlanePopup_Callback(hObject, eventdata, handles)   
    updatePCImages(handles); %update images on PC plot anytime we click on a new plane

% --- PLANE DROPDOWN - CREATE FUNCTION
function pcPlanePopup_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

    
% --- DATASET DROPDOWN - CALLBACK
function pcDatasetPopup_Callback(hObject, eventdata, handles)
    updatePCImages(handles); %update images on PC plot anytime we click on a new dataset

% --- DATASET DROPDOWN - CREATE FUNCTION
function pcDatasetPopup_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    
    
% --- DRAWROI BUTTON - CALLBACK
function drawROIbutton_Callback(hObject, eventdata, handles)
    axes(handles.pcPlanePlot); %make sure we do this on top PC plot 

    set(handles.loadCLpush,'Enable','off');
    set(handles.load2DPCbutton,'Enable','off');
    set(handles.pcPlanePopup,'Enable','off'); %make it so we can't select another plane
    set(handles.pcDatasetPopup,'Enable','off'); %make it so we can't select another plane
    set(handles.drawROIbutton,'Enable','off'); %make it so we can't draw a second ROI
    
    planeNum = get(handles.pcPlanePopup,'Value'); %get current plane of interest (eg AAo)
    mydlg = warndlg('Press enter when the ROI is set'); %open dialog warning 
    waitfor(mydlg); %MUST PRESS ENTER TO PROCEED
    circle = drawcircle('FaceAlpha',0.1,'Color','g','LineWidth',1,'Deletable',0); %draw circle on PC image
    while true
        w = waitforbuttonpress; %wait for enter push ...
        switch w 
            case 1 % if it was a keyboard press.
            key = get(gcf,'currentcharacter'); %get key that was pressed
                switch key
                    case 27 % escape key
                        set(handles.loadCLpush,'Enable','on');
                        set(handles.load2DPCbutton,'Enable','on');
                        set(handles.pcPlanePopup,'Enable','on'); %make it so we can't select another plane
                        set(handles.pcDatasetPopup,'Enable','on'); %make it so we can't select another plane
                        set(handles.drawROIbutton,'Enable','on');
                        broke = 1;
                        break % break out of the while loop
                    case 13 % 13 is the enter/return key 
                        circle.InteractionsAllowed = 'none'; %freeze circle
                        broke = 0;
                        break
                    otherwise 
                        %wait for a different command
                end
       end
    end
    
    if ~broke
        handles.pcDatasets(planeNum).ROI = circle; %temporarily hold circle (deleted once gone)
    else 
        handles.pcDatasets(planeNum).ROI = [];
    end 

    guidata(hObject,handles);
    updatePCImages(handles); 

    
% --- LOAD ROI BUTTON - CALLBACK
function loadROIbutton_Callback(hObject, eventdata, handles)
    planeNum = get(handles.pcPlanePopup,'Value'); %get current plane (eg AAo)
    handles.global.totalROIs = handles.global.totalROIs + 1; %add 1 to total ROI count

    v = handles.pcDatasets(planeNum).Images.v; %grab time-resolved velocity
    circle = handles.pcDatasets(planeNum).ROI; %pull circle data from handles
    radius = circle.Radius; %get radius of circle
    center = round(circle.Center); %get center coordinates

    [X,Y] = ndgrid(1:size(v,1),1:size(v,2));
    X = X-center(2); %shift coordinate grid
    Y = Y-center(1);
    mask = sqrt(X.^2 + Y.^2)<=radius; %anything outside radius is ignored

    %%% Create velocity/flow curves
    if isfield(handles.pcDatasets(planeNum).Info,'matrixx') %if radial data (pcvipr recon)
        matrixx = handles.pcDatasets(planeNum).Info.matrixx; %matrix size in x dimension
        fovx = handles.pcDatasets(planeNum).Info.fovx;  %field of view (mm)
        xres = fovx/matrixx; %resolution (mm). ASSUMED TO BE SAME IN Y DIMENSION
        frames = handles.pcDatasets(planeNum).Info.frames;
        timeres = handles.pcDatasets(planeNum).Info.timeres; %temporal resolution (ms)
    else 
        xres = handles.pcDatasets(planeNum).Info.PixelSpacing(1); %resolution (mm) ASSUMED SAME IN Y DIM
        rr_interval = handles.pcDatasets(planeNum).Info.NominalInterval; %average RR interval (ms)
        frames = handles.pcDatasets(planeNum).Info.CardiacNumberOfImages; %number of cardiac frames
        timeres = rr_interval/frames; %temporal resolution (ms)
    end 
    
    area = sum(mask(:))*(xres)^2; %ROI area (mm^2)
    for i=1:frames
        vTemp = v(:,:,i); %through-plane velocity in frame i
        roi_data_raw(:,i) = double(vTemp(mask)); %indexed velocities within mask
        mean_roi(i) = mean(roi_data_raw(:,i)); %mean velocity in frame i (mm/s)
        stdv_roi(i) = std(double(vTemp(mask))); %stdv of velocity in frame i (mm/s)
        flow_roi(i) = area.*mean_roi(i).*0.001; %flow in frame i (mm^3/s*0.001 = mL/s)
    end 
    times = double(timeres.*(0:(frames-1))); %original times
    
    % Unrap Code (use in for loop above as needed)
%     newPlane = vTemp.*mask;
%     newPlane = (newPlane/1500)*pi;
%     UW = Unwrap_TIE_DCT_Iter(newPlane);
%     UW = (UW/pi)*-1500;
%     figure; imshowpair(newPlane,UW,'montage');
%     plane = vTemp.*(~mask);
%     vTemp = -plane - UW;
    
    %%% Add ROI info and raw (uninterpolated) data to Raw structure
    Raw.radius = radius; 
    Raw.center = center;
    Raw.mask = mask; 
    Raw.raw_data = roi_data_raw;
    Raw.name = ''; %will get changed below
    Raw.roi_number = handles.global.totalROIs;
    Raw.times = times;
    Raw.frames = frames;
    Raw.timeres = timeres;
    Raw.mean_v = mean_roi;
    Raw.stdv_v = stdv_roi;
    Raw.flow = flow_roi;  

    %%% Linear interpolation (to get more points on flow curve)
    sample_density = 1000;
    tq = linspace(0,frames-1,sample_density); %interpolate time dimension
    times_interp = double(timeres.*tq); %interpolated times
    
    %set temporal resolution to 1ms, chop ends of waveforms
%     times_interp = double(0:1:timeres*(frames-1)); %interpolate to 1ms
%     if isnan(handles.global.minTime)
%         handles.global.minTime = int16(length(times_interp) - 0.05*length(times_interp));
%         times_interp = times_interp(1:handles.global.minTime);
%     else
%         times_interp = times_interp(1:handles.global.minTime);
%     end 

    mean_interp = interp1(times,mean_roi,times_interp,'linear');
    stdv_interp = interp1(times,stdv_roi,times_interp,'linear');
    flow_interp = interp1(times,flow_roi,times_interp,'linear');
    
    Interp.times = times_interp;
    Interp.sample_density = sample_density;
    Interp.mean_v = mean_interp; 
    Interp.stdv_v = stdv_interp;
    Interp.flow = flow_interp; 

    %%% Create Interpolated Curve with Gaussian Smoothing   
    if frames<=40
        smf = 80; %smooth factor (window)
    elseif frames<=80
        smf = 100;
    else
        smf = 120;
    end 
    mean_gaus = smoothdata(mean_interp,'gaussian',smf);
    stdv_gaus = smoothdata(stdv_interp,'gaussian',smf);
    flow_gaus = smoothdata(flow_interp,'gaussian',smf);

    Gaussian.times = times_interp;
    Gaussian.smf = smf;
    Gaussian.mean_v = mean_gaus; 
    Gaussian.stdv_v = stdv_gaus; 
    Gaussian.flow = flow_gaus; 
    
    %%% Shift Raw Data for PPG gating (then interpolate, then smooth)
    shift = round((handles.global.pgShift/100)*frames);  %turn % into # frames to shift on raw data
    mean_roi = circshift(Raw.mean_v,shift);
    stdv_roi = circshift(Raw.stdv_v,shift);
    flow_roi = circshift(Raw.flow,shift);
    mean_interp = interp1(times,mean_roi,times_interp,'linear');
    stdv_interp = interp1(times,stdv_roi,times_interp,'linear');
    flow_interp = interp1(times,flow_roi,times_interp,'linear');
    
    Shifted.times = times_interp;
    Shifted.mean_v = smoothdata(mean_interp,'gaussian',smf);
    Shifted.stdv_v = smoothdata(stdv_interp,'gaussian',smf);
    Shifted.flow = smoothdata(flow_interp,'gaussian',smf);
    

    %%% Save all data into handles
    if isstruct(handles.pcDatasets(planeNum).Interp)
        handles.pcDatasets(planeNum).Raw(end+1) = Raw;
        handles.pcDatasets(planeNum).Interp(end+1) = Interp;
        handles.pcDatasets(planeNum).Gaussian(end+1) = Gaussian;
        handles.pcDatasets(planeNum).Shifted(end+1) = Shifted;
    else
        handles.pcDatasets(planeNum).Raw = Raw;
        handles.pcDatasets(planeNum).Interp = Interp;
        handles.pcDatasets(planeNum).Gaussian = Gaussian;
        handles.pcDatasets(planeNum).Shifted = Shifted;
    end 

    dataDir = handles.global.homeDir; %directory in which plane data is located
    if dataDir(end)=='\' || dataDir(end)=='/' %kill the slash if it exists
        dataDir(end) = [];
    end 
    if ~exist([dataDir filesep 'ROIimages_' handles.global.dataType],'dir') %if the proposed directory doesn't exist
        mkdir([dataDir filesep 'ROIimages_' handles.global.dataType]); %make it
        cd([dataDir filesep 'ROIimages_' handles.global.dataType]); %move into it
        frame = getframe(handles.pcPlanePlot); %get a snapshot of the PC plane plot with ROI
        image = frame2im(frame); %make into image
        imwrite(image,[handles.pcDatasets(planeNum).Names '.png']) %write it out as PNG
    else
        cd([dataDir filesep 'ROIimages_' handles.global.dataType]); %if ROIimages already exists, move into it
        frame = getframe(handles.pcPlanePlot);
        image = frame2im(frame);
        imwrite(image,[handles.pcDatasets(planeNum).Names '.png'])
    end 
    cd(handles.global.homeDir); %lets go back home
    
    %%% Label each ROI w/ names (helpful because there may be 2 ROIs/plane)
    for i=1:numel(handles.pcDatasets)
        if isstruct(handles.pcDatasets(i).Raw) %if we've made ROI data for this dataset
            if length(handles.pcDatasets(i).Raw)==1
                handles.pcDatasets(i).Raw.name = handles.pcDatasets(planeNum).Names; %name ROI
            else
                for j=1:length(handles.pcDatasets(i).Raw)
                    name = handles.pcDatasets(i).Names; %get plane name
                    planeName = [name ' ROI ' num2str(j)]; %needed if more than one ROI/plane
                    handles.pcDatasets(i).Raw(j).name = planeName; %name ROI
                end 
            end 
        end 
    end 

    set(handles.interpolatePopup,'Enable','on'); %turn on interpolate button
    plotVelocity(handles); %plot flow curves
    
    set(handles.loadCLpush,'Enable','on');
    set(handles.load2DPCbutton,'Enable','on');
    set(handles.pcPlanePopup,'Enable','on'); %make it so we can't select another plane
    set(handles.pcDatasetPopup,'Enable','on'); %make it so we can't select another plane
    set(handles.drawROIbutton,'Enable','on');
    set(handles.interpolatePopup,'Enable','on'); %set all possible interpolation types
    set(handles.errorBarRadio,'Enable','on');

    guidata(hObject,handles);
    updatePCImages(handles); %update images (to remove green ROI circle)
    axes(handles.pcPlanePlot); %make sure we're still on PC plot
   

    % --- PLANE SLIDER - CALLBACK
function pcSlider_Callback(hObject, eventdata, handles)
    updatePCImages(handles); %update images if slider is moved

% --- PLANE SLIDER - CREATE FUNCTION
function pcSlider_CreateFcn(hObject, eventdata, handles)
    if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor',[.9 .9 .9]);
    end
    
    
% --- MINIMUM CONTRAST VALUE BOX - CALLBACK  
function minContrastBox_Callback(hObject, eventdata, handles)
    updatePCImages(handles)

% --- MINIMUM CONTRAST VALUE BOX - CREATE FUNCTION   
function minContrastBox_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


% --- MAXIMUM CONTRAST VALUE BOX - CALLBACK  
function maxContrastBox_Callback(hObject, eventdata, handles)
    updatePCImages(handles)
    
% --- MAXIMUM CONTRAST VALUE BOX - CREATE FUNCTION
function maxContrastBox_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
    
    
%%%%%%%%%%%% VELOCITY PLOT %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% --- VELOCITY PLOT - CREATE FUNCTION
function velocityPlot_CreateFcn(hObject, eventdata, handles)
    

% --- INTERPOLATE POPUP - CALLBACK
function interpolatePopup_Callback(hObject, eventdata, handles)
    interp = get(handles.interpolatePopup,'Value');
    switch interp
        case 1
            handles.global.interpType = 'None'; %set global flag
            set(handles.pgShift,'Enable','off');
        case 2
            handles.global.interpType = 'Gaussian';
            set(handles.pgShift,'Enable','off');
        case 3 
            handles.global.interpType = 'Shifted';
            set(handles.pgShift,'Enable','on');
    end 
 
    guidata(hObject, handles);
    plotVelocity(handles) %replot our velocity with interpolated data
    if handles.global.startAnalyzing %if we're already analyzing PWVs
        completeLoadingROI_Callback(hObject, eventdata, handles); %recompute PWVs with interpolated data
    end 

% --- INTERPOLATE POPUP - CREATE FUNCTION
function interpolatePopup_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end
    
    
% --- ERROR BAR RADIO - CALLBACK
function errorBarRadio_Callback(hObject, eventdata, handles)
    plotVelocity(handles) %replot flow curves with errors bars
    
    
% --- PG SHIFT RADIO - CALLBACK
function pgShift_Callback(hObject, eventdata, handles)
    pgShift = str2double(get(handles.pgShift,'String')); %percent shift
    handles.global.pgShift = pgShift; 
    
    frames = handles.pcDatasets(1).Raw.frames;
    shift = round((pgShift/100)*frames);  %turn % into # frames to shift on raw data
    for i=1:numel(handles.pcDatasets) %for all planes
        if isstruct(handles.pcDatasets(i).Raw) %if we've made ROI data for this dataset
            for j=1:length(handles.pcDatasets(i).Raw) %for each ROI
                roi = handles.pcDatasets(i).Raw(j);
                shifted = handles.pcDatasets(i).Shifted(j);
                times = roi.times;
                times_interp = handles.pcDatasets(i).Interp(j).times;
                smf = handles.pcDatasets(i).Gaussian(j).smf;

                mean_roi = circshift(roi.mean_v,shift);
                stdv_roi = circshift(roi.stdv_v,shift);
                flow_roi = circshift(roi.flow,shift);
                mean_interp = interp1(times,mean_roi,times_interp,'linear');
                stdv_interp = interp1(times,stdv_roi,times_interp,'linear');
                flow_interp = interp1(times,flow_roi,times_interp,'linear');
                
                shifted.mean_v = smoothdata(mean_interp,'gaussian',smf);
                shifted.stdv_v = smoothdata(stdv_interp,'gaussian',smf);
                shifted.flow = smoothdata(flow_interp,'gaussian',smf);

                %%% Save all data into handles
                handles.pcDatasets(i).Shifted(j) = shifted;
            end
        end
    end 
    guidata(hObject, handles);
    
    plotVelocity(handles) %replot flow curves with half cycle shift
    if handles.global.startAnalyzing %if we're already analyzing PWVs
        completeLoadingROI_Callback(hObject, eventdata, handles); %recompute PWVs with shifted curves
    end 
    guidata(hObject, handles);
    
% --- PG SHIFT RADIO - CREATE FUNCTION
function pgShift_CreateFcn(hObject, eventdata, handles)
    if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
        set(hObject,'BackgroundColor','white');
    end

   
    
% --- Executes on button press in completeLoadingROI.
function completeLoadingROI_Callback(hObject, eventdata, handles)
    handles.global.startAnalyzing = 1; %turn on flag to state that we are ready for PWV analysis
    handles.flow = organizeFlowInfo(handles);
    
    % COMPUTE TIME SHIFTS
    if ~isfield(handles.flow,'TTUpstroke')
        flow = computeTTs(handles.flow,handles.global);
        handles.flow = flow;
    end 
    
    % COMPUTE PWVs
    numCompares = numel(flow);
    distance = [0 cumsum(handles.centerline.PlaneDistances)];
    distance = distance(1:numCompares);
    TTPoint = [handles.flow.TTPoint]; 
    TTFoot = [handles.flow.TTFoot]; 
    TTUpstroke = [handles.flow.TTUpstroke]; 
    Xcorr = [handles.flow.Xcorr];

    numMethods = get(handles.ttpointRadio,'Value') + ...
        get(handles.ttfRadio,'Value') + ...
        get(handles.ttuRadio,'Value') + ...
        get(handles.xcorrRadio,'Value'); %add all PWV buttons turned on
    %numCompares = numel(distance); %get number of time shift methods  
    
    average = zeros(1,numCompares); %initialize average timeshift array
    
    cla(handles.TimeVsDistance,'reset'); %reset PWV plot
    axes(handles.TimeVsDistance); hold on; %force axes to PWV plot
    xlabel('Centerline Distance (mm)'); ylabel('Time Shift (ms)'); %label axes
    sz = 30; %create marker sizes (filled in circles) of 30 pixels
    legendSet = {}; %initialize legend cell array
    if get(handles.ttpointRadio,'Value')
        scatter(distance,TTPoint,sz,'filled','MarkerFaceColor',[0.850 0.325 0.098]); %orange 
        for i=1:numCompares
            average(i) = average(i)+TTPoint(i); %add all distances for each ROI location
        end 
        legendSet{end+1} = 'TTPoint';
    end 
    if get(handles.ttfRadio,'Value')
        scatter(distance,TTFoot,sz,'filled','MarkerFaceColor',[0.494 0.184 0.556]); %purple 
        for i=1:numCompares
            average(i) = average(i)+TTFoot(i); %keep adding distances
        end 
        legendSet{end+1} = 'TTFoot';
    end 
    if get(handles.ttuRadio,'Value')
        scatter(distance,TTUpstroke,sz,'filled','MarkerFaceColor',[0.4660 0.8740 0.1880]); %green
        for i=1:numCompares
            average(i) = average(i)+TTUpstroke(i);
        end 
        legendSet{end+1} = 'TTUpstroke';
    end 
    if get(handles.xcorrRadio,'Value')
        scatter(distance,Xcorr,sz,'filled','MarkerFaceColor',[0.8350 0.0780 0.1840]); %red 
        for i=1:numCompares 
            average(i) = average(i)+Xcorr(i);
        end 
        legendSet{end+1} = 'XCorr';
    end 

    D = max(distance);
    Tmin = min([TTPoint TTFoot TTUpstroke Xcorr]);
    Tmax = max([TTPoint TTFoot TTUpstroke Xcorr]);
    xlim([0 D+D*0.2]); ylim([Tmin-Tmin*0.02 Tmax+Tmax*0.02]);
        
    average = average./numMethods; %get average TT for each ROI
    scatter(distance,average,40,'black'); %open black circles (size 40 pixels)
    legendSet{end+1} = 'AVERAGE'; %add average to legend
    
    d = 0:round(D+D*0.2);
    % Here, we perform linear regression to find best fit line for each
    % time-to (TT) method. Note: lineFit(1)=slope; lineFit(2)=intercept

    %TTPoint Method
    PointFit = polyfit(distance,TTPoint,1); %fit line to TTPoint points
    linePoint = PointFit(1)*d + PointFit(2); %calculate line

    %TTFoot Method
    FootFit = polyfit(distance,TTFoot,1);
    lineFoot = FootFit(1)*d + FootFit(2);

    %TTUpstroke Method
    UpstrokeFit = polyfit(distance,TTUpstroke,1);
    lineUpstroke = UpstrokeFit(1)*d + UpstrokeFit(2);

    %Xcorr Method
    XcorrFit = polyfit(distance,Xcorr,1);
    lineXcorr = XcorrFit(1)*d + XcorrFit(2);

    %Average TT
    AverageFit = polyfit(distance,average,1);
    lineAverage= AverageFit(1)*d + AverageFit(2);

    PWVpoint = 1/PointFit(1); %PWV = 1/slope (mm/ms = m/s)
    PWVfoot = 1/FootFit(1);
    PWVupstroke = 1/UpstrokeFit(1);
    PWVxcorr = 1/XcorrFit(1);
    PWVaverage = 1/AverageFit(1);

    hold on; 
    if get(handles.ttpointRadio,'Value') %if our ttpoint button is on, plot average ttp
        plot(d,linePoint,':','LineWidth',0.2,'Color',[0.8500 0.3250 0.0980]); %orange 
        set(handles.ttpointData,'String',[sprintf('%0.2f',PWVpoint) ' m/s']); %write out PWV value in text field
    end 
    if get(handles.ttfRadio,'Value')
        plot(d,lineFoot,':','LineWidth',0.2,'Color',[0.4940 0.1840 0.5560]); %purple
        set(handles.ttfData,'String',[sprintf('%0.2f',PWVfoot) ' m/s']);
    end 
    if get(handles.ttuRadio,'Value')
        plot(d,lineUpstroke,':','LineWidth',0.2,'Color',[0.4660 0.8740 0.1880]); %green
        set(handles.ttuData,'String',[sprintf('%0.2f',PWVupstroke) ' m/s']);
    end 
    if get(handles.xcorrRadio,'Value')
        plot(d,lineXcorr,':','LineWidth',0.2,'Color',[0.8350 0.0780 0.1840]); %red
        set(handles.xcorrData,'String',[sprintf('%0.2f',PWVxcorr) ' m/s']);
    end 
    plot(d,lineAverage,'-k','LineWidth',0.2);
    legend(legendSet,'Location','northwest');
    hold off; 
    
    if numMethods>0 %if we have at least one ttbutton on
        set(handles.averageData,'String',[sprintf('%0.2f',PWVaverage) ' m/s']); %set PWV text field
    else %if have not methods selected (all ttbuttons are off)
        set(handles.averageData,'String','0 m/s'); %set PWV text field to 0 m/s
    end
    
    set(handles.ttpointRadio,'Enable','on');
    set(handles.ttuRadio,'Enable','on');
    set(handles.ttfRadio,'Enable','on');
    set(handles.xcorrRadio,'Enable','on');
    set(handles.exportAnalysisButton,'Enable','on');
    guidata(hObject, handles);
    
    

%%%%%%%%%%%% PWV PLOT %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% --- PWV PLOT - CREATE FUNCTION
function TimeVsDistance_CreateFcn(hObject, eventdata, handles)


% --- EXPORT ANALYSIS - CALLBACK
function exportAnalysisButton_Callback(hObject, eventdata, handles)
    date = datestr(now); %get current date/time
    chopDate = [date(1:2) '-' date(4:6) '-' date(10:11) '-' date(13:14) date(16:17)]; %chop date up
    globals = handles.global;
    
    for a=1:2
        flow = computeTTs(handles.flow,globals); %recalculate flow struct for each interp
        numCompares = numel(flow);
        handles.flow = flow;
        guidata(hObject, handles);

        distance = [0 cumsum(handles.centerline.PlaneDistances)];
        distance = distance(1:numCompares);
        ttpoint = [handles.flow.TTPoint]; 
        ttfoot = [handles.flow.TTFoot];
        ttupstroke = [handles.flow.TTUpstroke]; 
        xcorr = [handles.flow.Xcorr];
        ttaverage = mean([ttfoot; ttpoint; ttupstroke; xcorr],1); %get average time shift

        Distances = abs(diff([distance 0])); %add zero at end to get full distance
        TTPoint = [diff(ttpoint) sum(diff(ttpoint))];
        TTFoot = [diff(ttfoot) sum(diff(ttfoot))];
        TTUpstroke = [diff(ttupstroke) sum(diff(ttupstroke))];
        Xcorr = [diff(xcorr) sum(diff(xcorr))];
        TTaverage = [diff(ttaverage) sum(diff(ttaverage))];

        % Make first column for excel file
        %numCompares = length(Distances); %get number of PWV measurements
        for d=1:(numCompares-1)
            PLANES{d} = [flow(d).Name ' --> ' flow(d+1).Name]; %get names for Excel
        end 
        PLANES{numCompares} = [flow(1).Name ' --> ' flow(end).Name]; %get names for Excel
        PLANES{numCompares+1} = 'FIT_PWV'; %These rows be for PWV fit parameters
        PLANES{numCompares+2} = 'm (slope)';
        PLANES{numCompares+3} = 'b (y-intercept)';
        PLANES{numCompares+4} = 'R^2';

        PWV_Point = NaN(numCompares+4,1); %add dummy rows to match PLANES size
        PWV_Foot = NaN(numCompares+4,1);
        PWV_Upstroke = NaN(numCompares+4,1);
        PWV_Xcorr = NaN(numCompares+4,1);
        PWV_Average = zeros(1,numCompares);

        for i=1:numCompares
            PWV_Point(i) = Distances(i)/TTPoint(i); %calculate pointwise PWVs (not fits)
            PWV_Foot(i) = Distances(i)/TTFoot(i);
            PWV_Upstroke(i) = Distances(i)/TTUpstroke(i);
            PWV_Xcorr(i) = Distances(i)/Xcorr(i);
            PWV_Average(i) = (PWV_Point(i)+PWV_Foot(i)+PWV_Upstroke(i)+PWV_Xcorr(i))/4; %get simple average of all PWVs
        end 

        [linePointFit,S] = polyfit(distance,ttpoint,1); %get linear regression fit for all points 
        PWV_Point(numCompares+1) = 1/linePointFit(1); %calculate PWV (=1/slope)
        PWV_Point(numCompares+2) = linePointFit(1); %get slope
        PWV_Point(numCompares+3) = linePointFit(2); %get y-intercept
        PWV_Point(numCompares+4) = (1 - (S.normr/norm(ttpoint - mean(ttpoint)))^2); %R^2

        [lineFootFit,S] = polyfit(distance,ttfoot,1);
        PWV_Foot(numCompares+1) = 1/lineFootFit(1);
        PWV_Foot(numCompares+2) = lineFootFit(1);
        PWV_Foot(numCompares+3) = lineFootFit(2);
        PWV_Foot(numCompares+4) = (1 - (S.normr/norm(ttfoot - mean(ttfoot)))^2);

        [lineUpstrokeFit,~] = polyfit(distance,ttupstroke,1);
        PWV_Upstroke(numCompares+1) = 1/lineUpstrokeFit(1);
        PWV_Upstroke(numCompares+2) = lineUpstrokeFit(1);
        PWV_Upstroke(numCompares+3) = lineUpstrokeFit(2);
        PWV_Upstroke(numCompares+4) = (1 - (S.normr/norm(ttupstroke - mean(ttupstroke)))^2);

        [lineXcorrFit,~] = polyfit(distance,xcorr,1);
        PWV_Xcorr(numCompares+1) = 1/lineXcorrFit(1);
        PWV_Xcorr(numCompares+2) = lineXcorrFit(1);
        PWV_Xcorr(numCompares+3) = lineXcorrFit(2);
        PWV_Xcorr(numCompares+4) = (1 - (S.normr/norm(xcorr - mean(xcorr)))^2);

        [lineAverageFit,~] = polyfit(distance,ttaverage,1);
        PWV_Average(numCompares+1) = 1/lineAverageFit(1);
        PWV_Average(numCompares+2) = lineAverageFit(1);
        PWV_Average(numCompares+3) = lineAverageFit(2);
        PWV_Average(numCompares+4) = (1 - (S.normr/norm(ttaverage - mean(ttaverage)))^2);
        PWV_Average = PWV_Average';

        Distances(numCompares+1:numCompares+4) = NaN;
        TTPoint(numCompares+1:numCompares+4) = NaN;
        TTFoot(numCompares+1:numCompares+4) = NaN;
        TTUpstroke(numCompares+1:numCompares+4) = NaN;
        Xcorr(numCompares+1:numCompares+4) = NaN;
        TTaverage(numCompares+1:numCompares+4) = NaN;

        PLANES = PLANES'; %Needed for excel, won't save name if ' in table call
        Distances = Distances';
        TTPoint = TTPoint';
        TTFoot = TTFoot';
        TTUpstroke = TTUpstroke';
        Xcorr = Xcorr';
        TTaverage = TTaverage';

        % Make table for writing excel file
        pwvTable = table(PLANES,Distances,TTPoint,TTFoot,TTUpstroke,Xcorr,TTaverage,PWV_Point,PWV_Foot,PWV_Upstroke,PWV_Xcorr,PWV_Average);
        baseDir = globals.homeDir; %rejoin string to get name of folder one up from plane data
        dataDir = ['DataAnalysis_' globals.dataType];

        set(handles.exportDone,'String',['Saving Dataset ' num2str(a) ' ...']);
        if ~exist([baseDir filesep dataDir],'dir') %if directory doesn't exist
            mkdir([baseDir filesep dataDir]); %make it
        end 
        cd([baseDir filesep dataDir]); %go to it
        mkdir(globals.analysisType);
        cd(globals.analysisType);
        writetable(pwvTable,['Summary_' chopDate '.xlsx'],'FileType','spreadsheet'); %write excel sheet for each interp
        %if strcmp(globals.interpType,'Gaussian')
            saveTTplots(handles,flow);
        %end 
        %if ~exist('flow.mat','file')
            save('flow.mat','flow')
            save('pwvTable.mat','pwvTable');
            save('globals.mat','globals');
        %end 
        cd(globals.homeDir); %go back home  
        clear PLANES %need to do this because PLANES will keep getting transposed

        cd([baseDir filesep dataDir]); %go to it
        frame = getframe(handles.TimeVsDistance); %get snapshot of PWV plot 
        imwrite(frame2im(frame),'PWVanalysisPlot.png'); %write out to PNG
        pcDatasets = handles.pcDatasets;
        save('pcDatasets.mat','pcDatasets');
        cd(globals.homeDir); %go back home  
        if ~exist('PWV_2DPC_Analysis-v2','dir')
            mkdir('PWV_2DPC_Analysis-v2');
        end 
        movefile(dataDir,'PWV_2DPC_Analysis-v2');
        handles.globals.analysisType = 'Relative_Time';
        globals.analysisType = 'Relative_Time';
        guidata(hObject, handles);
    end 
    movefile(['ROIimages_' handles.global.dataType],'PWV_2DPC_Analysis-v2');
    set(handles.exportDone,'String','Export Completed!');

    guidata(hObject, handles);
    
   

% --- TTPoint READOUT - CREATE FUNCTION
function ttpointData_CreateFcn(hObject, eventdata, handles)
% --- TTPoint RADIO - CALLBACK
function ttpointRadio_Callback(hObject, eventdata, handles)
    if ~get(handles.ttpointRadio,'Value') %if we're turned off
        set(handles.ttpointData,'String',' '); %don't display PWV
    end 
    if handles.global.startAnalyzing %if we're analyzing PWVs
        completeLoadingROI_Callback(hObject, eventdata, handles); %reanalyze without TTpoint 
    end 

% --- TTUpstroke READOUT - CREATE FUNCTION
function ttuData_CreateFcn(hObject, eventdata, handles)
% --- TTUpstroke RADIO - CALLBACK
function ttuRadio_Callback(hObject, eventdata, handles)
    if ~get(handles.ttuRadio,'Value')
        set(handles.ttuData,'String',' ');
    end 
    if handles.global.startAnalyzing
        completeLoadingROI_Callback(hObject, eventdata, handles);
    end 

% --- TTFoot READOUT - CREATE FUNCTION
function ttfData_CreateFcn(hObject, eventdata, handles)
% --- TTFoot RADIO - CALLBACK
function ttfRadio_Callback(hObject, eventdata, handles)
    if ~get(handles.ttfRadio,'Value')
        set(handles.ttfData,'String',' ');
    end 
    if handles.global.startAnalyzing
        completeLoadingROI_Callback(hObject, eventdata, handles);
    end 

% --- Xcorr READOUT - CREATE FUNCTION
function xcorrData_CreateFcn(hObject, eventdata, handles)
% --- Xcorr RADIO - CALLBACK
function xcorrRadio_Callback(hObject, eventdata, handles)
    if ~get(handles.xcorrRadio,'Value')
        set(handles.xcorrData,'String',' ');
    end 
    if handles.global.startAnalyzing
        completeLoadingROI_Callback(hObject, eventdata, handles);
    end 

% --- AVERAGE READOUT - CREATE FUNCTION
function averageData_CreateFcn(hObject, eventdata, handles)





%%%%%%%%%%%% MY FUNCTIONS %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

% --- Update Images in PLANE PLOT
function updatePCImages(handles)
    axes(handles.pcPlanePlot); %force axes to PC plot
    
    planeNum = get(handles.pcPlanePopup,'Value'); % get current plane index (eg AAo)
    datasetNum = get(handles.pcDatasetPopup,'Value'); %get current dataset (eg MAG)
    
    dataset = handles.pcDatasets;
    if isstruct(dataset(planeNum).Images) %if Data is a structure
        imageSet = struct2cell(dataset(planeNum).Images); %convert from struct to cell
    else  %or if we already have an array
        imageSet = dataset(planeNum).Images; %just change the name
    end 

    if iscell(imageSet) %if our set of images are contained in a cell
        images = imageSet(datasetNum); %pull images for current dataset
        images = cell2mat(images); %turn to matrix
    else 
        images = imageSet; %otherwise, do nothing
    end 
    
    if ndims(images)<3 %if we are dealing with time-averaged images (ndim=2)
        maxSize = max(size(images,1),size(images,2));
        steps = [1 maxSize]; %set our slider as wide as possible so we can't slide
        set(handles.pcSlider,'SliderStep', steps);
        slice = images; %change name for consistency (see below)
    else %if we are dealing with time-resolved images (ndim=3)
        dim3size = size(images,3); %get the size of the third dimension
        steps = [1/(dim3size-1) 10/(dim3size-1)]; %set so one 'slide' moves to the next slice exactly
        set(handles.pcSlider,'SliderStep', steps);
        sliceNum = 1+round( get(handles.pcSlider,'Value').*(dim3size-1) ); %get slice number from slider
        slice = images(:,:,sliceNum); %pull slice from images
    end 
    
    minc = str2double(get(handles.minContrastBox,'String'));
    maxc = str2double(get(handles.maxContrastBox,'String'));

    if ~isempty(handles.pcDatasets(planeNum).ROI) %if we have an ROI placed
        hold on %make it so we can slide while keeping the ROI on the figure
        imshow(rescale(slice),[minc maxc]);
    else 
        cla(handles.pcPlanePlot,'reset') %otherwise, reset the plot
        imshow(rescale(slice),[minc maxc]) %then show the image
    end 
    
    
    
% --- "Time to" calculations (TTPoint, TTUpstroke, TTFoot, Xcorr)
function flow = computeTTs(flow,globals)
    numROIs = globals.totalROIs;
    for i=1:numROIs %for each ROI
        if globals.startAnalyzing %if we're analyzing (put here because we call computePWV in interpolatePopup)
            switch globals.interpType  %find the appropriate interp data
                case 'None'
                    flowTemp = flow(i).Interp.mean_v;
                case 'Gaussian'
                    flowTemp = flow(i).Gaussian.mean_v;  
                case 'Shifted'
                    flowTemp = flow(i).Shifted.mean_v;  
            end 
        else %else, we just use shifted data
            flowTemp = flow(i).Shifted.mean_v;  
        end 
        
        if mean(flowTemp)<0 %if our flow curve is mainly negative
            flowTemp = -1*flowTemp; %flip it upside down so we can find time shifts
        end 
        
        if strcmp(globals.analysisType,'Absolute_Time')
            times = flow(i).Interp.times; %get individual time frames (ms)
        elseif strcmp(globals.analysisType,'Relative_Time')
            times = flow(1).Interp.times; %use only 1st time frames (ms)
            flow(i).Raw.times = flow(1).Raw.times;
            flow(i).Interp.times = flow(1).Interp.times;
            flow(i).Gaussian.times = flow(1).Gaussian.times;
            flow(i).Shifted.times = flow(1).Shifted.times;
        else
            disp('Please select correct data analysis type (Absolute_Time or Relative_Time)');
        end 
        timeres = times(2)-times(1); %temporal resolution (ms)
        curvePoints(i).times = times;
        curvePoints(i).timeres = timeres;

        [maxPeakVel,maxPeakVelIdx] = max(flowTemp); %find max velocity value and its location
        upstroke = flowTemp(1:maxPeakVelIdx); %define 'upstroke' region of the flow curve
        curvePoints(i).maxPeakVelIdx = maxPeakVelIdx; %add max point to curvePoints struct
        curvePoints(i).maxPeakVel = maxPeakVel; %add max velocity to curvePoints struct
        
        [minPeakVel,minPeakVelIdx] = min(upstroke); %find max velocity value and its location
        curvePoints(i).minPeakVelIdx = minPeakVelIdx; %add max point to curvePoints struct
        curvePoints(i).minPeakVel = minPeakVel; %add max velocity to curvePoints struct
        
        upstroke = normalize(upstroke,'range'); %baseline correct
        
        [~,SeventyPointIdx] = min(abs(upstroke - 0.70)); %get point at 70% max peak
        curvePoints(i).SeventyPointIdx = SeventyPointIdx; %add to curvePoints struct
        curvePoints(i).SeventyPoint = flowTemp(SeventyPointIdx); %add 70% flow value to curvePoints
        
        [~,ThirtyPointIdx] = min(abs(upstroke - 0.30)); %get point at 30% max peak
        curvePoints(i).ThirtyPointIdx = ThirtyPointIdx; %add to curvePoints struct
        curvePoints(i).ThirtyPoint = flowTemp(ThirtyPointIdx);
        
        [~,FiftyPointIdx] = min(abs(upstroke - 0.50)); %get point at 50% max peak    
        curvePoints(i).FiftyPointIdx = FiftyPointIdx;
        curvePoints(i).FiftyPoint = flowTemp(FiftyPointIdx);

        flows(i,:) = normalize(flowTemp,'range'); %normalize curve from here on
        %flows(i,:) = flowTemp;
    end      
    
    %check if were missing the upslope (happens frequently with PG gating)
    peakAtStart = [curvePoints.maxPeakVelIdx]==1; 
    
    if sum(peakAtStart)
        for i=1:numROIs
            flow(i).TTPoint = NaN; %add to flow struct
            flow(i).TTUpstroke = NaN;
            flow(i).TTFoot = NaN;
            flow(i).Xcorr = NaN;
        end   
    else
        % TTPoint - time to point calculation
        for i=1:numROIs
            tp = curvePoints(i).times(curvePoints(i).FiftyPointIdx); %get time at 50% peak (ms) 
            flow(i).TTPoint = tp; %add to flow struct
        end    
        
        % TTFoot - time to foot calculation
        for i=1:numROIs
            leftIdx = curvePoints(i).ThirtyPointIdx;
            rightIdx = curvePoints(i).SeventyPointIdx;
            flowSeg = flows(i,leftIdx:rightIdx);
            timeSeg = curvePoints(i).times(leftIdx:rightIdx);
            offset = flows(i,curvePoints(1).minPeakVelIdx); %baseline offset
            p1 = polyfit(timeSeg,flowSeg,1);  %y=mx+b
            tp = (offset - p1(2))/p1(1);  %offset=mx+b -> x=(offset-b)/m
            flow(i).TTFoot = tp; %add to flow struct
            flow(i).P1 = p1;
        end  

        % TTCurvature - time to maximum curvature calculation
        % for i=1:numROIs
        %     leftIdx = curvePoints(i).minPeakVelIdx;
        %     rightIdx = curvePoints(i).maxPeakVelIdx;
        %     flowSeg = flows(i,leftIdx:rightIdx);
        %     timeSeg = curvePoints(i).times(leftIdx:rightIdx);
        %     dy = diff(flowSeg,1);  %dy/dt
        %     dy(end) = [];
        %     dx = diff(timeSeg,1);  %dx/dt = delta_t
        %     dx(end) = [];
        %     ddy = diff(flowSeg,2);  %ddy/dt (ddx/dt=0)
        %     curvature = ddy.*dx./(dx.^2 + dy.^2).^(3/2);
        %     [~,tIdx] = max(curvature);
        %     tp = timeSeg(tIdx-2);  %-2 accounts for dimension change after derivative
        %     flow(i).TTCurvature = tp; %max curvature of flow curve
        % end  
        
        % TTUpstroke - time to upstroke calculation (DOI: 10.1002/jmri.22570)
        for i=1:numROIs
            [sigmoid,sigTimes,tp,midpt] = sigFit(flows(i,:),curvePoints(i).times); %see sigFit function below
            flow(i).TTUpstroke = tp; %start of upstroke
            % flow(i).TTUpstroke = midpt; %midpoint of sigmoid
            flow(i).SigmoidFit = sigmoid;
            flow(i).SigmoidTimes = sigTimes;
        end 
        
        % XCorr - cross correlation calculation
        %we need 2 curves to get shift, set 1st time to arb. reference
        flow(1).Xcorr = curvePoints(1).ThirtyPointIdx;
        for i=2:numROIs
            flowSeg1 = flows(i-1,:);
            timeSeg1 = curvePoints(i-1).times;
            lastPoint1 = curvePoints(i-1).times(end); %get time endpoint on first curve
            flowSeg2 = flows(i,:); %chop second curve between 1 and right index
            timeSeg2 = curvePoints(i).times; %get times 
            lastPoint2 = curvePoints(i).times(end); %get time endpoint on second curve
            if lastPoint2>lastPoint1  %make sure times are aligned
                flowSeg2 = interp1(timeSeg2,flowSeg2,timeSeg1);
                timeres = curvePoints(i-1).timeres;
            else
                flowSeg1 = interp1(timeSeg1,flowSeg1,timeSeg2);
                timeres = curvePoints(i).timeres;
            end
            
            % Homemade cross-correlation (circshift, not FFT-based)
            % L = (1:length(flowSeg2))-(length(flowSeg2)/2);
            % for k=1:length(L)
            %     xc(i) = sum(flowSeg2 .* circshift(flowSeg1,L(k)));
            % end 
            % [~,idx] = max(xc);
            % shift = L(idx);
            
            %Built-in cross-correlation (fails if peaks are near edge)
            [Xcorrs,lags] = xcorr(flowSeg2,flowSeg1,'normalized'); %perform cross correlation between flow curves
            [~,maxLag] = max(Xcorrs);
            shift = lags(maxLag);
            flow(i).Xcorr = flow(i-1).Xcorr + timeres*shift; %cumulative sum from last time
        end  
    end 
    
    
% --- Plot Velocities    
function plotVelocity(handles)
    cla(handles.velocityPlot,'reset'); %reset axes
    axes(handles.velocityPlot); %make sure we plot on the right axis
    
    times = handles.pcDatasets(1).Interp(1).times;
    plot(times, zeros(1,length(times)) ,'Color','black','LineWidth',1.5); %line of y=0 (for visual reference)
    % Note that above we assume same time scale for each plane (MR scan)
    xlim([min(times) max(times)]);
    legendSet = {'Baseline'}; %add baseline to legend names
    for i=1:numel(handles.pcDatasets) %for all planes
        if isstruct(handles.pcDatasets(i).Interp) %if we've made ROI data for this dataset
            for j=1:length(handles.pcDatasets(i).Interp) %for each ROI
                legendSet{end+1} = handles.pcDatasets(i).Raw(j).name; %add name of ROI to list
                switch handles.global.interpType %check what interpolation we are using
                    case 'None'
                        velocity = handles.pcDatasets(i).Interp(j).mean_v; %grab mean velocity
                        if strcmp(handles.global.analysisType,'Absolute_Time')
                            times = handles.pcDatasets(i).Interp(j).times;
                        else
                            times = handles.pcDatasets(1).Interp(1).times;
                        end 
                        stdv = handles.pcDatasets(i).Interp(j).stdv_v; %grab stdv of velocity
                    case 'Gaussian'
                        velocity = handles.pcDatasets(i).Gaussian(j).mean_v;
                        if strcmp(handles.global.analysisType,'Absolute_Time')
                            times = handles.pcDatasets(i).Gaussian(j).times;
                        else
                            times = handles.pcDatasets(1).Gaussian(1).times;
                        end 
                        stdv = handles.pcDatasets(i).Gaussian(j).stdv_v;
                    case 'Shifted'
                        velocity = handles.pcDatasets(i).Shifted(j).mean_v;
                        if strcmp(handles.global.analysisType,'Absolute_Time')
                            times = handles.pcDatasets(i).Shifted(j).times;
                        else
                            times = handles.pcDatasets(1).Shifted(1).times;
                        end 
                        stdv = handles.pcDatasets(i).Shifted(j).stdv_v;
                end 

                if mean(velocity)<0 %if we are mainly negative velocities (as in descending aorta)
                    if get(handles.errorBarRadio,'Value') %and if we want to show error bars
                        hold on; errorbar(times,-1*velocity,stdv); %plot inverted velocity and errors
                    else
                        hold on; plot(times,-1*velocity); %else, plot inverted velocity
                    end 
                else %otherwise, don't invert velocity (as in ascending aorta)
                    if get(handles.errorBarRadio,'Value')
                        hold on; errorbar(times,velocity,stdv);
                    else
                        hold on; plot(times,velocity);
                    end
                end 
            end 
        end 
    end 
    legend(legendSet); hold off
    xlabel('Time (ms)'); ylabel('Mean Velocity in ROI (mm/s)'); %set axes labels
    %xlim([times(1),times(end)]); %chop limits to make curve full width
    

    
% --- Sigmoid Fit Function
function [sigmoid,t,t1,t2] = sigFit(mean_v,times)
%%% See the following article by Anas Dogui in JMRI:
% https://pubmed.ncbi.nlm.nih.gov/21591000/

    [~,MAX] = max(mean_v); %find max
    upslope = mean_v(1:MAX); %find upslope region of flow curve
    t0 = times(1:MAX);
    t = linspace(1,times(MAX),1000); %interpolate even more
    upslope = interp1(t0,upslope,t); %interpolate upslope
    upslope = rescale(upslope); %normalize from 0 to 1
    [~,MIN] = min(upslope);
    upslope = upslope(MIN:end);
    t = t(MIN:end);
    dt = t(2)-t(1); %new temporal resolution (=0.1)
    midpoint = round(length(upslope)/2);
    
    % c1 = b, c2 = a, c3 = x0, c4 = dx
    % Note that we could assume the equation e^t/(1+e^(t-t0)) since c1=1 and
    % c2=0. However, will keep the same as the Dogui paper.
    sigmoidModel = @(c) c(1) + ( (c(2)-c(1)) ) ./ ( 1+exp((t-c(3))./c(4)) ) - upslope;
    c0 = [0,1,t(midpoint),dt/2]; %initial params for upslope region
    opts = optimset('Display', 'off'); %turn off display output
    c = lsqnonlin(sigmoidModel,c0,[],[],opts); %get nonlinear LSQ solution
    sigmoid = c(1) + ( (c(2)-c(1)) ) ./ ( 1+exp((t-c(3))./c(4)) ); %calculate our sigmoid fit with our new params
    t2 = c(3); %midpoint of sigmoid
    
    dy = diff(sigmoid,1);
    dy(end) = [];
    dx = diff(t,1);
    dx(end) = [];
    ddy = diff(sigmoid,2);
    curvature = ddy.*dx./(dx.^2 + dy.^2).^(3/2);
    [~,tIdx] = max(curvature);
    t1 = t(tIdx); %point of max curvature of sigmoid

    
    
% --- Condense and organize flow data obtained from ROIs
function flow = organizeFlowInfo(handles)
% This function is designed to parse handles.pcDatasets. It is
% difficult to do analysis on the handles structure because some slices
% have 2 ROIs. It is much easier to have a structure that pulls out each
% ROI. It only adds a bit of memory since we aren't saving the raw images.
    count = 1; %overall iterator
    for i=1:numel(handles.pcDatasets) %for each PC dataset
        if isstruct(handles.pcDatasets(i).Interp) %do we have data?
            if length(handles.pcDatasets(i).Interp)==1 %if we just have one ROI
                flow(count).Name = handles.pcDatasets(i).Raw.name; %pull only relevant info for PWV calcs
                flow(count).Raw = handles.pcDatasets(i).Raw;
                flow(count).HeaderInfo = handles.pcDatasets(i).Info;
                flow(count).Interp = handles.pcDatasets(i).Interp;
                flow(count).Gaussian = handles.pcDatasets(i).Gaussian;
                flow(count).Shifted = handles.pcDatasets(i).Shifted;
                flow(count).pcDatasetREF = [i,1];
                count = count+1;
            else 
                for j=1:length(handles.pcDatasets(i).Interp) %if we have more than one ROI
                    flow(count).Name = handles.pcDatasets(i).Raw(j).name; %parse into individual ROIs in flow struct
                    flow(count).Raw = handles.pcDatasets(i).Raw(j);
                    flow(count).HeaderInfo = handles.pcDatasets(i).Info;
                    flow(count).Interp = handles.pcDatasets(i).Interp(j);
                    flow(count).Gaussian = handles.pcDatasets(i).Gaussian(j);
                    flow(count).Shifted = handles.pcDatasets(i).Shifted(j);
                    flow(count).pcDatasetREF = [i,j];
                    count = count+1;
                end 
            end 
        end 
    end 

% --- Save Time-To Plots for record-keeping and debugging
 function saveTTplots(handles,flow)   
    times = flow(1).Interp.times;
    figure('units','normalized','outerposition',[0 0 1 1]); 
    % Note that above we assume same time scale for each plane (MR scan)
    xlim([min(times) max(times)]);
    legendSet = {}; %add baseline to legend names
    count = 1;
    for i=1:numel(flow) %for all planes
        switch handles.global.interpType %check what interpolation we are using
            case 'None'
                vTemp = flow(i).Interp.mean_v; %grab mean velocity
                timeTemp = flow(i).Interp.times;
            case 'Gaussian'
                vTemp = flow(i).Gaussian.mean_v; %grab mean velocity
                timeTemp = flow(i).Gaussian.times;
            case 'Shifted'
                vTemp = flow(i).Shifted.mean_v; %grab mean velocity
                timeTemp = flow(i).Shifted.times;
        end 

        if mean(vTemp)<0
            vTemp = -1*vTemp;
        end
        vTemp = normalize(vTemp,'range');

        hold on; plot(timeTemp,vTemp);
        velocity(count,:) = vTemp;
        times(count,:) = timeTemp;
        count = count+1;
    end 
    wtime = 1.5;
    legend(legendSet,'Location','northeastoutside','AutoUpdate','off');
    xlabel('Time (ms)'); 
    ylim([0 1.05]);
    ylabel('Normalized Mean Velocity in ROI (mm/s)'); %set axes labels
    text(50,1.05,'Note: Waveforms may be distorted if HR different between scans!')

    ttp = [flow.TTPoint];
    for i=1:length(ttp)
        [~,idx] = min( abs(times(i,:)-ttp(i)) );
        h(i) = scatter(times(i,idx),velocity(i,idx),'k','filled');
        h2(i) = scatter(times(i,idx),0,'k');
        h3(i) = plot([times(i,idx) times(i,idx)],[velocity(i,idx),0],':k');
    end 
    title('Time-To-Point');
    frame = getframe(gcf);
    imwrite(frame2im(frame),[pwd filesep 'TTPoint.png']);
    pause(wtime)
    set(h,'Visible','off'); set(h2,'Visible','off'); set(h3,'Visible','off')
    
    ttf = [flow.TTFoot];
    for i=1:length(ttf)
        [~,idx] = min( abs(times(i,:)-ttf(i)) );
        [~,idx2] = max(velocity(i,:));
        P1 = flow(i).P1;
        x = times(i,idx:idx2);
        y = P1(1)*x + P1(2);
        h(i) = plot(x,y,'k');
        h2(i) = scatter(times(i,idx),0,'k');
    end 
    title('Time-To-Foot');
    frame = getframe(gcf);
    imwrite(frame2im(frame),[pwd filesep 'TTFoot.png']);
    pause(wtime)

    set(h,'Visible','off'); set(h2,'Visible','off');
        
    ttu = [flow.TTUpstroke];
    for i=1:length(ttu)
        [~,idx] = min( abs(times(i,:)-ttu(i)) );
        [~,maxIdx] = max(velocity(i,:));
        upstroke = velocity(i,1:maxIdx);
        [minVel,~] = min(upstroke);
        h(i) = scatter(times(i,idx),velocity(i,idx),'k','filled');
        h2(i) = scatter(times(i,idx),0,'k');
        h3(i) = plot([times(i,idx) times(i,idx)],[velocity(i,idx),0],':k');
        sigT = flow(i).SigmoidTimes;
        sig = flow(i).SigmoidFit;
        sigShift = sig - sig*minVel + minVel;
        h4(i) = plot(sigT,sigShift,'.k');
    end 
    title('Time-To-Upstroke');
    frame = getframe(gcf);
    imwrite(frame2im(frame),[pwd filesep 'TTUpstroke.png']);
    pause(wtime)
    set(h,'Visible','off'); set(h2,'Visible','off'); set(h3,'Visible','off'); set(h4,'Visible','off');
    
    xcorr = [flow.Xcorr];
    for i=1:length(xcorr)
        [~,idx] = min( abs(times(i,:)-xcorr(i)) );
        h(i) = scatter(times(i,idx),velocity(i,idx),'k','filled');
        h2(i) = scatter(times(i,idx),0,'k');
        h3(i) = plot([times(i,idx) times(i,idx)],[velocity(i,idx),0],':k');
    end 
    title('Cross Correlation Time Lag');
    frame = getframe(gcf);
    imwrite(frame2im(frame),[pwd filesep 'Xcorr.png']);
    pause(wtime)
    set(h,'Visible','off'); set(h2,'Visible','off'); set(h3,'Visible','off')
    
    hold off; close(gcf);   
    
    
    
% Load Dat files
function v = load_dat(name, res)
    [fid,errmsg]= fopen(name,'r');
    if fid < 0  %if name does not exist in directory
        disp(['Error Opening Data : ',errmsg]);
    end

    % Reads in as short, reshapes by image res.
    v = reshape(fread(fid,'short=>single'),res);
    fclose(fid);