0001 function [PRT,Phi] = prt_fs_modality(PRT,in, flag, addin)
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
0037
0038 prt_dir = regexprep(in.fname,'PRT.mat', '');
0039 def=prt_get_defaults('fs');
0040 fid = in.fid;
0041
0042
0043
0044 if nargin>=3 && flag
0045 if ~isfield(addin,'ID')
0046 ID= PRT.fs(fid).id_mat;
0047 n_mods=length(in.mod);
0048 mids=[];
0049 for i=1:n_mods
0050 if ~isempty(in.mod(i).mod_name)
0051 mids = [mids, i];
0052 end
0053 end
0054 else
0055 ID = addin.ID;
0056 mids = unique(ID(:,3));
0057 end
0058 n_mods=length(mids);
0059 if isfield(addin,'idvox_fas')
0060 idvox = addin.idvox_fas;
0061 else
0062 if in.tocomp(mids(1))
0063 idvox = 1:PRT.fas(mids(1)).dat.dim(2);
0064 else
0065 if ~isempty(PRT.fs(fid).modality(mids(1)).idfeat_fas)
0066 idvox = PRT.fs(fid).modality(mids(1)).idfeat_fas;
0067 else
0068 idvox = 1:PRT.fas(mids(1)).dat.dim(2);
0069 end
0070 end
0071 end
0072 n_vox = numel(idvox);
0073 nfa = 0;
0074 else
0075
0076 n_mods=length(in.mod);
0077 mids=[];
0078 for i=1:n_mods
0079 if ~isempty(in.mod(i).mod_name)
0080 mids = [mids, i];
0081 end
0082 end
0083 n_mods=length(mids);
0084 ID= PRT.fs(fid).id_mat;
0085 nfa = [];
0086 for m = 1:n_mods
0087 nfa = [nfa, PRT.fas(mids(m)).dat.dim(1)];
0088 if in.tocomp(mids(m))
0089 n_vox = PRT.fas(mids(m)).dat.dim(2);
0090 else
0091 if ~isempty(PRT.fs(fid).modality(m).idfeat_fas)
0092 idvox = PRT.fs(fid).modality(m).idfeat_fas;
0093 else
0094 idvox = 1:PRT.fas(mids(m)).dat.dim(2);
0095 end
0096 n_vox = numel(idvox);
0097 end
0098 end
0099 end
0100
0101
0102
0103
0104
0105
0106 n = size(ID,1);
0107 Phi = zeros(n);
0108
0109 mem = def.mem_limit;
0110 block_size = floor(mem/(8*3)/max([nfa, n]));
0111 n_block = ceil(n_vox/block_size);
0112
0113 bstart = 1; bend = min(block_size,n_vox);
0114 if nargin<3 || ~flag
0115 h = waitbar(0,'Please wait while preparing feature set');
0116 step=1;
0117 end
0118 for b = 1:n_block
0119 if nargin<3 || ~flag
0120 disp ([' > preparing block: ', num2str(b),' of ',num2str(n_block),' ...'])
0121 end
0122 vox_range = bstart:bend;
0123 block_size = length(vox_range);
0124 kern_vols = zeros(block_size,n);
0125 for m=1:n_mods
0126 mid=mids(m);
0127
0128
0129
0130
0131 ind_ddmask = PRT.fas(mid).idfeat_img(vox_range);
0132
0133
0134 if ~isempty(in.precmask{m})
0135 prec_mask = prt_load_blocks(in.precmask{m},ind_ddmask);
0136 else
0137 prec_mask = ones(block_size,1);
0138 end
0139
0140 indm = find(PRT.fs(fid).fas.im==mid);
0141 ifa = PRT.fs(fid).fas.ifa(indm);
0142 indm = find(ID(:,3) == mid);
0143
0144
0145
0146 if in.tocomp(mid)
0147 n_groups = length(PRT.group);
0148 sample_range=0;
0149 nfa=PRT.fas(mid).dat.dim(1);
0150 datapr=zeros(block_size,nfa);
0151
0152
0153 for gid = 1:n_groups
0154 for sid = 1:length(PRT.group(gid).subject)
0155 n_vols_s = size(PRT.group(gid).subject(sid).modality(mid).scans,1);
0156 sample_range = (1:n_vols_s)+max(sample_range);
0157 fname = PRT.group(gid).subject(sid).modality(mid).scans;
0158 datapr(:,sample_range) = prt_load_blocks(fname,ind_ddmask);
0159
0160 [inan,jnan] = find(isnan(datapr(:,sample_range)));
0161 if ~isempty(inan)
0162 disp('Warning: NaNs found in loaded data')
0163 disp('Consider updating 1st level mask for better performance')
0164 for inn=1:length(inan)
0165 datapr(inan(inn),sample_range(jnan(inn))) = 0;
0166 end
0167 end
0168
0169
0170 if in.mod(mid).detrend ~= 0
0171 if isfield(PRT.group(gid).subject(sid).modality(mid).design,'TR')
0172 TR = PRT.group(gid).subject(sid).modality(mid).design.TR;
0173 else
0174 try
0175 TR = PRT.group(gid).subject(sid).modality(mid).TR;
0176 catch
0177 error('detrend:TRnotfound','No TR in data, suggesting that detrend is not necessary')
0178 end
0179 end
0180 switch in.mod(mid).detrend
0181 case 1
0182 C = poly_regressor(length(sample_range), ...
0183 in.mod(mid).param_dt);
0184 case 2
0185 C = dct_regressor(length(sample_range), ...
0186 in.mod(mid).param_dt,TR);
0187 end
0188 R = eye(length(sample_range)) - C*pinv(C);
0189 datapr(:,sample_range) = datapr(:,sample_range)*R';
0190 end
0191 end
0192 end
0193
0194 if b==1
0195
0196 namedat=['Feature_set_',char(in.mod(mid).mod_name),'.dat'];
0197 fpd_clean(m) = fopen(fullfile(prt_dir,namedat), 'w','ieee-le');
0198 fwrite(fpd_clean(m), datapr', 'float32',0,'ieee-le');
0199 else
0200
0201 fwrite(fpd_clean(m), datapr', 'float32',0,'ieee-le');
0202 end
0203
0204
0205 kern_vols(:,indm) = datapr(:,ifa).* ...
0206 repmat(prec_mask~=0,1,length(ifa));
0207
0208 if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0209 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0210 repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0211 end
0212 clear datapr
0213 else
0214 kern_vols(:,indm) = (PRT.fas(mid).dat(ifa,idvox(vox_range)))';
0215
0216 if ~isempty(PRT.fs(fid).modality(m).normalise.scaling)
0217 kern_vols(:,indm) = kern_vols(:,indm)./ ...
0218 repmat(PRT.fs(fid).modality(m).normalise.scaling,block_size,1);
0219 end
0220
0221 end
0222 if nargin<3 || ~flag
0223 waitbar(step/ (n_block*n_mods),h);
0224 step=step+1;
0225 end
0226 end
0227 if size(kern_vols,2)>6e3
0228
0229
0230
0231 for ic=1:size(kern_vols,2)
0232 Phi(:,ic) = Phi(:,ic) + kern_vols' * kern_vols(:,ic);
0233 end
0234 else
0235 Phi = Phi + (kern_vols' * kern_vols);
0236 end
0237 bstart = bend+1; bend = min(bstart+block_size-1,n_vox);
0238 clear block_mask kern_vols
0239 end
0240 if nargin<3 || ~flag
0241 close(h)
0242 end
0243
0244
0245 if exist('fpd_clean','var')
0246 for ii=1:numel(fpd_clean)
0247 fclose(fpd_clean(ii));
0248 end
0249 end
0250
0251
0252
0253
0254
0255 function c = poly_regressor(n,order)
0256
0257
0258
0259 basis = repmat([1:n]',[1 order]);
0260 o = repmat([1:order],[n 1]);
0261 c = [ones(n,1) basis.^o];
0262 return
0263
0264 function c=dct_regressor(n,cut_off,TR)
0265
0266
0267
0268
0269 if cut_off<0
0270 error('cut off cannot be negative')
0271 end
0272
0273 T = n*TR;
0274 order = floor((T/cut_off)*2)+1;
0275 c = spm_dctmtx(n,order);
0276 return