Correct holes and islands in atlas label [Atlas] = vb_correct_atlas_label(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 [mm] Para.rate : number less than rate*(# of largest region) is deleted Para.Nmin : number less than Nmin is deleted --- Output Atlas : structure of atlas label --- Method extract connected region 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,result] = vb_correct_atlas_label(Atlas,brainfile,Para) 0002 % Correct holes and islands in atlas label 0003 % [Atlas] = vb_correct_atlas_label(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 [mm] 0013 % Para.rate : number less than rate*(# of largest region) is deleted 0014 % Para.Nmin : number less than Nmin is deleted 0015 % --- Output 0016 % Atlas : structure of atlas label 0017 % --- Method 0018 % extract connected region 0019 % fill holes by gaussian filter morphology 0020 % 0021 % 2006-11-12 M. Sato 0022 % 0023 % Copyright (C) 2011, ATR All Rights Reserved. 0024 % License : New BSD License(see VBMEG_LICENSE.txt) 0025 0026 0027 if ~exist('Para','var'), Para = []; end 0028 0029 % max radius for morphology 0030 if isfield(Para,'Rlabel') 0031 Rmax = Para.Rlabel/1000; % [mm] -> [m] 0032 else 0033 Rmax = 0.01;% [m] 0034 end 0035 if isfield(Para,'rate') 0036 rate = Para.rate; 0037 else 0038 rate = 0.1; 0039 end 0040 if isfield(Para,'Nmin') 0041 Nmin = Para.Nmin; 0042 else 0043 Nmin = 10; 0044 end 0045 0046 % load brain 0047 [V, F] = vb_load_cortex(brainfile); 0048 0049 F = F.F3; % patch index 0050 NV = size(V,1); % # of vertex 0051 0052 label = Atlas.label ; % label value; 0053 Narea = length(label); 0054 0055 % label value for each vertex [Nvertex x 1] 0056 xxP = Atlas.xxP; 0057 xxQ = zeros(NV,1); 0058 0059 Num0 = sum( xxP == 0 ); 0060 fprintf('--- Unlabeled vertex = %d\n',Num0) 0061 fprintf('--- Find connected region\n') 0062 0063 % --- Remove label for small disconnected region 0064 for n=1:Narea 0065 % vertex index for n-th label 0066 ix = find( xxP == label(n) ); 0067 0068 if isempty(ix), continue; end; 0069 0070 % find connected region 0071 [Vindx, Nall] = vb_connected_vertex(ix,F); 0072 0073 % set label to large connected region 0074 if ~isempty(Vindx) 0075 % set label for largest region 0076 xxQ([Vindx{1}]) = label(n); 0077 0078 jj = find( (Nall >= Nall(1)*rate) & (Nall > Nmin) ); 0079 0080 for m=1:length(jj) 0081 xxQ([Vindx{jj(m)}]) = label(n); 0082 end 0083 end 0084 0085 end 0086 0087 Num = sum( xxQ == 0 ); 0088 0089 fprintf('--- Unlabeled vertex after removal = %d\n',Num) 0090 fprintf('--- Label filling\n') 0091 0092 % unlabeled index 0093 ixfix = find( xxQ == 0 ); 0094 Nfix = length(ixfix); 0095 0096 % label weight for unlabeled points 0097 PP = zeros(Nfix,Narea); 0098 0099 % smoothing filter 0100 W = vb_spatial_gauss_filter(brainfile,Rmax,Rmax,[1:NV]); 0101 W = W(ixfix,:); 0102 0103 for n=1:Narea 0104 % vertex index for n-th label 0105 ix = find( xxQ == label(n) ); 0106 0107 if isempty(ix), continue; end; 0108 0109 % Gaussian filtering from n-th labeled points 0110 PP(:,n) = sum(W(:,ix), 2); 0111 0112 end 0113 0114 % --- Label filling 0115 for n=1:Nfix 0116 % find maximum weighted label 0117 [Pmax, id] = max(PP(n,:)); 0118 0119 if Pmax > 0, 0120 xxQ(ixfix(n)) = label(id); 0121 end 0122 end 0123 0124 %xxP = xxQ; 0125 %xxQ = zeros(NV,1); 0126 % 0127 %for n=1:Narea 0128 % % vertex index for n-th label 0129 % ix = find( xxP == label(n) ); 0130 % 0131 % if isempty(ix), continue; end; 0132 % 0133 % % find connected region 0134 % [Vindx, Nall] = vb_connected_vertex(ix,F); 0135 % 0136 % % set label to large connected region 0137 % if ~isempty(Vindx) 0138 % % set label for largest region 0139 % xxQ([Vindx{1}]) = label(n); 0140 % 0141 % jj = find( (Nall >= Nall(1)*rate) & (Nall > Nmin) ); 0142 % 0143 % for m=1:length(jj) 0144 % xxQ([Vindx{jj(m)}]) = label(n); 0145 % end 0146 % end 0147 % 0148 %end 0149 0150 Atlas.xxP = xxQ; 0151 0152 Num2 = sum( xxQ == 0 ); 0153 fprintf('--- Final unlabeled vertex = %d\n',Num2) 0154 0155 if Num2 >= Num, 0156 result = 0; 0157 else 0158 result = 1; 0159 end 0160 0161 %%% surface morphology 0162 % Iextract = vb_open_area(Jarea, R, nextIX, nextDD); 0163 % Iextract = vb_close_area(Jarea, R, nextIX, nextDD);