Home > . > prt_fs_modality.m

prt_fs_modality

PURPOSE ^

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

SYNOPSIS ^

function [PRT,Phi] = prt_fs_modality(PRT,in, flag, addin)

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

 in.fid:      index of feature set to create
 in.tocomp:   vector of booleans indicating whether to build the feature set
 in.precmask: cell array containing the names of the second-level mask for
           each modality to build

 flag:     set to 1 to compute one kernel per region as labelled in atlas
 addin:    additional inputs for this operation to optimize computation

 Outputs:
 --------
 Writes the kernel matrix to the path indicated by in.fs_name and the
 feature set in a file array if it needs to be computed
__________________________________________________________________________
 Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [PRT,Phi] = prt_fs_modality(PRT,in, flag, addin)
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 % in.fid:      index of feature set to create
0019 % in.tocomp:   vector of booleans indicating whether to build the feature set
0020 % in.precmask: cell array containing the names of the second-level mask for
0021 %           each modality to build
0022 %
0023 % flag:     set to 1 to compute one kernel per region as labelled in atlas
0024 % addin:    additional inputs for this operation to optimize computation
0025 %
0026 % Outputs:
0027 % --------
0028 % Writes the kernel matrix to the path indicated by in.fs_name and the
0029 % feature set in a file array if it needs to be computed
0030 %__________________________________________________________________________
0031 % Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory
0032 
0033 % Written by A. Marquand and J. Schrouff
0034 % $Id: prt_fs.m 650 2013-02-20 13:36:58Z amarquan $
0035 
0036 % Configure some variables and get defaults
0037 % -------------------------------------------------------------------------
0038 prt_dir  = regexprep(in.fname,'PRT.mat', ''); % or: fileparts(fname);
0039 def=prt_get_defaults('fs');
0040 fid = in.fid;
0041 
0042 % get the indexes of the samples and of the features to use if flag is set
0043 % to 1
0044 if nargin>=3 && flag
0045     if ~isfield(addin,'ID')
0046         ID= PRT.fs(fid).id_mat;
0047         n_mods=length(in.mod);
0048         mids=[];
0049         for i=1:n_mods
0050             if ~isempty(in.mod(i).mod_name)
0051                 mids = [mids, i];
0052             end
0053         end
0054     else
0055         ID = addin.ID;
0056         mids = unique(ID(:,3));
0057     end
0058     n_mods=length(mids);
0059     if isfield(addin,'idvox_fas')
0060         idvox = addin.idvox_fas;
0061     else
0062         if in.tocomp(mids(1)) % build everything
0063             idvox = 1:PRT.fas(mids(1)).dat.dim(2); %n_vox has to be the same for all concatenated modalities (version 1.1)
0064         else
0065             if ~isempty(PRT.fs(fid).modality(mids(1)).idfeat_fas)
0066                 idvox = PRT.fs(fid).modality(mids(1)).idfeat_fas;
0067             else
0068                 idvox = 1:PRT.fas(mids(1)).dat.dim(2);
0069             end
0070         end
0071     end
0072     n_vox = numel(idvox);
0073     nfa = 0;    
0074 else
0075     % get the index of the modalities for which the user wants a kernel/data
0076     n_mods=length(in.mod);
0077     mids=[];
0078     for i=1:n_mods
0079         if ~isempty(in.mod(i).mod_name)
0080             mids = [mids, i];
0081         end
0082     end
0083     n_mods=length(mids);
0084     ID= PRT.fs(fid).id_mat;
0085     nfa = [];
0086     for m = 1:n_mods
0087         nfa   = [nfa, PRT.fas(mids(m)).dat.dim(1)];
0088         if in.tocomp(mids(m))
0089             n_vox = PRT.fas(mids(m)).dat.dim(2); %n_vox has to be the same for all concatenated modalities (version 1.1)
0090         else
0091             if ~isempty(PRT.fs(fid).modality(m).idfeat_fas)
0092                 idvox = PRT.fs(fid).modality(m).idfeat_fas;
0093             else
0094                 idvox = 1:PRT.fas(mids(m)).dat.dim(2);
0095             end
0096             n_vox = numel(idvox);
0097         end
0098     end
0099 end
0100 
0101 
0102 
0103 % -------------------------------------------------------------------------
0104 % ---------------------Build file arrays and kernel------------------------
0105 % -------------------------------------------------------------------------
0106 n   = size(ID,1);
0107 Phi = zeros(n);
0108 % set memory limit
0109 mem         = def.mem_limit;
0110 block_size  = floor(mem/(8*3)/max([nfa, n])); % Block size (double = 8 bytes)
0111 n_block     = ceil(n_vox/block_size);
0112 
0113 bstart = 1; bend = min(block_size,n_vox);
0114 if nargin<3 || ~flag
0115     h = waitbar(0,'Please wait while preparing feature set');
0116     step=1;
0117 end
0118 for b = 1:n_block
0119     if nargin<3 || ~flag
0120         disp ([' > preparing block: ', num2str(b),' of ',num2str(n_block),' ...'])
0121     end
0122     vox_range  = bstart:bend;
0123     block_size = length(vox_range);
0124     kern_vols  = zeros(block_size,n);
0125     for m=1:n_mods
0126         mid=mids(m);
0127         %Parameters for the masks and indexes of the voxels
0128         %-------------------------------------------------------------------
0129         % get the indices of the voxels within the file array mask (data &
0130         % design step)
0131         ind_ddmask = PRT.fas(mid).idfeat_img(vox_range);
0132         
0133         %load the mask for that modality if another one was specified
0134         if ~isempty(in.precmask{m})
0135             prec_mask = prt_load_blocks(in.precmask{m},ind_ddmask);
0136         else
0137             prec_mask = ones(block_size,1);
0138         end
0139         %indexes to access the file array
0140         indm = find(PRT.fs(fid).fas.im==mid);
0141         ifa  = PRT.fs(fid).fas.ifa(indm);
0142         indm = find(ID(:,3) == mid);
0143         %get the data from each subject of each group and save its linear
0144         %detrended version in a file array
0145         %-------------------------------------------------------------------
0146         if in.tocomp(mid)  %need to build the file array corresponding to that modality
0147             n_groups = length(PRT.group);
0148             sample_range=0;
0149             nfa=PRT.fas(mid).dat.dim(1);
0150             datapr=zeros(block_size,nfa);
0151             
0152             %get the data for each subject of each group
0153             for gid = 1:n_groups
0154                 for sid = 1:length(PRT.group(gid).subject)
0155                     n_vols_s = size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0156                     sample_range = (1:n_vols_s)+max(sample_range);
0157                     fname = PRT.group(gid).subject(sid).modality(mid).scans;
0158                     datapr(:,sample_range) = prt_load_blocks(fname,ind_ddmask);
0159                     %check for NaNs, in case of beta maps
0160                     [inan,jnan] = find(isnan(datapr(:,sample_range)));
0161                     if ~isempty(inan)
0162                         disp('Warning: NaNs found in loaded data')
0163                         disp('Consider updating 1st level mask for better performance')
0164                         for inn=1:length(inan)
0165                             datapr(inan(inn),sample_range(jnan(inn))) = 0;
0166                         end
0167                     end
0168                     
0169                     %detrend if necessary
0170                     if in.mod(mid).detrend ~= 0
0171                         if  isfield(PRT.group(gid).subject(sid).modality(mid).design,'TR')
0172                             TR = PRT.group(gid).subject(sid).modality(mid).design.TR;
0173                         else
0174                             try
0175                                 TR = PRT.group(gid).subject(sid).modality(mid).TR;
0176                             catch
0177                                 error('detrend:TRnotfound','No TR in data, suggesting that detrend is not necessary')
0178                             end
0179                         end
0180                         switch in.mod(mid).detrend
0181                             case 1
0182                                 C = poly_regressor(length(sample_range), ...
0183                                         in.mod(mid).param_dt);
0184                             case 2
0185                                 C = dct_regressor(length(sample_range), ...
0186                                         in.mod(mid).param_dt,TR);
0187                         end
0188                         R = eye(length(sample_range)) - C*pinv(C);
0189                         datapr(:,sample_range) = datapr(:,sample_range)*R';
0190                     end
0191                 end
0192             end
0193             
0194             if b==1
0195                 % Write the detrended data into the file array .dat
0196                 namedat=['Feature_set_',char(in.mod(mid).mod_name),'.dat'];
0197                 fpd_clean(m) = fopen(fullfile(prt_dir,namedat), 'w','ieee-le'); %#ok<AGROW> % 'a' append
0198                 fwrite(fpd_clean(m), datapr', 'float32',0,'ieee-le');
0199             else
0200                 % Append the data in file .dat
0201                 fwrite(fpd_clean(m), datapr', 'float32',0,'ieee-le');
0202             end
0203             
0204             % get the data to build the kernel
0205             kern_vols(:,indm) = datapr(:,ifa).* ...
0206                 repmat(prec_mask~=0,1,length(ifa));
0207             % if a scaling was entered, apply it now
0208             if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0209                 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0210                     repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0211             end
0212             clear datapr
0213         else
0214             kern_vols(:,indm) = (PRT.fas(mid).dat(ifa,idvox(vox_range)))';
0215             % if a scaling was entered, apply it now
0216             if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0217                 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0218                     repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0219             end
0220             
0221         end
0222         if nargin<3 || ~flag
0223             waitbar(step/ (n_block*n_mods),h);
0224             step=step+1;
0225         end
0226     end
0227     if size(kern_vols,2)>6e3
0228         % Slower way of estimating kernel but using less memory.
0229         % size limit setup such that no more than ~1Gb of mem is required:
0230         % 1Gb/3(nr of matrices)/8(double)= ~40e6 -> sqrt -> 6e3 element
0231         for ic=1:size(kern_vols,2)
0232             Phi(:,ic) = Phi(:,ic) + kern_vols' * kern_vols(:,ic);
0233         end
0234     else
0235         Phi = Phi + (kern_vols' * kern_vols);
0236     end
0237     bstart = bend+1; bend = min(bstart+block_size-1,n_vox);
0238     clear block_mask kern_vols
0239 end
0240 if nargin<3 || ~flag
0241     close(h)
0242 end
0243 
0244 % closing feature file(s)
0245 if exist('fpd_clean','var')
0246     for ii=1:numel(fpd_clean)
0247         fclose(fpd_clean(ii));
0248     end
0249 end
0250 
0251 
0252 % -------------------------------------------------------------------------
0253 % Private functions
0254 % -------------------------------------------------------------------------
0255 function c = poly_regressor(n,order)
0256 % n:     length of the series
0257 % order: the order of polynomial function to fit the tend
0258 
0259 basis = repmat([1:n]',[1 order]);
0260 o = repmat([1:order],[n 1]);
0261 c = [ones(n,1) basis.^o];
0262 return
0263 
0264 function c=dct_regressor(n,cut_off,TR)
0265 % n:       length of the series
0266 % cut_off: the cut off perioed in second (1/ cut off frequency)
0267 % TR:      TR
0268 
0269 if cut_off<0
0270     error('cut off cannot be negative')
0271 end
0272 
0273 T = n*TR;
0274 order = floor((T/cut_off)*2)+1;
0275 c = spm_dctmtx(n,order);
0276 return

Generated on Tue 10-Feb-2015 18:16:33 by m2html © 2005