function ck_GetMRI_pRF(monkeys,models,output)

if nargin < 2
    fprintf('ERROR: Not enough arguments specified\n');
    return
end
clc;

%% Add nifti reading toolbox ==============================================
% location of this script
homefld = pwd;

% Get the root folder path for the hared data/code
cd ../../../
SHARED_ROOT_FLD = pwd;
cd(homefld)
tool_basepath = fullfile(SHARED_ROOT_FLD,'Toolboxes');
addpath(genpath(fullfile(tool_basepath, 'NIfTI')));

%% Set fMRI paths =========================================================
fitres_path = fullfile(SHARED_ROOT_FLD,'FitResults','MRI');
bids_path = fullfile(SHARED_ROOT_FLD,'Preprocessed_data','manual-masks');

%% collect the fit results ================================================
fprintf('Collecting MRI retinotopic maps...\n');

for m = 1:length(monkeys)
    fprintf(['Processing monkey: ' monkeys{m} '\n']);
    cMonkey = monkeys{m}; cMonkey(1)=upper(cMonkey(1));
    
    R(m).monkey = monkeys{m};
    R(m).mode = 'MRI';
    
    for mm = 1:length(models)
        fprintf(['Model: ' models{mm} '\n']);
        
        RES = load(fullfile(fitres_path,monkeys{m},models{mm},...
            ['pRF_Sess-' models{mm}]),'result');
        
        R(m).model(mm).prfmodel = models{mm};
        R(m).model(mm).hrf = RES.result.options.hrf;
        R(m).model(mm).voldim = size(RES.result.R2,4);
                
        for i = 1:size(RES.result.R2,4)
            R(m).model(mm).ang(i,:) = reshape(RES.result.ang(:,:,:,i),...
                [numel(RES.result.ang(:,:,:,i)),1]);
            R(m).model(mm).ecc(i,:) = reshape(RES.result.ecc(:,:,:,i),...
                [numel(RES.result.ang(:,:,:,i)),1]);
            R(m).model(mm).rfs(i,:) = reshape(RES.result.rfsize(:,:,:,i),...
                [numel(RES.result.ang(:,:,:,i)),1]);
            R(m).model(mm).gain(i,:) = reshape(RES.result.gain(:,:,:,i),...
                [numel(RES.result.ang(:,:,:,i)),1]);
            R(m).model(mm).R2(i,:) = reshape(RES.result.R2(:,:,:,i),...
                [numel(RES.result.ang(:,:,:,i)),1]);
            R(m).model(mm).fwhm(i,:) = reshape(RES.result.fwhm(:,:,:,i),...
                [numel(RES.result.ang(:,:,:,i)),1]);
            if models{mm}(1:3) == 'dog'
                R(m).model(mm).sdratio(i,:) = reshape(RES.result.sdratio(:,:,:,i),...
                    [numel(RES.result.ang(:,:,:,i)),1]);
                R(m).model(mm).normamp(i,:) = reshape(RES.result.normamp(:,:,:,i),...
                    [numel(RES.result.ang(:,:,:,i)),1]);
            else
                R(m).model(mm).sdratio = [];
                R(m).model(mm).normamp = [];
            end
            if models{mm}(1:3) == 'css'
                R(m).model(mm).expt(i,:) = reshape(RES.result.expt(:,:,:,i),...
                    [numel(RES.result.ang(:,:,:,i)),1]);
            else
                R(m).model(mm).expt = [];
            end
        end
    end
    clear RES
    
    %% Get ROI info =======================================================
    fprintf(['Processing atlas for monkey: ' monkeys{m} '\n']);
    
    labelfile = fullfile(bids_path,['sub-' monkeys{m}],'atlas','D99_labeltable.txt');
    D99 = ck_GetAtlasTable(labelfile);
    
    atlasfile = fullfile(bids_path,['sub-' monkeys{m}],'atlas',...
        ['D99_in_' cMonkey '_adj_inFunc.nii']);
    brainmaskfile = fullfile(bids_path,['sub-' monkeys{m}],'func',...
        ['sub-' monkeys{m} '_ref_func_mask_res-1x1x1.nii']);
    V1_maskfile = fullfile(bids_path,['sub-' monkeys{m}],'atlas',...
        'ROI_ELECTRODES', 'ELECTRODES_LH_V1.nii');
    V4_maskfile = fullfile(bids_path,['sub-' monkeys{m}],'atlas',...
        'ROI_ELECTRODES', 'ELECTRODES_LH_V4.nii');
    
    BRAINMASK = load_nii(brainmaskfile);
    R(m).BRAIN = reshape(BRAINMASK.img, [numel(BRAINMASK.img),1]);
    
    ROI = load_nii(atlasfile);
    R(m).ROI = reshape(ROI.img, [numel(ROI.img),1]);
    
    % additional ROI for approximate electrode locations
    ROI_V1 = load_nii(V1_maskfile);
    R(m).ELECTR_V1 = reshape(ROI_V1.img, [numel(ROI_V1.img),1]);
    ROI_V4 = load_nii(V4_maskfile);
    R(m).ELECTR_V4 = reshape(ROI_V4.img, [numel(ROI_V4.img),1]);
    
end

%% Save the combined results ==============================================
fprintf('Saving the combined MRI result-file\n');
[~,~,~] = mkdir(fullfile(fitres_path,'Combined'));
save(fullfile(fitres_path,'Combined',output),'R','D99')