function varargout = GTG_preprocess_GUI(varargin)

% Author: Jeffrey M. Spielberg (jspielb2@gmail.com)
% Version: Beta 0.43 (09.02.15)
% 
% 
% History:
% 02.27.14 - Beta 0.13 - initial public release
% 03.11.14 - Beta 0.20 - 1) small bugfixes, 2) major overhaul of user 
%                        interface into GUIs
% 03.17.14 - Beta 0.21 - small bugfixes
% 03.24.14 - Beta 0.22 - lots of small bugfixes
% 04.08.14 - Beta 0.23 - small bugfixes
% 04.23.14 - Beta 0.24 - no changes to this stage
% 05.07.14 - Beta 0.25 - no changes to this stage
% 06.10.14 - Beta 0.30 - 1) small bugfixes, 2) handles now used to pass 
%                        information between functions (rather than via 
%                        out, which was made global), allowing users to 
%                        launch processes from the same gui with less 
%                        chance of info from the previous process 
%                        interfering
% 06.25.14 - Beta 0.31 - 1) addition of Muschelli et al. (2014)'s procedure
%                        for partialing the first 5 principal components
%                        instead of the mean of white matter and
%                        ventricular signal (NOTE, if local white matter is
%                        selected, the mean of the local ROI will be
%                        computed, even if the PCA option is selected), 2) 
%                        addition of Patel et al. (2014)'s WaveletDespike 
%                        procedure
% 08.26.14 - Beta 0.32 - 1) made script compatible with both compressed &
%                        non-compressed nifti format, 2) fixed bug whereby
%                        wavelet despike was saving files in RAS format,
%                        but the script expected LAS; now the LR
%                        orientation is checked & changed to LAS if need be
% 09.22.14 - Beta 0.33 - no changes to this stage
% 10.06.14 - Beta 0.35 - 1) ROIlabel cell array must now contain 2 columns:
%                        the first contains the ROI labels and the second
%                        the numeric identifiers that correspond to the 3d
%                        ROI images. This change allows for the case that
%                        some participants may not have all ROIs, in which
%                        case NaNs will be used for the timeseries for this
%                        ROI, 2) added two new possible methods for
%                        extracting timeseries from each ROI (previously
%                        only the mean was available): median, largest
%                        principal component, 3) added the capability to
%                        erode white matter and ventricle masks
% 11.19.14 - Beta 0.36 - minor bugfixes
% 12.17.14 - Beta 0.37 - 1) minor bugfixes, 2) addition of option to
%                        extract timeseries from ROIs before all
%                        preprocessing is complete, then perform remaining
%                        preprocessing on mean extracted signal (to save
%                        time)
% 03.24.15 - Beta 0.38 - minor bugfixes
% 05.01.15 - Beta 0.39 - 1) fixed bug whereby the threshold for whether to
%                        assign NaNs during ROI extraction was >90% of the
%                        ROI was removed during masking instead of >50%
%                        (but logfile still said 50%), 2) added option to
%                        specify the minimum # of voxels required to be in
%                        an ROI (if >= that #, the timeseries is replaced
%                        with NaNs), previously this had been hard coded as
%                        5 (the current default value), 3) added option to
%                        specify the minimum % of voxels retained after
%                        masking by the functional mask required for an ROI
%                        (if >= that %, the timeseries is replaced with 
%                        NaNs), previously this had been hard coded as 50%
%                        (the current default value), 4) added option to
%                        use existing nii files (with correct name) instead
%                        of recomputing every file every time, 5) fixed a 
%                        bug whereby the mean was used to extract white 
%                        matter/ventricular signal after the second round 
%                        of processing when scrubbing (i.e., after removing
%                        bad time points) even if PCA was selected (PCA was
%                        used during the first round), 6) updated to use 
%                        the preferred version of matlab's pca 
%                        instantiation, 7) allowed PCA extraction of
%                        ventricular signal when using local white matter
%                        processing, 8) added safeguard whereby matlab will
%                        test whether text files (i.e., containing white
%                        matter, ventricular, or global signal or DVARS)
%                        already exist and, if so, add a number to the
%                        filename of newly created files in order to
%                        distinguish them and avoid confusion, 9) changed
%                        assignment of sign/direction of signal when using
%                        PCA to extract each ROI's timeseries; previously
%                        the sign/direction was determined by PCA alone;
%                        after conducting simulations, it was determined
%                        that PCA accurately reflected the sign/direction 
%                        of underlying signal 90-95% of the time (SVD was
%                        accurate only ~50% of the time), whereas the mean 
%                        always accurately captured the sign (PCA more
%                        accurately captured the underlying signal
%                        variance); therefore, the script now checks the 
%                        correlation between the mean and largest principal
%                        component, and if it is -1, the PC is multiplied 
%                        by -1; of note, this is the same procedure used by
%                        AFNI when using PCA to summarize ROI signal
% 06.30.15 - Beta 0.40 - minor bugfixes
% 07.20.15 - Beta 0.41 - added automatic naming if the user does not
%                        provide an output name
% 08.24.15 - Beta 0.42 - no changes
% 09.02.15 - Beta 0.43 - 1) bugfix with butterworth filter not filtering
%                        approprietly, 2) added calpability to do highpass
%                        or lowpass rather than bandpass, 3) changed to a
%                        matlab implementation to calculate DVARS, which
%                        should reduce computation time
%
%
% WARNING: This is a beta version. There no known bugs, but only limited 
% testing has been perfomed. This software comes with no warranty (even the
% implied warranty of merchantability or fitness for a particular purpose).
% Therefore, USE AT YOUR OWN RISK!!!
%
% Copyleft 2014-2015. Software can be modified and redistributed, but 
% modifed, redistributed versions must have the same rights




% Begin initialization code
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GTG_preprocess_GUI_OpeningFcn, ...
                   'gui_OutputFcn',  @GTG_preprocess_GUI_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

function GTG_preprocess_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
guidata(hObject, handles);

function varargout = GTG_preprocess_GUI_OutputFcn(hObject, eventdata, handles) 
varargout{1} = handles.output;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function IDs_edit_Callback(hObject, eventdata, handles) %#ok<*DEFNU,*INUSD>
set(handles.Start_pushbutton,'enable','on');
temp             = get(hObject,'String');
handles.out.subs = evalin('base',temp);
if size(handles.out.subs,1) < size(handles.out.subs,2)
    handles.out.subs = handles.out.subs';
end
if ~iscell(handles.out.subs)
    handles.out.subs = strtrim(cellstr(num2str(handles.out.subs)));
end
guidata(hObject,handles);

function IDs_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function IDs_edit_ButtonDownFcn(hObject, eventdata, handles)
set(hObject, 'Enable', 'On');
uicontrol(handles.IDs_edit);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Useavail_popupmenu_Callback(hObject, eventdata, handles)
contents = cellstr(get(hObject,'String'));
handles.out.useavail = contents{get(hObject,'Value')};
guidata(hObject,handles);

function Useavail_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ROIlab_edit_Callback(hObject, eventdata, handles)
temp1                      = get(hObject,'String');
temp2                      = evalin('base',temp1);
handles.out.ROI_labels     = temp2(:,1);
handles.out.ROI_num_labels = cell2mat(temp2(:,2));
handles.out.nROI           = length(handles.out.ROI_labels);
guidata(hObject,handles);

function ROIlab_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function ROIlab_edit_ButtonDownFcn(hObject, eventdata, handles)




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ROIextord_popupmenu_Callback(hObject, eventdata, handles)
contents = cellstr(get(hObject,'String'));
handles.out.ROIextord = contents{get(hObject,'Value')};
if strcmp(handles.out.ROIextord,'After only slice timing/motion correction/detrending')
    set(handles.Despike_checkbox,'enable','off');
    set(handles.Motcens_checkbox,'enable','off');
    set(handles.Motcens_FD_edit,'enable','off');
    set(handles.Motcens_DVARS_edit,'enable','off');
    set(handles.Globpart_GCOR_checkbox,'enable','off');
elseif strcmp(handles.out.ROIextord,'After all preprocessing')
    set(handles.Despike_checkbox,'enable','on');
    set(handles.Motcens_checkbox,'enable','on');
    set(handles.Motcens_FD_edit,'enable','on');
    set(handles.Motcens_DVARS_edit,'enable','on');
    set(handles.Globpart_GCOR_checkbox,'enable','on');
end
guidata(hObject,handles);

function ROIextord_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Outfile_edit_Callback(hObject, eventdata, handles)
handles.out.outname = get(hObject,'String');
guidata(hObject,handles);

function Outfile_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Outfile_pushbutton_Callback(hObject, eventdata, handles) %#ok<*INUSL>
handles.out.outname = [uigetdir(pwd,'Select the output directory'),'/'];
set(handles.Outfile_edit,'String',handles.out.outname);
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Saveintermed_popupmenu_Callback(hObject, eventdata, handles)
contents = cellstr(get(hObject,'String'));
handles.out.saveintermed = contents{get(hObject,'Value')};
guidata(hObject,handles);

function Saveintermed_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function nTP_edit_Callback(hObject, eventdata, handles)
handles.out.nTP = str2double(get(hObject,'String'));
guidata(hObject,handles);

function nTP_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function TR_edit_Callback(hObject, eventdata, handles)
%set(handles.Start_pushbutton,'enable','on');
handles.out.TR = str2double(get(hObject,'String'));
guidata(hObject,handles);

function TR_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MinROIvox_edit_Callback(hObject, eventdata, handles)
handles.out.MinROIvox = str2double(get(hObject,'String'));
guidata(hObject,handles);

function MinROIvox_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MinpercentROIvox_edit_Callback(hObject, eventdata, handles)
handles.out.MinpercentROIvox = str2double(get(hObject,'String'));
guidata(hObject,handles);

function MinpercentROIvox_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Funcfile_edit_Callback(hObject, eventdata, handles)
handles.out.func_first_filename = get(hObject,'String');
guidata(hObject,handles);

function Funcfile_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Funcfile_pushbutton_Callback(hObject, eventdata, handles)
[file,path]             = uigetfile({'*.nii.gz';'*.nii'},'Select the 4d EPI for the first participant');
handles.out.func_first_filename = [path,file];
set(handles.Funcfile_edit,'String',handles.out.func_first_filename);
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Funcmaskfile_pushbutton_Callback(hObject, eventdata, handles)
[file,path]                 = uigetfile({'*.nii.gz';'*.nii'},'Select the 3d functional brain mask (used to limit the voxels included) for the first participant');
handles.out.first_funcmask_filename = [path,file];
set(handles.Funcmaskfile_edit,'String',handles.out.first_funcmask_filename);
guidata(hObject,handles);

function Funcmaskfile_edit_Callback(hObject, eventdata, handles)
handles.out.first_funcmask_filename = get(hObject,'String');
guidata(hObject,handles);

function Funcmaskfile_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ROImaskfile_edit_Callback(hObject, eventdata, handles)
handles.out.first_parc_filename = get(hObject,'String');
guidata(hObject,handles);

function ROImaskfile_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function ROImaskfile_pushbutton_Callback(hObject, eventdata, handles)
[file,path]             = uigetfile({'*.nii.gz';'*.nii'},'Select the 3d image containing ROI masks for the first participant');
handles.out.first_parc_filename = [path,file];
set(handles.ROImaskfile_edit,'String',handles.out.first_parc_filename);
guidata(hObject,handles);

function ROImaskspace_popupmenu_Callback(hObject, eventdata, handles)
contents       = cellstr(get(hObject,'String'));
handles.out.parc_space = contents{get(hObject,'Value')};
guidata(hObject,handles);

function ROImaskspace_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function WMmaskfile_edit_Callback(hObject, eventdata, handles)
handles.out.first_WM_filename = get(hObject,'String');
guidata(hObject,handles);

function WMmaskfile_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function WMmaskfile_pushbutton_Callback(hObject, eventdata, handles)
[file,path]           = uigetfile({'*.nii.gz';'*.nii'},'Select the white matter mask for the first participant');
handles.out.first_WM_filename = [path,file];
set(handles.WMmaskfile_edit,'String',handles.out.first_WM_filename);
guidata(hObject,handles);

function WMmaskspace_popupmenu_Callback(hObject, eventdata, handles)
contents             = cellstr(get(hObject,'String'));
handles.out.WM_space = contents{get(hObject,'Value')};
guidata(hObject,handles);

function WMmaskspace_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Ventmaskfile_edit_Callback(hObject, eventdata, handles)
handles.out.first_vent_filename = get(hObject,'String');
guidata(hObject,handles);

function Ventmaskfile_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Ventmaskfile_pushbutton_Callback(hObject, eventdata, handles)
[file,path]             = uigetfile({'*.nii.gz';'*.nii'},'Select the ventricular mask for the first participant');
handles.out.first_vent_filename = [path,file];
set(handles.Ventmaskfile_edit,'String',handles.out.first_vent_filename);
guidata(hObject,handles);

function Ventmaskspace_popupmenu_Callback(hObject, eventdata, handles)
contents       = cellstr(get(hObject,'String'));
handles.out.vent_space = contents{get(hObject,'Value')};
guidata(hObject,handles);

function Ventmaskspace_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ST_checkbox_Callback(hObject, eventdata, handles)
handles.out.ST_correct = get(hObject,'Value');

if handles.out.ST_correct     == 1
    set(handles.STord_popupmenu,'enable','on');
elseif handles.out.ST_correct == 0
    set(handles.STord_popupmenu,'enable','off');
end
guidata(hObject,handles);

function STord_popupmenu_Callback(hObject, eventdata, handles)
contents   = cellstr(get(hObject,'String'));
handles.out.ST_ord = contents{get(hObject,'Value')};
guidata(hObject,handles);

function STord_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Detr_checkbox_Callback(hObject, eventdata, handles)
handles.out.detr = get(hObject,'Value');

if handles.out.detr     == 1
    set(handles.Detrord_popupmenu,'enable','on');
elseif handles.out.detr == 0
    set(handles.Detrord_popupmenu,'enable','off');
end
guidata(hObject,handles);

function Detrord_popupmenu_Callback(hObject, eventdata, handles)
contents     = cellstr(get(hObject,'String'));
handles.out.detr_ord = str2double(contents{get(hObject,'Value')});
guidata(hObject,handles);

function Detrord_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Despike_checkbox_Callback(hObject, eventdata, handles)
handles.out.despike = get(hObject,'Value');

if handles.out.despike     == 1
    set(handles.Motcens_checkbox,'enable','off');
    set(handles.Motcens_checkbox,'Value',0);
    handles.out.motcens = 0;
    set(handles.Motcens_FD_edit,'enable','off');
    set(handles.Motcens_DVARS_edit,'enable','off');
elseif handles.out.despike == 0
    set(handles.Motcens_checkbox,'enable','on');
end
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function BP_checkbox_Callback(hObject, eventdata, handles)
handles.out.BP = get(hObject,'Value');

if handles.out.BP     == 1
    set(handles.BPtype_popupmenu,'enable','on');
    set(handles.BP_HP_edit,'enable','on');
    set(handles.BP_LP_edit,'enable','on');
elseif handles.out.BP == 0
    set(handles.BPtype_popupmenu,'enable','off');
    set(handles.BP_HP_edit,'enable','off');
    set(handles.BP_LP_edit,'enable','off');
end
guidata(hObject,handles);

function BPtype_popupmenu_Callback(hObject, eventdata, handles)
contents    = cellstr(get(hObject,'String'));
handles.out.BP_type = contents{get(hObject,'Value')};

