Home > . > prt_init_fs.m

prt_init_fs

PURPOSE ^

function to initialise the kernel data structure

SYNOPSIS ^

function [fid,PRT,tocomp] = prt_init_fs(PRT, in, mids,mask,precmask,headers)

DESCRIPTION ^

 function to initialise the kernel data structure
 ------------------------------------------------

 FORMAT: Two modes are possible:
     fid = prt_init_fs(PRT, in)
     [fid, PRT, tocomp] = prt_init_fs(PRT, in)

 USAGE 1:
 -------------------------------------------------------------------------
 function will return the id of a feature set or an error if it doesn't
 exist in PRT.mat
 Input:
 ------
 in.fs_name: name for the feature set (string)

 Output:
 -------
 fid : is the identifier for the feature set in PRT.mat

 USAGE 2:
 -------------------------------------------------------------------------
 function will create the feature set in PRT.mat and overwrite it if it
 already exists.
 Input:
 ------
 in.fs_name: name for the feature set (string)
 in.fname:   name of PRT.mat

 in.mod(m).mod_name:  name of the modality
 in.mod(m).detrend:   type of detrending
 in.mod(m).mode:      'all_scans' or 'all_cond'
 in.mod(m).mask:       mask used to create the feature set
 in.mod(m).param_dt:  parameters used for detrending (if any)
 in.mod(m).normalise: scale the input scans or not
 in.mod(m).matnorm:   mat file used to scale the input scans

 Output:
 -------
 fid : is the identifier for the model constructed in PRT.mat

 Populates the following fields in PRT.mat (copied from above):
   PRT.fs(f).fs_name
   PRT.fs(f).fas
   PRT.fs(f).k_file
 Also computes the following fields:
   PRT.fs(f).id_mat:       Identifier matrix (useful later)
   PRT.fs(f).id_col_names: Columns in the id matrix

 Note: this function does not write PRT.mat. That should be done by the
       calling function
__________________________________________________________________________
 Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [fid,PRT,tocomp] = prt_init_fs(PRT, in, mids,mask,precmask,headers)
