function extract_roi_BOLD
flag='warf'; % Flag to select preprocessed fMRI data. 
extension='.nii';% Flag to select specific fMRI data format. 
data.fMRI_data_complete=[];
pre_post_paradigm=load('D:\Mohit\nicotine\spm_analsis\condition_pre_post.mat');% format similar to the SPM multiple condition  mat file. 
for ii=1:length(pre_post_paradigm.onsets{1,1})
    pre_base_idx(ii,:)=pre_post_paradigm.onsets{1,1}(ii):pre_post_paradigm.onsets{1,1}(ii)+pre_post_paradigm.durations{1,1}(ii)-1;
end
for ii=1:length(pre_post_paradigm.onsets{1,2})
    pre_reg_idx(ii,:)=pre_post_paradigm.onsets{1,2}(ii):pre_post_paradigm.onsets{1,2}(ii)+pre_post_paradigm.durations{1,1}(ii)-1;
end
training_paradigm=load('D:\Mohit\nicotine\spm_analsis\condition_training.mat');% format similar to the SPM multiple condition  mat file. 
for ii=1:length(training_paradigm.onsets{1,1})
    train_base_idx(ii,:)=pre_post_paradigm.onsets{1,1}(ii):pre_post_paradigm.onsets{1,1}(ii)+pre_post_paradigm.durations{1,1}(ii)-1;
end
for ii=1:length(training_paradigm.onsets{1,2})
    train_reg_idx(ii,:)=pre_post_paradigm.onsets{1,2}(ii):pre_post_paradigm.onsets{1,2}(ii)+pre_post_paradigm.durations{1,1}(ii)-1;
end
mni=[-35, 16, -8; 35 16 -8; 16,-24, 62];% MNI for the ROIs. Each row represents an ROI. 
rootpath=uigetdir(pwd,'Select the root data directory');% directory contianing the complete data. 
file_idx=0;
column=6;space=2;row=9;% for the visualization of the fMRI data
if sum(rootpath==0)>0
    disp('please select the root data directory')
    return
