Home > . > prt_build_region_weights.m

prt_build_region_weights

PURPOSE ^

SYNOPSIS ^

function [NW_roi,idfeatroi]=prt_build_region_weights(weight_fname,atlas_fname,build_im,comp_perm)

DESCRIPTION ^

function to compute the weights for each region as specified by the atlas
image (one value per region). Weights not in the atlas are comprised in an
additional region with name 'others'.
--------------------------------------------------------------------------
input: name of the weights image (weight_fname), name of the atlas image 
       (atlas_fname). build_im is set to 1 to build the resulting image.
       comp_perm is set to 1 to compute the ranking for the permutations.
output: -image file with the normalized weights in each region, which can
         be viewed in the results GUI as a weight image.
        -.mat file containing the weights of each ROI, in %, pH, the
        normalized weights for each region, in %, pHN and the list of
        region names.
--------------------------------------------------------------------------
Written by Jessica Schrouff, 19/09/2012
for PRoNTo

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [NW_roi,idfeatroi]=prt_build_region_weights(weight_fname,atlas_fname,build_im,comp_perm)
0002 %
0003 %function to compute the weights for each region as specified by the atlas
0004 %image (one value per region). Weights not in the atlas are comprised in an
0005 %additional region with name 'others'.
0006 %--------------------------------------------------------------------------
0007 %input: name of the weights image (weight_fname), name of the atlas image
0008 %       (atlas_fname). build_im is set to 1 to build the resulting image.
0009 %       comp_perm is set to 1 to compute the ranking for the permutations.
0010 %output: -image file with the normalized weights in each region, which can
0011 %         be viewed in the results GUI as a weight image.
0012 %        -.mat file containing the weights of each ROI, in %, pH, the
0013 %        normalized weights for each region, in %, pHN and the list of
0014 %        region names.
0015 %--------------------------------------------------------------------------
0016 %Written by Jessica Schrouff, 19/09/2012
0017 %for PRoNTo
0018 
0019 
0020 %select weight map
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 %select atlas
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 %set flag to 1 if not specified
0043 if nargin<3
0044     flag=1;
0045 else
0046     flag=build_im;
0047 end
0048 
0049 % in case the ranking has to be computed for the permutations, get the
0050 % names of the weight images for each permutation.
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         %get the names of the weight images
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 %resize atlas if needed
0079 %--------------------------------------------------------------------------
0080 
0081 %load images
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     %reslice
0091     fl_res = struct('mean',false,'interp',0,'which',1,'prefix','resized_');
0092     spm_reslice([dumb V1],fl_res);
0093 
0094     %build updated atlas
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     %put the files into the PRT directory
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 %compute histogram
0118 %--------------------------------------------------------------------------
0119 
0120 for ii=1:size(fperm,1)
0121     
0122     %Get the volumes into matrices
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 %when only the average across folds was computed
0128         w=VV(:);
0129     else
0130         for i=1:nfold+1  %number of folds + average
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      %Compute the volume of the 'others' region
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     % Compute the weights and normalized weights
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     % Put 0 if the fold has only NaNs
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     %Correct for the 'others' region (one time)
0163     if ii==1   
0164         %compute proportions as in PCA
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     %sum of weights in each region
0175     pH=H*100;
0176     oth_w=pH(1); %save the weight of the 'others' region
0177    
0178     %normalized sum of weights in each region
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 %save sorted H, HN and the list of corresponding ROIs for the 'true' image,
0188 %and the ranking of the average weights for each permutation
0189 W_roi=pH;
0190 NW_roi=pHN/100;
0191 SN=SN*100;
0192 [a,b,c]=fileparts(dumb.fname);
0193 % save(fullfile(a,['atlas_',b1,'_',b,'.mat']),'LR',...
0194 %     'W_roi','NW_roi','dwn','drwn','erwn','SN','ER','P_oth','oth_w');
0195 
0196 %build new image with the normalized weights and save values
0197 %--------------------------------------------------------------------------
0198 
0199 if flag
0200     disp('Building image of normalized weights')
0201     
0202     % if image exists, overwrite
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     %build image if flag
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     % Create weigths file
0220     %--------------------------------------------------------------------------
0221     disp('Creating image--------->>')
0222     No         = V(1).private;     % copy header
0223     No.dat     = img4d;            % change file_array
0224     No.descrip = 'Pronto weigths'; % description
0225     create(No);                    % write header
0226     disp('Done.')
0227 end
0228 
0229 
0230 
0231 % Previous functions that were saving all results in a separate .mat file
0232 
0233     %compute the rank of each region according to the weights
0234 %     [dub,ih]=sort(pH,1,'descend');
0235 %     [d1,dw]=sort(ih);
0236 %     %ranking distance between each fold and the "average" fold
0237 %     drw=zeros(1,nfold);
0238 %     for ifold=1:nfold
0239 %         drw(ifold)=prt_comp_ranking_dist(dw(:,ifold),dw(:,end));
0240 %     end
0241 %     %Expected value of the rank for each region
0242 %     rw=zeros(nr,1);
0243 %     for i=1:nr
0244 %         for j=1:nr
0245 %             tmp=length(find(dw(i,:)==j));
0246 %             rw(i)=rw(i)+j*tmp;
0247 %         end
0248 %     end
0249 %     rw=rw/nfold;

Generated on Tue 10-Feb-2015 18:16:33 by m2html © 2005