Home > . > prt_preproc.m

prt_preproc

PURPOSE ^

Function to preprocess the images, by loading each one of them (or the

SYNOPSIS ^

function prt_preproc(fname)

DESCRIPTION ^

 Function to preprocess the images, by loading each one of them (or the
 ones corresponding to the selected scans when a design was specified),
 applying the masks on them and, if asked, detrend along each voxel along
 the time series.
 
 INPUT: 
   fname   filename and path to PRT.mat

 OUTPUT:
   results are saved on disk.
_______________________________________________________________________
 Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function prt_preproc(fname)
0002 
0003 % Function to preprocess the images, by loading each one of them (or the
0004 % ones corresponding to the selected scans when a design was specified),
0005 % applying the masks on them and, if asked, detrend along each voxel along
0006 % the time series.
0007 %
0008 % INPUT:
0009 %   fname   filename and path to PRT.mat
0010 %
0011 % OUTPUT:
0012 %   results are saved on disk.
0013 %_______________________________________________________________________
0014 % Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory
0015 
0016 % Written by J.Rondina
0017 % $Id: prt_preproc.m 201 2011-10-25 14:44:41Z cphillip $
0018 
0019 %--------------------------------------------------------------------------
0020 % Load PRT.mat
0021 %--------------------------------------------------------------------------
0022 load(char(fname));    % Load data structure (PRT.mat)
0023 prt_dir=regexprep(char(fname),'PRT.mat', '');  % Get PRT.mat directory
0024 
0025 
0026 %-------------------------------------------------------------------------
0027 % Resizing masks
0028 %-------------------------------------------------------------------------
0029 
0030 % Get numbers of groups, modalities and subjects from PRT.mat
0031 n_groups = length(PRT.group);
0032 n_modalities = length(PRT.masks);
0033 tot_subj=sum(length([PRT.group(:).subject]));
0034 
0035 % Get the first scan from each modality as a sample to resize the
0036 % respective masks, if it is necessary
0037 sample_img=cell(n_modalities,1);
0038 sample_name=cell(n_modalities,1);
0039 mod_name=cell(n_modalities,1);
0040 for m=1:n_modalities
0041     sample_img{m} = nifti(char(PRT.group(1).subject(1).modality(m).scans(1,:)));
0042     sample_name{m} = char(PRT.group(1).subject(1).modality(m).scans(1,:));
0043     mod_name{m}=PRT.group(1).subject(1).modality(m).mod_name;
0044 end
0045 
0046 mnii=cell(n_modalities,1);
0047 mask=cell(n_modalities,1);
0048 for m=1:n_modalities
0049     maskmod=PRT.masks(m).mod_name;
0050     maskname =  regexprep(char(PRT.masks(strcmpi(mod_name, maskmod)).fname),',1','');   % get mask name
0051     mnii{m} = nifti(maskname);
0052     % Name of the resized mask to be saved
0053     newmask_name = [prt_dir 'updated_mask_m' PRT.group(1).subject(1).modality(m).mod_name '.img'];  
0054     % If dimensions of mask and scan don't match
0055     if size(mnii{m}.dat(:,:,:,1)) ~= size(sample_img{m}.dat(:,:,:,1)) 
0056        % Use spm_imcalc function for mask resizing
0057        img1_calc = spm_vol(sample_name{m}); 
0058        img2_calc = spm_vol(maskname);
0059        outputfile = img1_calc;
0060        outputfile.fname = newmask_name;
0061        spm_imcalc([img1_calc img2_calc],outputfile,'i2.*(i1>0)');
0062        mnii{m} = nifti(newmask_name); 
0063     end
0064     % Reshape masks (resized or not)for all modalities as vectors and
0065     % stores in the structure mask
0066     mask{m} = double(reshape(mnii{m}.dat,1,size(mnii{m}.dat,1)*size(mnii{m}.dat,2)*size(mnii{m}.dat,3)));
0067 end
0068 
0069 %-------------------------------------------------------------------------
0070 % Loading images, detrending, saving
0071 %-------------------------------------------------------------------------
0072 
0073 h = waitbar(0,'Please wait while images are pre-processed (details provided in the main window)');
0074 step=0;
0075 for g = 1:n_groups 
0076     disp(['Loading group ' int2str(g)]); 
0077     n_subj = length(PRT.group(g).subject);
0078 
0079     for s = 1:n_subj
0080         disp(['Loading subject ' int2str(s)]);
0081         n_mod = length(PRT.group(g).subject(s).modality);
0082         
0083         for m = 1:n_mod
0084             disp(['Loading modality ' int2str(m)]);                 
0085             n_scans = size(PRT.group(g).subject(s).modality(m).scans,1);
0086             m2=find(strcmpi(mod_name,PRT.group(g).subject(s).modality(m).mod_name));
0087             dim1 = size(sample_img{m2}.dat,1);
0088             dim2 = size(sample_img{m2}.dat,2);
0089             dim3 = size(sample_img{m2}.dat,3);
0090             n_voxels_scan = dim1*dim2*dim3; % get total number of voxels in the scan
0091             
0092             % Loading images
0093             img_allscans=zeros(n_scans,n_voxels_scan);
0094             for sc = 1:n_scans              
0095                 nii = nifti(PRT.group(g).subject(s).modality(m).scans(sc,:));
0096                 if size(nii.dat) ~= size(sample_img{m2}.dat)
0097                     % Check if some scan have different dimensionality
0098                     % (e.g. corrupted image)
0099                     disp (['Error - Scan ' int2str(sc) 'Modality ' int2str(m) 'Subject ' int2str(s) 'Group ' int2str(g) ' has wrong dimensions']);
0100                     break;
0101                 else
0102                     img = double(reshape(nii.dat,1,n_voxels_scan)); 
0103                     img(find(isnan(img))) = 0; % Replace NaNs for zeros
0104                     img_allscans(sc,:) = img.*mask{m2}; % Apply mask to scan and store it in a matrix formed for all the scans from the same subject and modality
0105                 end
0106             end
0107             
0108             step = step + 1;
0109             waitbar(step / (tot_subj*n_modalities*3));
0110           
0111             % Detrending
0112             if (PRT.group(g).subject(s).modality(m).detrend)
0113                 disp ('Detrending ....');
0114                 n_scans = length(PRT.group(g).subject(s).modality(m).scans);
0115                 mask_indices = find(mask{m2}>0);
0116                 n_voxels_mask = length(mask_indices); 
0117                 for v = 1:n_voxels_mask  % Detrend only voxels inside the mask
0118                     timeserie = img_allscans(:,mask_indices(v));
0119                     aux = detrend(timeserie);
0120                     img_allscans(:,mask_indices(v)) = aux;
0121                 end
0122             end
0123             
0124             step = step + 1;
0125             waitbar(step / (tot_subj*n_modalities*3));
0126             
0127             %Saving pre-processed images
0128             disp ('Saving pre-processed images ...');
0129             % Saving timeseries already detrended (in this case, the examples are extracted from the timeseries
0130             if isa(PRT.group(g).subject(s).modality(m).design,'struct') && PRT.group(g).subject(s).modality(m).detrend
0131                 % i.e. if there is a design and if the detrend was done in preprocessing
0132                 n_cond = length(PRT.group(g).subject(s).modality(m).design.conds);
0133                 for c = 1:n_cond 
0134                     filename = [prt_dir, prt_get_filename([g,s,m,c]),'.img']; %Get output file name
0135                     examples_list = PRT.group(g).subject(s).modality(m).design.conds(c).scans;
0136                     % Build 4d image with the extracted examples
0137                     img4d = file_array(filename,[dim1,dim2,dim3,length(examples_list)],'float64-le',0,1,0);
0138                     for i = 1:length(examples_list)
0139                         if ~(examples_list(i)>n_scans)
0140                             img1d = img_allscans(examples_list(i),:);
0141                             img3d = reshape(img1d,dim1,dim2,dim3);
0142                             img4d(:,:,:,i) = img3d;
0143                         else  %should not happen since checked in the data and design
0144                             disp('Warning - design exceeds timeseries')
0145                             disp('Corresponding scans were discarded')
0146                         end
0147                     end                
0148                     No         = sample_img{1}; % copy header
0149                     No.dat     = img4d;         % change file_array
0150                     No.descrip = 'Pronto data';
0151                     No.cal     = [0 1000];
0152                     create(No);                 % write header
0153                 end
0154             else % if design == 0 (i.e. structural) or if detrend == 0 (i.e., timeseries not detrended)
0155                 filename = [prt_dir, prt_get_filename([g,s,m]),'.img'];
0156                 img4d = file_array(filename,[dim1,dim2,dim3,n_scans],'FLOAT64-LE',0,1,0);
0157                 for i = 1:n_scans                 
0158                     img3d = reshape(img_allscans(i,:),dim1,dim2,dim3);  
0159                     img4d(:,:,:,i) = img3d;
0160                 end               
0161                 No         = sample_img{1}; % copy header
0162                 No.dat     = img4d;         % change file_array
0163                 No.descrip = 'Pronto data';
0164                 No.cal     = [0 1];
0165                 create(No);                 % write header
0166             end
0167             
0168             clear img_allscans;
0169             clear img4d;
0170             
0171             step = step+1;
0172             waitbar(step / (tot_subj*n_modalities*3));
0173         end        
0174     end
0175 end
0176 
0177 delete(h);

Generated on Sun 20-May-2012 13:24:48 by m2html © 2005