end
Groups=dir(rootpath);subj_idx=0;
for group=1:Groups
    for subj=3:length(subjects)
        if isdir([rootpath,filesep,subjects(subj).name])
            
            sessions= dir([rootpath,filesep,subjects(subj).name]);
            
            for sess=3:length(sessions)
                if isdir([rootpath,filesep,subjects(subj).name,filesep,sessions(sess).name]) && strcmp(sessions(sess).name(1),'t')
                    if strcmp(sessions(sess).name,'t1')
                        subj_idx=subj_idx+1;
                        total_runs=1;
                    end
                    
                    runs=dir([rootpath,filesep,subjects(subj).name,filesep,sessions(sess).name]);
                    for run=3:length(runs)
                        if isdir([rootpath,filesep,subjects(subj).name,filesep,sessions(sess).name,filesep,runs(run).name])&& strcmp(runs(run).name(1),'s')
                            files=dir([rootpath,filesep,subjects(subj).name,filesep,sessions(sess).name,filesep,runs(run).name,filesep,flag,'*',extension]);
                            if length(files)>0
                                disp(['processing: ' ,subjects(subj).name,' Session: ',sessions(sess).name,' Run: ',runs(run).name])
                                if subj_idx==1 && strcmp(sessions(sess).name,'t1') && strcmp(runs(run).name,'s1') % T1 = session 1 and S1= first run.
                                    aan=load_untouch_nii(strcat(rootpath,filesep,subjects(subj).name,filesep,sessions(sess).name,filesep,runs(run).name,filesep,files(1).name));
                                    [coordinate_voxels, mask_vox_ss] = mni2cor(mni , [aan.hdr.hist.srow_x;aan.hdr.hist.srow_y;aan.hdr.hist.srow_z; 0 0 0 1 ],[3 3 3], 'Sphere',aan.hdr.dime.dim(2:4));
                                    
                                end
                                iin2=0; x=2; y=2 ;
                                fmri_imag(2*(row+1)+row*aan.hdr.dime.dim(3),2*(column+1)+column*aan.hdr.dime.dim(2))=0;
                                mask_imag(2*(row+1)+row*aan.hdr.dime.dim(3),2*(column+1)+column*aan.hdr.dime.dim(2))=0;
                                for in2=1:aan.hdr.dime.dim(4)
                                    iin2=iin2+1;
                                    fmri_imag(x+1:x+aan.hdr.dime.dim(3),y+1:y+aan.hdr.dime.dim(2))= rot90(aan.img(:,:,in2),1);
                                    mask_imag(x+1:x+aan.hdr.dime.dim(3),y+1:y+aan.hdr.dime.dim(2))= rot90(coordinate_voxels(:,:,in2),1);
                                    %         figure, imshow(handles.b.pre.mask_imag,[0 3])
                                    if mod(iin2,column)==0
                                        y= space;x=floor(iin2/ column)*aan.hdr.dime.dim(3)+space*(floor(iin2/row)+1);
                                    else
                                        y=mod(iin2,column)*aan.hdr.dime.dim(2)+ space*(mod(iin2,column)+1);
                                    end
                                end
                                h=figure;
                                axes_h=axis;
                                fmri_imag=fmri_imag./max(max(fmri_imag));
                                hh=imshow(fmri_imag,[0 1]);
                                hold on
                                mask_rgb=ind2rgb(mask_imag+1,[0 0 0; 1 0 0; 0 1 0; 0 0 1; 1 1 0; 0 1 1]);
                                color_map_mask=[0 0 0; 1 0 0; 0 1 0; 0 0 1; 1 1 0; 0 1 1];
                                h=image(mask_rgb );
                                %                             colormap( );
                                a = 0.8;
                                set(h,'AlphaData',(mask_imag).*a)
                                %   colormap([0 0 0; 1 0 0; 0 1 0; 0 0 1; 1 1 0; 0 1 1])
                                for K = 1 :size(mni,1)
                                    %  hidden_h(K) = surf(uint8(K-[1 1;1 1]), 'edgecolor', 'none');
                                    if K~= size(mni,1)
                                        text_legend{1,K}=['ROI ', num2str(K)];
                                    else
                                        text_legend{1,K}='ROI Ref';
                                    end
                                end
                                L = line(ones(size(mni,1)),ones(size(mni,1)),'LineWidth',3);
                                MM=color_map_mask(2:size(mni,1)+1,:);
                                set(L,{'color'},mat2cell(MM ,ones(1,size(mni,1)),3));
                                legend(text_legend,'Location','southoutside','Orientation','horizontal');
                                
                                for file=1:length(files)
                                    clear aan
                                    aan=load_untouch_nii(strcat(rootpath,filesep,subjects(subj).name,filesep,sessions(sess).name,filesep,runs(run).name,filesep,files(file).name));
                                    data.fMRI_data_complete(:,file)=aan.img(:);
                                    roi_data= sum(mask_vox_ss.*repmat(data.fMRI_data_complete(:,file),1,size(mni,1)))./sum(mask_vox_ss);
                                    data.ROI_fMRI_data(:,file)=roi_data';
                                    clear roi_data
                                end
                                %   data.label(length(data.label)+1:length(data.label)+length(files))=0;
                                if strcmp(sessions(sess).name,'t1') || strcmp(sessions(sess).name,'t5')|| strcmp(sessions(sess).name,'t4')
                                    for ii=1:size(pre_reg_idx,1)
                                        baseline=mean(data.ROI_fMRI_data(:,pre_base_idx(ii,:)),2);
                                        regulation=mean(data.ROI_fMRI_data(:,pre_reg_idx(ii,:)),2);
                                        data.per_block(:,ii)=(regulation-baseline)./baseline;
                                    end
                                elseif strcmp(sessions(sess).name,'t2') || strcmp(sessions(sess).name,'t3')
                                    for ii=1:size(train_reg_idx,1)
                                        baseline=mean(data.ROI_fMRI_data(:,train_base_idx(ii,:)),2);                                        
                                        regulation=mean(data.ROI_fMRI_data(:,train_reg_idx(ii,:)),2);
                                        data.per_block(:,ii)=(regulation-baseline)./baseline;                                       
                                    end
                                end
                                % file_idx=file_idx+length(files);
                                save(strcat(rootpath,filesep,subjects(subj).name,filesep,sessions(sess).name,filesep,subjects(subj).name,'_',sessions(sess).name,'_',runs(run).name,'_analyzed.mat'),'data');
                            else
                                disp(['No file found specific to flag (' flag,'):', subjects(subj),' Session: ',sessions(sess).name,' Run: ',runs(run).name])
                                disp('Continuing with the rest of the data...');
                            end
                            if subj_idx==1 && strcmp(sessions(sess).name,'t1') && strcmp(runs(run).name,'s1')
                                data_session.data_fmri_roi=data.ROI_fMRI_data;
                                data_session.data_fmri_complete=data.fMRI_data_complete;
                                data_session.data_fmri_label=data.label;
                                data_session.per_block(:,total_runs,subj_idx)=mean(data.per_block,2);
                                total_runs=total_runs+1;
                            else
                                aa=size(data_session.data_fmri_roi,2);ab=size(data.ROI_fMRI_data,2);
                                data_session.data_fmri_roi(:,aa+1:aa+ab)=data.ROI_fMRI_data;
                                clear aa ab
                                aa=size(data_session.data_fmri_complete,2);ab=size(data.fMRI_data_complete,2);
                                data_session.data_fmri_complete(:,aa+1:aa+ab)=data.fMRI_data_complete;
                                clear aa ab
                                aa=size(data_session.data_fmri_label,2);ab=size(data.label,2);
                                data_session.data_fmri_label(:,aa+1:aa+ab)=data.label;
                                data_session.per_block(:,total_runs,subj_idx)=mean(data.per_block,2);
                                total_runs=total_runs+1;
                            end
                        end
                    end
                    
                end
                
            end
            
        end
        
    end
