0001 function [NW_roi,idfeatroi]=prt_build_region_weights(weight_fname,atlas_fname,build_im,comp_perm)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 if nargin<1
0022 f=spm_select(1,[],'Select weight map',[],pwd,'.img');
0023 if isempty(f)
0024 disp('No weight map selected, aborting')
0025 return
0026 end
0027 else
0028 f=weight_fname{1};
0029 end
0030
0031
0032 if nargin<2
0033 gi=spm_select(1,'image','Select atlas');
0034 if isempty(f)
0035 disp('No atlas selected, aborting')
0036 return
0037 end
0038 else
0039 gi=atlas_fname;
0040 end
0041
0042
0043 if nargin<3
0044 flag=1;
0045 else
0046 flag=build_im;
0047 end
0048
0049
0050
0051 if nargin<4
0052 comp_perm=0;
0053 end
0054
0055 if comp_perm
0056 [a,b]=fileparts(char(f));
0057 dirn=[a,filesep,'perm_',b];
0058 pth=pwd;
0059 if isdir(dirn)
0060
0061 cd(dirn)
0062 files=dir('*perm*.img');
0063 fp=char({files(:).name});
0064 fper={[repmat(dirn,length(files),1),repmat(filesep,length(files),1),fp]};
0065 fperm=char([fper;{f}]);
0066 cd(pth)
0067 else
0068 disp('No folder containing the weight images for permutations')
0069 disp('Computing normalized weights for the provided image only')
0070 comp_perm=0;
0071 fperm=f;
0072 end
0073 else
0074 fperm=f;
0075 end
0076
0077
0078
0079
0080
0081
0082 V=spm_vol(f);
0083 [xxx,bb]=fileparts(f);
0084 nfo=length(V);
0085 V1=spm_vol(gi);
0086 dumb=V(1);
0087
0088 if ~any(dumb.dim == V1.dim)
0089 disp('Resizing atlas--------->>')
0090
0091 fl_res = struct('mean',false,'interp',0,'which',1,'prefix','resized_');
0092 spm_reslice([dumb V1],fl_res);
0093
0094
0095 [V1_pth,V1_fn,V1_ext] = spm_fileparts(V1.fname);
0096 rV1_fn = [fl_res.prefix,V1_fn];
0097
0098 if strcmp(V1_ext,'.nii')
0099 V_in = spm_vol(fullfile(V1_pth,[rV1_fn,'.nii']));
0100 V_out = V_in; V_out.fname = fullfile(V1_pth,[rV1_fn,'.img']);
0101 spm_imcalc(V_in,V_out,'i1');
0102 end
0103
0104
0105 mfile_new=['resized_',V1_fn];
0106 pp=spm_fileparts(f);
0107 movefile(fullfile(V1_pth,[rV1_fn,'.img']), ...
0108 fullfile(pp,[mfile_new,'.img']));
0109 movefile(fullfile(V1_pth,[rV1_fn,'.hdr']), ...
0110 fullfile(pp,[mfile_new,'.hdr']));
0111 g=spm_vol(fullfile(pp,[mfile_new,'.img']));
0112 h=spm_read_vols(g);
0113 else
0114 h=spm_read_vols(V1);
0115 end
0116
0117
0118
0119
0120 for ii=1:size(fperm,1)
0121
0122
0123 V=spm_vol(fperm(ii,:));
0124 w=zeros(V(1).dim(1)*V(1).dim(2)*V(1).dim(3),length(V));
0125 VV=spm_read_vols(V);
0126 nfold=size(VV,4)-1;
0127 if nfold == 0
0128 w=VV(:);
0129 else
0130 for i=1:nfold+1
0131 tmp=VV(:,:,:,i);
0132 w(:,i)=tmp(:);
0133 end
0134 clear tmp
0135 end
0136 atlas=h(:);
0137 atlas(isnan(w(:,1)))=NaN;
0138
0139
0140 if ii==1
0141 N_other=length(find(atlas==0));
0142 P_oth=N_other/length(find(isnan(atlas)));
0143 disp(['Volume of the others region: ',num2str(P_oth)])
0144 end
0145
0146
0147 disp(['Computing weights in each ROI for image ',num2str(ii)])
0148 [H HN SN idfeatroi] = prt_region_histogram(w, atlas);
0149 nr=size(H,1);
0150
0151
0152 indnan = [];
0153 for i=1:size(HN,2)
0154 if length(find(isnan(HN(:,i))))==size(HN,1) || ...
0155 length(find(HN(:,i)==0))==size(HN,1)
0156 HN(:,i) = 0;
0157 H(:,i) = 0;
0158 indnan = [indnan i];
0159 end
0160 end
0161
0162
0163 if ii==1
0164
0165 r_min=min(atlas);
0166 R=max(atlas);
0167 if r_min==0
0168 corr=1;
0169 else
0170 corr=0;
0171 end
0172 end
0173
0174
0175 pH=H*100;
0176 oth_w=pH(1);
0177
0178
0179 inn= ~isnan(HN(:,1));
0180 shn=sum(HN(inn,:),1);
0181 pHN=(HN./repmat(shn,size(HN,1),1))*100;
0182 pHN(:,indnan)=0;
0183
0184
0185 end
0186
0187
0188
0189 W_roi=pH;
0190 NW_roi=pHN/100;
0191 SN=SN*100;
0192 [a,b,c]=fileparts(dumb.fname);
0193
0194
0195
0196
0197
0198
0199 if flag
0200 disp('Building image of normalized weights')
0201
0202
0203 if exist(fullfile( ...
0204 a,['ROI_',b,c]),'file')
0205 disp('Image of normalized weights per region already exists, overwriting...')
0206 end
0207
0208
0209 img_name=[a,filesep,'ROI_',b,c];
0210 img4d = file_array(img_name,size(VV),'float32-le',0,1,0);
0211 for km=1:size(w,2)
0212 for r=r_min:R
0213 w(idfeatroi{r+corr},km)=pHN(r+corr,km);
0214 end
0215 img4d(:,:,:,km)=reshape(w(:,km),[size(VV,1),size(VV,2), size(VV,3),1]);
0216 end
0217
0218
0219
0220
0221 disp('Creating image--------->>')
0222 No = V(1).private;
0223 No.dat = img4d;
0224 No.descrip = 'Pronto weigths';
0225 create(No);
0226 disp('Done.')
0227 end
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249