if strcmp(handles.out.BP_type,'Butterworth')
    set(handles.BP_buttord_popupmenu,'enable','on');
elseif ~strcmp(handles.out.BP_type,'Butterworth')
    set(handles.BP_buttord_popupmenu,'enable','off');
end
guidata(hObject,handles);

function BPtype_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function BP_HP_edit_Callback(hObject, eventdata, handles)
handles.out.BP_HP_cutoff = str2double(get(hObject,'String'));
guidata(hObject,handles);

function BP_HP_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function BP_LP_edit_Callback(hObject, eventdata, handles)
handles.out.BP_LP_cutoff = str2double(get(hObject,'String'));
guidata(hObject,handles);

function BP_LP_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function BP_buttord_popupmenu_Callback(hObject, eventdata, handles)
contents    = cellstr(get(hObject,'String'));
handles.out.butter_ord = str2double(contents{get(hObject,'Value')});
guidata(hObject,handles);

function BP_buttord_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function MC_checkbox_Callback(hObject, eventdata, handles)
handles.out.MC = get(hObject,'Value');
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Motcens_checkbox_Callback(hObject, eventdata, handles)
handles.out.motcens = get(hObject,'Value');

if handles.out.motcens     == 1
    set(handles.Motcens_FD_edit,'enable','on');
    set(handles.Motcens_DVARS_edit,'enable','on');
    set(handles.Despike_checkbox,'enable','off');
    set(handles.Despike_checkbox,'Value',0);
    handles.out.despike = 0;
elseif handles.out.motcens == 0
    set(handles.Motcens_FD_edit,'enable','off');
    set(handles.Motcens_DVARS_edit,'enable','off');
    set(handles.Despike_checkbox,'enable','on');
end
guidata(hObject,handles);

function Motcens_FD_edit_Callback(hObject, eventdata, handles)
handles.out.motcens_FD_cutoff = str2double(get(hObject,'String'));
guidata(hObject,handles);

function Motcens_FD_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Motcens_DVARS_edit_Callback(hObject, eventdata, handles)
handles.out.motcens_DVARS_cutoff = str2double(get(hObject,'String'));
guidata(hObject,handles);

function Motcens_DVARS_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Motparpart_checkbox_Callback(hObject, eventdata, handles)
handles.out.motpar_part = get(hObject,'Value');

if handles.out.motpar_part     == 1
    set(handles.Motparpart_t1_checkbox,'enable','on');
    set(handles.Motparpart_sqr_checkbox,'enable','on');
    set(handles.Motparpart_deriv_checkbox,'enable','on');
elseif handles.out.motpar_part == 0
    set(handles.Motparpart_t1_checkbox,'enable','off');
    set(handles.Motparpart_sqr_checkbox,'enable','off');
    set(handles.Motparpart_deriv_checkbox,'enable','off');
end
guidata(hObject,handles);

function Motparpart_t1_checkbox_Callback(hObject, eventdata, handles)
handles.out.motpart1_part    = get(hObject,'Value');
guidata(hObject,handles);

function Motparpart_sqr_checkbox_Callback(hObject, eventdata, handles)
handles.out.motparsqr_part   = get(hObject,'Value');
guidata(hObject,handles);

function Motparpart_deriv_checkbox_Callback(hObject, eventdata, handles)
handles.out.motparderiv_part = get(hObject,'Value');
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Globpart_checkbox_Callback(hObject, eventdata, handles)
handles.out.globsig_part = get(hObject,'Value');

if handles.out.globsig_part     == 1
    set(handles.Globpart_deriv_checkbox,'enable','on');
    set(handles.Globpart_GNI_checkbox,'enable','on');
    set(handles.Globpart_GCOR_checkbox,'enable','on');
elseif handles.out.globsig_part == 0
    set(handles.Globpart_deriv_checkbox,'enable','off');
    set(handles.Globpart_GNI_checkbox,'enable','off');
    set(handles.Globpart_GCOR_checkbox,'enable','off');
end
guidata(hObject,handles);

function Globpart_deriv_checkbox_Callback(hObject, eventdata, handles)
handles.out.globsigderiv_part = get(hObject,'Value');
guidata(hObject,handles);

function Globpart_GNI_checkbox_Callback(hObject, eventdata, handles)
handles.out.globsig_calcGNI = get(hObject,'Value');
guidata(hObject,handles);

function Globpart_GCOR_checkbox_Callback(hObject, eventdata, handles)
handles.out.globsig_calcGCOR = get(hObject,'Value');
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function WMpart_checkbox_Callback(hObject, eventdata, handles)
handles.out.WMsig_part = get(hObject,'Value');

if handles.out.WMsig_part == 1
    set(handles.WMpart_scope_popupmenu,'enable','on');
    set(handles.WMpart_deriv_checkbox,'enable','on');
    set(handles.WMmaskfile_edit,'enable','on');
    set(handles.WMmaskfile_pushbutton,'enable','on');
    set(handles.WMmaskspace_popupmenu,'enable','on');
    set(handles.PCA_checkbox,'enable','on');
    set(handles.Erodemasks_checkbox,'enable','on');
elseif handles.out.WMsig_part == 0
    set(handles.WMpart_scope_popupmenu,'enable','off');
    set(handles.WMpart_deriv_checkbox,'enable','off');
    set(handles.WMmaskfile_edit,'enable','off');
    set(handles.WMmaskfile_pushbutton,'enable','off');
    set(handles.WMmaskspace_popupmenu,'enable','off');
    set(handles.PCA_checkbox,'enable','off');
    set(handles.Erodemasks_checkbox,'enable','off');
end
guidata(hObject,handles);

function WMpart_scope_popupmenu_Callback(hObject, eventdata, handles)
contents = cellstr(get(hObject,'String'));
handles.out.WMmask_scope = contents{get(hObject,'Value')};
guidata(hObject,handles);

function WMpart_scope_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function WMpart_deriv_checkbox_Callback(hObject, eventdata, handles)
handles.out.WMsigderiv_part = get(hObject,'Value');
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Ventpart_checkbox_Callback(hObject, eventdata, handles)
handles.out.ventsig_part = get(hObject,'Value');

if handles.out.ventsig_part == 1
    set(handles.Ventpart_deriv_checkbox,'enable','on');
    set(handles.Ventmaskfile_edit,'enable','on');
    set(handles.Ventmaskfile_pushbutton,'enable','on');
    set(handles.Ventmaskspace_popupmenu,'enable','on');
    set(handles.PCA_checkbox,'enable','on');
    set(handles.Erodemasks_checkbox,'enable','on');
elseif handles.out.ventsig_part == 0
    set(handles.Ventpart_deriv_checkbox,'enable','off');
    set(handles.Ventmaskfile_edit,'enable','off');
    set(handles.Ventmaskfile_pushbutton,'enable','off');
    set(handles.Ventmaskspace_popupmenu,'enable','off');
    set(handles.PCA_checkbox,'enable','off');
    set(handles.Erodemasks_checkbox,'enable','off');
end
guidata(hObject,handles);

function Ventpart_deriv_checkbox_Callback(hObject, eventdata, handles)
handles.out.ventsigderiv_part = get(hObject,'Value');
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function PCA_checkbox_Callback(hObject, eventdata, handles)
handles.out.PCA_WMvent = get(hObject,'Value');
guidata(hObject,handles);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Erodemasks_checkbox_Callback(hObject, eventdata, handles)
handles.out.erode_masks = get(hObject,'Value');
if handles.out.erode_masks == 1
    set(handles.Eroderadius_edit,'enable','on');
elseif handles.out.erode_masks == 0
    set(handles.Eroderadius_edit,'enable','off');
end
guidata(hObject,handles);

function Eroderadius_edit_Callback(hObject, eventdata, handles)
handles.out.erode_radius = get(hObject,'String');
guidata(hObject,handles);

function Eroderadius_edit_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ROIsummeas_popupmenu_Callback(hObject, eventdata, handles)
contents                        = cellstr(get(hObject,'String'));
handles.out.ROI_summary_measure = contents{get(hObject,'Value')};
guidata(hObject,handles);

function ROIsummeas_popupmenu_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Start_pushbutton_Callback(hObject, eventdata, handles)
out = handles.out;

% Check whether inputs have been specified
if ~isfield(out,'subs')
    msgbox('Enter a cell array of participant IDs','Error','error')
    return
end

if any(~cell2mat(cellfun(@ischar,out.subs,'UniformOutput',false)))
    out.subs = cellfun(@num2str,out.subs,'UniformOutput',false);
end

if ~isfield(out,'useavail')
    out.saveintermed = 'No, run all processing steps';
end
if ~isfield(out,'ROI_labels')
    msgbox('Enter a cell array of ROI labels','Error','error')
    return
end
if ~isfield(out,'ROIextord')
    out.ROIextord = 'After all preprocessing';
end
if ~isfield(out,'saveintermed')
    out.saveintermed = 'No';
end
if ~isfield(out,'nTP')
    out.nTP = 120;
end
if ~isfield(out,'TR')
    out.TR = 2;
end
if ~isfield(out,'MinROIvox')
    out.MinROIvox = 5;
end
if ~isfield(out,'MinpercentROIvox')
    out.MinpercentROIvox = 50;
end
if ~isfield(out,'func_first_filename')
    msgbox('Enter the 4d EPI filename for the first participant','Error','error')
    return
end
if ~isfield(out,'first_funcmask_filename')
    msgbox('Enter the 3d brainmask filename for the first participant','Error','error')
    return
end
if ~isfield(out,'first_parc_filename')
    msgbox('Enter the filename for the 3d image containing ROI masks for the first participant','Error','error')
    return
end
if ~isfield(out,'parc_space')
    out.parc_space = 'Functional';
end
if ~isfield(out,'ST_correct')
    out.ST_correct = 0;
end
if out.ST_correct == 1 && ~isfield(out,'ST_ord')
    msgbox('Please select the slice aquisition order','Error','error')
    return
end
if ~isfield(out,'detr')
    out.detr = 0;
end
if out.detr == 1 && ~isfield(out,'detr_ord')
    out.detr_ord = 0;
elseif ~isfield(out,'detr_ord')
    out.detr_ord = 2;
end
if ~isfield(out,'BP')
    out.BP = 0;
end
if out.BP == 1 && ~isfield(out,'BP_type')
    msgbox('Please select the bandpass filter type','Error','error')
    return
end
if ~isfield(out,'despike')
    out.despike = 0;
end
if ~isfield(out,'BP_LP_cutoff')
    out.BP_LP_cutoff = 0.10;
elseif isempty(out.BP_LP_cutoff)
    out.BP_LP_cutoff = -1;
elseif out.BP_LP_cutoff<=0
    out.BP_LP_cutoff = -1;
end
if ~isfield(out,'BP_HP_cutoff')
    out.BP_HP_cutoff = 0.01;
elseif isempty(out.BP_HP_cutoff)
    out.BP_HP_cutoff = -1;
elseif out.BP_HP_cutoff<=0
    out.BP_HP_cutoff = -1;
end
if out.BP_HP_cutoff==-1 && out.BP_LP_cutoff==-1
    msgbox('Temporal filtering is selected but both lowpass and highpass values are invalid','Error','error')
    return
end
if isfield(out,'BP_type')
    if strcmp(out.BP_type,'Butterworth') && ~isfield(out,'butter_ord')
        out.butter_ord = 1;
    end
end
if ~isfield(out,'MC')
    out.MC = 0;
end
if ~isfield(out,'motcens')
    out.motcens = 0;
end
if ~isfield(out,'motcens_FD_cutoff')
    out.motcens_FD_cutoff = 0.3;
end
if ~isfield(out,'motcens_DVARS_cutoff')
    out.motcens_DVARS_cutoff = 2.5;
end
if ~isfield(out,'motpar_part')
    out.motpar_part = 0;
end
if ~isfield(out,'motpart1_part')
    out.motpart1_part = 0;
end
if ~isfield(out,'motparsqr_part')
    out.motparsqr_part = 0;
end
if ~isfield(out,'motparderiv_part')
    out.motparderiv_part = 0;
end
if ~isfield(out,'globsig_part')
    out.globsig_part = 0;
end
if ~isfield(out,'globsigderiv_part')
    out.globsigderiv_part = 0;
end
if ~isfield(out,'globsig_calcGCOR')
    out.globsig_calcGCOR = 0;
end
if ~isfield(out,'globsig_calcGNI')
    out.globsig_calcGNI = 0;
end
if ~isfield(out,'WMsig_part')
    out.WMsig_part = 0;
end
if ~isfield(out,'WMsigderiv_part')
    out.WMsigderiv_part = 0;
end
if out.WMsig_part == 1 && ~isfield(out,'first_WM_filename')
    msgbox('Enter the filename for the white matter mask for the first participant','Error','error')
    return
end
if ~isfield(out,'WM_space')
    out.WM_space = 'Functional';
end
if ~isfield(out,'WMmask_scope')
    out.WMmask_scope = 'Entire Mask';
end
if ~isfield(out,'ventsig_part')
    out.ventsig_part = 0;
end
if ~isfield(out,'ventsigderiv_part')
    out.ventsigderiv_part = 0;
end
if out.ventsig_part == 1 && ~isfield(out,'first_vent_filename')
    msgbox('Enter the filename for the ventricular mask for the first participant','Error','error')
    return
end
if ~isfield(out,'vent_space')
    out.vent_space = 'Functional';
end
if ~isfield(out,'PCA_WMvent')
    out.PCA_WMvent = 0;
end
if ~isfield(out,'erode_masks')
    out.erode_masks = 0;
end
if out.erode_masks == 1 && ~isfield(out,'erode_radius')
    out.erode_radius = 2;
end
if isfield(out,'erode_radius')
    out.erode_radius = round(out.erode_radius);
    if out.erode_radius < 2
        out.erode_radius = 2;
    end
end
if ~isfield(out,'ROI_summary_measure')
    out.ROI_summary_measure = 'Mean';
end
if strcmp(out.ROIextord,'After only slice timing/motion correction/detrending')
    if strcmp(out.BP_type,'FSL')
        out.BP_type = 'Ideal';
    end
    out.despike = 0;
    out.motcens = 0;
    out.globsig_calcGCOR = 0;
    out.WMmask_scope = 'Entire Mask';
end

set(handles.Start_pushbutton,'enable','off');

if isempty(strfind(out.outname,'/'))
    out.outname = [pwd,'/',out.outname];
elseif out.outname(end)=='/'
    out.outname = [out.outname,'out'];
end

if any(strfind(out.func_first_filename,'.nii.gz'))
    img_ext = '.nii.gz';
else
    img_ext = '.nii';
end

if out.MC == 0 && (out.motpar_part == 1 || out.motcens == 1)
    par_base_filename = strrep(out.func_first_filename,img_ext,'.par');
    if ~exist(par_base_filename,'file')
        [file,path]       = uigetfile('*.par','Select the .par file (containing motion correction parameters) for the first participant');
        par_base_filename = [path,file];
    end
    par_base_filename = strrep(par_base_filename,out.subs{1},'SUBNUM');
else
    par_base_filename = ' ';
end

% Set base filenames
func_base_filename     = strrep(out.func_first_filename,out.subs{1},'SUBNUM');
funcmask_base_filename = strrep(out.first_funcmask_filename,out.subs{1},'SUBNUM');
parc_base_filename     = strrep(out.first_parc_filename,out.subs{1},'SUBNUM');
if isfield(out,'first_WM_filename')
    WM_base_filename   = strrep(out.first_WM_filename,out.subs{1},'SUBNUM');
