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 ind_ddmask = PRT.fas(mid).idfeat_img(vox_range);
0082
0083
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
0090 indm = find(PRT.fs(fid).fas.im==mid);
0091 ifa = PRT.fs(fid).fas.ifa(indm);
0092
0093
0094
0095 if tocomp(mid)
0096 sample_range=0;
0097 nfa=PRT.fas(mid).dat.dim(1);
0098 datapr=zeros(block_size,nfa);
0099
0100
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
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
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
0138 namedat=['Feature_set_',char(in.mod(mid).mod_name),'.dat'];
0139 fpd_clean(m) = fopen(fullfile(prt_dir,namedat), 'a','ieee-le.l64');
0140 fwrite(fpd_clean(m), datapr', 'float64',0,'ieee-le.l64');
0141 else
0142
0143 fwrite(fpd_clean(m), datapr', 'float64',0,'ieee-le.l64');
0144 end
0145
0146
0147 kern_vols(:,indm) = datapr(:,ifa).* ...
0148 repmat(prec_mask,1,length(ifa));
0149
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
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
0170
0171
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
0184 if exist('fpd_clean','var')
0185 for ii=1:numel(fpd_clean)
0186 fclose(fpd_clean(ii));
0187 end
0188 end
0189
0190
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
0207
0208 function [mask, precmask, headers,PRT] = load_masks(PRT, prt_dir, in, mids)
0209
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
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
0228 mfile=in.mod(mid).mask;
0229 if ~isempty(mfile)
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
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
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
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))
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
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
0296
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
0305
0306
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