vb_remove_island_atlas

 Remove islands of atlas labels using patch information
 Note that original act.mat and area.mat files will be overwritten.

 ---- Syntax
 vb_remove_island_atlas(parm)

 ---- Input
 brain_file  : brain file path
 atlas_id    : ID of atlas ['aal' or 'brodmann']

 2016/5/27 Last modified by Y.Takeda

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

0001 function vb_remove_island_atlas(brain_file,atlas_id)
0002 % Remove islands of atlas labels using patch information
0003 % Note that original act.mat and area.mat files will be overwritten.
0004 %
0005 % ---- Syntax
0006 % vb_remove_island_atlas(parm)
0007 %
0008 % ---- Input
0009 % brain_file  : brain file path
0010 % atlas_id    : ID of atlas ['aal' or 'brodmann']
0011 %
0012 % 2016/5/27 Last modified by Y.Takeda
0013 %
0014 % Copyright (C) 2011, ATR All Rights Reserved.
0015 % License : New BSD License(see VBMEG_LICENSE.txt)
0016 
0017 if ~exist('atlas_id','var')
0018     atlas_id='aal';
0019 end
0020 
0021 % Set file names
0022 EXT_brain = '.brain.mat';
0023 brain_id  = brain_file(1:findstr(brain_file,EXT_brain)-1);
0024 area_file  = [brain_id '_' atlas_id '.area.mat'];
0025 atlas_file = [brain_id '_' atlas_id '.act.mat'];
0026 
0027 % Load atlas label to modify
0028 act = vb_get_act(atlas_file,atlas_id);
0029 area_list=unique(act.xxP(act.xxP>0));
0030 Narea=length(area_list);
0031 
0032 % Load patch information
0033 [tmp,F]=vb_load_cortex(brain_file,'Inflate');
0034 F0=F.F3;
0035 
0036 % Start to find and remove island clusters
0037 xxP=act.xxP;
0038 xxP_old=xxP+1;
0039 while 1
0040     
0041     % Find island clusters
0042     clear Area1
0043     for area1=1:Narea
0044         ix=find(xxP==area_list(area1));
0045         Nv1=length(ix);
0046         
0047         % Make clusters using patch information
0048         clear same
0049         for v1=1:Nv1
0050             a=find(F0(:,1)==ix(v1) | F0(:,2)==ix(v1) | F0(:,3)==ix(v1));
0051             f=F0(a,:);
0052             f=f(:);
0053             same(v1).ix=v1;
0054             for v2=v1+1:Nv1
0055                 if ~isempty(find(f==ix(v2)))
0056                     same(v1).ix=[same(v1).ix,v2];
0057                 end
0058             end
0059         end
0060         
0061         % Combine adjacent clusters
0062         while 1
0063             become_empty=0;
0064             for v1=1:Nv1-1
0065                 for v2=v1+1:Nv1
0066                     if ~isempty(same(v2).ix)
0067                         ix2=[same(v1).ix,same(v2).ix];
0068                         if length(unique(ix2))<length(ix2)
0069                             same(v1).ix=unique(ix2);
0070                             same(v2).ix=[];
0071                             become_empty=1;
0072                         end
0073                     end
0074                 end
0075             end
0076             if become_empty==0
0077                 break
0078             end
0079         end
0080         
0081         % Extract combined clusters
0082         v2=1;
0083         for v1=1:Nv1
0084             if ~isempty(same(v1).ix)
0085                 Area1(area1).cluster(v2).ix=ix(same(v1).ix);
0086                 Area1(area1).Nvertex(v2)=length(same(v1).ix);
0087                 v2=v2+1;
0088             end
0089         end
0090         
0091         % Define island clusters
0092         ix1=1:length(Area1(area1).Nvertex);
0093         [~,b]=max(Area1(area1).Nvertex);
0094         ix1(b)=[];
0095         Area1(area1).center=b;
0096         Area1(area1).island=ix1;
0097     end
0098     
0099     % Remove island clusters
0100     for area1=1:Narea
0101         Nisland=length(Area1(area1).island);
0102         if Nisland>0
0103             for island1=1:Nisland
0104                 
0105                 % Obtain surrounding labels of an island
0106                 ix_island_cluster=Area1(area1).island(island1);
0107                 vertices_in_an_island_cluster=Area1(area1).cluster(ix_island_cluster).ix;
0108                 Nn=zeros(Narea,1);
0109                 for v=1:length(vertices_in_an_island_cluster)
0110                     v2=vertices_in_an_island_cluster(v);
0111                     a=find(F0(:,1)==v2 | F0(:,2)==v2 | F0(:,3)==v2);
0112                     f=F0(a,:);
0113                     f=f(:);
0114                     ae=xxP(f);
0115                     ae=unique(ae(ae~=area_list(area1)));
0116                     for ae1=1:length(ae)
0117                         a=find(area_list==ae(ae1));
0118                         Nn(a,1)=Nn(a,1)+1;
0119                     end
0120                 end
0121                 
0122                 % Change the label of the island to the most common surrounding label
0123                 if sum(Nn)>0
0124                     [~,b]=max(Nn);
0125                     xxP(vertices_in_an_island_cluster)=area_list(b);
0126                 else
0127                     xxP(vertices_in_an_island_cluster)=0;
0128                 end
0129             end
0130         end
0131     end
0132     
0133     % Exit if xxP is not changed
0134     if norm(xxP-xxP_old)>0
0135         xxP_old=xxP;
0136     else
0137         break
0138     end
0139 end
0140 
0141 % Prepare for saving modified atlas label
0142 Atlas.key = atlas_id;
0143 Atlas.xxP = xxP;
0144 Atlas.label = act.label;
0145 Atlas.label_name = act.label_name;
0146 
0147 if exist(atlas_file,'file')
0148     delete(atlas_file);
0149 end
0150 if exist(area_file,'file')
0151     delete(area_file);
0152 end
0153 
0154 % Save modified label as actfile
0155 fprintf('Save modified atlas file\n');
0156 vb_add_act(atlas_file,Atlas,[],false);
0157 
0158 % Save modified area label into area file
0159 fprintf('Save modified area file\n');
0160 vb_save_atlas_label(area_file,Atlas.xxP,Atlas.label,Atlas.label_name);