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 571 2012-08-14 10:44:04Z schrouff $
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         %Parameters for the masks and indexes of the voxels
0078         %-------------------------------------------------------------------
0079         % get the indices of the voxels within the file array mask (data &
0080         % design step)
0081         ind_ddmask = PRT.fas(mid).idfeat_img(vox_range);
0082         
0083         %load the mask for that modality if another one was specified
0084         if ~isempty(precmask{m})
0085             prec_mask = prt_load_blocks(precmask{m},ind_ddmask);
0086         else
0087             prec_mask = ones(block_size,1);
0088         end
0089         %indexes to access the file array
0090         indm = find(PRT.fs(fid).fas.im==mid);
0091         ifa  = PRT.fs(fid).fas.ifa(indm);
0092         %get the data from each subject of each group and save its linear
0093         %detrended version in a file array
0094         %-------------------------------------------------------------------
0095         if tocomp(mid)  %need to build the file array corresponding to that modality
0096             sample_range=0;
0097             nfa=PRT.fas(mid).dat.dim(1);
0098             datapr=zeros(block_size,nfa);
0099             
0100             %get the data for each subject of each group
0101             for gid = 1:n_groups
0102                 for sid = 1:length(PRT.group(gid).subject)
0103                     n_vols_s = size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0104                     sample_range = (1:n_vols_s)+max(sample_range);
0105                     fname = PRT.group(gid).subject(sid).modality(mid).scans;
0106                     datapr(:,sample_range) = prt_load_blocks(fname,ind_ddmask);
0107                     %check for NaNs, in case of beta maps
0108                     [inan,jnan] = find(isnan(datapr(:,sample_range)));
0109                     if ~isempty(inan)
0110                         for inn=1:length(inan)
0111                             datapr(inan(inn),sample_range(jnan(inn))) = 0;
0112                         end
0113                     end
0114                     
0115                     %detrend if necessary
0116                     if in.mod(mid).detrend ~= 0
0117                         if  isfield(PRT.group(gid).subject(sid).modality(mid).design,'TR')
0118                             TR = PRT.group(gid).subject(sid).modality(mid).design.TR;
0119                         else
0120                             TR = PRT.group(gid).subject(sid).modality(mid).TR;
0121                         end
0122                         switch in.mod(mid).detrend
0123                             case 1
0124                                 C = poly_regressor(length(sample_range), ...
0125                                         in.mod(mid).param_dt);
0126                             case 2
0127                                 C = dct_regressor(length(sample_range), ...
0128                                         in.mod(mid).param_dt,TR);
0129                         end
0130                         R = eye(length(sample_range)) - C*pinv(C);
0131                         datapr(:,sample_range) = datapr(:,sample_range)*R';
0132                     end
0133                 end
0134             end
0135             
0136             if b==1
0137                 % Write the detrended data into the file array .dat
0138                 namedat=['Feature_set_',char(in.mod(mid).mod_name),'.dat'];
0139                 fpd_clean(m) = fopen(fullfile(prt_dir,namedat), 'a','ieee-le.l64'); %#ok<AGROW> % 'a' append
0140                 fwrite(fpd_clean(m), datapr', 'float64',0,'ieee-le.l64');
0141             else
0142                 % Append the data in file .dat
0143                 fwrite(fpd_clean(m), datapr', 'float64',0,'ieee-le.l64');
0144             end
0145             
0146             % get the data to build the kernel
0147             kern_vols(:,indm) = datapr(:,ifa).* ...
0148                 repmat(prec_mask,1,length(ifa));
0149             % if a scaling was entered, apply it now
0150             if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0151                 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0152                     repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0153             end
0154             clear datapr
0155         else
0156             kern_vols(:,indm) = (PRT.fas(mid).dat(ifa,vox_range))' .* ...
0157                 repmat(prec_mask,1,length(ifa));
0158             % if a scaling was entered, apply it now
0159             if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0160                 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0161                     repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0162             end
0163             
0164         end
0165         waitbar(step/ (n_block*n_mods),h);
0166         step=step+1;
0167     end
0168     if size(kern_vols,2)>6e3
0169         % Slower way of estimating kernel but using less memory.
0170         % size limit setup such that no more than ~1Gb of mem is required:
0171         % 1Gb/3(nr of matrices)/8(double)= ~40e6 -> sqrt -> 6e3 element
0172         for ic=1:size(kern_vols,2)
0173             Phi(:,ic) = Phi(:,ic) + kern_vols' * kern_vols(:,ic);
0174         end
0175     else
0176         Phi = Phi + (kern_vols' * kern_vols);
0177     end
0178     bstart = bend+1; bend = min(bstart+block_size-1,n_vox);
0179     clear block_mask kern_vols
0180 end
0181 close(h)
0182 
0183 % closing feature file(s)
0184 if exist('fpd_clean','var')
0185     for ii=1:numel(fpd_clean)
0186         fclose(fpd_clean(ii));
0187     end
0188 end
0189 
0190 % Save kernel and function output
0191 % -------------------------------------------------------------------------
0192 outfile = in.fname;
0193 disp('Saving feature set to: PRT.mat.......>>')
0194 disp(['Saving kernel to: ',in.fs_name,'.mat.......>>'])
0195 fs_file = [prt_dir,in.fs_name];
0196 if spm_matlab_version_chk('7') >= 0
0197     save(outfile,'-V6','PRT');
0198     save(fs_file,'-V6','Phi');
0199 else
0200     save(outfile,'PRT');
0201     save(fs_file,'Phi');
0202 end
0203 disp('Done.')
0204 
0205 % -------------------------------------------------------------------------
0206 % Private functions
0207 % -------------------------------------------------------------------------
0208 function [mask, precmask, headers,PRT] = load_masks(PRT, prt_dir, in, mids)
0209 % function to load the mask for each modality
0210 % -------------------------------------------
0211 n_mods=length(mids);
0212 mask  = cell(1,n_mods);
0213 precmask  = cell(1,n_mods);
0214 headers = cell(1,n_mods);
0215 for m = 1:n_mods
0216     mid=mids(m);
0217     
0218     %get mask for the within-brain voxels (from data and design)
0219     ddmask=PRT.masks(mid).fname;
0220     try
0221         M=nifti(ddmask);
0222     catch
0223         error('prt_prepare_data:CouldNotLoadFile',...
0224             'Could not load mask file');
0225     end
0226     
0227     %get mask for the kernel if one was specified
0228     mfile=in.mod(mid).mask;
0229     if ~isempty(mfile) %&&  mfile ~= 0
0230         try
0231             precM = spm_vol(char(mfile));
0232         catch
0233             error('prt_prepare_data:CouldNotLoadFile',...
0234                 'Could not load mask file for preprocessing');
0235         end
0236     end
0237     
0238     %get header of the first scan of that modality
0239     if isfield(PRT,'fas') && mid<=length(PRT.fas) && ~isempty(PRT.fas(mid).dat)
0240         N=PRT.fas(mid).hdr;
0241     else
0242         N=spm_vol(PRT.group(1).subject(1).modality(mid).scans(1,:));
0243     end
0244     headers{m}=N;
0245     
0246     
0247     % compute voxel dimensions and check for equality if n_mod > 1
0248     if m == 1
0249         n_vox = prod(N.dim(1:3));
0250     elseif n_mods > 1 && n_vox ~= prod(N.dim(1:3))
0251         error('prt_prepare_data:multipleModatlitiesVariableFeatures',...
0252             'Multiple modalities specified, but have variable numbers of features');
0253     end
0254     
0255     %resize the different masks if needed
0256     if any(size(M.dat(:,:,:,1)) ~= N.dim)
0257         warning('prt_prepare_data:maskAndImagesDifferentDim',...
0258             'Mask has different dimensions to the image files. Resizing...');
0259         
0260         V2 = spm_vol(char(ddmask));
0261         mfile_new       = N;
0262         mfile_new.fname = [prt_dir 'updated_1stlevel_mask_m',num2str(mid),'.img'];
0263         tmp             = spm_imcalc([N V2],mfile_new,'0.*i1+(i2>0)');
0264         mask{m}         = mfile_new.fname;
0265         PRT.masks(mid).fname = mfile_new.fname;
0266     else
0267         mask{m} = ddmask;
0268         
0269     end
0270     if ~isempty(mfile)  && any((precM.dim~= N.dim)) % && mfile ~= 0
0271         warning('prt_prepare_data:maskAndImagesDifferentDim',...
0272             'Preprocessing mask has different dimensions to the image files. Resizing...');
0273         
0274         V2 = spm_vol(char(mfile));
0275         mfile_new       = N;
0276         nummask=1;
0277         %if more than one 2nd level mask to resize
0278         while exist(fullfile( ...
0279                 prt_dir,['updated_2ndlevel_mask_m',num2str(mid),'_',...
0280                 num2str(nummask),'.img']))
0281                 nummask=nummask+1;
0282         end
0283         mfile_new.fname = [prt_dir 'updated_2ndlevel_mask_m',num2str(mid),...
0284             '_',num2str(nummask),'.img'];
0285         tmp             = spm_imcalc([N V2],mfile_new,'0.*i1+(i2>0)');
0286         precmask{m}     = mfile_new.fname;
0287     else
0288         precmask{m} = mfile;
0289     end
0290     clear M N precM V1 V2 mfile mfile_new
0291 end
0292 return
0293 
0294 function c=poly_regressor(n,order)
0295 %n: length of the series
0296 %order: the order of polynomial function to fit the tend
0297 
0298 basis=repmat([1:n]',[1 order]);
0299 o=repmat([1:order],[n 1]);
0300 c=[ones(n,1) basis.^o];
0301 return
0302 
0303 function c=dct_regressor(n,cut_off,TR)
0304 % n: length of the series
0305 %cut_off: the cut off perioed in second (1/ cut off frequency)
0306 %TR: TR
0307 
0308 if cut_off<0
0309     error('cut off cannot be negative')
0310 end
0311 
0312 T=n*TR;
0313 order=floor((T/cut_off)*2)+1;
0314 c=spm_dctmtx(n,order);
0315 return

Generated on Mon 03-Sep-2012 18:07:18 by m2html © 2005