0001 function [outfile] = prt_fs(PRT,in)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 prt_dir = regexprep(in.fname,'PRT.mat', '');
0032 n_groups = length(PRT.group);
0033
0034
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
0046 [mask,precmask,headers,PRT] = load_masks(PRT, prt_dir, in,mids);
0047
0048
0049 [fid,PRT,tocomp] = prt_init_fs(PRT,in,mids,mask,precmask,headers);
0050
0051
0052
0053
0054 n = length(PRT.fs(fid).id_mat);
0055 Phi = zeros(n);
0056
0057
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]));
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
0079
0080
0081
0082 ind_ddmask = PRT.fas(mid).idfeat_img(vox_range);
0083
0084
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
0091 indm=find(PRT.fs(fid).fas.im==mid);
0092 ifa=PRT.fs(fid).fas.ifa(indm);
0093
0094
0095
0096
0097 if tocomp(mid)
0098 sample_range=0;
0099 nfa=PRT.fas(mid).dat.dim(1);
0100 datapr=zeros(nfa,block_size);
0101
0102
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
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
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
0138 namedat=['Feature_set_',char(in.mod(mid).mod_name),'.dat'];
0139 fpd_clean = fopen(fullfile(prt_dir,namedat), 'a','ieee-le.l64');
0140 if b==1
0141
0142 fwrite(fpd_clean, datapr, 'float64',0,'ieee-le.l64');
0143 fclose(fpd_clean);
0144 else
0145
0146 fseek(fpd_clean,0,'eof');
0147 fwrite(fpd_clean, datapr, 'float64',0,'ieee-le.l64');
0148 fclose(fpd_clean);
0149 end
0150
0151
0152 kern_vols(:,indm) = datapr(ifa,:)'.* ...
0153 repmat(prec_mask,1,length(ifa));
0154
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
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
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
0196
0197 function [mask, precmask, headers,PRT] = load_masks(PRT, prt_dir, in, mids)
0198
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
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
0217 mfile=in.mod(mid).mask;
0218 if ~isempty(mfile)
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
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
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
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))
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
0277
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
0286
0287
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