function plot_fMRI_BSNIP(S_hat2,id_B, tr,mask)
%% Load results and true voxel locations (indices)
% Matlab function for plotting fMRI maps from the data. 
% Dependencies: SPM needed 

%Input:
% S_hat2 = fMRI data matrix (N*T) after flattaning, where N  = number of obervations(Subjects or timepoint or components), 
% T = Samples (voxels)
% id_b = ID of the maps, 1:N if we need to plot all the maps
% tr = Thresold for voxels
% mask = mask information needed for plotting

%Output:
% fMRI maps based on selected id_B



fileName=['test']; %filename for saving the figures
numOfCV=3; %total number of estimated components
%%
DX         = size(mask,1);  
DY         = size(mask,2); 
DZ         = size(mask,3); 
maskindx = mask(:);
[ind_f,col,v]  = find(maskindx);
%% Select components of interest and z-score threshold
% Sf=-sR123.signal{3,2}(:,1:48546); %{dataset,run}(voxels) %% 1:48546___48547:97092___97093:145638

if ~exist('monthresh') %set z-score threshold (higher threshold means fewer plotted voxels
    monthresh = tr
end

fprintf('Generating Data Display ... \n');
if ~exist('sele') %set the components to be plotted
    sele = [id_B];
end
store = ind_f;
ind_f = ind_f(1:size(S_hat2,2),1);
numOfCV = size(S_hat2,1); 
%% Create 4-D fMRI data and read in anatomical map
dat3f = zeros(numOfCV,DX, DY, DZ); 
dat3f(:,ind_f) = S_hat2(:,(1:length(ind_f)));
dat3f = reshape(dat3f,numOfCV,DX, DY, DZ);

fprintf('- Reading in Anatomical Data\n');
%anat = spm_read_vols(spm_vol('C:\Users\shume\Documents\MATLAB\MATLAB\GroupICATv4.0b\icatb\icatb_templates\nsingle_subj_T1_2_2_4.nii'));
anat = spm_read_vols(spm_vol('C:\Users\shume\Documents\MATLAB\Plotting Code\nsingle_subj_T1_2_2_3.img')); % change it based on the user
%anatFile   = 'C:\Users\shume\Documents\MATLAB\random_chunying\rch2bet.nii';
%anat      = spm_read_vols(icatb_spm_vol(anatFile));
% select slices to show ("1:1:46" will show all 46 slices)
z_choice = [2:5:46];
%%
%for individual images
for i=1:length(sele) %1:numOfCV
    tmp(:,:,:)=dat3f(abs(sele(i)),:,:,[8:5:52]);
    m=(sign(sele(i))*S_hat2/stdN(S_hat2)); % use if you want to normalize based on all voxels (as opposed to only the plotted ones (see line 43))

    figure;
    xlabel('fMRI','FontSize', 12);
    colorbar
    make_composite(sign(sele(i))*tmp/stdN(tmp),anat(:,:,z_choice),monthresh); %plot component
    clear tmp
end
%% save files as PDF
% 
% print('-depsc', fileName);
% hgsave(fileName);
% save(fileName);
% eps2pdf([fileName, '.eps'], 'C:\Program Files\gs\gs9.14\bin\gswin64.exe');   %for windows, the path of gs should be different
end