0002 % function to initialise the kernel data structure
0003 % ------------------------------------------------
0004 %
0005 % FORMAT: Two modes are possible:
0006 %     fid = prt_init_fs(PRT, in)
0007 %     [fid, PRT, tocomp] = prt_init_fs(PRT, in)
0008 %
0009 % USAGE 1:
0010 % -------------------------------------------------------------------------
0011 % function will return the id of a feature set or an error if it doesn't
0012 % exist in PRT.mat
0013 % Input:
0014 % ------
0015 % in.fs_name: name for the feature set (string)
0016 %
0017 % Output:
0018 % -------
0019 % fid : is the identifier for the feature set in PRT.mat
0020 %
0021 % USAGE 2:
0022 % -------------------------------------------------------------------------
0023 % function will create the feature set in PRT.mat and overwrite it if it
0024 % already exists.
0025 % Input:
0026 % ------
0027 % in.fs_name: name for the feature set (string)
0028 % in.fname:   name of PRT.mat
0029 %
0030 % in.mod(m).mod_name:  name of the modality
0031 % in.mod(m).detrend:   type of detrending
0032 % in.mod(m).mode:      'all_scans' or 'all_cond'
0033 % in.mod(m).mask:       mask used to create the feature set
0034 % in.mod(m).param_dt:  parameters used for detrending (if any)
0035 % in.mod(m).normalise: scale the input scans or not
0036 % in.mod(m).matnorm:   mat file used to scale the input scans
0037 %
0038 % Output:
0039 % -------
0040 % fid : is the identifier for the model constructed in PRT.mat
0041 %
0042 % Populates the following fields in PRT.mat (copied from above):
0043 %   PRT.fs(f).fs_name
0044 %   PRT.fs(f).fas
0045 %   PRT.fs(f).k_file
0046 % Also computes the following fields:
0047 %   PRT.fs(f).id_mat:       Identifier matrix (useful later)
0048 %   PRT.fs(f).id_col_names: Columns in the id matrix
0049 %
0050 % Note: this function does not write PRT.mat. That should be done by the
0051 %       calling function
0052 %__________________________________________________________________________
0053 % Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory
0054 
0055 % Written by A Marquand
0056 % $Id: prt_init_fs.m 542 2012-05-20 10:48:51Z cphillip $
0057 
0058 % find index for the new feature set
0059 fs_exists = false;
0060 if ~(prt_checkAlphaNumUnder(in.fs_name))
0061     beep
0062     disp('Feature set name should be entered in alphanumeric format only')
0063     disp('Please correct')
0064     return
0065 end
0066 if isfield(PRT,'fs')
0067     if any(strcmpi(in.fs_name,{PRT.fs(:).fs_name}))
0068         fid = find(strcmpi(in.fs_name,{PRT.fs(:).fs_name}));
0069         fs_exists = true;
0070     else
0071         fid = length(PRT.fs)+1;
0072     end
0073 else
0074     fid = 1;
0075 end
0076 
0077 % do we want to initialise the feature set?
0078 if nargout == 1
0079     if fs_exists
0080         % just display message and exit (returning id)
0081         disp(['Feature set ''',in.fs_name,''' found in PRT.mat.']);
0082     else
0083         error('prt_init_fs:fsNotFoundinPRT',...
0084             ['Feature set ''',in.fs_name,''' not found in PRT.mat.']);
0085     end
0086 else
0087     % initialise
0088     
0089     [pathName fileName]=fileparts(in.fname);
0090     if fs_exists
0091         warning('prt_init_fs:overwriteFsInPRT',...
0092             ['Feature set ''',in.fs_name,''' found in PRT.mat. Overwriting ...']);
0093     else
0094         % doesn't exist. initialise the structure
0095         disp(['Feature set ''',in.fs_name,''' not found in PRT.mat. Creating...'])
0096     end
0097     
0098     PRT.fs(fid).fs_name = in.fs_name;
0099     PRT.fs(fid).k_file = in.fs_name;
0100     PRT.fs(fid).id_col_names = {'group','subject','modality','condition','block','scan'};
0101     PRT.fs(fid).fas=struct('im',[],'ifa',[]);
0102     n_vox=0;
0103     n_mods=length(mids);
0104     for m = 1:n_mods
0105         PRT.fs(fid).modality(m).mod_name = in.mod(mids(m)).mod_name;
0106         PRT.fs(fid).modality(m).detrend  = in.mod(mids(m)).detrend;
0107         PRT.fs(fid).modality(m).param_dt = in.mod(mids(m)).param_dt;
0108         PRT.fs(fid).modality(m).mode     = in.mod(mids(m)).mode;
0109         %get indexes from mask specified in the data and design step
0110         vm = spm_vol(mask{m});
0111         vm = spm_read_vols(vm);
0112         PRT.fs(fid).modality(m).feat_idx_img = find(vm>0);
0113         if isempty(find(vm>0))
0114             error('prt_prepare_data:NoVoxelinMask',...
0115                 ['Mask of modality ',num2str(m),' does not contain any voxel >0'])
0116         end
0117         mid=mids(m);
0118         if m==1
0119             n_vox=length(find(vm>0));
0120         end
0121         if n_vox~=length(find(vm>0))
0122             error('prt_prepare_data:MasksNotConsistent',...
0123                 'Masks access areas of different sizes across modalities')
0124         end
0125         %get subindexes from mask specified in the data prepare step
0126         if ~isempty(precmask{m})
0127             vm = spm_vol(precmask{m});
0128             vm = spm_read_vols(vm);
0129             if isempty(find(vm>0))
0130                 error('prt_prepare_data:NoVoxelinMask',...
0131                     ['2nd level mask of modality ',num2str(m),' does not contain any voxel >0'])
0132             end
0133             [d,PRT.fs(fid).modality(m).idfeat_fas] = intersect(PRT.fs(fid).modality(m).feat_idx_img, find(vm>0));
0134         else
0135             PRT.fs(fid).modality(m).idfeat_fas=[];
0136         end
0137         PRT.fs(fid).modality(m).normalise=struct('type',[],'scaling',[]);
0138     end
0139     
0140     indm = zeros(n_mods,1);
0141     szm = zeros(n_mods,1);
0142     
0143     % First count the total number of samples. Loops are needed since each
0144     % subject may have a variable number of scans
0145     n = 0;
0146     for gid = 1:length(PRT.group) % group
0147         for sid = 1:length(PRT.group(gid).subject);  % subject
0148             for m = 1:n_mods
0149                 mid = mids(m);
0150                 if strcmpi(in.mod(mid).mode,'all_scans');
0151                     n = n + size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0152                 elseif strcmpi(in.mod(mid).mode,'all_cond')
0153                     if ~isfield(PRT.group(gid).subject(sid).modality(mid).design,'conds')
0154                         error('prt_model:fsIsAllCondModelisAllScans',...
0155                             ['''All conditions'' selected for modality ', num2str(m)...
0156                             ' but no design was specified. This syntax is invalid, '...
0157                             'Please use ''All Scans'' instead.']);
0158                     end
0159                     for cid = 1:length(PRT.group(gid).subject(sid).modality(mid).design.conds)    % condition
0160                         n = n + length(PRT.group(gid).subject(sid).modality(mid).design.conds(cid).scans);
0161                     end
0162                 end
0163             end  % modality
0164         end  % subject
0165     end  % group
0166     PRT.fs(fid).id_mat = zeros(n,length(PRT.fs(fid).id_col_names));
0167     PRT.fs(fid).fas.im = zeros(n,1);
0168     PRT.fs(fid).fas.ifa= zeros(n,1);
0169     
0170     % Count the total number of samples and set sample ids for the kernel
0171     % Set fas for the file arrays
0172     sample_range = 0;
0173     for gid = 1:length(PRT.group) % group
0174         for sid = 1:length(PRT.group(gid).subject);  % subject
0175             for m = 1:n_mods
0176                 mid = mids(m);
0177                 
0178                 if strcmpi(in.mod(mid).mode,'all_scans')
0179                     n_vols_s  = size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0180                     all_scans = 1:n_vols_s;
0181                     
0182                     % configure indices
0183                     sample_range = (1:n_vols_s)+max(sample_range);
0184                     PRT.fs(fid).id_mat(sample_range,1) = gid;
0185                     PRT.fs(fid).id_mat(sample_range,2) = sid;
0186                     PRT.fs(fid).id_mat(sample_range,3) = mid;
0187                     
0188                     if isfield(PRT.group(gid).subject(sid).modality(mid).design,'conds')
0189                         conds = PRT.group(gid).subject(sid).modality(mid).design.conds;
0190                         for cid = 1:length(conds)
0191                             scans  = PRT.group(gid).subject(sid).modality(mid).design.conds(cid).scans;
0192                             blocks = PRT.group(gid).subject(sid).modality(mid).design.conds(cid).blocks;
0193                             
0194                             PRT.fs(fid).id_mat(sample_range(scans),4) = cid;
0195                             PRT.fs(fid).id_mat(sample_range(scans),5) = blocks;
0196                             %PRT.fs(fid).id_mat(sample_range(scans),6) = 1:length(scans);
0197                         end
0198                         
0199                         PRT.fs(fid).id_mat(sample_range,6) = 1:length(all_scans);
0200                     else
0201                         scans  = 1:size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0202                         PRT.fs(fid).id_mat(sample_range,6) = scans;
0203                     end
0204                     
0205                     sctoadd=(1:n_vols_s)+indm(m);
0206                     PRT.fs(fid).fas.ifa(sample_range)=sctoadd';
0207                     PRT.fs(fid).fas.im(sample_range)=mid*ones(n_vols_s,1);
0208                     %configure indices for the file array
0209                     indm(m)=n_vols_s+max(indm(m));
0210                 elseif strcmpi(in.mod(mid).mode,'all_cond')
0211                     conds     = PRT.group(gid).subject(sid).modality(mid).design.conds;
0212                     n_vols_s  = size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0213                     
0214                     % now loop over conditions
0215                     for cid = 1:length(conds)    % condition
0216                         scans     = PRT.group(gid).subject(sid).modality(mid).design.conds(cid).scans;
0217                         blocks    = PRT.group(gid).subject(sid).modality(mid).design.conds(cid).blocks;
0218                         n_vol_s_c = length(scans);
0219                         
0220                         % configure indices
0221                         sample_range = (1:n_vol_s_c)+max(sample_range);
0222                         PRT.fs(fid).id_mat(sample_range,1) = gid;
0223                         PRT.fs(fid).id_mat(sample_range,2) = sid;
0224                         PRT.fs(fid).id_mat(sample_range,3) = mid;
0225                         PRT.fs(fid).id_mat(sample_range,4) = cid;
0226                         PRT.fs(fid).id_mat(sample_range,5) = blocks;
0227                         %PRT.fs(fid).id_mat(sample_range,6) = 1:length(sample_range);
0228                         PRT.fs(fid).id_mat(sample_range,6) = scans;
0229                         
0230                         %configure indices for the file array
0231                         sctoadd=scans+indm(m);
0232                         PRT.fs(fid).fas.ifa(sample_range)=sctoadd';
0233                         PRT.fs(fid).fas.im(sample_range)=mid*ones(n_vol_s_c,1);
0234                     end
0235                     %configure indices for the file array
0236                     indm(m)=n_vols_s+max(indm(m));
0237                 end
0238                 szm(m)=szm(m)+size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0239             end  % modality
0240         end  % subject
0241     end  % group
0242     
0243     %initialize the file arrays if they do not exist already or if the
0244     %detrending parameters were modified
0245     if ~isfield(PRT,'fas');
0246         % initialise all modalities (not just those we're working on)
0247         for m = 1:length(PRT.masks)
0248             
0249             PRT.fas(m)=struct('mod_name',[],'dat',[],'detrend',[],'param_dt',[],'hdr',[]);
0250             PRT.fas(m).mod_name = PRT.masks(m).mod_name;
0251         end
0252     end
0253     tocomp=zeros(1,length(in.mod));
0254     prt_dir=fileparts(in.fname);
0255     for i=1:n_mods
0256         % check whether we need to recreate the file array
0257         if mids(i)>length(PRT.fas) ||...
0258                 isempty(PRT.fas(mids(i)).dat) || exist(PRT.fas(mids(i)).dat.fname)==0 ||...
0259                 PRT.fas(mids(i)).detrend ~= in.mod(mids(i)).detrend  || ...
0260                 (isempty(PRT.fas(mids(i)).param_dt) && ~isempty(in.mod(mids(i)).param_dt)) || ...
0261                 (~isempty(PRT.fas(mids(i)).param_dt) && isempty(in.mod(mids(i)).param_dt)) || ...
0262                 ((~isempty(PRT.fas(mids(i)).param_dt) && ~isempty(in.mod(mids(i)).param_dt)) && ...
0263                 PRT.fas(mids(i)).param_dt~=in.mod(mids(i)).param_dt)
0264             
0265             if mids(i)>length(PRT.fas) || isempty(PRT.fas(mids(i)).dat)
0266                 disp(['File array does not exist for modality ''',...
0267                     char(in.mod(mids(i)).mod_name),'''. Creating...'])
0268             elseif PRT.fas(mids(i)).detrend ~= in.mod(mids(i)).detrend ...
0269                     && any(strcmpi(fieldnames(PRT.fas(mids(i)).dat),'fname')) ...
0270                     && exist(PRT.fas(mids(i)).dat.fname,'file')
0271                 
0272                 warning('prt_init_fs:overwriteFileArray',...
0273                     ['File array already exists for modality ''',...
0274                     char(in.mod(mids(i)).mod_name),''', but parameters ',...
0275                     'have changed. Re-creating ...']);
0276                 
0277                 delete(PRT.fas(mids(i)).dat.fname);
0278             end
0279             
0280             tocomp(mids(i))=1;
0281             %PRT.fas(mids(i)).mod_name = in.mod(mids(i)).mod_name;
0282             PRT.fas(mids(i)).detrend = in.mod(mids(i)).detrend;
0283             PRT.fas(mids(i)).param_dt = in.mod(mids(i)).param_dt;
0284             PRT.fas(mids(i)).hdr = headers{i};
0285             PRT.fas(mids(i)).idfeat_img = PRT.fs(fid).modality(i).feat_idx_img;                % index of voxels in the full image (nifti)
0286             datname=[prt_dir,filesep,'Feature_set_',char(in.mod(mids(i)).mod_name),'.dat'];
0287             PRT.fas(mids(i)).dat = file_array(...
0288                 datname, ...                 % fname     - filename
0289                 [szm(i),n_vox],...           % dim       - dimensions (default = [0 0] )
0290                 spm_type('float64'), ...  % dtype     - datatype   (default = 'float')
0291                 0, ...                       % offset    - offset into file (default = 0)
0292                 1);                          % scl_slope - scalefactor (default = 1)
0293         else
0294             disp(['Using existing file array for modality ''', ...
0295                 char(in.mod(mids(i)).mod_name),'''.'])
0296         end
0297         
0298         % check that the input .mat for the scaling have the right size
0299         if in.mod(mids(i)).normalise==2
0300             try
0301                 load(in.mod(mids(i)).matnorm);
0302             catch
0303                 error('prt_prepare_data:ScalingMatUnloadable',...
0304                     'Could not load the .mat file containing the scaling')
0305             end
0306             try
0307                 szin=length(scaling);
0308             catch
0309                 error('prt_prepare_data:ScalingNotinFile',...
0310                     'This file does not contain the "scaling" field required')
0311             end
0312             if szin~=szm(i)
0313                 error('prt_prepare_data:Scalingdimensionwrong',...
0314                     'The dimension of the .mat file does not correspond to the number of scans in that modality')
0315             end
0316             PRT.fs(fid).modality(i).normalise.type=2;
0317             PRT.fs(fid).modality(i).normalise.scaling=reshape(scaling,1,szm(i));
0318         elseif in.mod(mids(i)).normalise==1
0319             PRT.fs(fid).modality(i).normalise.type=1;
0320         else
0321             PRT.fs(fid).modality(i).normalise.type=0;
0322         end
0323     end
0324     
0325     PRT.fs(fid).modality=rmfield(PRT.fs(fid).modality,'feat_idx_img');
0326 end

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