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
0032
0033
0034
0035
0036 prt_dir = regexprep(in.fname,'PRT.mat', '');
0037
0038
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
0049
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;
0054 end
0055 end
0056 end
0057
0058 [mask,precmask,headers,PRT,ratl] = load_masks(PRT, prt_dir, in, mids);
0059
0060
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
0068
0069
0070 Phi = [];
0071 igd = [];
0072 PRT.fs(fid).multkernelROI = 0;
0073 PRT.fs(fid).multkernel = 0;
0074 nroi = 0;
0075
0076 if in.flag_mm
0077 for i = 1:n_mods
0078
0079 idtk = PRT.fs(fid).id_mat(:,3) == mids(i);
0080 nimm = length(unique(PRT.fs(fid).id_mat(:,3) == mids(i)));
0081
0082
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
0091 if isfield(in.mod(mids(i)),'multroi') ...
0092 && in.mod(mids(i)).multroi
0093 atl=spm_vol(ratl{i});
0094
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
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
0120
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
0129
0130
0131
0132
0133
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
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
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
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
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
0213
0214
0215 function [mask, precmask, headers,PRT, ratl] = load_masks(PRT, prt_dir, in, mids)
0216
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
0227 ddmask = PRT.masks(mid).fname;
0228 try
0229 M = nifti(ddmask);
0230 catch
0231 error('prt_fs:CouldNotLoadFile',...
0232 'Could not load mask file');
0233 end
0234
0235
0236 mfile = in.mod(mid).mask;
0237 if ~isempty(mfile)
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
0247 if isfield(in.mod(mid),'atlasroi')
0248 alfile = in.mod(mid).atlasroi;
0249 if ~isempty(alfile)
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
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
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
0279 if N.dim(3)==1, Npdim = N.dim(1:2); else Npdim = N.dim; end
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
0286 fl_res = struct('mean',false,'interp',0,'which',1,'prefix','tmp_');
0287 spm_reslice([N V2],fl_res)
0288
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
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))
0308 warning('prt_fs:maskAndImagesDifferentDim',...
0309 'Preprocessing mask has different dimensions to the image files. Resizing...');
0310 V2 = spm_vol(char(mfile));
0311
0312 fl_res = struct('mean',false,'interp',0,'which',1,'prefix','tmp_');
0313 spm_reslice([N V2],fl_res)
0314
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
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
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
0345 fl_res = struct('mean',false,'interp',0,'which',1,'prefix','tmp_');
0346 spm_reslice([N V2],fl_res)
0347
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
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
0371
0372
0373
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;
0377 else
0378 idt = 1:length(PRT.fas(mids).idfeat_img);
0379 end
0380 idm1 = PRT.fas(mids).idfeat_img(idt);
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
0387
0388 PRT.fs(fid).modality(ind_mod).idfeat_img = cell(nroi,1);
0389 igd = [];
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
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
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