Home > . > prt_fs.m

prt_fs

PURPOSE ^

Function to build file arrays containing the (linearly detrended) data

SYNOPSIS ^

function [outfile] = prt_fs(PRT,in)

DESCRIPTION ^

 Function to build file arrays containing the (linearly detrended) data
 and compute a linear (dot product) kernel from them

 Inputs:
 -------
 in.fname:      filename for the PRT.mat (string)
 in.fs_name:    name of fs and relative path filename for the kernel matrix

 in.mod(m).mod_name:  name of modality to include in this kernel (string)
 in.mod(m).detrend:   detrend (scalar: 0 = none, 1 = linear)
 in.mod(m).param_dt:  parameters for the kernel detrend (e.g. DCT bases)
 in.mod(m).mode:      'all_cond' or 'all_scans' (string)
 in.mod(m).mask:      mask file used to create the kernel
 in.mod(m).normalise: 0 = none, 1 = normalise_kernel, 2 = scale modality
 in.mod(m).matnorm:   filename for scaling matrix

 Outputs:
 --------
 Calls prt_init_fs to populate basic fields in PRT.fs(f)...
 Writes PRT.mat
 Writes the kernel matrix to the path indicated by in.fs_name
__________________________________________________________________________
 Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [outfile] = prt_fs(PRT,in)