else
    WM_base_filename   = ' ';
end
if isfield(out,'first_vent_filename')
    vent_base_filename = strrep(out.first_vent_filename,out.subs{1},'SUBNUM');
else
    vent_base_filename = ' ';    
end
if strcmp(out.parc_space,'Structural') || strcmp(out.WM_space,'Structural') || strcmp(out.vent_space,'Structural')
    [file,path]                          = uigetfile('*.mat','Select the structural-to-functional transformation matrix for the first participant');
    filename                             = [path,file];
    struct2func_regmat_base_filename = strrep(filename,out.subs{1},'SUBNUM');
else
    struct2func_regmat_base_filename = ' ';
end
if strcmp(out.parc_space,'Standard') || strcmp(out.WM_space,'Standard') || strcmp(out.vent_space,'Standard')
    [path,file]                            = uigetfile({'*.nii.gz';'*.nii'},'Select the standard-to-functional warp file for the first participant');
    filename                               = [path,file];
    standard2func_regmat_base_filename = strrep(filename,out.subs{1},'SUBNUM');
else
    standard2func_regmat_base_filename = ' ';
end

toolboxes = ver;
use_parfor = any(strcmpi({toolboxes.Name},'Parallel Computing Toolbox'));
if use_parfor
    if isempty(gcp('nocreate'))
        num_par_workers = str2double(inputdlg(sprintf('The Parallel Computing Toolbox was found on your system.\n\nEnter the number of workers you want to use (enter 1 to not use the PCT).\n\nNote: this must be <= the number of cores'),'PCT Workers',2));
        if num_par_workers > 12
            num_par_workers = 12;
        end
        if num_par_workers > feature('numCores')
            num_par_workers = feature('numCores');
        end
        if num_par_workers > 1
            try
                parpool('open',num_par_workers);
            catch
                matlabpool('open',num_par_workers);
            end
        else
            use_parfor = false;
        end
    end
end

out.preproc_outname = [out.outname '_preproc'];

% Write non-participant specific info to logfile
preproc_logfile_outname = [out.outname '_preproc_logfile.txt'];
logfile_fid = fopen(preproc_logfile_outname,'w');
fprintf(logfile_fid,'Output filename = %s\n',out.preproc_outname);
fprintf(logfile_fid,'# of ROIs = %u\n',out.nROI);
fprintf(logfile_fid,'# of timepoints = %u\n',out.nTP);
fprintf(logfile_fid,'(User specified) TR = %u\n',out.TR);
fprintf(logfile_fid,'Minimum # of voxels required in masked ROI = %u\n',out.MinROIvox);
fprintf(logfile_fid,'Minimum percentage of ROI required to be retained after masking = %u\n',out.MinpercentROIvox);
if strcmp(out.ROIextord,'After all preprocessing')
    fprintf(logfile_fid,'The ROI timeseries were extracted after all preprocessing had occurred\n');
elseif strcmp(out.ROIextord,'After only slice timing/motion correction/detrending')
    fprintf(logfile_fid,'The ROI timeseries were extracted after only slice timing\motion correction/detrending (if these steps were performed), after which all preprocessing was performed on the extracted timeseries\n');
end
if strcmp(out.saveintermed,'Yes, do not run step if output file exists')
    fprintf(logfile_fid,'Existing files were used if available for a participant\n');
end
if out.ST_correct == 1
    fprintf(logfile_fid,'Slicetiming correction was performed with %s (user specified) slice order\n',out.ST_ord);
end
if out.MC == 1
    fprintf(logfile_fid,'Motion correction was performed\n');
end
if out.detr == 1
    fprintf(logfile_fid,'The data were detrended for up to polynomials of order %u\n',out.detr_ord);
end
if out.despike == 1
    fprintf(logfile_fid,'The data were wavelet despiked\n');
end
if out.BP == 1 && strcmp(out.BP_type,'Ideal')
    if out.BP_LP_cutoff==-1
        fprintf(logfile_fid,'The data were highpass filtered with the "Ideal" filter (cutoff = %6.4f)\n',out.BP_HP_cutoff);
    elseif out.BP_HP_cutoff==-1
        fprintf(logfile_fid,'The data were lowpass filtered with the "Ideal" filter (cutoff = %6.4f)\n',out.BP_LP_cutoff);
    else
        fprintf(logfile_fid,'The data were bandpass filtered with the "Ideal" filter (lowpass cutoff = %6.4f, highpass cutoff = %6.4f)\n',out.BP_LP_cutoff,out.BP_HP_cutoff);
    end
elseif out.BP == 1 && strcmp(out.BP_type,'Butterworth')
    if out.BP_LP_cutoff==-1
        fprintf(logfile_fid,'The data were highpass filtered with a Butterworth filter of order %u (cutoff = %6.4f)\n',out.butter_ord,out.BP_HP_cutoff);
    elseif out.BP_HP_cutoff==-1
        fprintf(logfile_fid,'The data were lowpass filtered with a Butterworth filter of order %u (cutoff = %6.4f)\n',out.butter_ord,out.BP_LP_cutoff);
    else
        fprintf(logfile_fid,'The data were bandpass filtered with a Butterworth filter of order %u (lowpass cutoff = %6.4f, highpass cutoff = %6.4f)\n',out.butter_ord,out.BP_HP_cutoff,out.BP_LP_cutoff);
    end
elseif out.BP == 1 && strcmp(out.BP_type,'FSL')
    if out.BP_LP_cutoff==-1
        fprintf(logfile_fid,'The data were filtered with FSL''s non-linear highpass filter (cutoff = %6.4f)\n',out.BP_HP_cutoff);
    elseif out.BP_HP_cutoff==-1
        fprintf(logfile_fid,'The data were filtered with FSL''s Gaussian linear lowpass filter (cutoff = %6.4f)\n',out.BP_LP_cutoff);
    else
        fprintf(logfile_fid,'The data were filtered with FSL''s non-linear highpass and Gaussian linear lowpass filter (lowpass cutoff = %6.4f, highpass cutoff = %6.4f)\n',out.BP_LP_cutoff,out.BP_HP_cutoff);
    end
end
if out.motcens == 1
    fprintf(logfile_fid,'Motion-censoring was performed with FD cutoff = %5.3fmm and DVARS cutoff = %6.3f\n',out.motcens_FD_cutoff,out.motcens_DVARS_cutoff);
end
fprintf(logfile_fid,'\nBase filenames:\n');
fprintf(logfile_fid,'Base filename for functional data = %s\n',func_base_filename);
fprintf(logfile_fid,'Base filename for functional mask = %s\n',funcmask_base_filename);
fprintf(logfile_fid,'Base filename for ROI mask = %s, which was input in %s space\n',parc_base_filename,out.parc_space);
if out.WMsig_part == 1
    fprintf(logfile_fid,'Base filename for white matter mask = %s, which was input in %s space\n',WM_base_filename,out.WM_space);
end
if out.ventsig_part == 1
    fprintf(logfile_fid,'Base filename for ventricle mask = %s, which was input in %s space\n',vent_base_filename,out.vent_space);
end
if strcmp(out.parc_space,'Structural') || strcmp(out.WM_space,'Structural') || strcmp(out.vent_space,'Structural')
    fprintf(logfile_fid,'Base filename for structural to functional warp = %s\n',struct2func_regmat_base_filename);
end
if strcmp(out.parc_space,'Standard') || strcmp(out.WM_space,'Standard') || strcmp(out.vent_space,'Standard')
    fprintf(logfile_fid,'Base filename for standard space to functional warp = %s\n',standard2func_regmat_base_filename);
end
fprintf(logfile_fid,'\nData partialing:\n');
if out.motpar_part == 1
    fprintf(logfile_fid,'Motion paramaters were partialed from the timeseries data\n');
    if out.motpart1_part == 1
        fprintf(logfile_fid,'The t - 1 motion paramaters were partialed from the timeseries data\n');
    end
    if out.motparderiv_part == 1
        fprintf(logfile_fid,'The 1st derivatives of the motion paramaters were partialed from the timeseries data\n');
    end
    if out.motparsqr_part == 1
        fprintf(logfile_fid,'The squares of the motion paramaters were partialed from the timeseries data\n');
    end
end
if out.globsig_part == 1 && out.globsig_calcGNI == 0
    fprintf(logfile_fid,'Global signal was partialed from the timeseries data for all participants\n');
    if out.globsigderiv_part == 1
        fprintf(logfile_fid,'The 1st derivative of the global signal was partialed from the timeseries data\n');
    end
elseif out.globsig_part == 1 && out.globsig_calcGNI == 1
    fprintf(logfile_fid,'Global signal was partialed from the timeseries data for participants with GNI < 3\n');
    if out.globsigderiv_part == 1
        fprintf(logfile_fid,'The 1st derivative of the global signal was partialed from the timeseries data for participants with GNI < 3\n');
    end
end
if out.erode_masks == 1 && out.WMsig_part == 1
    fprintf(logfile_fid,'The white matter mask was (3d) eroded with a sphere of diameter %3.0f voxels\n',out.erode_radius);
end
if out.WMsig_part == 1 && out.PCA_WMvent == 1
    fprintf(logfile_fid,'First 5 principal components of white matter signal (from entire WM mask) were partialed from the timeseries data\n');
elseif out.WMsig_part == 1
    if strcmp(out.WMmask_scope,'Entire Mask')
        fprintf(logfile_fid,'White matter signal (from entire WM mask) was partialed from the timeseries data\n');
    elseif out.WMsig_part == 1 && strcmp(out.WMmask_scope,'Local mask')
        fprintf(logfile_fid,'White matter signal (from WM within a ~45mm radius sphere around each voxel) was partialed from the timeseries data\n');
    end
end
if out.WMsig_part == 1 && out.WMsigderiv_part == 1
    fprintf(logfile_fid,'The 1st derivative of the white matter signal(s) was partialed from the timeseries data\n');
end
if out.erode_masks == 1 && out.ventsig_part == 1
    fprintf(logfile_fid,'The ventricle mask was (3d) eroded with a sphere of diameter %3.0f voxels\n',out.erode_radius);
end
if out.ventsig_part == 1
    if out.PCA_WMvent == 1
        fprintf(logfile_fid,'First 5 principal components of ventricular signal were partialed from the timeseries data\n');
    else
        fprintf(logfile_fid,'Ventricular signal was partialed from the timeseries data\n');
    end
    if out.ventsigderiv_part == 1
        fprintf(logfile_fid,'The 1st derivative of the ventricular signal(s) was partialed from the timeseries data\n');
    end
end
fprintf(logfile_fid,'The summary measure used to extract the timeseries from each ROI was the %s\n',out.ROI_summary_measure);
fprintf(logfile_fid,'\nParticipant specific information:\n');

if use_parfor
    subs        = out.subs;
    MC          = out.MC;
    motpar_part = out.motpar_part;
    motcens     = out.motcens;
    parc_space  = out.parc_space;
    WM_space    = out.WM_space;
    vent_space  = out.vent_space;
    loc_out     = out;
    
    parfor currsub = 1:length(subs)
        sub = subs{currsub};
        func_filename          = strrep(func_base_filename,'SUBNUM',sub);
        parc_filename          = strrep(parc_base_filename,'SUBNUM',sub);
        funcmask_filename      = strrep(funcmask_base_filename,'SUBNUM',sub);
        
        if strcmp('WM_base_filename',' ')
            WM_filename        = ' ';
        else
            WM_filename        = strrep(WM_base_filename,'SUBNUM',sub);
        end
        if strcmp('vent_base_filename',' ')
            vent_filename      = ' ';
        else
            vent_filename      = strrep(vent_base_filename,'SUBNUM',sub);
        end
        if MC == 0 && (motpar_part == 1 || motcens == 1) && ~strcmp('par_base_filename',' ')
            par_filename       = strrep(par_base_filename,'SUBNUM',sub);
        else
            par_filename       = ' ';
        end
        if strcmp(parc_space,'Structural')
            parc_warp_filename = strrep(struct2func_regmat_base_filename,'SUBNUM',sub);
        elseif strcmp(parc_space,'Standard')
            parc_warp_filename = strrep(standard2func_regmat_base_filename,'SUBNUM',sub);
        else
            parc_warp_filename = ' ';
        end
        if strcmp(WM_space,'Structural')
            WM_warp_filename   = strrep(struct2func_regmat_base_filename,'SUBNUM',sub);
        elseif strcmp(WM_space,'Standard')
            WM_warp_filename   = strrep(standard2func_regmat_base_filename,'SUBNUM',sub);
        else
            WM_warp_filename   = ' ';
        end
        if strcmp(vent_space,'Structural')
            vent_warp_filename = strrep(struct2func_regmat_base_filename,'SUBNUM',sub);
        elseif strcmp(vent_space,'Standard')
            vent_warp_filename = strrep(standard2func_regmat_base_filename,'SUBNUM',sub);
        else
            vent_warp_filename = ' ';
        end
        
        if exist(func_filename,'file') && exist(parc_filename,'file') && exist(funcmask_filename,'file')
            [ts{currsub,1},num_censored_vols(currsub,1),num_censored_vols_FD(currsub,1),num_censored_vols_DVARS(currsub,1),GCOR(currsub,1),mean_FD(currsub,1),mean_DVARS(currsub,1)] = preprocess_for_graph(loc_out,sub,logfile_fid,img_ext,func_filename,parc_filename,funcmask_filename,WM_filename,vent_filename,parc_warp_filename,WM_warp_filename,vent_warp_filename,par_filename); %#ok<*PFOUS>
        else
            fprintf(logfile_fid,'Missing files for %s\n',sub);
            ts{currsub,1}                      = NaN;
            num_censored_vols(currsub,1)       = NaN;
            num_censored_vols_FD(currsub,1)    = NaN;
            num_censored_vols_DVARS(currsub,1) = NaN;
            GCOR(currsub,1)                    = NaN;
            mean_FD(currsub,1)                 = NaN;
            mean_DVARS(currsub,1)              = NaN;
        end
    end
    out.ts = ts;
    out.num_censored_vols = num_censored_vols;
    out.num_censored_vols_FD = num_censored_vols_FD;
    out.num_censored_vols_DVARS = num_censored_vols_DVARS;
    out.GCOR = GCOR;
    out.mean_FD = mean_FD;
    out.mean_DVARS = mean_DVARS;
