vb_atlas_label

 Create an area indices based on 
   the Automated anatomical labeling 'aal'
   or 'Brodmann' atlas file
 
 ---- Syntax
 vb_atlas_label(brainfile,atlas_id,Rmax)

 ---- Input (Field names in 'parm')
 brainfile  : brain file path
 ---- Optional
 atlas_id   : ID of atlas ['aal']
                 : 'aal' or 'brodmann'
 Rmax: Box size (radius of sphere) (mm)  [ 7 ]

  Label is asigned to the label with maximum counting number 
  inside the box (sphere) around the vertex by counting number of each label

 2015/12/8 Made by M.Sato 

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    vb_atlas_label(brainfile, atlas_id, Rmax)
0002 % Create an area indices based on
0003 %   the Automated anatomical labeling 'aal'
0004 %   or 'Brodmann' atlas file
0005 %
0006 % ---- Syntax
0007 % vb_atlas_label(brainfile,atlas_id,Rmax)
0008 %
0009 % ---- Input (Field names in 'parm')
0010 % brainfile  : brain file path
0011 % ---- Optional
0012 % atlas_id   : ID of atlas ['aal']
0013 %                 : 'aal' or 'brodmann'
0014 % Rmax: Box size (radius of sphere) (mm)  [ 7 ]
0015 %
0016 %  Label is asigned to the label with maximum counting number
0017 %  inside the box (sphere) around the vertex by counting number of each label
0018 %
0019 % 2015/12/8 Made by M.Sato
0020 %
0021 % Copyright (C) 2011, ATR All Rights Reserved.
0022 % License : New BSD License(see VBMEG_LICENSE.txt)
0023 
0024 % ----- DEBUG
0025 % parm=[];
0026 % braindir  = 'D:\MatData\MNI-ICBM152\vbmeg\';
0027 % brainfile = [braindir 'mni_icbm152_t1_tal_nlin_asym_09c.brain.mat'];
0028 %% ----- DEBUG
0029 
0030 
0031 if nargin < 2,    atlas_id = 'aal'; end;
0032 
0033 % Box size for neighbor: 7 mm
0034 if nargin < 3,    Rmax = 7; end;
0035 
0036 %
0037 % ---- Define atlas
0038 %
0039 atlas_dir  = 'atlas2vb_dir/MNI_atlas_templates/';
0040 atlas_text = which([atlas_dir atlas_id '.txt']);
0041 atlas_file = which([atlas_dir atlas_id '.img']);
0042 
0043 % value corresponding to Cerebrum cortex region
0044 %  AAL: 1-90: Cerebrum, 91-108: Cerebellum, 109-116: Vermis
0045 %       Odd  : Left  &  Even : Right
0046 % Brodmann: 1-47; brodmann area , 100- subcortical area
0047 MaxCortex = 90;
0048 
0049 EXT_brain = '.brain.mat';
0050 brain_id  = brainfile(1:findstr(brainfile,EXT_brain)-1);
0051 
0052 save_areafile  = [brain_id '_' atlas_id '.area.mat'];
0053 save_atlasfile = [brain_id '_' atlas_id '.act.mat'];
0054 
0055 % ----- Load MNI coordinate [m] of brain
0056 Vmni = vb_load_cortex(brainfile, 'MNI')*1000;
0057 Vmni = floor(Vmni);
0058 Npoint = size(Vmni,1);
0059 [Ndipole, NdipoleL] = vb_load_cortex_info(brainfile);
0060 
0061 if Npoint ~= Ndipole,
0062     error('Brain file info error');
0063 end
0064 
0065 Vlabel = zeros(Npoint,1);
0066 Vcount = zeros(Npoint,1);
0067 
0068 % Cortex area index of 'V'
0069 % Vinfo.cortexL =
0070 % Vinfo.cortexR =
0071 fprintf('Load Cortex area index\n')
0072 
0073 load(brainfile, 'Vinfo');
0074 
0075 Vindx  = [Vinfo.cortexL; Vinfo.cortexR];
0076 Ncortex = length(Vindx);
0077 
0078 %
0079 % ---- Load Atlas template file
0080 %      Template file is neurological format (right handed spm).
0081 %
0082 
0083 % ----- Read Template image-header
0084 [dim, Trans] = vb_get_image_info(atlas_file);
0085 invTrans = inv(Trans');
0086 
0087 %[dim, Trans, XYZmm] = vb_get_image_info(atlas_file);
0088 % dim(1:3) - the x, y and z dimensions of the volume
0089 % Trans    - a 4x4 affine transformation matrix mapping from
0090 %            voxel coordinates to MNI-mm-coordinates (RAS)
0091 % XYZmm    = Trans * [XYZvox ; ones(1,Npoint)] [ 3 x Npoint]
0092 %          : mm coordinate in MNI coordinate
0093 %
0094 %   [x  y  z  1] = [i  j  k  1] * Trans'
0095 % XYZmni = Vox  * Trans'
0096 % Vox    = XYZmni * inv(Trans')
0097 
0098 % ----- Read atlas-template image
0099 % Zlabel - 3D image data (label value)
0100 avw = avw_img_read(atlas_file);
0101 Zlabel = avw.img;
0102 
0103 Zlabel = Zlabel(:);
0104 Nlabel = max(Zlabel);
0105 
0106 % Make box index
0107 x_indx = -Rmax:Rmax;
0108 [X,Y,Z] = ndgrid(x_indx,x_indx,x_indx);
0109 Box_ix = [X(:),Y(:),Z(:)];
0110 
0111 % Select inside of sphere
0112 ix = find( sum(Box_ix.^2,2) <= Rmax^2);
0113 Box_ix = Box_ix(ix,:);
0114 
0115 Nbox = size(Box_ix,1);
0116 ilabel = (1:Nbox)';
0117 
0118 for n = 1:Ncortex
0119     % MNI coordinate of np-th point
0120     np = Vindx(n);
0121     V  = Vmni(np,:);
0122     
0123     % Box coordinate around V
0124     XYZ = vb_repadd(Box_ix, V);
0125     
0126     % ----- Transform MNI-mm-coordinates to voxel coordinates
0127     Vox = fix( [XYZ ones(Nbox,1)] * invTrans );
0128     
0129     ix = find(  (Vox(:,1) > 0) ...
0130               & (Vox(:,2) > 0) ...
0131               & (Vox(:,3) > 0) ...
0132               & (Vox(:,1) <= dim(1) ) ...
0133               & (Vox(:,2) <= dim(2) ) ...
0134               & (Vox(:,3) <= dim(3) ) ...
0135               );
0136     
0137     % ----- Transform to 1D-index from [X,Y,Z] voxel coordinates
0138     Vindex = sub2ind(dim,Vox(ix,1),Vox(ix,2),Vox(ix,3));
0139     
0140     % Map label to vertex in the brain
0141     vlabel = Zlabel(Vindex) ;
0142     
0143     % if vlabel > MaxCortex = 90, it is in subcortical area
0144     jj = find( vlabel > MaxCortex );
0145     vlabel(jj) = 0;
0146     
0147     if strcmp( atlas_id, 'aal')
0148         if np > NdipoleL, 
0149         %  Right :  label = Even
0150             jj = find( mod(vlabel,2) == 1);
0151             vlabel(jj) = 0;
0152         else
0153         %  Left : label = Odd
0154             jj = find( mod(vlabel,2) == 0);
0155             vlabel(jj) = 0;
0156         end
0157     end
0158     
0159     % Make vlabel >= 1 by adding 1
0160     vlabel = vlabel + 1;
0161     
0162     % ---- Make SPARSE table to count label number in the box
0163     % S = SPARSE(i,j,s,m,n)
0164     S = sparse(ix , vlabel, 1, Nbox, Nlabel+1);
0165     % S(i, j ) = 1 for j = vlabel(i)
0166     % count_label(j) : counting number of j-th label inside the box
0167     count_label = full(sum(S,1));
0168     
0169     % Find Max counting number (remove label=0 component)
0170     [count,imax] = max(count_label(2:Nlabel+1));
0171     
0172     if count > 0
0173         Vlabel(np) = imax(1);
0174         Vcount(np) = count;
0175         
0176         if strcmp( atlas_id, 'brodmann')
0177             if np > NdipoleL, 
0178                 % Add 1000 to label in Right area
0179                 Vlabel(np) = Vlabel(np) + 1000;
0180             end
0181         end
0182     end
0183 end
0184 
0185 %
0186 % ---- Save labeled area
0187 %
0188 
0189 % Read atlas text file and extract area keys.
0190 [label, label_name] = vb_read_atlas_label(atlas_text);
0191 
0192 % Add Left-Right to label_name for 'brodmann'
0193 if strcmp( atlas_id, 'brodmann')
0194     % Add 1000 to label in Right hemisphere
0195     label = [label(:) ; label(:) + 1000];
0196     
0197     % Add Left-Right to label_name
0198     N = length(label_name);
0199     for n = 1:N
0200         %name = strrep(label_name{n},'brodmann area ','brodmann_');
0201         name = label_name{n};
0202         label_name{n}   = [name '_L'];
0203         label_name{n+N} = [name '_R'];
0204     end    
0205     
0206 end
0207 
0208 %% Add Corpus label
0209 %Narea = length(label);
0210 %label(Narea+1) = corpus_label;
0211 %label_name{Narea+1} = corpus_key;
0212 
0213 % Intensity of 'Act.xxP' is the label for each vertex
0214 Atlas.key = [atlas_id];
0215 Atlas.xxP = Vlabel;
0216 Atlas.label      = label;
0217 Atlas.label_name = label_name;
0218 
0219 if exist(save_atlasfile,'file')
0220     delete(save_atlasfile);
0221 end
0222 if exist(save_areafile,'file')
0223     delete(save_areafile);
0224 end
0225 
0226 % Save label as actfile
0227 fprintf('Save the atlas file as "%s"  \n', save_atlasfile);
0228 vb_add_act(save_atlasfile,Atlas,[],false);
0229 
0230 % Save area label into area file
0231 fprintf('Save the area file as "%s" \n', save_areafile);
0232 
0233 vb_save_atlas_label(save_areafile,Vlabel,label,label_name);
0234 
0235 NLABEL = sum( Vlabel > 0 );    % # of labeled points
0236 fprintf('# of all vertex     = %d\n',Npoint)
0237 fprintf('# of cortex points  = %d\n',Ncortex)
0238 fprintf('# of labeled points = %d\n',NLABEL)
0239