0002 % Function to build file arrays containing the (linearly detrended) data
0003 % and compute a linear (dot product) kernel from them
0004 %
0005 % Inputs:
0006 % -------
0007 % in.fname:      filename for the PRT.mat (string)
0008 % in.fs_name:    name of fs and relative path filename for the kernel matrix
0009 %
0010 % in.mod(m).mod_name:  name of modality to include in this kernel (string)
0011 % in.mod(m).detrend:   detrend (scalar: 0 = none, 1 = linear)
0012 % in.mod(m).param_dt:  parameters for the kernel detrend (e.g. DCT bases)
0013 % in.mod(m).mode:      'all_cond' or 'all_scans' (string)
0014 % in.mod(m).mask:      mask file used to create the kernel
0015 % in.mod(m).normalise: 0 = none, 1 = normalise_kernel, 2 = scale modality
0016 % in.mod(m).matnorm:   filename for scaling matrix
0017 %
0018 % Outputs:
0019 % --------
0020 % Calls prt_init_fs to populate basic fields in PRT.fs(f)...
0021 % Writes PRT.mat
0022 % Writes the kernel matrix to the path indicated by in.fs_name
0023 %__________________________________________________________________________
0024 % Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory
0025 
0026 % Written by A. Marquand and J. Schrouff
0027 % $Id: prt_fs.m 542 2012-05-20 10:48:51Z cphillip $
0028 
0029 % Configure some variables and get defaults
0030 % -------------------------------------------------------------------------
0031 prt_dir  = regexprep(in.fname,'PRT.mat', ''); % or: fileparts(fname);
0032 n_groups = length(PRT.group);
0033 
0034 % get the index of the modalities for which the user wants a kernel/data
0035 n_mods=length(in.mod);
0036 mids=[];
0037 for i=1:n_mods
0038     if ~isempty(in.mod(i).mod_name)
0039         mids = [mids, i];
0040     end
0041 end
0042 n_mods=length(mids);
0043 def=prt_get_defaults('fs');
0044 
0045 % Load mask(s) and resize if necessary
0046 [mask,precmask,headers,PRT] = load_masks(PRT, prt_dir, in,mids);
0047 
0048 % Initialize the file arrays, kernel and feature set parameters
0049 [fid,PRT,tocomp] = prt_init_fs(PRT,in,mids,mask,precmask,headers);
0050 
0051 % -------------------------------------------------------------------------
0052 % ---------------------Build file arrays and kernel------------------------
0053 % -------------------------------------------------------------------------
0054 n   = length(PRT.fs(fid).id_mat);
0055 Phi = zeros(n);
0056 
0057 % set memory limit
0058 nfa = [];
0059 for m=1:n_mods
0060     nfa   = [nfa, PRT.fas(mids(m)).dat.dim(1)];
0061     n_vox = PRT.fas(mids(m)).dat.dim(2);
0062 end
0063 mem         = def.mem_limit;
0064 block_size  = floor(mem/8/max([nfa, n])); % Block size (double = 8 bytes)
0065 n_block     = ceil(n_vox/block_size);
0066 
0067 bstart=1; bend=min(block_size,n_vox);
0068 h = waitbar(0,'Please wait while preparing feature set');
0069 step=1;
0070 for b = 1:n_block
0071     disp ([' > preparing block: ', num2str(b),' of ',num2str(n_block),' ...'])
0072     vox_range = bstart:bend;
0073     block_size=length(vox_range);
0074     kern_vols=zeros(block_size,n);
0075     for m=1:n_mods
0076         mid=mids(m);
0077         
0078         %Parameters for the masks and indexes of the voxels
0079         %-------------------------------------------------------------------
0080         % get the indices of the voxels within the file array mask (data &
0081         % design step)
0082         ind_ddmask = PRT.fas(mid).idfeat_img(vox_range);
0083         
0084         %load the mask for that modality if another one was specified
0085         if ~isempty(precmask{m})
0086             prec_mask = prt_load_blocks(precmask{m},ind_ddmask);
0087         else
0088             prec_mask = ones(block_size,1);
0089         end
0090         %indexes to access the file array
0091         indm=find(PRT.fs(fid).fas.im==mid);
0092         ifa=PRT.fs(fid).fas.ifa(indm);
0093         
0094         %get the data from each subject of each group and save its linear
0095         %detrended version in a file array
0096         %-------------------------------------------------------------------
0097         if tocomp(mid)  %need to build the file array corresponding to that modality
0098             sample_range=0;
0099             nfa=PRT.fas(mid).dat.dim(1);
0100             datapr=zeros(nfa,block_size);
0101             
0102             %get the data for each subject of each group
0103             for gid = 1:n_groups
0104                 for sid = 1:length(PRT.group(gid).subject)
0105                     n_vols_s = size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0106                     sample_range=(1:n_vols_s)+max(sample_range);
0107                     fname = PRT.group(gid).subject(sid).modality(mid).scans;
0108                     datapr(sample_range,:) = (prt_load_blocks(fname,ind_ddmask))';
0109                     %check for NaNs, in case of beta maps
0110                     [inan,jnan]=find(isnan(datapr(sample_range,:)));
0111                     if ~isempty(inan)
0112                         for inn=1:length(inan)
0113                             datapr(sample_range(inan(inn)),jnan(inn))=0;
0114                         end
0115                     end
0116                     
0117                     %detrend if necessary
0118                     if in.mod(mid).detrend ~= 0
0119                         if  isfield(PRT.group(gid).subject(sid).modality(mid).design,'TR')
0120                             TR=PRT.group(gid).subject(sid).modality(mid).design.TR;
0121                         else
0122                             TR=PRT.group(gid).subject(sid).modality(mid).TR;
0123                         end
0124                         switch in.mod(mid).detrend
0125                             case 1
0126                                 c= poly_regressor(length(sample_range),in.mod(mid).param_dt);
0127                             case 2
0128                                 c=dct_regressor(length(sample_range),in.mod(mid).param_dt,TR);
0129                         end
0130                         C = [c];
0131                         R = eye(length(sample_range)) - C*pinv(C);
0132                         datapr(sample_range,:) = R*datapr(sample_range,:);
0133                     end
0134                 end
0135             end
0136             
0137             % Write the detrended data into the file array
0138             namedat=['Feature_set_',char(in.mod(mid).mod_name),'.dat'];
0139             fpd_clean = fopen(fullfile(prt_dir,namedat), 'a','ieee-le.l64'); % 'a' append
0140             if b==1
0141                 % write the data in file .dat
0142                 fwrite(fpd_clean, datapr, 'float64',0,'ieee-le.l64');
0143                 fclose(fpd_clean);
0144             else
0145                 % Append the data in file .dat
0146                 fseek(fpd_clean,0,'eof');
0147                 fwrite(fpd_clean, datapr, 'float64',0,'ieee-le.l64');
0148                 fclose(fpd_clean);
0149             end
0150             
0151             % get the data to build the kernel
0152             kern_vols(:,indm) = datapr(ifa,:)'.* ...
0153                 repmat(prec_mask,1,length(ifa));
0154             % if a scaling was entered, apply it now
0155             if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0156                 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0157                     repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0158             end
0159             clear datapr
0160         else
0161             kern_vols(:,indm) = (PRT.fas(mid).dat(ifa,vox_range))' .* ...
0162                 repmat(prec_mask,1,length(ifa));
0163             % if a scaling was entered, apply it now
0164             if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0165                 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0166                     repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0167             end
0168             
0169         end
0170         waitbar(step/ (n_block*n_mods),h);
0171         step=step+1;
0172     end
0173     Phi = Phi + (kern_vols' * kern_vols);
0174     bstart=bend+1; bend=min(bstart+block_size-1,n_vox);
0175     clear block_mask kern_vols
0176 end
0177 close(h)
0178 
0179 % Save kernel and function output
0180 % -------------------------------------------------------------------------
0181 outfile = in.fname;
0182 disp('Saving feature set to: PRT.mat.......>>')
0183 disp(['Saving kernel to: ',in.fs_name,'.mat.......>>'])
0184 fs_file = [prt_dir,in.fs_name];
0185 if spm_matlab_version_chk('7') >= 0
0186     save(outfile,'-V6','PRT');
0187     save(fs_file,'-V6','Phi');
0188 else
0189     save(outfile,'PRT');
0190     save(fs_file,'Phi');
0191 end
0192 disp('Done.')
0193 
0194 % -------------------------------------------------------------------------
0195 % Private functions
0196 % -------------------------------------------------------------------------
0197 function [mask, precmask, headers,PRT] = load_masks(PRT, prt_dir, in, mids)
0198 % function to load the mask for each modality
0199 % -------------------------------------------
0200 n_mods=length(mids);
0201 mask  = cell(1,n_mods);
0202 precmask  = cell(1,n_mods);
0203 headers = cell(1,n_mods);
0204 for m = 1:n_mods
0205     mid=mids(m);
0206     
0207     %get mask for the within-brain voxels (from data and design)
0208     ddmask=PRT.masks(mid).fname;
0209     try
0210         M=nifti(ddmask);
0211     catch
0212         error('prt_prepare_data:CouldNotLoadFile',...
0213             'Could not load mask file');
0214     end
0215     
0216     %get mask for the kernel if one was specified
0217     mfile=in.mod(mid).mask;
0218     if ~isempty(mfile) %&&  mfile ~= 0
0219         try
0220             precM = spm_vol(char(mfile));
0221         catch
0222             error('prt_prepare_data:CouldNotLoadFile',...
0223                 'Could not load mask file for preprocessing');
0224         end
0225     end
0226     
0227     %get header of the first scan of that modality
0228     if isfield(PRT,'fas') && mid<=length(PRT.fas) && ~isempty(PRT.fas(mid).dat)
0229         N=PRT.fas(mid).hdr;
0230     else
0231         N=spm_vol(PRT.group(1).subject(1).modality(mid).scans(1,:));
0232     end
0233     headers{m}=N;
0234     
0235     
0236     % compute voxel dimensions and check for equality if n_mod > 1
0237     if m == 1
0238         n_vox = prod(N.dim(1:3));
0239     elseif n_mods > 1 && n_vox ~= prod(N.dim(1:3))
0240         error('prt_prepare_data:multipleModatlitiesVariableFeatures',...
0241             'Multiple modalities specified, but have variable numbers of features');
0242     end
0243     
0244     %resize the different masks if needed
0245     if any(size(M.dat(:,:,:,1)) ~= N.dim)
0246         warning('prt_prepare_data:maskAndImagesDifferentDim',...
0247             'Mask has different dimensions to the image files. Resizing...');
0248         
0249         V2 = spm_vol(char(ddmask));
0250         mfile_new       = N;
0251         mfile_new.fname = [prt_dir 'updated_1stlevel_mask_m',num2str(mid),'.img'];
0252         tmp             = spm_imcalc([N V2],mfile_new,'0.*i1+(i2>0)');
0253         mask{m}         = mfile_new.fname;
0254         PRT.masks(mid).fname = mfile_new.fname;
0255     else
0256         mask{m} = ddmask;
0257         
0258     end
0259     if ~isempty(mfile)  && any((precM.dim~= N.dim)) % && mfile ~= 0
0260         warning('prt_prepare_data:maskAndImagesDifferentDim',...
0261             'Preprocessing mask has different dimensions to the image files. Resizing...');
0262         
0263         V2 = spm_vol(char(mfile));
0264         mfile_new       = N;
0265         mfile_new.fname = [prt_dir 'updated_2ndlevel_mask_m',num2str(mid),'.img'];
0266         tmp             = spm_imcalc([N V2],mfile_new,'0.*i1+(i2>0)');
0267         precmask{m}     = mfile_new.fname;
0268     else
0269         precmask{m} = mfile;
0270     end
0271     clear M N precM V1 V2 mfile mfile_new
0272 end
0273 return
0274 
0275 function c=poly_regressor(n,order)
0276 %n: length of the series
0277 %order: the order of polynomial function to fit the tend
0278 
0279 basis=repmat([1:n]',[1 order]);
0280 o=repmat([1:order],[n 1]);
0281 c=[ones(n,1) basis.^o];
0282 return
0283 
0284 function c=dct_regressor(n,cut_off,TR)
0285 % n: length of the series
0286 %cut_off: the cut off perioed in second (1/ cut off frequency)
0287 %TR: TR
0288 
0289 if cut_off<0
0290     error('cut off cannot be negative')
0291 end
0292 
0293 T=n*TR;
0294 order=floor((T/cut_off)*2)+1;
0295 c=spm_dctmtx(n,order);
0296 return

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