else
    progressbar('Total Progress')
    
    for currsub = length(out.subs):-1:1
        sub = out.subs{currsub};
        func_filename          = strrep(func_base_filename,'SUBNUM',sub);
        parc_filename          = strrep(parc_base_filename,'SUBNUM',sub);
        funcmask_filename      = strrep(funcmask_base_filename,'SUBNUM',sub);
        if exist('WM_base_filename','var')
            WM_filename        = strrep(WM_base_filename,'SUBNUM',sub);
        else
            WM_filename        = ' ';
        end
        if exist('vent_base_filename','var')
            vent_filename      = strrep(vent_base_filename,'SUBNUM',sub);
        else
            vent_filename      = ' ';
        end
        if out.MC == 0 && (out.motpar_part == 1 || out.motcens == 1)
            par_filename       = strrep(par_base_filename,'SUBNUM',sub);
        else
            par_filename       = ' ';
        end
        if strcmp(out.parc_space,'Structural')
            parc_warp_filename = strrep(struct2func_regmat_base_filename,'SUBNUM',sub);
        elseif strcmp(out.parc_space,'Standard')
            parc_warp_filename = strrep(standard2func_regmat_base_filename,'SUBNUM',sub);
        else
            parc_warp_filename = ' ';
        end
        if strcmp(out.WM_space,'Structural')
            WM_warp_filename   = strrep(struct2func_regmat_base_filename,'SUBNUM',sub);
        elseif strcmp(out.WM_space,'Standard')
            WM_warp_filename   = strrep(standard2func_regmat_base_filename,'SUBNUM',sub);
        else
            WM_warp_filename   = ' ';
        end
        if strcmp(out.vent_space,'Structural')
            vent_warp_filename = strrep(struct2func_regmat_base_filename,'SUBNUM',sub);
        elseif strcmp(out.vent_space,'Standard')
            vent_warp_filename = strrep(standard2func_regmat_base_filename,'SUBNUM',sub);
        else
            vent_warp_filename = ' ';
        end
        if exist(func_filename,'file') && exist(parc_filename,'file') && exist(funcmask_filename,'file')
            [out.ts{currsub,1},out.num_censored_vols(currsub,1),out.num_censored_vols_FD(currsub,1),out.num_censored_vols_DVARS(currsub,1),out.GCOR(currsub,1),out.mean_FD(currsub,1),out.mean_DVARS(currsub,1)] = preprocess_for_graph(out,sub,logfile_fid,img_ext,func_filename,parc_filename,funcmask_filename,WM_filename,vent_filename,parc_warp_filename,WM_warp_filename,vent_warp_filename,par_filename);
        else
            fprintf(logfile_fid,'Missing files for %s\n',sub);
            out.ts{currsub,1}                      = NaN;
            out.num_censored_vols(currsub,1)       = NaN;
            out.num_censored_vols_FD(currsub,1)    = NaN;
            out.num_censored_vols_DVARS(currsub,1) = NaN;
            out.GCOR(currsub,1)                    = NaN;
            out.mean_FD(currsub,1)                 = NaN;
            out.mean_DVARS(currsub,1)              = NaN;
        end
        
        prog = 1-((currsub-1)/length(out.subs));
        progressbar(prog)
    end
end

if use_parfor
    try
        parpool close
    catch
        matlabpool close
    end
end

fclose(logfile_fid);
save(out.preproc_outname,'out');
set(handles.Start_pushbutton,'enable','on');




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% Conduct preprocessing for a participant
function [ts,num_censored_vols,num_censored_vols_FD,num_censored_vols_DVARS,GCOR,mean_FD,mean_DVARS] = preprocess_for_graph(subfunc_out,sub,logfile_fid,img_ext,orig_func_filename,parc_filename,funcmask_filename,WM_filename,vent_filename,parc_warp_filename,WM_warp_filename,vent_warp_filename,par_filename)

if subfunc_out.BP == 1
    if strcmp(subfunc_out.BP_type,'Butterworth')
        if subfunc_out.BP_HP_cutoff==-1
            [butter_zeds,butter_poles,butter_gain] = butter(subfunc_out.butter_ord,subfunc_out.BP_LP_cutoff/((1/subfunc_out.TR)/2),'low');
            [butter_sos,butter_gain] = zp2sos(butter_zeds,butter_poles,butter_gain);
            butter_filt = dfilt.df2sos(butter_sos,butter_gain);
        elseif subfunc_out.BP_LP_cutoff==-1
            [butter_zeds,butter_poles,butter_gain] = butter(subfunc_out.butter_ord,subfunc_out.BP_HP_cutoff/((1/subfunc_out.TR)/2),'high');
            [butter_sos,butter_gain] = zp2sos(butter_zeds,butter_poles,butter_gain);
            butter_filt = dfilt.df2sos(butter_sos,butter_gain);
        else
            [butter_zeds,butter_poles,butter_gain] = butter(subfunc_out.butter_ord,[subfunc_out.BP_HP_cutoff,subfunc_out.BP_LP_cutoff]/((1/subfunc_out.TR)/2),'bandpass');
            [butter_sos,butter_gain] = zp2sos(butter_zeds,butter_poles,butter_gain);
            butter_filt = dfilt.df2sos(butter_sos,butter_gain);
        end
    elseif strcmp(subfunc_out.BP_type,'Ideal')
        filt_ind = calc_IdealFilter(subfunc_out.nTP,subfunc_out.TR,subfunc_out.BP_HP_cutoff,subfunc_out.BP_LP_cutoff);
    elseif strcmp(subfunc_out.BP_type,'FSL')
        if subfunc_out.BP_HP_cutoff==-1
            HP_sigma = -1;
        else
            HP_sigma = (1./subfunc_out.BP_HP_cutoff)./subfunc_out.TR;
        end
        if subfunc_out.BP_LP_cutoff==-1
            LP_sigma = -1;
        else
            LP_sigma = (1./subfunc_out.BP_LP_cutoff)./subfunc_out.TR;
        end
    end
end

if subfunc_out.ST_correct == 1
    switch subfunc_out.ST_ord
        case 'Sequential Up'
            slice_order = '';
        case 'Sequential Down'
            slice_order = '--down';
        case 'Interleaved'
            slice_order = '--odd';
    end
    if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
        if ~exist(strrep(orig_func_filename,'.nii','_st.nii'),'file')
            eval(['!slicetimer -i ',strrep(orig_func_filename,img_ext,''),' -r ',num2str(subfunc_out.TR),' ',slice_order]);
        end
    else
        eval(['!slicetimer -i ',strrep(orig_func_filename,img_ext,''),' -r ',num2str(subfunc_out.TR),' ',slice_order]);
    end
    orig_func_filename = strrep(orig_func_filename,'.nii','_st.nii');
end

if subfunc_out.MC == 1
    if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
        if ~exist(strrep(orig_func_filename,'.nii','_mcf.nii'),'file')
            eval(['!mcflirt -in ',orig_func_filename,' -plots']);
        end
    else
        eval(['!mcflirt -in ',orig_func_filename,' -plots']);
    end
    orig_func_filename = strrep(orig_func_filename,'.nii','_mcf.nii');
    par_filename = strrep(orig_func_filename,img_ext,'.par');
end

func_filename = orig_func_filename;

if subfunc_out.detr == 1
    if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
        if ~exist(strrep(func_filename,'.nii',['_detr',num2str(subfunc_out.detr_ord),'.nii']),'file')
            detrend_image(func_filename,funcmask_filename,subfunc_out.detr_ord);
        end
    else
        detrend_image(func_filename,funcmask_filename,subfunc_out.detr_ord);
    end
    
    func_filename = strrep(func_filename,'.nii',['_detr',num2str(subfunc_out.detr_ord),'.nii']);
end