end
save(strcat(rootpath,filesep,'total_data.mat'),'data_session');
function [coordinate_voxels, mask_vox_ss] = mni2cor(mni , T,roi_size, type,dim)
% function coordinate = mni2cor(mni, T)
% convert mni coordinate to matrix coordinate
%
% mni: a Nx3 matrix of mni coordinate
% T: (optional) transform matrix
% coordinate is the returned coordinate in matrix
%
% caution: if T is not specified, we use:
% T = ...
%     [-4     0     0    84;...
%      0     4     0  -116;...
%      0     0     4   -56;...
%      0     0     0     1];
%
%

% if isempty(mni )
%     coordinate_voxels = [];
%     return;
%  end

% if isempty(T)
%     T = ...
%         [-4     0     0    84;...
%          0     4     0  -116;...
%          0     0     4   -56;...
%          0     0     0     1];
% end
% mni=[0 0 0; 30 30 0]
nn=floor(roi_size(1)/2);
nn1=floor(roi_size(2)/2);
nn2=floor(roi_size(3)/2);
coordinate =[mni(:,1) mni(:,2) mni(:,3) ones(size(mni,1),1)]*(inv(T))';
coordinate(:,4) = [];
coordinate = round(coordinate);
coordinate_voxels=zeros( dim(1)* dim(2)* dim(3),1);
mask_vox_ss= zeros( dim(1)* dim(2)* dim(3),size((mni),1));
if strcmp(type,'Sphere')
    for iin=1:size(coordinate,1)
        [xx ,yy, zz] = meshgrid(1:dim(1),1:dim(2),1:dim(3));
        vv = sqrt((xx-coordinate(iin,1)).^2+(yy-coordinate(iin,2)).^2+(zz-coordinate(iin,3)).^2)<=(roi_size(1)/2);
        for ini=1:dim(3)
            voxels(:,:,ini)=rot90(vv(:,:,ini),1);
            %           if sum(sum(vv(:,:,ini)))>0
            %               ini
            %               figure, imshow(vv(:,:,ini),[0 4])
            %               figure, imshow(voxels(:,:,ini),[0 4])
            %           end
        end
        coordinate_voxels (voxels(:)>0)=iin;
        mask_vox_ss(:,iin)=(voxels(:)>0);
    end
elseif strcmp(type,'Cube')
    for iin=1:size(coordinate,1)
        voxels(dim(1),dim(2),dim(3))=0;
        xx= coordinate(1,1)-nn:1:coordinate(1,1)+nn ;
        yy= coordinate(1,2)-nn1:1:coordinate(1,2)+nn1 ;
        zz= coordinate(1,3)-nn2:1:coordinate(1,3)+nn2 ;
        comb_voxel=combvec(xx,yy,zz);
        voxels(comb_voxel(1,:),comb_voxel(2,:),comb_voxel(3,:))=1;
        coordinate_voxels (voxels >0)=iin;
        mask_vox_ss(:,iin)=(voxels(:)>0);
    end
end
coordinate_voxels=reshape(coordinate_voxels,[dim(1),dim(2),dim(3)]);