vb_correct_atlas_label_LR

 Correct holes and islands in 'brodmann'-atlas label by morphology
   [Atlas] = vb_correct_atlas_label_LR(Atlas,brainfile)
 --- Input
 Atlas : structure of atlas label
 Atlas.xxP        : label value for each vertex
 Atlas.label      : label value
 Atlas.label_name : label_name corresponding to 'label'
 
 brainfile

 Para.Rlabel : max radius for morphology  [m]
 Para.rate   : number less than rate*(# of largest region) is deleted
 --- Output
 Atlas : structure of atlas label
 --- Method
 extract connected region in left/right cortex sepalately
 fill holes by gaussian filter morphology

 2006-11-12 M. Sato

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

0001 function    [Atlas] = vb_correct_atlas_label_LR(Atlas,brainfile,Para)
0002 % Correct holes and islands in 'brodmann'-atlas label by morphology
0003 %   [Atlas] = vb_correct_atlas_label_LR(Atlas,brainfile)
0004 % --- Input
0005 % Atlas : structure of atlas label
0006 % Atlas.xxP        : label value for each vertex
0007 % Atlas.label      : label value
0008 % Atlas.label_name : label_name corresponding to 'label'
0009 %
0010 % brainfile
0011 %
0012 % Para.Rlabel : max radius for morphology  [m]
0013 % Para.rate   : number less than rate*(# of largest region) is deleted
0014 % --- Output
0015 % Atlas : structure of atlas label
0016 % --- Method
0017 % extract connected region in left/right cortex sepalately
0018 % fill holes by gaussian filter morphology
0019 %
0020 % 2006-11-12 M. Sato
0021 %
0022 % Copyright (C) 2011, ATR All Rights Reserved.
0023 % License : New BSD License(see VBMEG_LICENSE.txt)
0024 
0025 if ~exist('Para','var'), Para = []; end
0026 
0027 % max radius for morphology
0028 if isfield(Para,'Rlabel')
0029     Rmax  = Para.Rlabel;
0030 else
0031     Rmax  = 0.012;% [m]
0032 end
0033 
0034 if isfield(Para,'rate')
0035     rate  = Para.rate;
0036 else
0037     rate  = 0.1;
0038 end
0039 
0040 % load brain
0041 [V, F] = vb_load_cortex(brainfile);
0042 
0043 % Left/Right cortical patch index
0044 NL = F.NdipoleL;
0045 FL = F.F3L;
0046 FR = F.F3R;
0047 NV = size(V,1);    % # of vertex
0048 
0049 name  = Atlas.label_name ; % label_name corresponding to 'label'
0050 label = Atlas.label  ; % label value;
0051 Narea = length(label);
0052 
0053 % New label to separate L-R
0054 Label = [];
0055 Name  = [];
0056 Lid   = 0;
0057 Rid   = max(label)+10;
0058 Lext  = '_L';
0059 Rext  = '_R';
0060 
0061 % label value for each vertex [Nvertex x 1]
0062 xxP   = Atlas.xxP; 
0063 xxQ   = zeros(NV,1);
0064 
0065 Num = sum( xxP == 0 );
0066 fprintf('--- Unlabeled vertex = %d\n',Num)
0067 fprintf('--- Find connected region\n')
0068 
0069 id = 0;
0070 
0071 for n=1:Narea
0072     % vertex index for n-th label
0073     ix = find( xxP == label(n) );
0074     
0075     if isempty(ix), continue; end;
0076     
0077     %
0078     % --- Left cortex
0079     %
0080     % find connected region
0081     ixL = ix( ix <= NL );
0082     [Vindx, Nall] = vb_connected_vertex(ixL,FL);
0083 
0084     % set label to largest connected region
0085     if ~isempty(Vindx)
0086         jj = find( Nall >= Nall(1)*rate );
0087         xxQ([Vindx{jj'}]) = label(n) + Lid;
0088     end
0089     id = id + 1;
0090     Label(id) = label(n) + Lid;
0091     Name{id}  = [name{n} Lext];
0092     
0093     %
0094     % --- Right cortex
0095     %
0096     % find connected region
0097     ixR = ix( ix > NL );
0098     [Vindx, Nall] = vb_connected_vertex(ixR,FR);
0099 
0100     % set label to largest connected region
0101     if ~isempty(Vindx)
0102         jj = find( Nall >= Nall(1)*rate );
0103         xxQ([Vindx{jj'}]) = label(n) + Rid;
0104     end
0105     id = id + 1;
0106     Label(id) = label(n) + Rid;
0107     Name{id}  = [name{n} Rext];
0108     
0109 end
0110 
0111 fprintf('--- Morphological smoothing\n')
0112 
0113 % New label
0114 Narea = length(Label);
0115 
0116 % unlabeled index
0117 ixfix = find( xxQ == 0 );
0118 Nfix  = length(ixfix);
0119 
0120 % label weight
0121 PP = zeros(Nfix,Narea);
0122 
0123 % smoothing filter
0124 W  = vb_spatial_gauss_filter(brainfile,Rmax,Rmax,[1:NV]);
0125 W  = W(ixfix,:);
0126 
0127 for n=1:Narea
0128     % vertex index for n-th label
0129     ix = find( xxQ == Label(n) );
0130     
0131     if isempty(ix), continue; end;
0132     
0133     PP(:,n) = sum(W(:,ix), 2);
0134     
0135 end
0136 
0137 for n=1:Nfix
0138     [Pmax, id] = max(PP(n,:));
0139     
0140     if Pmax > 0,
0141         xxQ(ixfix(n)) = Label(id);
0142     end
0143 end
0144 
0145 Atlas.xxP   = xxQ; % label value for each vertex
0146 Atlas.label =  Label; % label value
0147 Atlas.label_name =  Name ; % label_name corresponding to 'label'
0148 
0149 Num = sum( xxQ == 0 );
0150 fprintf('--- Final unlabeled vertex = %d\n',Num)
0151 
0152 %%% surface morphology
0153 %    Iextract = vb_open_area(Jarea, R, nextIX, nextDD);
0154 %    Iextract = vb_close_area(Jarea, R, nextIX, nextDD);