if strcmp(subfunc_out.ROIextord,'After all preprocessing')
    if subfunc_out.despike == 1
        if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
            if ~exist(strrep(func_filename,'.nii','_wds.nii'),'file')
                WaveletDespike(func_filename,strrep(func_filename,img_ext,''),'verbose',0,'sp',0);
                func_filename = strrep(func_filename,'.nii','_wds.nii');
                orient = evalc(['!fslorient -getorient ',func_filename]);
                if strcmp(orient(1:12),'NEUROLOGICAL')
                    eval(['!fslswapdim ',func_filename,' -x y z ',func_filename]);
                    eval(['!fslorient -swaporient ',func_filename]);
                end
            else
                func_filename = strrep(func_filename,'.nii','_wds.nii');
                orient = evalc(['!fslorient -getorient ',func_filename]);
                if strcmp(orient(1:12),'NEUROLOGICAL')
                    eval(['!fslswapdim ',func_filename,' -x y z ',func_filename]);
                    eval(['!fslorient -swaporient ',func_filename]);
                end
            end
        else
            WaveletDespike(func_filename,strrep(func_filename,img_ext,''),'verbose',0,'sp',0);
            func_filename = strrep(func_filename,'.nii','_wds.nii');
            orient = evalc(['!fslorient -getorient ',func_filename]);
            if strcmp(orient(1:12),'NEUROLOGICAL')
                eval(['!fslswapdim ',func_filename,' -x y z ',func_filename]);
                eval(['!fslorient -swaporient ',func_filename]);
            end
        end
    end
    
    if ~strcmp(subfunc_out.parc_space,'Functional')
        out_parc_filename = strrep(parc_filename,'.nii','_warp2func.nii');
        if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
            if ~exist(out_parc_filename,'file')
                if any(strfind(parc_warp_filename,'.mat'))
                    eval(['!flirt -interp nearestneighbour -in ',parc_filename,' -ref ',funcmask_filename,' -out ',out_parc_filename,' -init ',parc_warp_filename,' -applyxfm']);
                elseif any(strfind(parc_warp_filename,'.nii'))
                    eval(['!applywarp --ref=',funcmask_filename,' --in=',parc_filename,' --warp=',parc_warp_filename,' --out=',out_parc_filename,' --interp=nn']);
                end
            end
        else
            if any(strfind(parc_warp_filename,'.mat'))
                eval(['!flirt -interp nearestneighbour -in ',parc_filename,' -ref ',funcmask_filename,' -out ',out_parc_filename,' -init ',parc_warp_filename,' -applyxfm']);
            elseif any(strfind(parc_warp_filename,'.nii'))
                eval(['!applywarp --ref=',funcmask_filename,' --in=',parc_filename,' --warp=',parc_warp_filename,' --out=',out_parc_filename,' --interp=nn']);
            end
        end
        parc_filename = out_parc_filename;
    end
    
    if ~strcmp(subfunc_out.WM_space,'Functional')
        out_WM_filename = strrep(WM_filename,'.nii','_warp2func.nii');
        if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
            if ~exist(out_WM_filename,'file')
                if any(strfind(WM_warp_filename,'.mat'))
                    eval(['!flirt -interp nearestneighbour -in ',WM_filename,' -ref ',funcmask_filename,' -out ',out_WM_filename,' -init ',WM_warp_filename,' -applyxfm']);
                elseif any(strfind(vent_warp_filename,'.nii'))
                    eval(['!applywarp --ref=',funcmask_filename,' --in=',WM_filename,' --warp=',WM_warp_filename,' --out=',out_WM_filename,' --interp=nn']);
                end
            end
        else
            if any(strfind(WM_warp_filename,'.mat'))
                eval(['!flirt -interp nearestneighbour -in ',WM_filename,' -ref ',funcmask_filename,' -out ',out_WM_filename,' -init ',WM_warp_filename,' -applyxfm']);
            elseif any(strfind(vent_warp_filename,'.nii'))
                eval(['!applywarp --ref=',funcmask_filename,' --in=',WM_filename,' --warp=',WM_warp_filename,' --out=',out_WM_filename,' --interp=nn']);
            end
        end
        WM_filename = out_WM_filename;
    end

    if ~strcmp(subfunc_out.vent_space,'Functional')
        out_vent_filename = strrep(vent_filename,'.nii','_warp2func.nii');
        if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
            if ~exist(out_vent_filename,'file')
                if any(strfind(vent_warp_filename,'.mat'))
                    eval(['!flirt -interp nearestneighbour -in ',vent_filename,' -ref ',funcmask_filename,' -out ',out_vent_filename,' -init ',vent_warp_filename,' -applyxfm']);
                elseif any(strfind(vent_warp_filename,'.nii'))
                    eval(['!applywarp --ref=',funcmask_filename,' --in=',vent_filename,' --warp=',vent_warp_filename,' --out=',out_vent_filename,' --interp=nn']);
                end
            end
        else
            if any(strfind(vent_warp_filename,'.mat'))
                eval(['!flirt -interp nearestneighbour -in ',vent_filename,' -ref ',funcmask_filename,' -out ',out_vent_filename,' -init ',vent_warp_filename,' -applyxfm']);
            elseif any(strfind(vent_warp_filename,'.nii'))
                eval(['!applywarp --ref=',funcmask_filename,' --in=',vent_filename,' --warp=',vent_warp_filename,' --out=',out_vent_filename,' --interp=nn']);
            end
        end
        vent_filename = out_vent_filename;
    end
    
    desmat = [];
    
    if subfunc_out.WMsig_part == 1
        if subfunc_out.erode_masks == 1
            WM_eroded_filename = strrep(WM_filename,img_ext,['_eroded_',num2str(subfunc_out.erode_radius),img_ext]);
            if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
                if ~exist(WM_eroded_filename,'file')
                    tempmask = load_nii_gz(WM_filename);
                    [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                    save_nii_gz(tempmask,WM_eroded_filename);
                end
            else
                tempmask = load_nii_gz(WM_filename);
                [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                save_nii_gz(tempmask,WM_eroded_filename);
            end
            WM_filename = WM_eroded_filename;
        end
        if strcmp(subfunc_out.WMmask_scope,'Entire Mask')
            if subfunc_out.PCA_WMvent == 1
                WM_extract_filename = strrep(WM_filename,img_ext,'_EIG.txt');
                if ~exist(WM_extract_filename,'file')
                    eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename,' --eig --order=5']);
                else
                    vernum = 0;
                    while 1
                        if exist(strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                            vernum = vernum+1;
                        else
                            WM_extract_filename = strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                            break
                        end
                    end
                    eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename,' --eig --order=5']);
                end
                WM_preds = dlmread(WM_extract_filename);
                WM_preds = spm_detrend(WM_preds,subfunc_out.detr_ord);
                desmat = [desmat,WM_preds];
                
                if subfunc_out.WMsigderiv_part == 1
                    WM1stderiv_preds = WM_preds(3:end,:) - WM_preds(1:(end-2),:);
                    WM1stderiv_preds = [WM1stderiv_preds(1,:);WM1stderiv_preds;WM1stderiv_preds(end,:)];
                    WM1stderiv_preds = spm_detrend(WM1stderiv_preds,subfunc_out.detr_ord);
                    desmat = [desmat,WM1stderiv_preds];
                end
            else
                WM_extract_filename = strrep(WM_filename,img_ext,'.txt');
                if ~exist(WM_extract_filename,'file')
                    eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename]);
                else
                    vernum = 0;
                    while 1
                        if exist(strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                            vernum = vernum+1;
                        else
                            WM_extract_filename = strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                            break
                        end
                    end
                    eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename]);
                end
                WM_pred = dlmread(WM_extract_filename);
                WM_pred = spm_detrend(WM_pred,subfunc_out.detr_ord);
                desmat = [desmat,WM_pred];
                
                if subfunc_out.WMsigderiv_part == 1
                    WM1stderiv_pred = WM_pred(3:end) - WM_pred(1:(end-2));
                    WM1stderiv_pred = [WM1stderiv_pred(1);WM1stderiv_pred;WM1stderiv_pred(end)];
                    WM1stderiv_pred = spm_detrend(WM1stderiv_pred,subfunc_out.detr_ord);
                    desmat = [desmat,WM1stderiv_pred];
                end
            end
        else
            full_WM_mask = load_nii_gz(WM_filename);
            WM_fmask_inds = find(full_WM_mask.img == 1);
            WM_rad = round(45/mean(full_WM_mask.hdr.dime.pixdim(2:4)));
            WM_fmask_dims = size(full_WM_mask.img);
            clear full_WM_mask
            [vx,vy,vz] = meshgrid(-WM_rad:WM_rad);
            V = sqrt((vx.^2)+(vy.^2)+(vz.^2));
            V(V <= WM_rad) = 1;
            V(V > WM_rad) = 0;
            [WM_sp_x,WM_sp_y,WM_sp_z] = ind2sub(size(V),find(V == 1));
            WM_sp_x = WM_sp_x - WM_rad - 1;
            WM_sp_y = WM_sp_y - WM_rad - 1;
            WM_sp_z = WM_sp_z - WM_rad - 1;
        end
    end
    
    if subfunc_out.ventsig_part == 1
        if subfunc_out.erode_masks == 1
            vent_eroded_filename = strrep(vent_filename,img_ext,['_eroded_',num2str(subfunc_out.erode_radius),img_ext]);
            if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
                if ~exist(vent_eroded_filename,'file')
                    tempmask = load_nii_gz(vent_filename);
                    [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                    save_nii_gz(tempmask,vent_eroded_filename);
                end
            else
                tempmask = load_nii_gz(vent_filename);
                [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                save_nii_gz(tempmask,vent_eroded_filename);
            end
            vent_filename = vent_eroded_filename;
        end
        if subfunc_out.PCA_WMvent == 1
            vent_extract_filename = strrep(vent_filename,img_ext,'_EIG.txt');
            if ~exist(vent_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename,' --eig --order=5']);
            else
                vernum = 0;
                while 1
                    if exist(strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        vent_extract_filename = strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename,' --eig --order=5']);
            end
            vent_preds = dlmread(vent_extract_filename);
            vent_preds = spm_detrend(vent_preds,subfunc_out.detr_ord);
            desmat = [desmat,vent_preds];
            
            if subfunc_out.ventsigderiv_part == 1
                vent1stderiv_preds = vent_preds(3:end,:) - vent_preds(1:(end-2),:);
                vent1stderiv_preds = [vent1stderiv_preds(1,:);vent1stderiv_preds;vent1stderiv_preds(end,:)];
                vent1stderiv_preds = spm_detrend(vent1stderiv_preds,subfunc_out.detr_ord);
                desmat = [desmat,vent1stderiv_preds];
            end
        else
            vent_extract_filename = strrep(vent_filename,img_ext,'.txt');
            if ~exist(vent_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename]);
            else
                vernum = 0;
                while 1
                    if exist(strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        vent_extract_filename = strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename]);
            end
            vent_pred = dlmread(vent_extract_filename);
            vent_pred = spm_detrend(vent_pred,subfunc_out.detr_ord);
            desmat = [desmat,vent_pred];
            
            if subfunc_out.ventsigderiv_part == 1
                vent1stderiv_pred = vent_pred(3:end) - vent_pred(1:(end-2));
                vent1stderiv_pred = [vent1stderiv_pred(1);vent1stderiv_pred;vent1stderiv_pred(end)];
                vent1stderiv_pred = spm_detrend(vent1stderiv_pred,subfunc_out.detr_ord);
                desmat = [desmat,vent1stderiv_pred];
            end
        end
    end
    
    if subfunc_out.globsig_part == 1
        if subfunc_out.globsig_calcGNI == 1
            GNI = Rest2GNI_edit(func_filename,funcmask_filename,100);
            fprintf(logfile_fid,'Original GNI for %s = %6.3f\n',sub,GNI);
        else
            GNI = 0;
        end
        if GNI <= 3
            funcmask_extract_filename = strrep(funcmask_filename,img_ext,'.txt');
            if ~exist(funcmask_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',funcmask_extract_filename,' -m ',funcmask_filename]);
            else
                vernum = 0;
                while 1
                    if exist(strrep(funcmask_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        funcmask_extract_filename = strrep(funcmask_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',funcmask_extract_filename,' -m ',funcmask_filename]);
            end
            global_pred = dlmread(funcmask_extract_filename);
            global_pred = spm_detrend(global_pred,subfunc_out.detr_ord);
            desmat = [desmat,global_pred];
            
            if subfunc_out.globsigderiv_part == 1
                global1stderiv_pred = global_pred(3:end) - global_pred(1:(end-2));
                global1stderiv_pred = [global1stderiv_pred(1);global1stderiv_pred;global1stderiv_pred(end)];
                global1stderiv_pred = spm_detrend(global1stderiv_pred,subfunc_out.detr_ord);
                desmat = [desmat,global1stderiv_pred];
            end
        end
    end
    
    if subfunc_out.motpar_part == 1
        mot_pred = dlmread(par_filename);
        for par = 1:size(mot_pred,2)
            mot_pred(:,par) = spm_detrend(mot_pred(:,par),subfunc_out.detr_ord);
        end
        desmat = [desmat,mot_pred];
        if subfunc_out.motpart1_part == 1
            mot1back_pred = [mot_pred(2:end,:);mot_pred(end,:)];
            for par = 1:size(mot_pred,2)
                mot1back_pred(:,par) = spm_detrend(mot1back_pred(:,par),subfunc_out.detr_ord);
            end
            desmat = [desmat,mot1back_pred];
        end
        if subfunc_out.motparderiv_part == 1
            mot1stderiv_pred = mot_pred(3:end,:) - mot_pred(1:(end-2),:);
            mot1stderiv_pred = [mot1stderiv_pred(1,:);mot1stderiv_pred;mot1stderiv_pred(end,:)];
            for par = 1:size(mot_pred,2)
                mot1stderiv_pred(:,par) = spm_detrend(mot1stderiv_pred(:,par),subfunc_out.detr_ord);
            end
            desmat = [desmat,mot1stderiv_pred];
        end
        if subfunc_out.motparsqr_part == 1
            motsquare_pred = mot_pred.^2;
            for par = 1:size(mot_pred,2)
                motsquare_pred(:,par) = spm_detrend(motsquare_pred(:,par),subfunc_out.detr_ord);
            end
            desmat = [desmat,motsquare_pred];
            if subfunc_out.motpart1_part == 1
                mot1backsquare_pred = mot1back_pred.^2;
                for par = 1:size(mot_pred,2)
                    mot1backsquare_pred(:,par) = spm_detrend(mot1backsquare_pred(:,par),subfunc_out.detr_ord);
                end
                desmat = [desmat,mot1backsquare_pred];
            end
        end
    end
    
    func_data = load_nii_gz(func_filename);
    funcmask = load_nii_gz(funcmask_filename);
    func_data.img = double(func_data.img);
    funcmask.img = double(funcmask.img);
    
    out_func_data.img = zeros(size(func_data.img));
    
    if subfunc_out.BP == 1 && (subfunc_out.WMsig_part == 1 || subfunc_out.ventsig_part == 1 || (subfunc_out.globsig_part == 1 && GNI <= 3) || subfunc_out.motpar_part == 1)
        desmat = [ones(size(desmat,1),1),desmat];
        orig_desmat = desmat;
        if rank(desmat) ~= min(size(desmat))
            error(['Initial desmat to partial out unwanted variance for ',sub,' is not full rank']);
        end
        func_filename = strrep(func_filename,'.nii','_remnuisance_bpfilt.nii');
        Q = desmat*pinv(full(desmat));
        for x = 1:size(func_data.img,1)
            for y = 1:size(func_data.img,2)
                for z = 1:size(func_data.img,3)
                    if funcmask.img(x,y,z) == 1
                        if strcmp(subfunc_out.WMmask_scope,'Local Mask')
                            curr_WM_sp_x = WM_sp_x + x;
                            curr_WM_sp_y = WM_sp_y + y;
                            curr_WM_sp_z = WM_sp_z + z;
                            inrangeinds = curr_WM_sp_x > 0 & curr_WM_sp_x <= size(func_data.img,1) & curr_WM_sp_y > 0 & curr_WM_sp_y <= size(func_data.img,2) & curr_WM_sp_z > 0 & curr_WM_sp_z <= size(func_data.img,3);
                            WM_sp_inds = sub2ind(WM_fmask_dims,curr_WM_sp_x(inrangeinds),curr_WM_sp_y(inrangeinds),curr_WM_sp_z(inrangeinds));
                            loc_WM_inds = intersect(WM_fmask_inds,WM_sp_inds);
                            if ~isempty(loc_WM_inds)
                                for tp = 1:size(func_data.img,4)
                                    tp_func = func_data.img(:,:,:,tp);
                                    loc_WM_pred(tp,1) = sum(tp_func(loc_WM_inds))/length(loc_WM_inds); %#ok<*AGROW>
                                end
                                loc_WM_pred = spm_detrend(loc_WM_pred,subfunc_out.detr_ord);
                                desmat = [orig_desmat,loc_WM_pred];
                                if subfunc_out.WMsigderiv_part == 1
                                    loc_WM1stderiv_pred = loc_WM_pred(3:end) - loc_WM_pred(1:(end-2));
                                    loc_WM1stderiv_pred = [loc_WM1stderiv_pred(1);loc_WM1stderiv_pred;loc_WM1stderiv_pred(end)];
                                    loc_WM1stderiv_pred = spm_detrend(loc_WM1stderiv_pred,subfunc_out.detr_ord);
                                    desmat = [desmat,loc_WM1stderiv_pred];
                                end
                            end
                            Q = desmat*pinv(full(desmat));
                        end
                        tempdata = squeeze(func_data.img(x,y,z,:))-(Q*squeeze(func_data.img(x,y,z,:)));
                        tempdata = tempdata-mean(tempdata);
                        
                        if exist('butter_filt','var')
                            tempdata_padded = [zeros(15,size(tempdata,2));tempdata;zeros(15,size(tempdata,2))];
                            temp_butter = filter(butter_filt,tempdata_padded);
                            temp_butter(end:-1:1) = filter(butter_filt,temp_butter(end:-1:1));
                            out_func_data.img(x,y,z,:) = temp_butter(16:(end-15));
                        elseif exist('filt_ind','var')
                            padlength = rest_nextpow2_one35(subfunc_out.nTP);
                            tempdata = tempdata-repmat(mean(tempdata),size(tempdata,1),1);
                            tempdata = [tempdata;zeros(padlength-subfunc_out.nTP,size(tempdata,2))];
                            freq = fft(tempdata);
                            freq(filt_ind,:) = 0;
                            tempdata = ifft(freq);
                            out_func_data.img(x,y,z,:) = tempdata(1:subfunc_out.nTP,:);
                        else
                            out_func_data.img(x,y,z,:) = tempdata;
                        end
                    end
                end
            end
        end
        if exist('HP_sigma','var')
            out_func_data.hdr = func_data.hdr;
            out_func_data.hdr.dime.dim(2:5) = size(out_func_data.img);
            out_func_data.hdr.dime.cal_max = max(out_func_data.img(:));
            out_func_data.hdr.dime.cal_min = min(out_func_data.img(out_func_data.img~=0));
            out_func_data.hdr.dime.datatype = 64;
            out_func_data.hdr.dime.bitpix = 64;
            save_nii_gz(out_func_data,func_filename);
            eval(['!fslmaths ',func_filename,' -bptf ',num2str(HP_sigma),' ',num2str(LP_sigma),' ',func_filename]);
            out_func_data = load_nii_gz(func_filename);
        end
    elseif subfunc_out.BP == 1
        func_filename = strrep(func_filename,'.nii','_bpfilt.nii');
        if exist('HP_sigma','var')
            out_func_data.img = func_data.img;
            out_func_data.hdr = func_data.hdr;
            out_func_data.hdr.dime.dim(2:5) = size(out_func_data.img);
            out_func_data.hdr.dime.cal_max = max(out_func_data.img(:));
            out_func_data.hdr.dime.cal_min = min(out_func_data.img(out_func_data.img~=0));
            out_func_data.hdr.dime.datatype = 64;
            out_func_data.hdr.dime.bitpix = 64;
            save_nii_gz(out_func_data,func_filename);
            eval(['!fslmaths ',func_filename,' -bptf ',num2str(HP_sigma),' ',num2str(LP_sigma),' ',func_filename]);
            out_func_data = load_nii_gz(func_filename);
        else
            for x = 1:size(func_data.img,1)
                for y = 1:size(func_data.img,2)
                    for z = 1:size(func_data.img,3)
                        if funcmask.img(x,y,z) == 1
                            tempdata = squeeze(func_data.img(x,y,z,:));
                            tempdata = tempdata-mean(tempdata);
                            if exist('butter_filt','var')
                                tempdata_padded = [zeros(15,size(tempdata,2));tempdata;zeros(15,size(tempdata,2))];
                                temp_butter = filter(butter_filt,tempdata_padded);
                                temp_butter(end:-1:1) = filter(butter_filt,temp_butter(end:-1:1));
                                out_func_data.img(x,y,z,:) = temp_butter(16:(end-15));
                            else
                                padlength = rest_nextpow2_one35(subfunc_out.nTP);
                                tempdata = [tempdata;zeros(padlength-subfunc_out.nTP,size(tempdata,2))];
                                freq = fft(tempdata);
                                freq(filt_ind,:) = 0;
                                tempdata = ifft(freq);
                                out_func_data.img(x,y,z,:) = tempdata(1:subfunc_out.nTP,:);
                            end
                        end
                    end
                end
            end
        end
    elseif subfunc_out.WMsig_part == 1 || subfunc_out.ventsig_part == 1 || (subfunc_out.globsig_part == 1 && GNI <= 3) || subfunc_out.motpar_part == 1
        desmat = [ones(size(desmat,1),1),desmat];
        orig_desmat = desmat;
        if rank(desmat) ~= min(size(desmat))
            error(['Initial desmat to partial out unwanted variance for ' sub ' is not full rank']);
        end
        func_filename = strrep(func_filename,'.nii','_remnuisance.nii');
        Q = desmat*pinv(full(desmat));
        for x = 1:size(func_data.img,1)
            for y = 1:size(func_data.img,2)
                for z = 1:size(func_data.img,3)
                    if funcmask.img(x,y,z) == 1
                        if strcmp(subfunc_out.WMmask_scope,'Local Mask')
                            curr_WM_sp_x = WM_sp_x + x;
                            curr_WM_sp_y = WM_sp_y + y;
                            curr_WM_sp_z = WM_sp_z + z;
                            inrangeinds = curr_WM_sp_x > 0 & curr_WM_sp_x <= size(func_data.img,1) & curr_WM_sp_y > 0 & curr_WM_sp_y <= size(func_data.img,2) & curr_WM_sp_z > 0 & curr_WM_sp_z <= size(func_data.img,3);
                            WM_sp_inds = sub2ind(WM_fmask_dims,curr_WM_sp_x(inrangeinds),curr_WM_sp_y(inrangeinds),curr_WM_sp_z(inrangeinds));
                            loc_WM_inds = intersect(WM_fmask_inds,WM_sp_inds);
                            if ~isempty(loc_WM_inds)
                                for tp = 1:size(func_data.img,4)
                                    tp_func = func_data.img(:,:,:,tp);
                                    loc_WM_pred(tp,1) = sum(tp_func(loc_WM_inds))/length(loc_WM_inds);
                                end
                                loc_WM_pred = spm_detrend(loc_WM_pred,subfunc_out.detr_ord);
                                desmat = [orig_desmat,loc_WM_pred];
                                if subfunc_out.WMsigderiv_part == 1
                                    loc_WM1stderiv_pred = loc_WM_pred(3:end) - loc_WM_pred(1:(end-2));
                                    loc_WM1stderiv_pred = [loc_WM1stderiv_pred(1);loc_WM1stderiv_pred;loc_WM1stderiv_pred(end)];
                                    loc_WM1stderiv_pred = spm_detrend(loc_WM1stderiv_pred,subfunc_out.detr_ord);
                                    desmat = [desmat,loc_WM1stderiv_pred];
                                end
                            end
                            Q = desmat*pinv(full(desmat));
                        end
                        out_func_data.img(x,y,z,:) = squeeze(func_data.img(x,y,z,:))-(Q*squeeze(func_data.img(x,y,z,:)));
                    end
                end
            end
        end
    else
        out_func_data.img = func_data.img;
    end
    
    if strcmp(subfunc_out.saveintermed,'Yes')
        out_func_data.hdr = func_data.hdr;
        out_func_data.hdr.dime.dim(2:5) = size(out_func_data.img);
        out_func_data.hdr.dime.cal_max = max(out_func_data.img(:));
        out_func_data.hdr.dime.cal_min = min(out_func_data.img(out_func_data.img~=0));
        out_func_data.hdr.dime.datatype = 64;
        out_func_data.hdr.dime.bitpix = 64;
        save_nii_gz(out_func_data,func_filename);
    end
    
    if subfunc_out.motcens == 1
        mot_pred = dlmread(par_filename);
        FD = cat(1,0,sum(abs(diff([(5000-5000*cos(mot_pred(:,1:3))),mot_pred(:,4:6)])),2));
        mean_FD = mean(FD);
        
        [~,DVARS] = clct_DVARS(out_func_data.img,funcmask.img);
        DVARS = cat(1,0,DVARS);
        mean_DVARS = mean(DVARS);
        
        num_censored_vols_FD = sum(FD >= subfunc_out.motcens_FD_cutoff);
        num_censored_vols_DVARS = sum(DVARS >= subfunc_out.motcens_DVARS_cutoff);
        
        temporal_mask = logical(FD < subfunc_out.motcens_FD_cutoff & DVARS < subfunc_out.motcens_DVARS_cutoff);
        censored_vols = [0;find(temporal_mask == 0)];
        for c_vol = 1:(length(censored_vols)-1)
            seg_length = censored_vols(c_vol+1)-censored_vols(c_vol);
            if seg_length <= 5
                temporal_mask(censored_vols(c_vol)+1:censored_vols(c_vol+1)) = 0;
                temporal_mask(censored_vols(c_vol)+1:censored_vols(c_vol+1)) = 0;
            end
        end
        num_censored_vols = sum(temporal_mask == 0);
        if num_censored_vols > 0
            if (subfunc_out.nTP-num_censored_vols >= 50) && (num_censored_vols/subfunc_out.nTP <= .75)
                func_filename = orig_func_filename;
                
                if subfunc_out.detr == 1
                    detrend_image(func_filename,funcmask_filename,subfunc_out.detr_ord,temporal_mask);
                    func_filename = strrep(orig_func_filename,'.nii',['_motcensor_detr',num2str(subfunc_out.detr_ord),'.nii']);
                end
                
                desmat = [];
                
                if subfunc_out.WMsig_part == 1 && strcmp(subfunc_out.WMmask_scope,'Entire Mask')
                    if subfunc_out.PCA_WMvent == 1
                        WM_extract_filename = strrep(WM_filename,img_ext,'_EIG_motcensor.txt');
                        if ~exist(WM_extract_filename,'file')
                            eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename,' --eig --order=5']);
                        else
                            vernum = 0;
                            while 1
                                if exist(strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                                    vernum = vernum+1;
                                else
                                    WM_extract_filename = strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                                    break
                                end
                            end
                            eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename,' --eig --order=5']);
                        end
                        WM_preds = dlmread(WM_extract_filename);
                        WM_preds = spm_detrend(WM_preds,subfunc_out.detr_ord);
                        desmat = [desmat,WM_preds];
                        
                        if subfunc_out.WMsigderiv_part == 1
                            WM1stderiv_preds = WM_preds(3:end,:) - WM_preds(1:(end-2),:);
                            WM1stderiv_preds = [WM1stderiv_preds(1,:);WM1stderiv_preds;WM1stderiv_preds(end,:)];
                            WM1stderiv_preds = spm_detrend(WM1stderiv_preds,subfunc_out.detr_ord);
                            desmat = [desmat,WM1stderiv_preds];
                        end
                    else
                        WM_extract_filename = strrep(WM_filename,img_ext,'_motcensor.txt');
                        if ~exist(WM_extract_filename,'file')
                            eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename]);
                        else
                            vernum = 0;
                            while 1
                                if exist(strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                                    vernum = vernum+1;
                                else
                                    WM_extract_filename = strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                                    break
                                end
                            end
                            eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename]);
                        end
                        WM_pred = dlmread(WM_extract_filename);
                        WM_pred = spm_detrend(WM_pred,subfunc_out.detr_ord);
                        desmat = [desmat,WM_pred];
                        if subfunc_out.WMsigderiv_part == 1
                            WM1stderiv_pred = WM_pred(3:end) - WM_pred(1:(end-2));
                            WM1stderiv_pred = [WM1stderiv_pred(1);WM1stderiv_pred;WM1stderiv_pred(end)];
                            WM1stderiv_pred = spm_detrend(WM1stderiv_pred,subfunc_out.detr_ord);
                            desmat = [desmat,WM1stderiv_pred];
                        end
                    end
                end
                
                if subfunc_out.ventsig_part == 1
                    if subfunc_out.PCA_WMvent == 1
                        vent_extract_filename = strrep(vent_filename,img_ext,'_EIG_motcensor.txt');
                        if ~exist(vent_extract_filename,'file')
                            eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename,' --eig --order=5']);
                        else
                            vernum = 0;
                            while 1
                                if exist(strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                                    vernum = vernum+1;
                                else
                                    vent_extract_filename = strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                                    break
                                end
                            end
                            eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename,' --eig --order=5']);
                        end
                        vent_preds = dlmread(vent_extract_filename);
                        vent_preds = spm_detrend(vent_preds,subfunc_out.detr_ord);
                        desmat = [desmat,vent_preds];
                        
                        if subfunc_out.ventsigderiv_part == 1
                            vent1stderiv_preds = vent_preds(3:end,:) - vent_preds(1:(end-2),:);
                            vent1stderiv_preds = [vent1stderiv_preds(1,:);vent1stderiv_preds;vent1stderiv_preds(end,:)];
                            vent1stderiv_preds = spm_detrend(vent1stderiv_preds,subfunc_out.detr_ord);
                            desmat = [desmat,vent1stderiv_preds];
                        end
                    else
                        vent_extract_filename = strrep(vent_filename,img_ext,'_motcensor.txt');
                        if ~exist(vent_extract_filename,'file')
                            eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename]);
                        else
                            vernum = 0;
                            while 1
                                if exist(strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                                    vernum = vernum+1;
                                else
                                    vent_extract_filename = strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                                    break
                                end
                            end
                            eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename]);
                        end
                        vent_pred = dlmread(vent_extract_filename);
                        vent_pred = spm_detrend(vent_pred,subfunc_out.detr_ord);
                        desmat = [desmat,vent_pred];
                        
                        if subfunc_out.ventsigderiv_part == 1
                            vent1stderiv_pred = vent_pred(3:end) - vent_pred(1:(end-2));
                            vent1stderiv_pred = [vent1stderiv_pred(1);vent1stderiv_pred;vent1stderiv_pred(end)];
                            vent1stderiv_pred = spm_detrend(vent1stderiv_pred,subfunc_out.detr_ord);
                            desmat = [desmat,vent1stderiv_pred];
                        end
                    end
                end
                
                if subfunc_out.globsig_part == 1
                    if subfunc_out.globsig_calcGNI == 1
                        GNI = Rest2GNI_edit(func_filename,funcmask_filename,1000); % If your machine freezes/slows down tremendously at this step b/c you have low RAM, you can change the 1000 value to 100 (or less) and it will use less RAM
                        fprintf(logfile_fid,'Motion-censored GNI for %s = %6.3f\n',sub,GNI);
                    else
                        GNI = 0;
                    end
                    if GNI <= 3
                        funcmask_extract_filename = strrep(funcmask_filename,img_ext,'_motcensor.txt');
                        
                        if ~exist(funcmask_extract_filename,'file')
                            eval(['!fslmeants -i ',func_filename,' -o ',funcmask_extract_filename,' -m ',funcmask_filename]);
                        else
                            vernum = 0;
                            while 1
                                if exist(strrep(funcmask_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                                    vernum = vernum+1;
                                else
                                    funcmask_extract_filename = strrep(funcmask_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                                    break
                                end
                            end
                            eval(['!fslmeants -i ',func_filename,' -o ',funcmask_extract_filename,' -m ',funcmask_filename]);
                        end
                        global_pred = dlmread(funcmask_extract_filename);
                        global_pred = spm_detrend(global_pred,subfunc_out.detr_ord);
                        desmat = [desmat,global_pred];
                        if subfunc_out.globsigderiv_part == 1
                            global1stderiv_pred = global_pred(3:end) - global_pred(1:(end-2));
                            global1stderiv_pred = [global1stderiv_pred(1);global1stderiv_pred;global1stderiv_pred(end)];
                            global1stderiv_pred = spm_detrend(global1stderiv_pred,subfunc_out.detr_ord);
                            desmat = [desmat,global1stderiv_pred];
                        end
                    end
                end
                
                if subfunc_out.motpar_part == 1
                    mot_pred = dlmread(par_filename);
                    if subfunc_out.motpart1_part == 1
                        mot1back_pred = [mot_pred(2:end,:);mot_pred(end,:)];
                        mot1back_pred = mot1back_pred(temporal_mask,:);
                        for par = 1:size(mot_pred,2)
                            mot1back_pred(:,par) = spm_detrend(mot1back_pred(:,par),subfunc_out.detr_ord);
                        end
                        desmat = [desmat,mot1back_pred];
                    end
                    if subfunc_out.motparderiv_part == 1
                        mot1stderiv_pred = mot_pred(3:end,:) - mot_pred(1:(end-2),:);
                        mot1stderiv_pred = [mot1stderiv_pred(1,:);mot1stderiv_pred;mot1stderiv_pred(end,:)];
                        mot1stderiv_pred = mot1stderiv_pred(temporal_mask,:);
                        for par = 1:size(mot_pred,2)
                            mot1stderiv_pred(:,par) = spm_detrend(mot1stderiv_pred(:,par),subfunc_out.detr_ord);
                        end
                        desmat = [desmat,mot1stderiv_pred];
                    end
                    if subfunc_out.motparsqr_part == 1
                        motsquare_pred = mot_pred.^2;
                        motsquare_pred = motsquare_pred(temporal_mask,:);
                        for par = 1:size(mot_pred,2)
                            motsquare_pred(:,par) = spm_detrend(motsquare_pred(:,par),subfunc_out.detr_ord);
                        end
                        desmat = [desmat,motsquare_pred];
                    end
                    mot_pred = mot_pred(temporal_mask,:);
                    for par = 1:size(mot_pred,2)
                        mot_pred(:,par) = spm_detrend(mot_pred(:,par),subfunc_out.detr_ord);
                    end
                    desmat = [desmat,mot_pred];
                end
                
                func_data = load_nii_gz(func_filename);
                func_data.img = double(func_data.img);
                
                out_func_data.img = zeros(size(func_data.img));
                if subfunc_out.WMsig_part == 1 || subfunc_out.ventsig_part == 1 || (subfunc_out.globsig_part == 1 && GNI <= 3) || subfunc_out.motpar_part == 1
                    desmat = [ones(size(desmat,1),1),desmat];
                    orig_desmat = desmat;
                    if rank(desmat) ~= min(size(desmat))
                        error(['Motion censoring desmat to partial out unwanted variance for ' sub ' is not full rank']);
                    end
                    func_filename = strrep(func_filename,'.nii','_remnuisance.nii');
                    Q = desmat*pinv(full(desmat));
                    for x = 1:size(func_data.img,1)
                        for y = 1:size(func_data.img,2)
                            for z = 1:size(func_data.img,3)
                                if funcmask.img(x,y,z) == 1
                                    if strcmp(subfunc_out.WMmask_scope,'Local Mask')
                                        curr_WM_sp_x = WM_sp_x + x;
                                        curr_WM_sp_y = WM_sp_y + y;
                                        curr_WM_sp_z = WM_sp_z + z;
                                        WM_sp_inds = sub2ind(WM_fmask_dims,[curr_WM_sp_x,curr_WM_sp_y,curr_WM_sp_z]);
                                        loc_WM_inds = intersect(WM_fmask_inds,WM_sp_inds);
                                        if ~isempty(loc_WM_inds)
                                            for tp = 1:size(func_data.img,4)
                                                tp_func = func_data.img(:,:,:,tp);
                                                loc_WM_pred(tp,1) = mean(tp_func(loc_WM_inds));
                                            end
                                            loc_WM_pred = spm_detrend(loc_WM_pred,subfunc_out.detr_ord);
                                            desmat = [orig_desmat,loc_WM_pred];
                                            if subfunc_out.WMsigderiv_part == 1
                                                loc_WM1stderiv_pred = loc_WM_pred(3:end) - loc_WM_pred(1:(end-2));
                                                loc_WM1stderiv_pred = [loc_WM1stderiv_pred(1);loc_WM1stderiv_pred;loc_WM1stderiv_pred(end)];
                                                loc_WM1stderiv_pred = spm_detrend(loc_WM1stderiv_pred,subfunc_out.detr_ord);
                                                desmat = [desmat,loc_WM1stderiv_pred];
                                            end
                                        end
                                        Q = desmat*pinv(full(desmat));
                                    end
                                    out_func_data.img(x,y,z,:) = squeeze(func_data.img(x,y,z,:))-(Q*squeeze(func_data.img(x,y,z,:)));
                                end
                            end
                        end
                    end
                    out_func_data.hdr = func_data.hdr;
                    out_func_data.hdr.dime.dim(2:5) = size(out_func_data.img);
                    out_func_data.hdr.dime.cal_max = max(out_func_data.img(:));
                    out_func_data.hdr.dime.cal_min = min(out_func_data.img(out_func_data.img~=0));
                    out_func_data.hdr.dime.datatype = 64;
                    out_func_data.hdr.dime.bitpix = 64;
                    if strcmp(subfunc_out.saveintermed,'Yes')
                        save_nii_gz(out_func_data,func_filename);
                    end
                else
                    out_func_data = func_data;
                end
                
                if subfunc_out.BP == 1
                    func_filename = strrep(func_filename,'.nii','_bpfilt.nii');
                    orig_times = (1:subfunc_out.nTP);
                    uncensored_times = orig_times(temporal_mask);
                    freq_bins = 1:1:subfunc_out.nTP;
                    for bin = length(freq_bins):-1:1
                        theta = (1/(2*freq_bins(bin)))*atan(sum(sin(2*freq_bins(bin)*uncensored_times))/sum(cos(2*freq_bins(bin)*uncensored_times)));
                        cosmult(:,bin) = cos(freq_bins(bin).*(uncensored_times-theta))';
                        cosdiv(bin) = sum(cos(freq_bins(bin).*(uncensored_times-theta)).^2);
                        sinmult(:,bin) = sin(freq_bins(bin).*(uncensored_times-theta))';
                        sindiv(bin) = sum(sin(freq_bins(bin).*(uncensored_times-theta)).^2);
                        interpmult(:,:,bin) = [cos(freq_bins(bin).*(orig_times-theta));sin(freq_bins(bin).*(orig_times-theta))];
                    end
                    out_func_data.img = out_func_data.img - repmat(mean(out_func_data.img,4),[1,1,1,size(out_func_data.img,4)]);
                    
                    if exist('HP_sigma','var')
                        temp_data = zeros([size(func_data.img,1),size(func_data.img,2),size(func_data.img,3),subfunc_out.nTP]);
                        for x = 1:size(func_data.img,1)
                            for y = 1:size(func_data.img,2)
                                for z = 1:size(func_data.img,3)
                                    if funcmask.img(x,y,z) == 1
                                        for bin = length(freq_bins):-1:1
                                            cosval = sum(squeeze(out_func_data.img(x,y,z,:)).*cosmult(:,bin))./cosdiv(bin);
                                            sinval = sum(squeeze(out_func_data.img(x,y,z,:)).*sinmult(:,bin))./sindiv(bin);
                                            interp_data(:,bin) = [cosval,sinval]*interpmult(:,:,bin);
                                        end
                                        temp_data(x,y,z,:) = sum(interp_data,2);
                                    end
                                end
                            end
                        end
                        out_func_data.img = temp_data;
                        out_func_data.hdr = func_data.hdr;
                        out_func_data.hdr.dime.dim(2:5) = size(out_func_data.img);
                        out_func_data.hdr.dime.cal_max = max(out_func_data.img(:));
                        out_func_data.hdr.dime.cal_min = min(out_func_data.img(out_func_data.img~=0));
                        out_func_data.hdr.dime.datatype = 64;
                        out_func_data.hdr.dime.bitpix = 64;
                        save_nii_gz(out_func_data,func_filename);
                        eval(['!fslmaths ',func_filename,' -bptf ',num2str(HP_sigma),' ',num2str(LP_sigma),' ',func_filename]);
                        out_func_data = load_nii_gz(func_filename);
                        out_func_data.img = out_func_data.img(:,:,:,temporal_mask);
                        out_func_data.hdr.dime.dim(2:5) = size(out_func_data.img);
                        out_func_data.hdr.dime.cal_max = max(out_func_data.img(:));
                        out_func_data.hdr.dime.cal_min = min(out_func_data.img(out_func_data.img~=0));
                        out_func_data.hdr.dime.datatype = 64;
                        out_func_data.hdr.dime.bitpix = 64;
                        save_nii_gz(out_func_data,func_filename);
                    else
                        for x = 1:size(func_data.img,1)
                            for y = 1:size(func_data.img,2)
                                for z = 1:size(func_data.img,3)
                                    if funcmask.img(x,y,z) == 1
                                        for bin = length(freq_bins):-1:1
                                            cosval = sum(squeeze(out_func_data.img(x,y,z,:)).*cosmult(:,bin))./cosdiv(bin);
                                            sinval = sum(squeeze(out_func_data.img(x,y,z,:)).*sinmult(:,bin))./sindiv(bin);
                                            interp_data(:,bin) = [cosval,sinval]*interpmult(:,:,bin);
                                        end
                                        interp_data = sum(interp_data,2);
                                        if exist('butter_filt','var')
                                            interp_data_padded = [zeros(15,size(interp_data,2));interp_data;zeros(15,size(interp_data,2))];
                                            temp_butter = filter(butter_filt,interp_data_padded);
                                            temp_butter(end:-1:1) = filter(butter_filt,temp_butter(end:-1:1));
                                            temp_butter = temp_butter(16:(end-15));
                                            out_func_data.img(x,y,z,:) = temp_butter(temporal_mask,:);
                                        elseif exist('filt_ind','var')
                                            padlength = rest_nextpow2_one35(subfunc_out.nTP);
                                            interp_data = [interp_data;zeros(padlength-subfunc_out.nTP,size(interp_data,2))];
                                            freq = fft(interp_data);
                                            freq(filt_ind,:) = 0;
                                            interp_data = ifft(freq);
                                            interp_data = interp_data(1:subfunc_out.nTP,:);
                                            out_func_data.img(x,y,z,:) = interp_data(temporal_mask,:);
                                        end
                                    end
                                end
                            end
                        end
                    end
                    out_func_data.hdr = func_data.hdr;
                    out_func_data.hdr.dime.dim(2:5) = size(out_func_data.img);
                    out_func_data.hdr.dime.cal_max = max(out_func_data.img(:));
                    out_func_data.hdr.dime.cal_min = min(out_func_data.img(out_func_data.img~=0));
                    out_func_data.hdr.dime.datatype = 64;
                    out_func_data.hdr.dime.bitpix = 64;
                    if strcmp(subfunc_out.saveintermed,'Yes')
                        save_nii_gz(out_func_data,func_filename);
                    end
                end
                data_accept = 1;
            else
                data_accept = 0;
            end
        else
            data_accept = 1;
        end
    else
        mean_FD = NaN;
        mean_DVARS = NaN;
        num_censored_vols = 0;
        num_censored_vols_FD = 0;
        num_censored_vols_DVARS = 0;
        data_accept = 1;
    end
    
    GCOR = NaN;
    if data_accept == 1
        if subfunc_out.globsig_calcGCOR == 1
            GCOR_data = reshape(out_func_data.img,[numel(funcmask.img),size(out_func_data.img,4)]);
            GCOR_data = GCOR_data(logical(reshape(funcmask.img,[numel(funcmask.img),1])),:);
            GCOR_data = GCOR_data - (sum(GCOR_data,2)./size(GCOR_data,2))*ones(1,size(GCOR_data,2));
            GCOR = norm(sum(GCOR_data./(sqrt(sum(GCOR_data.^2,2))*ones(1,size(GCOR_data,2))))./size(GCOR_data,1)).^2;
        end
        
        parc_data = load_nii_gz(parc_filename);
        nredu_TP = subfunc_out.nTP-num_censored_vols;
        ts = zeros(subfunc_out.nROI,nredu_TP);
        
        for ROI = 1:subfunc_out.nROI
            curr_label = subfunc_out.ROI_num_labels(ROI);
            if iscell(curr_label)
                curr_label = curr_label{:};
            end
            ROImask = zeros(size(parc_data.img));
            if ~isempty(find(parc_data.img == curr_label,1))
                ROImask(logical(parc_data.img == curr_label)) = 1;
                nvox_orig = sum(ROImask(:));
                ROImask = ROImask.*funcmask.img;
                nvox_masked = sum(ROImask(:));
                
                if nvox_masked >= subfunc_out.MinROIvox
                    if (nvox_masked/nvox_orig) >= (subfunc_out.MinpercentROIvox/100)
                        if strcmp(subfunc_out.ROI_summary_measure,'Mean')
                            for t = 1:nredu_TP
                                ts(ROI,t) = sum(sum(sum(out_func_data.img(:,:,:,t).*ROImask)))./nvox_masked;
                            end
                        elseif strcmp(subfunc_out.ROI_summary_measure,'Median')
                            for t = 1:nredu_TP
                                temp = out_func_data.img(:,:,:,t);
                                ts(ROI,t) = median(temp(logical(ROImask)));
                            end
                        elseif strcmp(subfunc_out.ROI_summary_measure,'Largest Principal Component')
                            for t = 1:nredu_TP
                                tempmeansig(t,1) = sum(sum(sum(out_func_data.img(:,:,:,t).*ROImask)))./nvox_masked;
                            end
                            ROImask = logical(reshape(ROImask,[(size(ROImask,1)*size(ROImask,2)*size(ROImask,3)),1]));
                            temp = reshape(out_func_data.img,[(size(ROImask,1)*size(ROImask,2)*size(ROImask,3)),nredu_TP]);
                            temp(~ROImask,:) = [];
                            try
                                [~,PCA_scores] = pca(temp','Algorithm','svd');
                            catch
                                [~,PCA_scores] = princomp(temp');
                            end
                            if sign(corr(tempmeansig,PCA_scores(:,1)))==(-1)
                                PCA_scores(:,1) = PCA_scores(:,1)*(-1);
                            end
                            ts(ROI,:) = PCA_scores(:,1);
                        end
                        if sum(isnan(ts(ROI,:))) > 0
                            fprintf(logfile_fid,'There is at least one NaN in the timeseries of ROI # %u for %s\n',curr_label,sub);
                        end
                    else
                        fprintf(logfile_fid,'Less than 50%% of the original voxels are in the masked version of ROI # %u for %s\n',curr_label,sub);
                        ts(ROI,:) = NaN;
                    end
                else
                    fprintf(logfile_fid,'Less than 5 voxels in masked ROI # %u for %s\n',curr_label,sub);
                    ts(ROI,:) = NaN;
                end
            else
                fprintf(logfile_fid,'No mask for ROI # %u for %s\n',curr_label,sub);
                ts(ROI,:) = NaN;
            end
        end
    else
        ts(:,:) = NaN;
    end
elseif strcmp(subfunc_out.ROIextord,'After only slice timing/motion correction/detrending')
    if ~strcmp(subfunc_out.parc_space,'Functional')
        out_parc_filename = strrep(parc_filename,'.nii','_warp2func.nii');
        if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
            if ~exist(out_parc_filename,'file')
                if any(strfind(parc_warp_filename,'.mat'))
                    eval(['!flirt -interp nearestneighbour -in ',parc_filename,' -ref ',funcmask_filename,' -out ',out_parc_filename,' -init ',parc_warp_filename,' -applyxfm']);
                elseif any(strfind(parc_warp_filename,'.nii'))
                    eval(['!applywarp --ref=',funcmask_filename,' --in=',parc_filename,' --warp=',parc_warp_filename,' --out=',out_parc_filename,' --interp=nn']);
                end
            end
        else
            if any(strfind(parc_warp_filename,'.mat'))
                eval(['!flirt -interp nearestneighbour -in ',parc_filename,' -ref ',funcmask_filename,' -out ',out_parc_filename,' -init ',parc_warp_filename,' -applyxfm']);
            elseif any(strfind(parc_warp_filename,'.nii'))
                eval(['!applywarp --ref=',funcmask_filename,' --in=',parc_filename,' --warp=',parc_warp_filename,' --out=',out_parc_filename,' --interp=nn']);
            end
        end
        parc_filename = out_parc_filename;
    end
    
    if ~strcmp(subfunc_out.WM_space,'Functional')
        out_WM_filename = strrep(WM_filename,'.nii','_warp2func.nii');
        if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
            if ~exist(out_WM_filename,'file')
                if any(strfind(WM_warp_filename,'.mat'))
                    eval(['!flirt -interp nearestneighbour -in ',WM_filename,' -ref ',funcmask_filename,' -out ',out_WM_filename,' -init ',WM_warp_filename,' -applyxfm']);
                elseif any(strfind(vent_warp_filename,'.nii'))
                    eval(['!applywarp --ref=',funcmask_filename,' --in=',WM_filename,' --warp=',WM_warp_filename,' --out=',out_WM_filename,' --interp=nn']);
                end
            end
        else
            if any(strfind(WM_warp_filename,'.mat'))
                eval(['!flirt -interp nearestneighbour -in ',WM_filename,' -ref ',funcmask_filename,' -out ',out_WM_filename,' -init ',WM_warp_filename,' -applyxfm']);
            elseif any(strfind(vent_warp_filename,'.nii'))
                eval(['!applywarp --ref=',funcmask_filename,' --in=',WM_filename,' --warp=',WM_warp_filename,' --out=',out_WM_filename,' --interp=nn']);
            end
        end
        WM_filename = out_WM_filename;
    end
    
    if ~strcmp(subfunc_out.vent_space,'Functional')
        out_vent_filename = strrep(vent_filename,'.nii','_warp2func.nii');
        if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
            if ~exist(out_vent_filename,'file')
                if any(strfind(vent_warp_filename,'.mat'))
                    eval(['!flirt -interp nearestneighbour -in ',vent_filename,' -ref ',funcmask_filename,' -out ',out_vent_filename,' -init ',vent_warp_filename,' -applyxfm']);
                elseif any(strfind(vent_warp_filename,'.nii'))
                    eval(['!applywarp --ref=',funcmask_filename,' --in=',vent_filename,' --warp=',vent_warp_filename,' --out=',out_vent_filename,' --interp=nn']);
                end
            end
        else
            if any(strfind(vent_warp_filename,'.mat'))
                eval(['!flirt -interp nearestneighbour -in ',vent_filename,' -ref ',funcmask_filename,' -out ',out_vent_filename,' -init ',vent_warp_filename,' -applyxfm']);
            elseif any(strfind(vent_warp_filename,'.nii'))
                eval(['!applywarp --ref=',funcmask_filename,' --in=',vent_filename,' --warp=',vent_warp_filename,' --out=',out_vent_filename,' --interp=nn']);
            end
        end
        vent_filename = out_vent_filename;
    end
    
    parc_data = load_nii_gz(parc_filename);
    func_data = load_nii_gz(func_filename);
    funcmask = load_nii_gz(funcmask_filename);
    func_data.img = double(func_data.img);
    funcmask.img = double(funcmask.img);
    ts = zeros(subfunc_out.nROI,subfunc_out.nTP);
    
    for ROI = 1:subfunc_out.nROI
        curr_label = subfunc_out.ROI_num_labels(ROI);
        if iscell(curr_label)
            curr_label = curr_label{:};
        end
        ROImask = zeros(size(parc_data.img));
        if ~isempty(find(parc_data.img == curr_label,1))
            ROImask(logical(parc_data.img == curr_label)) = 1;
            nvox_orig = sum(ROImask(:));
            ROImask = ROImask.*funcmask.img;
            nvox_masked = sum(ROImask(:));
            
            if nvox_masked >= subfunc_out.MinROIvox
                if (nvox_masked/nvox_orig) >= (subfunc_out.MinpercentROIvox/100)
                    if strcmp(subfunc_out.ROI_summary_measure,'Mean')
                        for t = 1:subfunc_out.nTP
                            ts(ROI,t) = sum(sum(sum(func_data.img(:,:,:,t).*ROImask)))./nvox_masked;
                        end
                    elseif strcmp(subfunc_out.ROI_summary_measure,'Median')
                        for t = 1:subfunc_out.nTP
                            temp = func_data.img(:,:,:,t);
                            ts(ROI,t) = median(temp(logical(ROImask)));
                        end
                    elseif strcmp(subfunc_out.ROI_summary_measure,'Largest Principal Component')
                        for t = 1:subfunc_out.nTP
                            tempmeansig(t,1) = sum(sum(sum(func_data.img(:,:,:,t).*ROImask)))./nvox_masked;
                        end
                        ROImask = logical(reshape(ROImask,[(size(ROImask,1)*size(ROImask,2)*size(ROImask,3)),1]));
                        temp = reshape(func_data.img,[(size(ROImask,1)*size(ROImask,2)*size(ROImask,3)),subfunc_out.nTP]);
                        temp(~ROImask,:) = [];
                        try
                            [~,PCA_scores] = pca(temp','Algorithm','svd');
                        catch
                            [~,PCA_scores] = princomp(temp');
                        end
                        if sign(corr(tempmeansig,PCA_scores(:,1)))==(-1)
                            PCA_scores(:,1) = PCA_scores(:,1)*(-1);
                        end
                        ts(ROI,:) = PCA_scores(:,1);
                    end
                    if sum(isnan(ts(ROI,:))) > 0
                        fprintf(logfile_fid,'There is at least one NaN in the timeseries of ROI # %u for %s\n',curr_label,sub);
                    end
                else
                    fprintf(logfile_fid,'Less than 50%% of the original voxels are in the masked version of ROI # %u for %s\n',curr_label,sub);
                    ts(ROI,:) = NaN;
                end
            else
                fprintf(logfile_fid,'Less than 5 voxels in masked ROI # %u for %s\n',curr_label,sub);
                ts(ROI,:) = NaN;
            end
        else
            fprintf(logfile_fid,'No mask for ROI # %u for %s\n',curr_label,sub);
            ts(ROI,:) = NaN;
        end
    end
    
    desmat = [];
    
    if subfunc_out.WMsig_part == 1
        if subfunc_out.erode_masks == 1
            WM_eroded_filename = strrep(WM_filename,img_ext,['_eroded_',num2str(subfunc_out.erode_radius),img_ext]);
            if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
                if ~exist(WM_eroded_filename,'file')
                    tempmask = load_nii_gz(WM_filename);
                    [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                    save_nii_gz(tempmask,WM_eroded_filename);
                end
            else
                tempmask = load_nii_gz(WM_filename);
                [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                save_nii_gz(tempmask,WM_eroded_filename);
            end
            WM_filename = WM_eroded_filename;
        end
        if subfunc_out.PCA_WMvent == 1
            WM_extract_filename = strrep(WM_filename,img_ext,'_EIG.txt');
            if ~exist(WM_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename,' --eig --order=5']);
            else
                vernum = 0;
                while 1
                    if exist(strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        WM_extract_filename = strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename,' --eig --order=5']);
            end
            WM_preds = dlmread(WM_extract_filename);
            WM_preds = spm_detrend(WM_preds,subfunc_out.detr_ord);
            desmat = [desmat,WM_preds];
            
            if subfunc_out.WMsigderiv_part == 1
                WM1stderiv_preds = WM_preds(3:end,:) - WM_preds(1:(end-2),:);
                WM1stderiv_preds = [WM1stderiv_preds(1,:);WM1stderiv_preds;WM1stderiv_preds(end,:)];
                WM1stderiv_preds = spm_detrend(WM1stderiv_preds,subfunc_out.detr_ord);
                desmat = [desmat,WM1stderiv_preds];
            end
        else
            WM_extract_filename = strrep(WM_filename,img_ext,'.txt');
            if ~exist(WM_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename]);
            else
                vernum = 0;
                while 1
                    if exist(strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        WM_extract_filename = strrep(WM_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',WM_extract_filename,' -m ',WM_filename]);
            end
            WM_pred = dlmread(WM_extract_filename);
            WM_pred = spm_detrend(WM_pred,subfunc_out.detr_ord);
            desmat = [desmat,WM_pred];
            
            if subfunc_out.WMsigderiv_part == 1
                WM1stderiv_pred = WM_pred(3:end) - WM_pred(1:(end-2));
                WM1stderiv_pred = [WM1stderiv_pred(1);WM1stderiv_pred;WM1stderiv_pred(end)];
                WM1stderiv_pred = spm_detrend(WM1stderiv_pred,subfunc_out.detr_ord);
                desmat = [desmat,WM1stderiv_pred];
            end
        end
    end
    
    if subfunc_out.ventsig_part == 1
        if subfunc_out.erode_masks == 1
            vent_eroded_filename = strrep(vent_filename,img_ext,['_eroded_',num2str(subfunc_out.erode_radius),img_ext]);
            if strcmp(subfunc_out.saveintermed,'Yes, do not run step if output file exists')
                if ~exist(vent_eroded_filename,'file')
                    tempmask = load_nii_gz(vent_filename);
                    [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                    save_nii_gz(tempmask,vent_eroded_filename);
                end
            else
                tempmask = load_nii_gz(vent_filename);
                [tempmask.img] = erode_mask(tempmask.img,subfunc_out.erode_radius);
                save_nii_gz(tempmask,vent_eroded_filename);
            end
            vent_filename = vent_eroded_filename;
        end
        if subfunc_out.PCA_WMvent == 1
            vent_extract_filename = strrep(vent_filename,img_ext,'_EIG.txt');
            if ~exist(vent_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename,' --eig --order=5']);
            else
                vernum = 0;
                while 1
                    if exist(strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        vent_extract_filename = strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename,' --eig --order=5']);
            end
            vent_preds = dlmread(vent_extract_filename);
            vent_preds = spm_detrend(vent_preds,subfunc_out.detr_ord);
            desmat = [desmat,vent_preds];
            
            if subfunc_out.ventsigderiv_part == 1
                vent1stderiv_preds = vent_preds(3:end,:) - vent_preds(1:(end-2),:);
                vent1stderiv_preds = [vent1stderiv_preds(1,:);vent1stderiv_preds;vent1stderiv_preds(end,:)];
                vent1stderiv_preds = spm_detrend(vent1stderiv_preds,subfunc_out.detr_ord);
                desmat = [desmat,vent1stderiv_preds];
            end
        else
            vent_extract_filename = strrep(vent_filename,img_ext,'.txt');
            if ~exist(vent_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename]);
            else
                vernum = 0;
                while 1
                    if exist(strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        vent_extract_filename = strrep(vent_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',vent_extract_filename,' -m ',vent_filename]);
            end
            vent_pred = dlmread(vent_extract_filename);
            vent_pred = spm_detrend(vent_pred,subfunc_out.detr_ord);
            desmat = [desmat,vent_pred];
            
            if subfunc_out.ventsigderiv_part == 1
                vent1stderiv_pred = vent_pred(3:end) - vent_pred(1:(end-2));
                vent1stderiv_pred = [vent1stderiv_pred(1);vent1stderiv_pred;vent1stderiv_pred(end)];
                vent1stderiv_pred = spm_detrend(vent1stderiv_pred,subfunc_out.detr_ord);
                desmat = [desmat,vent1stderiv_pred];
            end
        end
    end
    
    if subfunc_out.globsig_part == 1
        if subfunc_out.globsig_calcGNI == 1
            GNI = Rest2GNI_edit(func_filename,funcmask_filename,100);
            fprintf(logfile_fid,'Original GNI for %s = %6.3f\n',sub,GNI);
        else
            GNI = 0;
        end
        if GNI <= 3
            funcmask_extract_filename = strrep(funcmask_filename,img_ext,'.txt');
            if ~exist(funcmask_extract_filename,'file')
                eval(['!fslmeants -i ',func_filename,' -o ',funcmask_extract_filename,' -m ',funcmask_filename]);
            else
                vernum = 0;
                while 1
                    if exist(strrep(funcmask_extract_filename,'.txt',['_',num2str(vernum),'.txt']),'file')
                        vernum = vernum+1;
                    else
                        funcmask_extract_filename = strrep(funcmask_extract_filename,'.txt',['_',num2str(vernum),'.txt']);
                        break
                    end
                end
                eval(['!fslmeants -i ',func_filename,' -o ',funcmask_extract_filename,' -m ',funcmask_filename]);
            end
            global_pred = dlmread(funcmask_extract_filename);
            global_pred = spm_detrend(global_pred,subfunc_out.detr_ord);
            desmat = [desmat,global_pred];
            if subfunc_out.globsigderiv_part == 1
                global1stderiv_pred = global_pred(3:end) - global_pred(1:(end-2));
                global1stderiv_pred = [global1stderiv_pred(1);global1stderiv_pred;global1stderiv_pred(end)];
                global1stderiv_pred = spm_detrend(global1stderiv_pred,subfunc_out.detr_ord);
                desmat = [desmat,global1stderiv_pred];
            end
        end
    end
    
    if subfunc_out.motpar_part == 1
        mot_pred = dlmread(par_filename);
        for par = 1:size(mot_pred,2)
            mot_pred(:,par) = spm_detrend(mot_pred(:,par),subfunc_out.detr_ord);
        end
        desmat = [desmat,mot_pred];
        if subfunc_out.motpart1_part == 1
            mot1back_pred = [mot_pred(2:end,:);mot_pred(end,:)];
            for par = 1:size(mot_pred,2)
                mot1back_pred(:,par) = spm_detrend(mot1back_pred(:,par),subfunc_out.detr_ord);
            end
            desmat = [desmat,mot1back_pred];
        end
        if subfunc_out.motparderiv_part == 1
            mot1stderiv_pred = mot_pred(3:end,:) - mot_pred(1:(end-2),:);
            mot1stderiv_pred = [mot1stderiv_pred(1,:);mot1stderiv_pred;mot1stderiv_pred(end,:)];
            for par = 1:size(mot_pred,2)
                mot1stderiv_pred(:,par) = spm_detrend(mot1stderiv_pred(:,par),subfunc_out.detr_ord);
            end
            desmat = [desmat,mot1stderiv_pred];
        end
        if subfunc_out.motparsqr_part == 1
            motsquare_pred = mot_pred.^2;
            for par = 1:size(mot_pred,2)
                motsquare_pred(:,par) = spm_detrend(motsquare_pred(:,par),subfunc_out.detr_ord);
            end
            desmat = [desmat,motsquare_pred];
            if subfunc_out.motpart1_part == 1
                mot1backsquare_pred = mot1back_pred.^2;
                for par = 1:size(mot_pred,2)
                    mot1backsquare_pred(:,par) = spm_detrend(mot1backsquare_pred(:,par),subfunc_out.detr_ord);
                end
                desmat = [desmat,mot1backsquare_pred];
            end
        end
    end
    
    if subfunc_out.BP == 1 && (subfunc_out.WMsig_part == 1 || subfunc_out.ventsig_part == 1 || (subfunc_out.globsig_part == 1 && GNI <= 3) || subfunc_out.motpar_part == 1)
        desmat = [ones(size(desmat,1),1),desmat];
        if rank(desmat) ~= min(size(desmat))
            error(['Initial desmat to partial out unwanted variance for ',sub,' is not full rank']);
        end
        Q = desmat*pinv(full(desmat));
        for ROI = 1:subfunc_out.nROI
            if sum(isnan(ts(ROI,:))) == 0
                tempdata = squeeze(ts(ROI,:))-(Q*squeeze(ts(ROI,:)));
                tempdata = tempdata-mean(tempdata);
                if exist('butter_filt','var')
                    tempdata_padded = [zeros(15,size(tempdata,2));tempdata;zeros(15,size(tempdata,2))];
                    temp_butter = filter(butter_filt,tempdata_padded);
                    temp_butter(end:-1:1) = filter(butter_filt,temp_butter(end:-1:1));
                    ts(ROI,:) = temp_butter(16:(end-15));
                elseif exist('filt_ind','var')
                    padlength = rest_nextpow2_one35(subfunc_out.nTP);
                    tempdata = [tempdata;zeros(padlength-subfunc_out.nTP,size(tempdata,2))];
                    freq = fft(tempdata);
                    freq(filt_ind,:) = 0;
                    tempdata = ifft(freq);
                    ts(ROI,:) = tempdata(1:subfunc_out.nTP,:);
                end
            end
        end
    elseif subfunc_out.BP == 1
        for ROI = 1:subfunc_out.nROI
            if sum(isnan(ts(ROI,:))) == 0
                tempdata = squeeze(ts(ROI,:));
                tempdata = tempdata-mean(tempdata);
                if exist('butter_filt','var')
                    tempdata_padded = [zeros(15,size(tempdata,2));tempdata;zeros(15,size(tempdata,2))];
                    temp_butter = filter(butter_filt,tempdata_padded);
                    temp_butter(end:-1:1) = filter(butter_filt,temp_butter(end:-1:1));
                    ts(ROI,:) = temp_butter(16:(end-15));
                else
                    padlength = rest_nextpow2_one35(subfunc_out.nTP);
                    tempdata = [tempdata;zeros(padlength-subfunc_out.nTP,size(tempdata,2))];
                    freq = fft(tempdata);
                    freq(filt_ind,:) = 0;
                    tempdata = ifft(freq);
                    ts(ROI,:) = tempdata(1:subfunc_out.nTP,:);
                end
            end
        end
    elseif subfunc_out.WMsig_part == 1 || subfunc_out.ventsig_part == 1 || (subfunc_out.globsig_part == 1 && GNI <= 3) || subfunc_out.motpar_part == 1
        desmat = [ones(size(desmat,1),1),desmat];
        if rank(desmat) ~= min(size(desmat))
            error(['Initial desmat to partial out unwanted variance for ' sub ' is not full rank']);
        end
        Q = desmat*pinv(full(desmat));
        for ROI = 1:subfunc_out.nROI
            if sum(isnan(ts(ROI,:))) == 0
                ts(ROI,:) = squeeze(ts(ROI,:))-(Q*squeeze(ts(ROI,:)));
            end
        end
    end
    
    mean_FD = NaN;
    mean_DVARS = NaN;
    num_censored_vols = 0;
    num_censored_vols_FD = 0;
    num_censored_vols_DVARS = 0;
end



%%%% Calculate the "ideal" bandpass filter
function [filt_ind] = calc_IdealFilter(nTP,TR,LowCutoff_HighPass,HighCutoff_LowPass)
% Created using the REST toolbox's rest_IdealFilter.m function

if LowCutoff_HighPass<=0
    LowCutoff_HighPass = 0;
elseif HighCutoff_LowPass<=0
    HighCutoff_LowPass = 0;
end
    
sampleFreq 	 = 1/TR;
paddedLength = rest_nextpow2_one35(nTP); %2^nextpow2(sampleLength);

% Get the frequency index
if (LowCutoff_HighPass >= sampleFreq/2) % All high stop
    idxLowCutoff_HighPass = paddedLength/2+1;
else % high pass, such as freq > 0.01 Hz
    idxLowCutoff_HighPass = ceil(LowCutoff_HighPass * paddedLength * TR+1);
end

if (HighCutoff_LowPass>=sampleFreq/2)||(HighCutoff_LowPass==0) % All low pass
    idxHighCutoff_LowPass = paddedLength/2+1;
else % Low pass, such as freq < 0.1 Hz
    idxHighCutoff_LowPass = fix(HighCutoff_LowPass * paddedLength * TR+1);
end

FrequencyMask = zeros(paddedLength,1);
FrequencyMask(idxLowCutoff_HighPass:idxHighCutoff_LowPass,1) = 1;
FrequencyMask(paddedLength-idxLowCutoff_HighPass+2:-1:paddedLength-idxHighCutoff_LowPass+2,1) = 1;

filt_ind = find(FrequencyMask==0);



%%%% Detrend a timeseries for a defined order
function detrend_image(func_filename,funcmask_filename,detr_poly_ord,varargin)

orig = load_nii_gz(func_filename);
mask = load_nii_gz(funcmask_filename);
if ~isa(orig.img,'double')
    orig.img = double(orig.img);
end

if ~isempty(varargin)
    temporal_mask = varargin{1};
    out_filename = strrep(func_filename,'.nii',['_motcensor_detr',num2str(detr_poly_ord),'.nii']);
else
    temporal_mask = true(size(orig.img,4),1);
    out_filename = strrep(func_filename,'.nii',['_detr',num2str(detr_poly_ord),'.nii']);
end

orig.img = orig.img(:,:,:,temporal_mask);
outdata.img = zeros(size(orig.img));

for x = 1:size(orig.img,1)
    for y = 1:size(orig.img,2)
        for z = 1:size(orig.img,3)
            if mask.img(x,y,z) == 1
                outdata.img(x,y,z,:) = spm_detrend(squeeze(orig.img(x,y,z,:)),detr_poly_ord);
            end
        end
    end
end

outdata.hdr = orig.hdr;
outdata.hdr.dime.dim(2:5) = size(outdata.img);
outdata.hdr.dime.datatype = 64;
outdata.hdr.dime.bitpix = 64;
outdata.hdr.dime.cal_max = max(outdata.img(:));
outdata.hdr.dime.cal_min = min(outdata.img(outdata.img~=0));
save_nii_gz(outdata,out_filename);