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
 in.mod(m).multroi    1 if one kernel per region required
 in.mod(m).atlasroi   name of the atlas to build one kernel per region

 in.flag_mm:   Perform multi-kernel learning (1) or not (0)? If yes, the
 kernel is saved as a cell vector, with one kernel per modality

 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 % in.mod(m).multroi    1 if one kernel per region required
0018 % in.mod(m).atlasroi   name of the atlas to build one kernel per region
0019 %
0020 % in.flag_mm:   Perform multi-kernel learning (1) or not (0)? If yes, the
0021 % kernel is saved as a cell vector, with one kernel per modality
0022 %
0023 % Outputs:
0024 % --------
0025 % Calls prt_init_fs to populate basic fields in PRT.fs(f)...
0026 % Writes PRT.mat
0027 % Writes the kernel matrix to the path indicated by in.fs_name
0028 %__________________________________________________________________________
0029 % Copyright (C) 2011 Machine Learning & Neuroimaging Laboratory
0030 
0031 % Written by A. Marquand and J. Schrouff
0032 % $Id$
0033 
0034 % Configure some variables and get defaults
0035 % -------------------------------------------------------------------------
0036 prt_dir  = regexprep(in.fname,'PRT.mat', ''); % or: fileparts(fname);
0037 
0038 % get the index of the modalities for which the user wants a kernel/data
0039 n_mods=length(in.mod);
0040 mids=[];
0041 for i=1:n_mods
0042     if ~isempty(in.mod(i).mod_name)
0043         mids = [mids, i];
0044     end
0045 end
0046 n_mods=length(mids);
0047 
0048 % Load mask(s) and resize if necessary
0049 %force second-level masking by atlas if provided
0050 for i=1:length(in.mod)
0051     if isfield(in.mod(i),'atlasroi')
0052         if ~isempty(in.mod(i).atlasroi) && isempty(in.mod(i).mask)
0053             in.mod(i).mask = in.mod(i).atlasroi;    %option only available for one modality
0054         end
0055     end
0056 end
0057 
0058 [mask,precmask,headers,PRT,ratl] = load_masks(PRT, prt_dir, in, mids);
0059 
0060 % Initialize the file arrays, kernel and feature set parameters
0061 [fid,PRT,tocomp] = prt_init_fs(PRT,in,mids,mask,precmask,headers);
0062 
0063 in.tocomp = tocomp;
0064 in.precmask = precmask;
0065 in.fid = fid;
0066 
0067 % Build the feature set and kernel
0068 %--------------------------------------------------------------------------
0069 
0070 Phi = [];
0071 igd = [];
0072 PRT.fs(fid).multkernelROI = 0; % Multiple kernels with an atlas
0073 PRT.fs(fid).multkernel = 0;    % Multiple kernels from different modalities
0074 nroi = 0;
0075 
0076 if in.flag_mm   % One kernel per modality so need to treat them independently
0077     for i = 1:n_mods  % multiple modalities
0078         % For each modality, get the corresponding ID mat and sample index
0079         idtk = PRT.fs(fid).id_mat(:,3) == mids(i);
0080         nimm = length(unique(PRT.fs(fid).id_mat(:,3) == mids(i)));
0081         
0082         %check that modalities have the same dimensions in terms of samples
0083         nim1 =length(unique(PRT.fs(fid).id_mat(:,3) == mids(1)));
0084         if nimm~= nim1
0085             error('prt_fs:MultKernMod_DifIm',...
0086                 'Modalities should have the same number of samples to be considered for MKL')
0087         end
0088         addin.ID = PRT.fs(fid).id_mat(idtk,:);
0089         
0090         % If second-level, i.e. atlas-based kernels as well
0091         if isfield(in.mod(mids(i)),'multroi') ...
0092                 && in.mod(mids(i)).multroi
0093             atl=spm_vol(ratl{i});
0094             %Initialize all fields and compute the feature sets if needed
0095             if any(in.tocomp)
0096                 [PRT] = prt_fs_modality(PRT,in,1,addin);
0097             end
0098             nroiprev = nroi;
0099             [PRT,Phim,igk,nroi] = prt_compute_ROI_kernels(PRT,in,fid,mids(i),atl,addin,i);
0100             PRT.fs(fid).atlas_name = ratl;            
0101             igd = [igd,igk+nroiprev];
0102         else
0103             [PRT,Phim] = prt_fs_modality(PRT,in,1,addin);
0104             [d1,idmax] = max(Phim);
0105             [d1,idmin] = min(Phim);
0106             min_max = find(idmax==idmin);
0107             if isempty(min_max) || unique(Phim(:,min_max))~=0 %Kernel does not contain a whole line of zeros
0108                 igd = [igd,i];
0109             else
0110                 beep
0111                 disp('No overlap between data and mask/atlas for at least one sample')
0112                 disp(['Kernel ',num2str(i),' will be removed from further analysis'])
0113             end
0114             Phim = {Phim};
0115             PRT.fs(fid).atlas_name = {};
0116         end
0117         Phi = [Phi, Phim];
0118     end
0119     %post-hoc: the ID mat should be the same for all modalities involved,
0120     %so only the first one will be saved
0121     if ~isempty(igd) && isempty(PRT.fs(fid).atlas_name)
0122         PRT.fs(fid).modality = PRT.fs(fid).modality(igd);
0123     end
0124     indm=PRT.fs(fid).fas.im==1;
0125     PRT.fs(fid).id_mat=PRT.fs(fid).id_mat(indm,:);
0126     PRT.fs(fid).multkernel = 1;
0127 else
0128     % Concatenate the modalities in samples or only one modality
0129     % Does the same as before without playing with the ID matrix
0130     % If second-level, i.e. atlas-based kernels as well
0131     
0132     % First check that all concatenated modalitlies have the same flag and
0133     % atlas, or no atlas
0134     if n_mods>1
0135         mult = zeros(n_mods,1);
0136         for i = 1:n_mods
0137             if isfield(in.mod(mids(i)),'multroi') ...
0138                     && in.mod(mids(i)).multroi
0139                 mult(i) = 1;
0140                 atl = ratl{i};
0141                 if i ==1
0142                     atlmod = atl;
0143                 end
0144                 if ~strcmpi(atl,atlmod)
0145                     error('prt_fs:ConcatenatingWithDifferentAtlases',...
0146                         'Concatenated multiple modalities should have the same atlas');
0147                 end
0148             end
0149         end
0150         if length(unique(mult))~=1
0151             error('prt_fs:ConcatenateModalitiesWithandWithoutAtlas',...
0152                 'Modalities cannot be concatenated unless they all have no or the same atlas')
0153         end
0154     end
0155     
0156     atl=spm_vol(ratl{1});
0157     addin = struct();
0158     %Initialize all fields and compute the feature sets if needed
0159     if isfield(in.mod(mids(1)),'multroi') ...
0160             && in.mod(mids(1)).multroi
0161         if any(in.tocomp)
0162             [PRT] = prt_fs_modality(PRT,in,0,[]);
0163         end
0164         [PRT,Phim,igd] = prt_compute_ROI_kernels(PRT,in,fid,mids(1),atl,addin,1);
0165         PRT.fs(fid).atlas_name{1} = ratl{1};
0166         Phi = Phim;
0167     else % Simply concatenate the modalities in samples
0168         [PRT,Phim] = prt_fs_modality(PRT,in,0,[]);
0169         PRT.fs(fid).multkernel = 0;
0170         PRT.fs(fid).atlas_name = {};
0171         [d1,idmax] = max(Phim);
0172         [d1,idmin] = min(Phim);
0173         min_max = find(idmax==idmin);
0174         if isempty(min_max) || unique(Phim(:,min_max))~=0 %Kernel does not contain a whole line of zeros
0175             igd = 1;
0176         else
0177             error('prt_fs:NoDataInMask',...
0178                 'No overlap between data and mask/atlas for at least one sample, cannot create kernel')
0179         end
0180         Phi{1}=Phim;
0181     end
0182 end
0183 
0184 
0185 clear Phim
0186 
0187 if isempty(igd)
0188     error('prt_fs:NoDataInMask',...
0189         'No overlap between data and mask/atlas for at least one sample, cannot create kernel')
0190 else
0191     Phi = Phi(igd);
0192     PRT.fs(fid).igood_kerns = igd;
0193 end
0194 
0195 
0196 % Save kernel and function output
0197 % -------------------------------------------------------------------------
0198 outfile = in.fname;
0199 disp('Saving feature set to: PRT.mat.......>>')
0200 disp(['Saving kernel to: ',in.fs_name,'.mat.......>>'])
0201 fs_file = [prt_dir,in.fs_name];
0202 if spm_check_version('MATLAB','7') < 0
0203     save(outfile,'-V6','PRT');
0204     save(fs_file,'-V6','Phi');
0205 else
0206     save(outfile,'PRT');
0207     save(fs_file,'Phi');
0208 end
0209 disp('Done.')
0210 
0211 %--------------------------------------------------------------------------
0212 %------------------------- Private function -------------------------------
0213 %--------------------------------------------------------------------------
0214 
0215 function [mask, precmask, headers,PRT, ratl] = load_masks(PRT, prt_dir, in, mids)
0216 % function to load the mask for each modality
0217 % -------------------------------------------
0218 n_mods   = length(mids);
0219 mask     = cell(1,n_mods);
0220 precmask = cell(1,n_mods);
0221 headers  = cell(1,n_mods);
0222 ratl     = cell(1,n_mods);
0223 for m = 1:n_mods
0224     mid = mids(m);
0225     
0226     % get mask for the within-brain voxels (from data and design)
0227     ddmask = PRT.masks(mid).fname;
0228     try
0229         M = nifti(ddmask);
0230     catch %#ok<*CTCH>
0231         error('prt_fs:CouldNotLoadFile',...
0232             'Could not load mask file');
0233     end
0234     
0235     % get mask for the kernel if one was specified
0236     mfile = in.mod(mid).mask;
0237     if ~isempty(mfile) %&&  mfile ~= 0
0238         try
0239             precM = spm_vol(char(mfile));
0240         catch
0241             error('prt_fs:CouldNotLoadFile',...
0242                 'Could not load mask file for preprocessing');
0243         end
0244     end
0245     
0246     % get atlas for the ROI based kernel if one was specified
0247     if isfield(in.mod(mid),'atlasroi')
0248         alfile = in.mod(mid).atlasroi;
0249         if ~isempty(alfile) %&&  mfile ~= 0
0250             try
0251                 precA = spm_vol(char(alfile));
0252             catch
0253                 error('prt_fs:CouldNotLoadFile',...
0254                     'Could not load mask file for preprocessing');
0255             end
0256         end
0257     else
0258         alfile=[];
0259     end
0260     
0261     % get header of the first scan of that modality
0262     if isfield(PRT,'fas') && mid<=length(PRT.fas) && ...
0263             ~isempty(PRT.fas(mid).dat)
0264         N = PRT.fas(mid).hdr;
0265     else
0266         N = spm_vol(PRT.group(1).subject(1).modality(mid).scans(1,:));
0267     end
0268     headers{m}=N;
0269     
0270     % compute voxel dimensions and check for equality if n_mod > 1
0271     if m == 1
0272         n_vox = prod(N.dim(1:3));
0273     elseif n_mods > 1 && n_vox ~= prod(N.dim(1:3))
0274         error('prt_fs:multipleModatlitiesVariableFeatures',...
0275             'Multiple modalities specified, but have variable numbers of features');
0276     end
0277     
0278     % resize the different masks if needed
0279     if N.dim(3)==1, Npdim = N.dim(1:2); else Npdim = N.dim; end % handling case of 2D images
0280     if any(size(M.dat(:,:,:,1)) ~= Npdim)
0281         warning('prt_fs:maskAndImagesDifferentDim',...
0282             'Mask has different dimensions to the image files. Resizing...');
0283         
0284         V2 = spm_vol(char(ddmask));
0285         % reslicing V2
0286         fl_res = struct('mean',false,'interp',0,'which',1,'prefix','tmp_');
0287         spm_reslice([N V2],fl_res)
0288         % now renaming the file
0289         [V2_pth,V2_fn,V2_ext] = spm_fileparts(V2.fname);
0290         rV2_fn = [fl_res.prefix,V2_fn];
0291         if strcmp(V2_ext,'.nii')
0292             % turn .nii into .img/.hdr image!
0293             V_in = spm_vol(fullfile(V2_pth,[rV2_fn,'.nii']));
0294             V_out = V_in; V_out.fname = fullfile(V2_pth,[rV2_fn,'.img']);
0295             spm_imcalc(V_in,V_out,'i1');
0296         end
0297         mfile_new = ['updated_1stlevel_mask_m',num2str(mid)];
0298         movefile(fullfile(V2_pth,[rV2_fn,'.img']), ...
0299             fullfile(prt_dir,[mfile_new,'.img']));
0300         movefile(fullfile(V2_pth,[rV2_fn,'.hdr']), ...
0301             fullfile(prt_dir,[mfile_new,'.hdr']));
0302         PRT.masks(mid).fname = fullfile(prt_dir,[mfile_new,'.img']);
0303         mask{m} = PRT.masks(mid).fname;
0304     else
0305         mask{m} = ddmask;
0306     end
0307     if ~isempty(mfile) && any((precM.dim~= N.dim)) % && mfile ~= 0
0308         warning('prt_fs:maskAndImagesDifferentDim',...
0309             'Preprocessing mask has different dimensions to the image files. Resizing...');
0310         V2 = spm_vol(char(mfile));
0311         % reslicing V2
0312         fl_res = struct('mean',false,'interp',0,'which',1,'prefix','tmp_');
0313         spm_reslice([N V2],fl_res)
0314         % now renaming the file
0315         [V2_pth,V2_fn,V2_ext] = spm_fileparts(V2.fname);
0316         rV2_fn = [fl_res.prefix,V2_fn];
0317         if strcmp(V2_ext,'.nii')
0318             % turn .nii into .img/.hdr image!
0319             V_in = spm_vol(fullfile(V2_pth,[rV2_fn,'.nii']));
0320             V_out = V_in; V_out.fname = fullfile(V2_pth,[rV2_fn,'.img']);
0321             spm_imcalc(V_in,V_out,'i1');
0322         end
0323         % if more than one 2nd level mask to resize
0324         nummask = 1;
0325         while exist(fullfile( ...
0326                 prt_dir,['updated_2ndlevel_mask_m',num2str(mid),'_',...
0327                 num2str(nummask),'.img']),'file')
0328             nummask = nummask+1;
0329         end
0330         mfile_new = ['updated_2ndlevel_mask_m',num2str(mid),...
0331             '_',num2str(nummask)];
0332         movefile(fullfile(V2_pth,[rV2_fn,'.img']), ...
0333             fullfile(prt_dir,[mfile_new,'.img']));
0334         movefile(fullfile(V2_pth,[rV2_fn,'.hdr']), ...
0335             fullfile(prt_dir,[mfile_new,'.hdr']));
0336         precmask{m} = fullfile(prt_dir,[mfile_new,'.img']);
0337     else
0338         precmask{m} = mfile;
0339     end
0340     if ~isempty(alfile) && any((precA.dim~= N.dim))
0341         warning('prt_fs:atlasAndImagesDifferentDim',...
0342             'Atlas has different dimensions to the image files. Resizing...');
0343         V2 = spm_vol(char(alfile));
0344         % reslicing V2
0345         fl_res = struct('mean',false,'interp',0,'which',1,'prefix','tmp_');
0346         spm_reslice([N V2],fl_res)
0347         % now renaming the file
0348         [V2_pth,V2_fn,V2_ext] = spm_fileparts(V2.fname);
0349         rV2_fn = [fl_res.prefix,V2_fn];
0350         if strcmp(V2_ext,'.nii')
0351             % turn .nii into .img/.hdr image!
0352             V_in = spm_vol(fullfile(V2_pth,[rV2_fn,'.nii']));
0353             V_out = V_in; V_out.fname = fullfile(V2_pth,[rV2_fn,'.img']);
0354             spm_imcalc(V_in,V_out,'i1');
0355         end
0356         alfile_new = ['updated_atlas_',V2_fn];
0357         movefile(fullfile(V2_pth,[rV2_fn,'.img']), ...
0358             fullfile(prt_dir,[alfile_new,'.img']));
0359         movefile(fullfile(V2_pth,[rV2_fn,'.hdr']), ...
0360             fullfile(prt_dir,[alfile_new,'.hdr']));
0361         ratl{m} = fullfile(prt_dir,[alfile_new,'.img']);
0362     else
0363         ratl{m} = alfile;
0364     end
0365     clear M N precM V1 V2 mfile mfile_new
0366 end
0367 
0368 
0369 function [PRT,Phi,igd,nroi] = prt_compute_ROI_kernels(PRT,in,fid,mids,atl,addin,ind_mod)
0370 % function to load the mask for each modality
0371 % -------------------------------------------
0372 
0373 %For each region, get the indexes of the voxels in the 2nd level mask
0374 h=spm_read_vols(atl);
0375 if ~isempty(PRT.fs(fid).modality(ind_mod).idfeat_fas)
0376     idt = PRT.fs(fid).modality(ind_mod).idfeat_fas; %indexes of voxels in second-level mask
0377 else
0378     idt = 1:length(PRT.fas(mids).idfeat_img);
0379 end
0380 idm1 = PRT.fas(mids).idfeat_img(idt); %indexes of voxels in first level mask
0381 interh = h(idm1);
0382 roi = unique(interh(interh>0));
0383 nroi = length(roi);
0384 Phi=cell(nroi,1);
0385 in.tocomp = zeros(1,length(in.tocomp));
0386 %For each region, compute kernel and save the indexes in the image for
0387 %further computation of the weights
0388 PRT.fs(fid).modality(ind_mod).idfeat_img = cell(nroi,1);
0389 igd = []; %indexes of non 0 kernels
0390 for i=1:nroi
0391     disp ([' > Computing kernel: ', num2str(i),' of ',num2str(nroi),' ...'])
0392     addin.idvox_fas = idt(interh == roi(i));
0393     [PRT,Phim] = prt_fs_modality(PRT,in,1,addin);
0394     [d1,idmax] = max(Phim);
0395     [d1,idmin] = min(Phim);
0396     min_max = find(idmax==idmin);
0397     if isempty(min_max) || unique(Phim(:,min_max))~=0 %Kernel does not contain a whole line of zeros
0398         igd = [igd,i];
0399     else
0400         beep
0401         disp('No overlap between data and mask/atlas for at least one sample')
0402         disp(['Region ',num2str(i),' will be removed from further analysis'])
0403     end
0404     Phi{i}=Phim;
0405     %         idts = idt(interh == roi(i));
0406     PRT.fs(fid).modality(ind_mod).idfeat_img{i} = find(interh == roi(i)) ;
0407 end
0408 PRT.fs(fid).multkernelROI = 1;
0409 if ~isempty(igd)
0410     PRT.fs(fid).modality(ind_mod).idfeat_img = PRT.fs(fid).modality(ind_mod).idfeat_img(igd);
0411     PRT.fs(fid).modality(ind_mod).num_ROI = roi(igd);
0412 end
0413 
0414 
0415 return

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